Extending Schrage Multiplication

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