Programming Praxis

One version is shown in its entirety on the next page and uses trial division, Pollard’s rho algorithm, Pollard’s p-1 algorithm for a single stage only, and the two-stage modern elliptic–curve algorithm; we do not use wheel factorization, Fermat’s method, or the second stage of Pollard’s p-1 method.

Processing is controlled by the factors function, which we will look at in pieces; you can see the entire function on the next page. We begin with the internal definition of the factor? function. The rho, p-1 and elliptic curve functions either return a divisor or #f to indicate failure, and the divisor is not necessarily prime. Factor? handles the various return values. If the factorization fails, factor? returns #f to its caller. If the divisor is prime, factor? updates the current result (stored in the variables n and facts), exits if the remaining co-factor is prime, or returns #t so the factorization can continue. If the divisor is not prime, factor? calls factors recursively to factor the non-prime divisor, exits with failure if the non-prime divisor can’t be factored, exits with success if the factorization of the non-prime divisor completes the overall factorization, or returns #t so the overall factorization can continue.

The body of the factors function first checks if the input number is prime, then performs trial division, rho, p-1 and elliptic curve factorization according to the input parameters, and finally reports failure if the factorization is incomplete. The exit function, provided by call-with-current-continuation, exits as soon as the factorization is complete.

The key to factors is the list of parameters at the top of the function, which may be tuned as desired; the given list is only a starting point. We perform trial divisions to 100000, then five trials of Pollard’s rho method with 100000 iterations each (trials that find a factor don’t count against the limit of five), then Pollard’s p-1 to a first-stage limit of 500000, then two-stage elliptic curve factorization with up to 100 curves, with the bounds increasing for each curve. There is nothing particularly magical about these parameters, which were determined by watching the progress of multiple factorizations, sometimes trying alternate parameters to see the effect.

Pollard’s rho algorithm is the only one that changed from our original implementation. The original kept working until it found a factor or fell into a cycle, but our new implementation continues until it reaches an iteration limit.

