AIM: to calculate the following operation:
switch (X % p) {
case 0: break;
case 1: X *= p * s; break;
default: undefined behaviour
}
Here, X is a variable, s is positive, p is a Mersenne prime 2**n-1.

IMPLEMENTATION: Define e ≜ 2**n, a ≜ initial X = p*q+r,
t ≜ r ? a*s : a/p, u = t-r*s.
X starts as a, Y starts at 0. We're aiming for final X = t*p.
We have a = 0 or 1 mod p (thus r = 0 or 1).

        a 0
        a a*p           [Y+=X, p times]
        a p             [Y/=X]
        a q             [Y/=X]
        a q*e           [Y+=Y, n times]
      q-r q*e           [X-=Y; q*e-a = q*e-q*(e-1)-r = q-r]
  e*(q-r) q*e           [X+=X, n times]
  e*(q-r) r*e           [Y-=X]
      q*e r*e           [X+=Y]
      q*e r*q*e*e       [Y*=X]
      q*e (r*e+1)*q*e   [Y+=X]
      q*e 1-r           [Y=X/Y; q*e/((r*e+1)*q*e) = 1/(r*e+1) =
                                when r = 0: 1/1 = 1 = 1-r
                                when r = 1: 1/(e+1) = 0 = 1-r]
      q*e (1-r)*q*e     [Y*=X]
      q*e -r*q*e        [Y-=X]
      e*u -r*q*e        [X-=Y, s*p-1 times; q*e-(s*p-1)*-r*q*e =
                                when r = 0: q*e = a/p*e
                                          = t*e = e*(t-r*s) = e*u
                                when r = 1: q*e+(q*e*(s*p-1))
                                          = q*e*s*p = e*p*q*s
                                          = e*(a-r)*s = e*(a*s-r*s)
                                          = e*(t-r*s) = e*u]
      e*u -r*q*e*e*u    [Y*=X]
      e*u (1-r*q*e)*e*u [Y+=X]
      e*u 1-r           [Y=X/Y; (e*u)/(1-r*q*e)*e*u = 1/(1-r*q*e) =
                                when r = 0: 1/(1-0) = 1 = 1-r
                                when r = 1: 0 = 1-r]
      e*u (1-r)*e*u     [Y*=X]
      e*u -r*e*u        [Y-=X]
      e*u -r            [Y/=X]
      e*t -r            [X-=Y, e*s times; e*u+r*e*s
                                        = e*(t-r*s)+e*r*s = e*t]
      e*t e*t-r         [Y+=X]
      e*t 1-r           [Y/=X]
      e*t 0             [Y/=X]
      e*t e*t           [Y+=X]
      e*t 1             [Y/=X]
      e*t e             [Y+=Y, n times]
        t e             [X/=Y]
        t 0             [Y-=Y]
        t t             [Y+=X]
      e*t t             [X+=X, n times]
      p*t t             [X-=Y]
      p*t 0             [Y/=X]

p*t is our target, so we're done.

WHAT THIS GAINS US:

Suppose we have a set of Mersenne primes. For each such prime, we can find
a prime power of it that's congruent to 1 modulo all the other Mersenne primes;
this gives us a set of prime powers.

We can now compile a program from The Waterfall Model to 2-variable Blindfolded
Arithmetic; associate a different prime power from the set with each waterclock,
and encode a state from The Waterfall Model via raising each waterclock's prime
power to the power of its value, and multiplying all the results together,
storing that in X.

With respect to each of the original Mersenne primes, X will equal 1 (mod the
prime) if the corresponding waterclock is at 0 (because it's a product of
values that are 1 mod that prime), or 0 (mod the prime) if the corresponding
waterclock is nonzero (because we multiplied a power of that prime into X). So
the operation above (multitply X by p*s if X mod p is 1, do nothing if X mod p
is 0, undefined behaviour otherwise) lets us implement a zeroing trigger from
The Waterfall Model, by choosing p*s as an encoding of the values added to the
various waterclocks (note that this must contain a factor of p, because each
waterclock has to increase its own value from 0 when zeroed). Only one
waterclock is allowed to be at 0 at a time, and the operation is safe to run
even after previous operations have run earlier in the cycle (it'll just make
no net change), so we can write one of these operations per waterclock in order
to implement all the zeroing triggers.

Finally, we just have to implement the steady decrement. We can multiply
together all the prime powers to produce a number d, and the decrement is
accomplished simply by dividing X by d, as follows:

        a 0
        a a*d           [Y+=X, d times]
        a d             [Y/=X]
      a/d d             [X/=Y]
      a/d 0             [Y-=Y]

Combining the zeroing triggers and the steady decrement, we have a working
implementation of The Waterfall Model. (Note that I/O is not supported in
this construction, nor is starting or halting the program. You can add extra
zeroing triggers on Mersenne primes that are not "naturally" decremented in
order to implement start and halt. Assuming that X starts at 1, the first
zeroing trigger in the program will run unconditionally on the first cycle,
so you can use it to initialise X, including its prime power in X; if it
isn't part of the steady decrement, it will never run again. Meanwhile,
halting can be implemented via doing X/=Y at the point in the zeroing
trigger where X is e*u and Y is 1-r; if Y is 0 (i.e. the halt waterclock is
nonzero, thus r is 0) this is a no-op that divides X by 1; if Y is 1 (i.e.
the halt waterclock is zero, thus r is 1) this exits the program.)