HP Forum Archive 13

[ Return to Index | Top of Index ]

Miller-Rabin Primality Test for HP-42S Message #1 Posted by Erik Ehrling (Sweden) on 12 Nov 2003, 4:16 p.m. Hi! This is a short note to tell that I have just added an implementation of the Miller-Rabin Primality Test for the HP-42S to my homepage. (By a partial implementation of addition, subtraction and multiplication of double size integers it can handle primality tests for any integer up to 999 999 999 999.) Best regards, Erik Ehrling (Sweden) Homepage:http://w1.322.telia.com/~u32220482/index.html

Re: Miller-Rabin Primality Test for HP-42S Message #2 Posted by hugh on 13 Nov 2003, 1:07 p.m., in response to message #1 by Erik Ehrling (Sweden) top program. i like it! 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.

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 Quote: 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. Hi! 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 :-) Best regards, Erik Ehrling (Sweden)

Re: Miller-Rabin Primality Test for HP-42S Message #4 Posted by hugh on 13 Nov 2003, 8:46 p.m., in response to message #3 by Erik Ehrling (Sweden) i was on the same road. the 42s is saved by the mod function being able to handle large numbers. without a mod function (that works), i was stuck trying to find a way to calculate, ahi*bhi*B^2 (mod n) and the cross-product term * B (mod n). maybe i could call a*b (mod n) recursively? in the end i left it using its own precision, but of course this means i can only verify numbers up to 6 digits (or so). consequently i dropped rounds 7 and 11. here is my original effort: http://www.voidware.com/calcs/hp15c.htm

Re: Miller-Rabin Primality Test for HP-42S Message #5 Posted by Erik Ehrling (Sweden) on 14 Nov 2003, 11:30 a.m., in response to message #4 by hugh Quote: i was on the same road. the 42s is saved by the mod function being able to handle large numbers. without a mod function (that works), i was stuck trying to find a way to calculate, ahi*bhi*B^2 (mod n) and the cross-product term * B (mod n). maybe i could call a*b (mod n) recursively? I think there is a way around this. Actually I implemented the a*b mod n slightly different first, not relying on the mod function being able to handle large numbers. The idea that I came across in an excellent paper ("http://home.pipeline.com/~hbaker1/AB-mod-N.html") was to calculate a*b/n in normal precision giving a floating point approximation. By then calculating a*b - int(a*b/n) * n in double length using the floating point approximation of a*b/n then this is max +/-n from a*b mod n so the result can be checked and n added/subtracted accordingly. Voila! Regards, Erik Ehrling (Sweden)

Re: Miller-Rabin Primality Test for HP-42S Message #6 Posted by hugh on 15 Nov 2003, 1:12 p.m., in response to message #5 by Erik Ehrling (Sweden) i see. you are still stuck with doing a*b - q*n in double prevision though. unless there's some way to get the low bits of a*b from a calculator.

[ Return to Index | Top of Index ]

Go back to the main exhibit hall