We are lucky to have an example factorization that uses all five factorization methods — trial division, Pollard rho, Pollard p-1, elliptic curve first stage, and elliptic curve second stage — and even finds a non-prime factor that must be resolved recursively; it takes about twenty minutes on a recent-vintage personal computer to factor the 180-digit number consisting of the digit 1 repeated 180 times (verbose? was set to #t and the random-number seed was set to 87654321 immediately before this computation):

> (factors (/ (- (expt 10 180) 1) 9))
Trial division: bound=100000
  Found factors (3 3 7 11 13 19 31 37 41 61 101 181 211 241 271 2161 3541 9091 9901 27961 29611 52579)
  Remaining co-factor 29956560507825300590108042955447329957647469004935061390885594928146755845912192010432318854917916514410327624616947653196954237
Pollard rho: bound=100000, constant=18962011251229359009156890919192824170473724536426022746062529389468884409874447055866434987889860434844346440554063158257045492
  Found non-prime factor 79639973227
Trial division: bound=100000
Pollard rho: bound=100000, constant=35873039617
  Found factor 238681, remaining co-factor 333667
  Factorization complete
Pollard rho: bound=100000, constant=121263926659793047018967494638620147561553996947786243821468765881240485063608641306099219373338487279701814068269130
  Found factor 2906161, remaining co-factor 129431854036878665789788503606102934914555794443624857905233676070610306084009413816127915624631961962838855271
Pollard rho: bound=100000, constant=45293981637973688716048178266534217044761458238954448479313897289960291213667743367473724822037170518030739204
  Found factor 3762091, remaining co-factor 34404232655956133381619026122999931398404715474353187603711254212248004124304652337258167233230658684981
Pollard rho: bound=100000, constant=20517728786518624670946330699278063974297482243611702891096060311818021462529351004291224402194861041299
  Found factor 4188901, remaining co-factor 8213188293530005455277894159589813986629121928246379564403946097615580822823134835905209321784081
Pollard rho: bound=100000, constant=4560197729896101667659568704690844968251309845849091469566069791335054974684537648861925647476044
  Found factor 39526741, remaining co-factor 207788147612018037492084008635819836161780247155878081737220533754998440747319260039809741
Pollard rho: bound=100000, constant=176169430932959721732008089947128227769657917130383865383836754714219292359112710057268097
Pollard rho: bound=100000, constant=28291279336105116145302486655691410751803760689614970472004961213557588504503755195273851
Pollard rho: bound=100000, constant=24635869495669727118364795744405870824391004297581436443775082216192183581337777679578324
Pollard rho: bound=100000, constant=63899881778541156357940625601259544218421993173612956606370101816278815575887160524056613
Pollard rho: bound=100000, constant=165011410118268861689665393474767919953612916472803118126628391610433127415018909589055692
Pollard p-1: bound=500000
  Found factor 999999000001, remaining co-factor 207788355400165649302333145319662822679283263616462414936019007359069780809741
Pollard p-1: bound=500000
Elliptic curve 0: b1=1000, b2=50000, s=52439241705992280661387058759290444743990610173959283911521868385057817701317
Elliptic curve 1: b1=2000, b2=100000, s=23279137057405283816111948906736436185210031613947531580991006682668895057278
Elliptic curve 2: b1=3000, b2=150000, s=170650334841400385219838918873772923071036693721758574028579640868776517209654
Elliptic curve 3: b1=4000, b2=200000, s=154469146893428290169400461559393736575108289910998024730405156527921313203208
  In stage 2, found factor 8985695684401, remaining co-factor 23124348152684756076934883891491014602536740744674119535584829341
Elliptic curve 0: b1=1000, b2=50000, s=10724031799023089051485614020101272571865279149801267113289353397
Elliptic curve 1: b1=2000, b2=100000, s=21987089378435962857828548879935051743999728564922906702686651801
Elliptic curve 2: b1=3000, b2=150000, s=18779529741658528711391999864415567667144476253190937469340609872
Elliptic curve 3: b1=4000, b2=200000, s=625799164448954968839754305761837948037219419865316836890359848
Elliptic curve 4: b1=5000, b2=250000, s=9601320745076068107561857962702051154691572899832371366704240144
Elliptic curve 5: b1=6000, b2=300000, s=17713059863961516663524905085782656646223057263386564535063386410
Elliptic curve 6: b1=7000, b2=350000, s=564045493311964263014188284960340817355223965403678820931420862
Elliptic curve 7: b1=8000, b2=400000, s=5209572398908800472884639156441147206366332979213368755593510447
Elliptic curve 8: b1=9000, b2=450000, s=6307401793881717891393948563658230301312361022350588116976507873
Elliptic curve 9: b1=10000, b2=500000, s=8524110460679756182810787173637228353735134988639595955294912671
Elliptic curve 10: b1=11000, b2=550000, s=7611028230787730241892988661437446486650365910515305841415972452
Elliptic curve 11: b1=12000, b2=600000, s=8451781005093146397208812652830098074879002269679339743065284434
Elliptic curve 12: b1=13000, b2=650000, s=4401980878820585104498142243941323133807753933049912166929981401
Elliptic curve 13: b1=14000, b2=700000, s=15428311655116387373460352792307706205615019845103236794302196947
Elliptic curve 14: b1=15000, b2=750000, s=10733305756211278361058223902781933360833403226963411353464785884
Elliptic curve 15: b1=16000, b2=800000, s=1618526775041751795897969972151616161592612106350877053331465249
Elliptic curve 16: b1=17000, b2=850000, s=6038479190818225190402508992550492901950152933519090367463163435
Elliptic curve 17: b1=18000, b2=900000, s=10328547593983626097575695652363785246285886089801292101455335562
Elliptic curve 18: b1=19000, b2=950000, s=969530092681983488970960097823639285099162511725242630615028003
Elliptic curve 19: b1=20000, b2=1000000, s=18925079058060887009587570182124654975185812216920632984172068540
Elliptic curve 20: b1=21000, b2=1050000, s=18397527807253943118964978667521059809173201222005435160180753058
  In stage 1, found factor 4185502830133110721, remaining co-factor 5524867403314917121546955801099447513812160221
Elliptic curve 0: b1=1000, b2=50000, s=442355151974767168222445376029473660327237320
Elliptic curve 1: b1=2000, b2=100000, s=2206279433446644464853552159469250192349741531
Elliptic curve 2: b1=3000, b2=150000, s=2470411758512542972916861813911461170745896364
Elliptic curve 3: b1=4000, b2=200000, s=3152426015192178269083280351032804556110091009
Elliptic curve 4: b1=5000, b2=250000, s=3241421810107419967692728817194174483252688627
Elliptic curve 5: b1=6000, b2=300000, s=1331808138068966876196591255633098781998212710
Elliptic curve 6: b1=7000, b2=350000, s=5144443227054066442787727013859364843439638809
Elliptic curve 7: b1=8000, b2=400000, s=3249522815228253985639188745265464291959947644
Elliptic curve 8: b1=9000, b2=450000, s=3153533460463054869926598296859021979036971800
Elliptic curve 9: b1=10000, b2=500000, s=1621407799784796766200384939744535080488780543
Elliptic curve 10: b1=11000, b2=550000, s=4095419245851841507994433224204801483498665212
Elliptic curve 11: b1=12000, b2=600000, s=4045167994206024037786407710332902241905850355
Elliptic curve 12: b1=13000, b2=650000, s=5220455565588998999268694929566131028591130436
Elliptic curve 13: b1=14000, b2=700000, s=1485685910296159279101857529830913283433978911
Elliptic curve 14: b1=15000, b2=750000, s=809250966885343234362051280556440746715262276
Elliptic curve 15: b1=16000, b2=800000, s=807343000822671616249769698232061009672059568
Elliptic curve 16: b1=17000, b2=850000, s=1753554154384882060230527140268505636478347650
Elliptic curve 17: b1=18000, b2=900000, s=5342735723053689337091910951882225804429510497
Elliptic curve 18: b1=19000, b2=950000, s=4905400540969127434345203024220955942738137241
Elliptic curve 19: b1=20000, b2=1000000, s=1668116958645439977782361779537025192551810464
Elliptic curve 20: b1=21000, b2=1050000, s=894484358747944959126286325382761795292590044
Elliptic curve 21: b1=22000, b2=1100000, s=64856520117754637564184176885118286475090339
Elliptic curve 22: b1=23000, b2=1150000, s=1435951695790406616090404332137979816537009685
Elliptic curve 23: b1=24000, b2=1200000, s=3023084460639970139947507816841710893340471000
Elliptic curve 24: b1=25000, b2=1250000, s=5207140187890936637194397654617304189641192072
Elliptic curve 25: b1=26000, b2=1300000, s=4654249288339770969381418821317541178947644723
Elliptic curve 26: b1=27000, b2=1350000, s=1551363541784361022573400405379561003905533469
Elliptic curve 27: b1=28000, b2=1400000, s=2914203338229303561756297496111623956245252208
Elliptic curve 28: b1=29000, b2=1450000, s=4871356015837380049242176282335393946088671421
Elliptic curve 29: b1=30000, b2=1500000, s=2332005390602449057041184340592681493680091666
  In stage 2, found factor 4999437541453012143121, remaining co-factor 1105097795002994798105101
  Factorization complete
(3 3 7 11 13 19 31 37 41 61 101 181 211 241 271 2161 3541
9091 9901 27961 29611 52579 238681 333667 2906161 3762091
4188901 39526741 999999000001 8985695684401
4185502830133110721 4999437541453012143121
1105097795002994798105101)

The complete program is a “greatest hits” of our prime-number exercises. From the Standard Prelude we have expm, isqrt, ilog, digits, rand, randint, sort, and when. We save the primes to a billion (depending on your parameters, you may only need primes to five or ten million) for use by next-prime, which requires both the regular sieve of Eratosthenes and also the segmented sieve to compute and save the primes the first time it is run. We use the Baillie-Wagstaff primality checker. And we use the trial division, Pollard rho, Pollard p-1 and two-stage elliptic curve factorization algorithms. All of this code is collected on the next page and at http://programmingpraxis.codepad.org/wcKg7nEw.

Pages: 1 2 3