(* * LANGUAGE : ANS Forth * PROJECT : Forth Environments * DESCRIPTION : Factorize huge numbers, using GINT * CATEGORY : Infinite precision package * AUTHOR : Marcel Hendrix * LAST CHANGE : Friday, January 24, 2003 8:49 PM, Marcel Hendrix; No known bugs anymore. * LAST CHANGE : Sunday, January 19, 2003 7:55 PM, Marcel Hendrix; occasional access violation in .FACTOR * LAST CHANGE : January 12, 2003, Marcel Hendrix *) NEEDS -miscutil NEEDS -gint REVISION -factor "ÄÄÄ Factor with GINT Version 1.01 ÄÄÄ" PRIVATES DOC (* The sources of the full GiantInt package can be downloaded from: http://www.perfsci.com/free/giantint/giantint.html (Release 22 OCT 2001) "You may download giantint at any time, free of charge. PSI makes no warranty for this free software and we disclaim any and all implied warranties or conditions, including any implied warranty of title, of noninfringement, of merchantability, or of fitness for a particular purpose. If these terms are acceptable to you and you wish to download the free software, click here..." The sources can be compiled to gint.dll, with some extra work. What to do exactly will be come clear when I release the Linux variant. BUGS ---- S" 99999999999999999999999991111111111111111173212028727617181111" .factor 33107 * 560929 * 5384833477785736418660700040798137594423688327731437 [composite?] S" 5384833477785736418660700040798137594423688327731437" .factor ( should give: 548469301596654329 * 9817930487832035819402468332988053 or time out ) This times out. However, this could be simply because Lenstra's method is better. From: jerry@binkley.cs.mcgill.ca (Jerry Kuch) Newsgroups: sci.crypt,sci.math,sci.math.symbolic Subject: Re: Looking for an article (factoring large numbers) Date: 28 Jul 1993 18:00:24 -0400 Organization: SOCS - McGill University, Montreal, Canada Lines: 87 Message-ID: <236sto$amq@binkley.cs.mcgill.ca> References: NNTP-Posting-Host: binkley.cs.mcgill.ca Xref: msuinfo sci.crypt:18216 sci.math:48477 sci.math.symbolic:8362 In article freitag@elrond.toppoint.de (Claus Schoenleber) writes: >Someone told me about a method for factoring large numbers >using elliptic equations. I didn't find anything about that. >Could someone give me a pointer? I think you're talking about Lenstra's work. You can find it in: Lenstra, H.W. Jr. "Factoring Integers with elliptic curves." Annals of Mathematics, 126 (1987) 649-673. Essentially the same text exists in a University of Amsterdam tech report and also, I believe, in one from MSRI. Lenstra's algorithm is based on the Pollard p-1 method for integer factorization. In this scheme, given an integer, n, to be factored, you: 1. Choose an integer k such that k is a multiple of all (or most of the) integers less than some bound B where B can be as simple as lcm(all ints <= B) or B!. 2. Choose an integer, a from the set {2, n-2} 3. Compute a^k mod n efficiently (there are easy tricks for this that are better than multiplying a by itself k times) 4. Compute d = gcd(a^k - 1, n) using the Euclidean algorithm and the result from step 3. 5. If d isn't a nontrivial divisor of n then choose new values for a and/or k. Then start over. If p is an unknown prime factor of n, and p-1 has no large prime divisors then the above method is quite likely to find the factor, p. To get some idea of how it works, assume that k is divisible by all of the positive integers less than or equal to the bound B, above. Suppose as well that p is a prime divisor of n such that p-1 is the product of small prime powers all less than B. If so, then k must be a multiple of p-1 since it is a multiple of all the prime powers whose product is p-1. By Fermat's Little Theorem a^k is congruent to 1 modulo p. Because of this p divides gcd(a^k-1,n) and we only fail to produce nontrivial factors of n (in step 4 of the above algo) if a^k happens to be congruent to 1 mod n. Lenstra's method remedies the following defect. The above method suffers when all of n's prime divsiors, p are such that p-1 is divisible by a relatively large prime, or, by a relatively large prime power. The problem in essence is that we're relying on the group (Z/pZ)^* as p runs through all of n's prime divisors. For a given n, these groups aare all fixed, and if all of these groups have order divisible by a large prime they're next to useless. The elliptic curve method overcomes the problem by using elliptic curves over F_p = Z/pZ. These curves provide an easy source of different groups to choose from while performing the factorization, enough in fact that in practice we can typically expect to find one whose order is not divisible by a large prime or prime power. If we end up choosing a group over F_p that suffers from the problem described above, we toss it out by randomly choosing the group corresponding to another elliptic curve over F_p. Maple's integer factorization package provides the option of using Lenstra's elliptic curve method (as does Mathematica). It is interesting to note that it seems to used a pre-fixed sequence of elliptic curves and thus, the performance of the package on factoring a sequence of integers may vary extremely widely depending on what order they are input in. I don't know of a good way to intelligently feed Maple to avoid this situation; even knowing in advance what sequence Maple tries, it's still probably a big grind to test the suitability of a given group for a given problem (since it amounts to trying to solve the problem anyway). If you're really interested in integer factorization you might also want to look into a couple of other algorithms: - multiple polynomial quadratic sieve - the number field sieve The asymptotic performance of Lenstra's algo is of the same order as the conjectured time estimates for some of today's best factoring algorithms. It does have a couple of advantages though. One of concern to implementors is the fact that is has a pretty small storage requirement. Also, it is exceptionally fast when n is divisible by a prime much less than n^(1/2). Also if the factors of n are relatively small, it tends to beat its competitors whose runtimes tend to be moreorless independent of the size of individual factors. -- Jerry Kuch jerry@cs.mcgill.ca *) ENDDOC 1 =: FERMAT -1 =: MERSENNE 0 =: NOT-SPECIAL NOT-SPECIAL VALUE mode #11 VALUE /passes ( 11 ) -- mode: +1=Fermat, 0, -1=Mersenne -- FACTOR does not need to initialize anything. : (factor) ( a b mode -- ) /passes 4 GINT-LIB: factor FOREIGN @+ TYPE 1 ?GINT-WARNING ;P : .FACTOR ( c-addr u -- ) NOT-SPECIAL (factor) ; : .FERMAT ( u -- ) 0 SWAP FERMAT (factor) ; : .MERSENNE ( u -- ) 0 SWAP MERSENNE (factor) ; -- Below words must initialize everything ---------------------------------------------------------------------------------- : INIT-GIANT ( maxbits -- ) 1 GINT-LIB: giantinit FOREIGN DROP ; : 'GIANT ( -- addr ) 0 GINT-LIB: popg FOREIGN ;P : GIANT CREATE 'GIANT , ( compile: sel -- ) DOES> @ ; ( runtime: -- addr ) #2048 ( #240000 ) =: SIZE-IN-BITS SIZE-IN-BITS INIT-GIANT ( all giants will be SIZE-IN-BITS [256] bits maximum ) GIANT s0 GIANT s1 GIANT s2 GIANT s3 GIANT s4 GIANT s5 GIANT s6 GIANT s7 : >GIANT ( c-addr u giant -- ) 3 GINT-LIB: stringtog FOREIGN ?GINT-WARNING ; : GIANT> ( giant -- c-addr u ) 1 GINT-LIB: gianttos FOREIGN @+ ; : .GIANT ( giant -- ) GIANT> TYPE ; : #GBITS ( giant -- #bits ) 1 GINT-LIB: bitlen FOREIGN ; : GBIT? ( giant u -- 2^val ) 2 GINT-LIB: bitval FOREIGN ; : G->G ( g1 g2 -- ) 2 GINT-LIB: gtog VFOREIGN ; : GALLOC ( -- ganon ) 0 GINT-LIB: popg FOREIGN ; : GFREE ( u#giants -- ) 1 GINT-LIB: pushg VFOREIGN ; : GS@ ( giant -- n ) 1 GINT-LIB: gtoi FOREIGN ; : G1= ( giant -- true=1 ) 1 GINT-LIB: isone FOREIGN 0<> ; : G0= ( giant -- true=0 ) 1 GINT-LIB: isZero FOREIGN 0<> ; : GSIGN ( giant -- -1|0|1 ) 1 GINT-LIB: gsign FOREIGN ; : G0< ( giant -- true=0< ) GSIGN 0< ; : G0> ( giant -- true=0> ) GSIGN 0> ; : GCOMPARE ( g1 g2 -- ) 2 GINT-LIB: gcompg FOREIGN ; \ Returns 1, 0, -1 when a>b, a=b, a> ( giant u -- ) SWAP 2 GINT-LIB: gshiftright VFOREIGN ; : GGCD ( giant_x giant_n -- ) SWAP 2 GINT-LIB: gcdg VFOREIGN ; \ General GCD, x:= GCD(n,x). : Gm^n ( m n giant -- ) 3 GINT-LIB: powerg VFOREIGN ; \ g1 <- m^n : GS^MOD ( g1 n g2 -- ) 3 GINT-LIB: powermod VFOREIGN ; \ g1 <- g1^n MOD g2 : GG^MOD ( g1 g2 g3 -- ) 3 GINT-LIB: powermodg VFOREIGN ; \ g1 <- g1^g2 MOD g3 : !INVERSE ( g1 gr -- ) 2 GINT-LIB: make_recip VFOREIGN ; \ gr = 2^(2b)/gd, b = #gbits(d-1) : GG/_R ( g1 gr g2 -- ) SWAP ROT 3 GINT-LIB: divg_via_recip VFOREIGN ; \ g1 <- g1 / g2 using stored gr = 1/d : GGMOD_R ( g1 gr g2 -- ) SWAP ROT 3 GINT-LIB: modg_via_recip VFOREIGN ; \ g1 <- g1 % g2 using stored gr = 1/d : FERMAT^MOD ( g1 n q -- ) 3 GINT-LIB: fermatpowermod VFOREIGN ; \ g1 <- g1^n MOD 2^q+1 : FERMATG^MOD ( g1 gn q -- ) 3 GINT-LIB: fermatpowermodg VFOREIGN ; \ g1 <- g1^gn MOD 2^q+1 : FERMATMOD ( g1 n -- ) SWAP 2 GINT-LIB: fermatmod VFOREIGN ; \ g1 <- g1 MOD 2^n+1 : MERSENNE^MOD ( g1 n q -- ) 3 GINT-LIB: mersennepowermod VFOREIGN ; \ g1 <- g1^n MOD 2^q-1 : MERSENNEG^MOD ( g1 gn q -- ) 3 GINT-LIB: mersennepowermodg VFOREIGN ; \ g1 <- g1^gn MOD 2^q-1 : MERSENNEMOD ( g1 n -- ) SWAP 2 GINT-LIB: mersennemod VFOREIGN ; \ g1 <- g1 MOD 2^n-1 : #TRAILING0S ( g1 -- exponent ) DUP 1 GINT-LIB: numtrailzeros FOREIGN DUP >S G>> S> ; ( g1 = g1' * 2^exponent ) -- not limited to 16 bits : GS! ( u giant -- ) OVER #16 RSHIFT DUP IF OVER 2 GINT-LIB: itog VFOREIGN DUP #16 G<< SWAP $FFFF AND GS+ EXIT ENDIF DROP 2 GINT-LIB: itog VFOREIGN ; CREATE smtab PRIVATE 0 C, 0 C, 1 C, 1 C, 0 C, 1 C, 0 C, 1 C, -- Because bigprimeq needs argument > 5 : GPrime? ( giant -- true=PRIME ) DUP #GBITS 4 < IF GS@ smtab + C@ EXIT ENDIF 1 GINT-LIB: bigprimeq FOREIGN ; \ giant should be > 5 : .FAC ( u -- ) DUP 0< ABORT" FAC :: argument should be positive" GALLOC LOCAL tmp 1 tmp GS! DUP 2 < IF DROP ELSE 1+ 1 ?DO tmp I GS* LOOP ENDIF tmp .GIANT 1 GFREE ; -- Calculate the largest number whose square is smaller than g1, put it -- in g1. Uses Newton-Raphson : sqrt(v) = lim n->ì An+1=(An+v/An)/2 : GSQRT ( g1 -- ) LOCAL g1 g1 G0= ?EXIT GALLOC GALLOC GALLOC LOCALS| u1 u2 u3 | g1 u1 G->G g1 g1 #GBITS 1- 2/ G>> \ Use l.s. half of the bits u1 #GBITS \ bits per word = too much! 0 ?DO u1 u2 G->G \ restore v, original number g1 u3 G->G \ save An, previous approximation u2 g1 GG/ \ v/An g1 u2 GG+ g1 1 G>> \ ( v/An + An ) / 2 g1 u3 GG= ?LEAVE LOOP 3 GFREE ; : CHOOSE<>0 ( u1 -- u2 ) BEGIN DUP CHOOSE DUP 0= WHILE DROP REPEAT NIP ;P : RANDOM<>0 ( -- u2 ) BEGIN RANDOM $FFFF AND DUP 0= WHILE DROP REPEAT ;P -- Rabin? tells us whether g1 **might be** prime; FALSE: certainly not, TRUE: probably 0 VALUE save : RABIN? ( g1 -- flag ) GALLOC GALLOC GALLOC 0 LOCALS| two's d probe ss g1 | g1 ss G->G ss 1 GS- ss #TRAILING0S TO two's \ find power of two g1 #GBITS #16 < IF g1 GS@ CHOOSE<>0 ELSE RANDOM<>0 ENDIF probe GS! \ start PROBE with a random number < g1 probe ss g1 GG^MOD probe G1= IF 3 GFREE TRUE EXIT ENDIF \ probe == 1 => certainly prime g1 d !INVERSE BEGIN probe 1 GS+ probe g1 GG= IF 3 GFREE TRUE EXIT ENDIF \ probably prime two's WHILE probe 1 GS- \ restore probe probe G^2 probe d g1 GGMOD_R \ probe^2 mod g1; result in probe -1 +TO two's REPEAT 3 GFREE FALSE ;P \ not prime -- Doing RABIN? 30 times makes the result almost certain. : IsPrimeMR? ( g1 -- bool ) >S TRUE #30 0 DO S RABIN? 0= IF DROP FALSE LEAVE ENDIF LOOP -S ; -- S" 59649589127497217" s1 >giant -- : test_old CR TIMER-RESET 0 1000 0 DO s1 GPrime? + LOOP .ELAPSED SPACE . ; -- test_old 0.151 seconds elapsed. 1000 ok -- : test_new CR TIMER-RESET 0 1000 0 DO s1 IsPrimeMR? + LOOP .ELAPSED SPACE . ; -- test_new 4.886 seconds elapsed. -1000 ( probe G^2 probe d g1 GGMOD_R ) DOC (* Lucas-Lehmer ============ The Lucas-Lehmer test is Let p > 3 be prime. Define a(1) = 4 and a(n+1) = a(n)^2 - 2 for n >= 1. Then 2^p - 1 is prime if and only if 2^p - 1 divides a(p-1). The sequence begins 4, 14, 194, 37634, .. 14 is divisible by 7 and 37634 by 31, so 7 and 31 are prime. Proofs can observe a(n) = alpha^(2^n) + beta^(2^n) where alpha = sqrt(3/2) + sqrt(1/2) and beta = sqrt(3/2) - sqrt(1/2). The hypothesis ensures 2 is a quadratic residue modulo 2^p - 1 but 3 is a non-residue. The values of a(n) can be reduced modulo 2^p - 1. *) ENDDOC -- Test Mersenne(p) = 2^p - 1, p < 32 : .LL(31) ( p -- ) DUP 3 < ABORT" p should be >= 3" DUP #31 > ABORT" p should be < 32" LOCAL p 4 LOCAL A(n) p 2^x 1- LOCAL modulo p 2 ?DO A(n) DUP UM* 2. D- modulo UM/MOD DROP TO A(n) LOOP ." 2^" p 0 .R ." -1 is " A(n) ( A[p-1] ) modulo MOD IF ." not " ENDIF ." prime, check: " modulo s1 GS! s1 Gprime? 0= IF ." not " ENDIF ." prime." ; : TEST-Mn(31) ( -- ) #32 3 DO CR I .LL(31) LOOP ; -- Test Mersenne(p) = 2^p - 1 : .LL ( p -- ) DUP 3 < ABORT" p should be >= 3" DUP SIZE-IN-BITS > ABORT" p too large" ( squaring is allowed! ) LOCAL p GALLOC LOCAL A(n) GALLOC LOCAL modulo 4 A(n) GS! 1 modulo GS! modulo p G<< modulo 1 GS- p 2 ?DO A(n) G^2 A(n) 2 GS- A(n) modulo GGMOD LOOP ." 2^" p 0 .R ." -1 is " A(n) modulo GGMOD A(n) G0= 0= IF ." not " ENDIF ." prime, check: " modulo Gprime? 0= IF ." not " ENDIF ." prime." 2 GFREE ; : TEST-Mn ( uhigh ulow -- ) 3 UMAX ?DO CR I .LL LOOP ; -- Classical gcd : GCGCD ( v a -- ) OVER GABS OVER G0= IF 2DROP EXIT ENDIF GALLOC GALLOC LOCALS| u r a v | a u G->G u GABS BEGIN v G0= 0= WHILE u r G->G r v GGMOD v u G->G r v G->G REPEAT u v G->G 2 GFREE ; [DEFINED] testing [IF] 2 32 s1 Gm^n s1 cr .giant ( 4294967296 ) 2 64 s1 Gm^n s1 cr .giant ( 18446744073709551616 ) 2 128 s1 Gm^n s1 cr .giant ( 340282366920938463463374607431768211456 ) 2 256 s1 Gm^n s1 cr .giant ( 115792089237316195423570985008687907853269984665640564039457584007913129639936 ) 2 32 s2 Gm^n s2 cr .giant ( 4294967296 ) s1 s2 GGMOD s1 cr .giant ( 0 ) 2 256 s1 Gm^n 16 2^x 1- s2 GS! s1 s2 GGMOD s1 cr .giant ( 1 ) 2 256 s1 Gm^n s1 16 2^x 1- GS/MOD cr . ( 1 ) s1 cr .giant ( 1766874025136433896750911497805568136923323180981774075523881651146915841 ) s1 16 2^x 1- GS* s1 cr .giant ( 115792089237316195423570985008687907853269984665640564039457584007913129639935 ) 2 64 s1 Gm^n s1 1 GS+ s1 cr .giant ( 18446744073709551617 ) S" 18446744073709551617" s1 >GIANT S" 67280421310721" s2 >GIANT s1 s2 GG/ s1 cr .giant ( 274177 ) cr S" 2771" .factor ( 17 * 163 ) cr S" 9221" .factor ( 9221 ) cr S" 314159265" .factor ( 3 * 3 * 5 * 7 * 127 * 7853 ) cr S" 536870911" .factor ( 233 * 1103 * 2089 ) cr S" 1234567891" .factor ( 1234567891 ) cr S" 4225910033" .factor ( 65003 * 65011 ) cr S" 5704689200685129054721" s1 >giant s1 Gprime? . ( 1 ) cr S" 37866809061660057264219253397" s1 >giant s1 Gprime? . ( 1 ) cr S" 59649589127497217" s1 >giant s1 Gprime? . ( 1 ) cr S" 93461639715357977769163558199606896584051237541638188580280321" s1 >giant s1 Gprime? . ( 1 ) cr S" 124912135" s1 >giant s1 Gprime? . ( 0 ) cr S" 124912135" .factor ( 5 * 29 * 47 * 18329 ) CR 171 .FAC ( 1,241,018,070,217,667,823,424,840,524,103,103,992,616,605,577,501,693,185,388,951,803,611,996,075,22 1,691,752,992,751,978,120,487,585,576,464,959,501,670,387,052,809,889,858,690,710,767,331,242,032,21 8,484,364,310,473,577,889,968,548,278,290,754,541,561,964,852,153,468,318,044,293,239,598,173,696,89 9,657,235,903,947,616,152,278,558,180,061,176,365,108,428,800,000,000,000,000,000,000,000,000,000,00 0,000,000,000 ) S" 2000000000000000000000000000000" s4 >giant s4 GSQRT s4 cr .giant ( 1414213562373095 ) 10 3 TEST-Mn ( 2^3-1 is prime, check: prime. 2^4-1 is not prime, check: not prime. 2^5-1 is prime, check: prime. 2^6-1 is not prime, check: not prime. 2^7-1 is prime, check: prime. 2^8-1 is not prime, check: not prime. 2^9-1 is not prime, check: not prime. ) 129 120 TEST-Mn ( 2^120-1 is not prime, check: not prime. 2^121-1 is not prime, check: not prime. 2^122-1 is not prime, check: not prime. 2^123-1 is not prime, check: not prime. 2^124-1 is not prime, check: not prime. 2^125-1 is not prime, check: not prime. 2^126-1 is not prime, check: not prime. 2^127-1 is prime, check: prime. 2^128-1 is not prime, check: not prime. ) DOC (* TEST* on p54c ------------- 10 times 10x20 longmuls: 0.000 seconds elapsed. 10 times 100x200 longmuls: 0.001 seconds elapsed. 10 times 1000x2000 longmuls: 0.251 seconds elapsed. 10 times 2000x4000 longmuls: 0.444 seconds elapsed. 10 times 4000x8000 longmuls: 0.951 seconds elapsed. 10 times 8000x16000 longmuls: 2.256 seconds elapsed. TEST* on Athlon --------------- 10 times 10x20 longmuls: 0.000 seconds elapsed. 10 times 100x200 longmuls: 0.000 seconds elapsed. 10 times 1000x2000 longmuls: 0.023 seconds elapsed. 10 times 2000x4000 longmuls: 0.044 seconds elapsed. 10 times 4000x8000 longmuls: 0.093 seconds elapsed. 10 times 8000x16000 longmuls: 0.271 seconds elapsed. *) ENDDOC : TEST* ( msize -- ) 4 MAX DUP 8 * LOCALS| bitsize msize | bitsize 3 * SIZE-IN-BITS U> ABORT" test* :: increase SIZE-IN-BITS" 1 s1 GS! s1 bitsize G<< 1 s3 GS! s3 bitsize 2* G<< s3 s2 G->G CR ." 10 times " msize 0 .R &x EMIT msize 2* DEC. ." longmuls: " TIMER-RESET #10 0 ?DO s3 s2 G->G LOOP diff0 @ TIMER-PRESET #10 0 ?DO s3 s2 G->G s2 s1 GG* LOOP .ELAPSED ; [THEN] :ABOUT CR .~ Try: S" num" .FACTOR~ CR ." .FERMAT" CR ." .MERSENNE" ; .ABOUT -factor CR DEPRIVE (* End of Source *)