|Re: Miller-Rabin Primality Test for HP-42S|
Message #3 Posted by Erik Ehrling (Sweden) on 13 Nov 2003, 2:40 p.m.,
in response to message #2 by hugh
i have tried to implement this on the 15c (and others). the difficult bit (as you've noticed) is a*b mod n. do you use the complex number just as a convenient holder or can complex be used to help the precision (eg single precision c1*c2 requires double precision working. in practice the machines use all their internals if they can).
you can do the a*b mod n in binary. i tried it, but calculators take all day (literally) so i abandoned the idea.
Interesting to hear that someone else has been struggling with the same problem!
First I didn't realise the problems with calculating a*b mod n, headed for a straight implementation and ran into trouble, but also got frustrated enough to not be able to put it away...
I use complex numbers mainly as a place holder, but for both addition and subtraction they come in handy for the calculation as well. The real part of the complex number is used to represent the higher (i.e. larger than 1e12) portion of the integer, while the imaginary part represented the lower portion (i.e. a normal precision integer).
As I only needed to deal with positive integers I realised that I could subtract 1e12 from one of the lower parts (i.e. the imaginary part of one of the terms) before adding the two long (complex) numbers together. Then I could check if the lower part of their sum was negative. If it was negative then I just added 1e12 to it again, otherwise I added 1 to the higher (real) part of the sum. (This is done on line 297-315 in the code).
Subtraction uses a variant of the same technique (line 284-296 in the code).
Multiplication was worse. However I only needed to multiply two normal length integers giving a double length product. For this purpose I divided a in two "small" halves a_hi and a_lo and similarly b into b_hi and b_lo. The higher portion here holds the part of the number that is larger than 1e6 and the lower portion holds the part that is smaller than 1e6. Having split a and b in higher and smaller part their double length product can be calculated. The higher (double length) part of a*b can be calculated as a_hi*b_hi + int((a_hi*b_lo + a_lo*b_hi / 1e6)) and the lower part can be calculated as a_lo * b_lo + (a_hi*b_lo + a_lo*b_hi mod 1e6). (The multiplication is made on line 229-283 in the code)
The last step, when a*b is finally calculated (represented as a double length integer), is to determine what the product is mod n. Here I realised that the HP-42S actually calculates very large numbers mod n absolutely correctly (e.g. 1e20 mod 999 999 999 961). So I simply used the rule that ( i + j ) mod n = ((i mod n) + (j mod n)) mod n. In this case this translates into (if i is the larger part of a*b and j the smaller part) --> a*b mod n = ( (( i * 1e12 ) mod n ) + j mod n )) mod n. (This is done on line 195-228 in the code).
I guess the same technique could be used on the HP-15C. Program length might be an issue, though (even if some hundred bytes of PRM? are just for program encapsulation as to preserve the stack). On the other hand this means that implementing Miller-Rabin on the HP-15C would be much more of a challenge than on the HP-42S :-)
Erik Ehrling (Sweden)