Extending Schrage Multiplication

[Lew]

Actually, the problem is with division. If x < 0, y > 0, should x/y be
floor(x/y) or ceil(x/y)? If it is floor(x/y), then everything works out
fine. x - y * floor(x/y) will be in [0,y) for y > 0, and everyone would
be happy.

What most CPUs and programming languages do is, for y > 0, take x/y as
floor(x/y) if x > 0, and ceil(x/y) if x < 0, which I think you'll find
is the option least preferred by mathematicians.

True. I was, of course, referring to division as implemented by "most CPUs
and programming languages", particularly Java, as that was the language of
interest for the OP.

There is a reason why computer engineers prefer the way they do it, which
reason escapes my memory at the moment. Something about the consistency of
|x/y| being the same regardless of operand sign.

 

[Tim Smith]

There is a reason why computer engineers prefer the way they do it, which
reason escapes my memory at the moment. Something about the consistency of
|x/y| being the same regardless of operand sign.

That's what I've always assumed. They want (-x)/y = x/(-y) = -(x/y),
probably because it lets them design dividors that can ignore sign, and
then they can set the sign of the answer at the end.

It's just that mathematically, that's not a particularly interesting
property to have for integer division, so I'd prefer that integer
division be done in a way that preserves x=y*(x/y) + x%y and makes x%y
come out in [0,y). A lot of things in mathematics count on that
relationship, and it gets annoying to have to take extra steps in most
programming languages to deal with that. I'd much rather have that just
work, and if a case ever arose where I needed (-x)/y = -(x/y), I'd
rather that be the case that requires extra code.

 

[Colin Barker]

/**
* 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."); } ...
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:

That seems overly complicated. Note that

a*b = a*(b%m) mod m

That is an excellent suggestion, and I have incorporated it.

So, you can start with

a = a%m;
b = b%m;

If I was writing a shell round the basic method then I would probably
do that. I am reluctant to add an initial divide to a recursive
algorithm where it will only be potentially useful for the top level
of recursion.

to get both your inputs into [0,m-1]. If you really need b < m-1,

I do, Schrage can give incorrect results if b == m-1.

you can then special case the b=m-1 case. Note that in that case b%m
= -1,
so your answer is simply (m-a)%m.

Again an excellent suggesition which I have incorporated.

I wrote that as (m-a)%m there, rather than (-a)%m, to take into account
the fact that most programming languages don't do something sensible
when asked for x%y and x < 0, y > 0. By doing (m-a)%m, I make sure the
answer is the one that would make sense to a mathematician, rather than
the answer that for some inexplicable reason seems to make sense to most
programming language designers and most CPU designers.

Agreed, that was one of my motives in avoiding negative inputs in this
discussion.

I ran some timing tests on my original algorithm (Schrage), an
improved version (Fast Schrage) with your suggestions incorporated and
a basic non-overflowing mulMod function similar to the addMod I gave
earlier:

Over 3000000 random tests, there were 0 mismatched answers.
Schrage took 9.485 seconds.
Fast Schrage took 9.266 seconds.
mulMod took 15.734 seconds.

I suspect that the Schrage method is faster because it can operate on
larger numbers and so needs less depth of recursion to reach a point
where it can calculate the answer without overflowing.

Thankyou for your suggestions,

rossum


// --- Code ---

static int schrageFast(int a, int b, int m) {
if (a < 2) { return (a * b) % m; }

int result;
// Check for b >= m-1
if (b > m - 1) {
result = schrageFast(a, b % m, m);
} else if (b == m - 1) {
result = (m - a) % m;
} else {
// Check for rem >= quot
int quot = m / a;
int rem = m % a;
if (rem >= quot) {
result = schrageFast(a/2, b, m);
result = addMod(result, result, m);
if (a % 2 == 1) { result = addMod(result, b, m); }
} else {
// Schrage method
result = a * (b % quot) - rem * (b / quot);
if (result < 0) { result += m; }
} // end if
} // end if
return result;
} // end schrageFast()

// --- End Code ---

 

[Colin Barker]

/**

* 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;
}

Assuming 0 <= a,b < m wouldn't the following be more efficient?

int addMod(int a, int b, int m) {
return (a < m - b) ? (a + b) : (a - (m - b));
}

Sorry about the null reply I made a couple of minutes ago.

 

[rossum]

Assuming 0 <= a,b < m wouldn't the following be more efficient?

int addMod(int a, int b, int m) {
return (a < m - b) ? (a + b) : (a - (m - b));
}

Yes it would, it roughly halved the timings of my tests:

Old values:
Schrage took 9.485 seconds.
Fast Schrage took 9.266 seconds.

New values (with your addMod() method)
Schrage took 5.907 seconds.
Fast Schrage took 5.813 seconds.


I also tested against my old version of addMod() and there were no
differences for positive inputs.

Thank you,

rossum