George Woltman wrote:
>The bug was in the routine mulmod which was supposed to compute a * b
mod c
>The implementation was a rather sloppy:
> unsigned long tmp;
> tmp = a * (b & 0x3FF);
> tmp += ((a << 10) % c) * (b >> 10);
> return (tmp % c);
>You can see in the first computation of tmp there will be an overflow
if
>a exceeds 2^22 or 4194304. The new implementation is:
> double tmp;
> unsigned long q;
> tmp = (double) a * (double) b;
> q = (unsigned long) (tmp / (double) c);
> return ((unsigned long) (tmp - (double) q * (double) c));
>This implementation isn't perfect, but will work as long as a * b does
>not exceed 53 bits.
This got me thinking, since I've had a use for modular multiplication in
some of the little programs that I've bee writing. Based on your code, I
wrote the following:
unsigned long mulmod(unsigned long a, unsigned long b, unsigned long n)
{
long double tmp, q;
tmp = a * (long double)b;
q = floorl(a * (long double)b / n);
return tmp - q * n;
}
In my testing, this appears to give the correct result for any 32-bit
values of a, b, and n, thanks to the 64-bit precision of Intel's 80-bit
floating point type. Is there anything wrong with this approach? It
depends on having a long double version of the floor function, which is
available on my Borland compiler.
Get Your Private, Free Email at http://www.hotmail.com
________________________________________________________________
Unsubscribe & list info -- http://www.scruz.net/~luke/signup.htm
