jurjen@cwi.nl (Jurjen NE Bos) (01/31/91)
This is a program that can factor large integers (up to 63 bits) on the
HP48. It is inspired by the postings of Klaus Kalb of last week. In fact, I
used two of Klaus' routines (see below).
The program has the following features:
- very fast (average 4 seconds for numbers of 12 digits,
1 minute for 18 digits.)
- uses two different algorithms for factoring (trial division and Pollard rho)
and proves all factors prime (using Selfridge's theorem).
- combines RPL (the language in the manual) with forth (what's in the ROM) and
machine code. This gives speed and servicability at the same time.
- The recursive calls are nicely shown during the computation.
- A separate prime testing routine.
There's also a bug:
- It doesn't quite work with wordsizes that are a multiple of four. This is
due to the carries. If you set the wordsize to 64 bits, it will be
modified to 63 to make it work.
How to use:
Download the object at the end of this posting. Use ASC-> to convert it into
a directory. Save it under the name FACTOR. Enter the directory. Press
DEMO. See what happens.
Remark: factoring numbers with large factors (more than 7 digits) can take
up to ten minutes.
Mind the bug above: if you set the wordsize to a multiple of four, factoring
might take forever.
Some interesting numbers to try:
3^37+1 (to get this, type 63 STWS #3 #37 #1 NEG POWMOD 1 +). The program
thinks it found a large prime factor, but quickly finds out it isn't.
100264053529. This number looks even more prime.
(Intermezzo:)
Some puzzles for experienced HP48 users:
(For information, send Email. Please realize these are PUZZLES, so
they aren't trivial, and explanation is complicated.)
- Compute 'SIN(180_\^o)' . You do not get 0. Why?
- The time to compute '253!' is much longer than the time to compute '253.1!'.
- To compute the factorial for large negative x, use '\pi*x/(SIN(\pi*x)*(-x)!)'
instead of 'x!' for faster results. (Don't forget RAD.)
Happy Hacking,
Jurjen Bos
-------------------Listing of the FACTOR directory (do not download)------
For the nosy people among you, here the listing of the routines:
(Don't bother downloading these, it's all in the ASC string at the end.)
DEMO @ Factor a random integer and time the computation
\<< DEC TICKS RAND 1000000 * R\->B RAND 1,E12 * + RAND 1,E18 * + FACTOR
TICKS ROT - B\->R 1,220703125E-4_s * SWAP
\>>
PRIME? @Prime test
\<< # 0d +
CASE DUP # 2d <
THEN DROP 0
END DUP # 210d BGCD # 1d ==
THEN RCWS 63 MIN STWS 0 'L' STO LPRT 'L' PURGE
END DUP # 121d <
THEN B\->R { 2 3 5 7 } SWAP POS 0 \=/
END DROP 0
END
\>>
FACTOR @The factoring routine
\<< RCWS 63 MIN STWS # 0d + 0 'L' STO LFAC 'L' PURGE
\>>
MULMOD @Modulo multiplication of binaries
(In machine code.
a b n -> a*b mod n
)
POWMOD @Modulo powering of binaries
(In machine code.
a e n -> a^e mod n
)
BGCD @Greatest commond divisor of binaries
(In machine code. Written by Klaus Kalb.
)
TRIAL @Trial division routine
(In machine code and forth. Machine code part written by Klaus Kalb.
)
MILRAB @Miller-Rabin test
(In forth.
n base -> 1 if n is pseudoprime
)
LFAC @Internal factoring routine
\<< DUP 'L' INCR DISP "Trial division" 'L' INCR DISP # 2000d TRIAL
IF DUP # 1d ==
THEN DROP { }
ELSE 1 \->LIST @List of composite factors
END
WHILE DUP SIZE
REPEAT SWAP OVER 1 GET
CASE DUP # 4012009d <
THEN B\->R +
END DUP L 1 - DISP "Primetest" L DISP DUP LPRT
THEN
IF DUP # 1000000000000d <
THEN B\->R
END +
END "Pollard \Gr" L DISP
WHILE P\Gr DUP # 1d ==
REPEAT DROP
END 2 \->LIST ROT SWAP + SWAP
END SWAP 2 OVER SIZE SUB
END DROP 'L' "
" OVER DECR DISP 1 STO-
\>>
LPRT @Internal primtesting routine
\<<
CASE DUP # 3d MILRAB NOT
THEN 0
END DUP # 25000000000d >
THEN SRT DUP
END DUP # 2d MILRAB NOT
THEN 0
END DUP # 1373653d <
THEN 1
END DUP # 5d MILRAB NOT
THEN 0
END DUP # 25326001d <
THEN 1
END DUP # 7d MILRAB NOT
THEN 0
END DUP # 3215031751d \=/
END SWAP DROP
\>>
SRT @Test primality of large numbers with Selfridge test.
\<< DUP 1 - LFAC # 3d \-> n l a
\<< "Selfridge" L DISP 9 CF 1 1 l SIZE
FOR k
IF l k GET SWAP OVER \=/ 8 FC?C OR
THEN
CASE a n 3 PICK / n POWMOD # 1d DUP2 \=/
THEN SWAP 3 PICK # 0d + n POWMOD \=/ 8 + SF
END + 'a' STO+ RAND ,1 \>=
THEN
END n a MILRAB
THEN
END 9 SF
END
END 8 FS? 9 FS? -
STEP DROP 9 FC?C
\>>
\>>
P\Gr @Pollard rho factoring algorithm.
\<< \-> n
\<< RAND n B\->R * R\->B DUP NEWOB
WHILE n Code n BGCD DUP # 1d ==
@ Code is a machine code routine that does 50 rounds
REPEAT DROP
END ROT ROT DROP2 n
\>> OVER /
\>>
-------------------FACTOR directory in downloadable form------------------
This is the entire program. Download in your machine, use ASC->,
and save under the name FACTOR. You'll need about 7K free to download it,
it is about 2430 bytes.
%%HP: T(3)A(D)F(,);
"69A20FF77E1100000020057920D9D20E16321C432D6E2010E6E16329B1C1D6E2
010E6BB691EEDA1B969178BF1CB2A133032D6E2010E6CCD208E0008FD8F357C8
0AF2E610A156110815711091923A98118AF77380B16A977970B16108119A9779
60B1610911099250A9EB12A9711A7B4010A18013A13896CFA781011015011111
5111128DD69508F2D76014713416917414713517901A9082281F832D0A1A9945
0B10A1699150B1991FCDA9601D6E2010E684E20402474344478BF1E4A2060000
1279E1D50328DBF149632E0CF1E0CF13FBF1D6E2010E6EF53292CF150FA19363
2B21304B1003035254530D9D20E163278BF19C2A290DA184E2040C4641434E4A
206000031C432D6E2010E6D6E2010C6D6E201016E1632C2A20710003556C6662
79646765684E2010C4485A1173A25D2C19C2A29C2A2D6E2010C68B9C10A132D6
E2010B63CE22D6E2010C6D6E2010B66C7D1DBBF192CF1D9AE1C53A2025C1908E
1AFE22D9D20D8732D9D20D6E201016D6E2010E63F2A2A9CF150FA1D6E2010E68
4E206005F475D4F444E4A205100010000000000000002ABF1D9AE18A732D9D20
DBBF13F2A2A9CF1E4A2051000000000000000000076BA1D6E2010E684E206005
F475D4F444D9AE1C53A276BA1472C1B21305DF2276BA145632D6E20101697632
B44029B1C1339209990000000000010B9DE18A732D9D20B21305DF22D6E2010E
6D6E20101684E2060D494C42514248A732D9D20B21305DF22173A2472C1B2130
5DF22B21305DF22C53A2313C1173A2313C190DA1083328DBF1173A2025C1EF53
293632B21300C20040C405254540D9D20E1632D8732D9D2078BF1E4A20510003
00000000000000084E2060D494C4251424F88E18A7324B2A25DF2278BF1E4A20
5100000ABD12D50000000D5CE18A732D9D2084E203035254578BF1B21305DF22
78BF1E4A2051000200000000000000084E2060D494C4251424F88E18A7324B2A
25DF2278BF1E4A20510005D5F410000000000EBBE18A7329C2A25DF2278BF1E4
A2051000500000000000000084E2060D494C4251424F88E18A7324B2A25DF227
8BF1E4A20510001B17281000000000EBBE18A7329C2A25DF2278BF1E4A205100
0700000000000000084E2060D494C4251424F88E18A7324B2A25DF2278BF1E4A
20510007CD71AFB00000000D9AE1B21305DF22DBBF18DBF193632B2130A22004
0C464143440D9D20E163278BF14563284E2010C4976324F802485A1C2A201200
045279616C602469667963796F6E64563284E2010C4976324F802485A1E4A205
10000D7000000000000084E205045259414C43CE2278BF1E4A20510001000000
000000000279E1AFE22D9D208DBF147A20B2130B21305BF22D9D209C2A2387C1
B21305DF223303278BF18B9C1D5032D9D20DBBF192CF19C2A26C7D1D8732D9D2
078BF1E4A20510009E73D30000000000EBBE18A732D9D20BB69176BA1B21305D
F2278BF184E2010C49C2A290DA1485A1C2A2071000052796D6564756374784E2
010C4485A178BF184E2040C40525458A732D9D203CE2278BF1E4A20510000001
5A4D8E000000EBBE1AFE22BB6915DF2276BA1B21305DF22C2A207100005F6C6C
6162746027984E2010C4485A13303284E2020057978BF1E4A205100010000000
00000000279E1D50328DBF149632ED2A2387C1E0CF1DBBF176BA1DBBF1B21305
DF22DBBF1ED2A292CF18B9C1C58C1B2130496328DBF14563284E2010C497632C
2A2070000A092CF1AA902485A19C2A28350293632B2130F230060D494C425142
460D9D20FDE8111920BB000D9D202C230E4A206000010BE3588130CCD2042000
8FD8F3582281C832AFA14B148DD69505AC26A321684E206005F475D4F4448813
0E4A206000019D445FC7A239916D9D20E7F069C2A2B21302A17088130C12169D
445B67A2EF116A32169D445B67A264B30EE170D9D2088130A321684E2060D455
C4D4F44432230E5D3532230B21305E17044230CE445B9F06B2130B2130741005
045259414C450D9D20FDE8111920BB000D9D20CCD20B61008F77F3510110AAF2
10810B2081B58082444000C21366570242462642466264264684242486462462
664246264242A2A021224011BD2BF6BF6BF6BF6BF6BF61088F2D7608F735608F
B97601112F8DD6950AF015A097CB12081B580824C7FFFC21366EDF118C24A910
81181129F2C8111118AF3AF19F262A76B779FE7F81EA75A7F9F280B7AB7597F9
E11A9F180AF910A97CD3AF4101113132AF226301A7A103208F2D7608F735608F
B9760113130657F1606E3F47A20D6E2010B6D6E201027B21300D470D6E2010B6
6AC30A2170D9D20D6E2010B63C370FBD81D6E2010B65233043370B2130D6E201
0B695450D6E20102779470B2130B213012200402474344440D9D20FDE81ECD46
CCD20F70008F77F35AF397845AFE978C48087091AFE80870F081E81CB7765EF8
08B09081E65FF9F650AFEB7A9780181C808608F65EFAFA97BC0A74A7F64FF8F5
9E35B21304A0006005F475D4F44460D9D200FE8111920BBB00CCD20D80008FD8
F35101208F2D7608F77F35121A98A90B6482281E108832901197C10121A96721
012111891E8D8DD6950A97A96A9082281F832D0A1A99450B10A1699150B1991F
CD01B2130BB00060D455C4D4F44460D9D200FE8111920BBB00CCD20E50008FD8
F35101208F2D7608F77F35A97119A9182281F83201A1847099550B11A1447099
250B1A91F6DA948DD6950B2130C80006064143445F42560D9D20E1632EF5C133
92010000000000003603ECB15C5C1E4A2051000000000000000000076BA14B2A
24563284E2010C497632DCC0284E2040C46414344563284E2010C497632EFE02
93632B2130BA00060052594D454F360D9D20E1632E4A20510000000000000000
00076BA1D8732D9D2078BF1E4A20510002000000000000000EBBE18A732D9D20
8DBF14B2A2B21305DF2278BF1E4A20510002D0000000000000084E2040247434
44E4A20510001000000000000000279E18A732D9D20EF5C13392010000000000
003603ECB15C5C14B2A24563284E2010C497632DCC0284E2040C405254545632
84E2010C497632EFE02B21305DF2278BF1E4A20510009700000000000000EBBE
18A732D9D20BB69147A20ED2A23F2A2D13A2743A2B2130DBBF14BAC14B2A2D9A
E1B21305DF228DBF14B2A2B21305DF2293632B2130BD100404454D4F440D9D20
E1632475C1D28919B1C1339206000000000000010EEDA1B96919B1C133920210
0000000000010EEDA176BA19B1C1339208100000000000010EEDA176BA184E20
6064143445F425D2891E0CF190DA1BB691ADA20339206990052130702210C2A2
0700003768B01B2130EEDA1DBBF193632B21306F87"
--- end of posting ---