Re: Prime factorization in the 48GX
Message #9 Posted by Joe Horn on 16 July 2008, 12:53 p.m.,
in response to message #8 by Gerson W. Barbosa
Mark Adler's FACNUM (on Goodies Disk #2) is twice as fast as my HP48 factorizer. Its listing follows, with Mark's comments.
%%HP: T(3)A(R)F(.);
@ by Mark Adler
@ FACNUM (BYTES = 277.5, #74FFh)
@ Given an integer, return its prime factorization as a list.
@ For example, 16353 FACNUM returns { 3 3 23 79 }.
\<<
@ initialize factor list
{ } SWAP
@ For this part of the program the stack is: n factorlist, where the two are
@ kept so that the product of the list times n is the original integer.
@ factor out 2's
WHILE DUP 2 MOD NOT REPEAT
2 / SWAP 2 + SWAP END
@ factor out 3's
WHILE DUP 3 MOD NOT REPEAT
3 / SWAP 3 + SWAP END
@ start factor search at 5
5
@ The stack is now: k n factorlist, where the second two are maintained as
@ before, and k is the largest 5 mod 6 integer less than or equal to the
@ last factor found (execpt initially when it is set to 5).
@ search from k to sqrt(n) for factors---k must be 5 mod 6
WHILE OVER 1 \=/ REPEAT @ do while n is not one
OVER \v/ FLOOR @ go up to the floor of the square root
IFERR @ (divide by zero used as a loop breaker)
FOR i @ look at factors that are 1 and 5 mod 6
IF DUP i MOD NOT THEN
i 0 / END @ if 5 mod 6 divides n, then cause error
IF DUP i 2 + MOD NOT THEN
i 2 + 0 / END @ if 1 mod 6 divides n, then cause error
6 STEP
THEN DROP @ got an error---trash the zero
ELSE DUP @ end of loop---n divides n
END @ after this, stack is: factor n factorlist
ROT OVER + ROT ROT @ add factor to list (factor n factorlist')
SWAP OVER / SWAP @ divide out divisor (factor n' factorlist')
DUP 6 MOD 5 - 2 / + @ set k' to largest 5 mod 6 <= divisor
END @ find next factor (k' n' factorlist')
@ return list, dropping n and k
DROP2
\>>
If you want REALLY fast factorizing on the HP48, you'll have to go beyond User RPL. The best one ever written (I think) is Klaus Kalb's "FCTR" on Goodies Disk #8, in the MATH directory. It is very fast, though not as fast as FACTOR in the HP 50g, of course.
-Joe-
| 
Go back to the main exhibit hall