The Schrage Method and its Limitations

Schrage's algorithm ("A More Portable Fortran Random Number
Generator." ACM Trans. Math. Software 5, 132-138, 1979.) allows a
multiplication with modulus calculation without overflow, given
certain limitations on the numbers used:

 /**
  * Calculates (a * b) % m without overflow
  * given certain conditions
  */
 int schrage1(int a, int b, int m) {
   // Check limitations
   if (b >= m - 1) { throw new RuntimeException(
                        "Schrage: b >= m-1."); }
   if (a == 0) { return 0; }
   int quot = m / a;
   int rem = m % a;
   if (rem >= quot) { throw new RuntimeException(
                        "Schrage: rem >= quot."); }

   // Schrage's algorithm
   int result = a * (b % quot) - rem * (b / quot);
   if (result < 0) { result += m; }

   return result;
 }

The check for a == 0 is required because of the division m / a; we
need to avoid a divide by zero.

Schrage's algorithm works well, but only if the two given conditions
are met: b < m-1 and rem < quot.  If either of these two is not met
then the algorithm will fail.  

Since I wanted a way to multiply and mod any numbers without
overflowing I looked for a way to extend Schrage's algorithm to deal
with any positive inputs.  This meant finding a way to remove the two
limitations on the parameters.

Removing the First, b < m-1, Limitation

If we have b >= m-1 then we can try to get round it by reducing the
value of b.  If b is even then a * b = (a * b/2) + (a * b/2), if b is
odd then a * b = (a * b/2) + (a * b/2) + a  [here b/2 is an integer
division with fractional parts ignored].  We can use this to recurse
with a reduced value of b and remove the first limitation:

 /**
  * Calculates (a * b) % m without overflow
  * given certain conditions
  */
 int schrage2(int a, int b, int m) {
   if (a == 0) { return 0; }
   if (b >= m - 1) {
     int result = schrage2(a, b/2, m);
     result = addMod(result, result, m);
     if (b % 2 == 1) { result = addMod(result, a, m); }
     return result;
   }

   // Check second limitation
   int quot = m / a;
   int rem = m % a;
   if (rem >= quot) { throw new RuntimeException(
                        "Schrage: rem >= quot."); }

   // Schrage's algorithm
   int result = a * (b % quot) - rem * (b / quot);
   if (result < 0) { result += m; }

   return result;
 }

I have moved the a == 0 check to avoid unneccessary recursion.  By
halving b for each recursive call we ensure that eventually b < m
which meets the first limitation of the basic Schrage algorithm.  From
that we can build up the final result as the recursion unwinds.

The addMod() function is a non-overflowing modular addition:

 /**
  * Calculates (a + b) mod m without overflow.
  */
 int addMod(int a, int b, int m) {
   int result;
   if (b <= Integer.MAX_VALUE - a) {
     result = (a + b) % m;
   } else {
     result = addMod(a, b/2, m);
     result = addMod(result, b/2, m);
     if (b % 2 == 1) { result = addMod(result, 1, m); }
   }
   return result;
 }

This uses a similar method of recursing with a half size operand to
avoid overflow.

Removing the Second, rem < quot, Limitation

Since quot = m / a, by reducing a we increase the value of quot.  We
also reduce the maximum possible value of rem, since rem = m % a.
Hence by reducing a we can bring quot and rem closer to allowed
values.  If a is even then a * b = (a/2 * b) + (a/2 * b), if a is odd
then a * b = (a/2 * b) + (a/2 * b) + b  [here a/2 is an integer
division with fractional parts ignored].  We can use this to recurse
with a reduced value of a and remove the second limitation:

 /**
  * Calculates (a * b) % m without overflow
  */
 int schrage3(int a, int b, int m) {
   if (a == 0) { return 0; }
   if (b >= m - 1) {
     int result = schrage3(a, b/2, m);
     result = addMod(result, result, m);
     if (b % 2 == 1) { result = addMod(result, a, m); }
     return result;
   }
   int quot = m / a;
   int rem = m % a;
   if (rem >= quot) {
     int result = schrage3(a/2, b, m);
     result = addMod(result, result, m);
     if (a % 2 == 1) { result = addMod(result, b, m); }
     return result;
   }

   // Schrage's algorithm
   int result = a * (b % quot) - rem * (b / quot);
   if (result < 0) { result += m; }

   return result;
 }

I have tested this revised algorithm and it works with random positive
inputs, a >= 0, b >= 0, m >= 1.  It does not work with negative inputs
- for my application I only need to deal with positive numbers.  It
would not be difficult to add a shell round the algorithm to handle
signed inputs.

I have presented this example with integer inputs, which allows
testing against the same calculations in longs.  For my real
application I have it using longs in order to avoid going into Java
BigIntegers and slowing things down.

My apologies to those of you to whom I am teaching ovisuction.  I hope
that some others among you might find it useful.

rossum