diff --git a/changelog b/changelog index b219341..ec2debe 100644 --- a/changelog +++ b/changelog @@ -1,3 +1,4 @@ +20080312 tpd src/algebra/intfact.spad speed BasicSieve, prime, add docs 20080305 tpd src/hyper/bookvol11 add additional hyperdoc page translations 20080304 tpd src/hyper/bookvol11 add additional hyperdoc page translations 20080303 tpd src/hyper/bookvol11 add additional hyperdoc page translations diff --git a/src/algebra/intfact.spad.pamphlet b/src/algebra/intfact.spad.pamphlet index 996b923..fd43be1 100644 --- a/src/algebra/intfact.spad.pamphlet +++ b/src/algebra/intfact.spad.pamphlet @@ -10,6 +10,7 @@ \tableofcontents \eject \section{package PRIMES IntegerPrimesPackage} +We've expanded the list of small primes to include those between 1 and 10000. <>= )abbrev package PRIMES IntegerPrimesPackage ++ Author: Michael Monagan @@ -59,19 +60,221 @@ IntegerPrimesPackage(I:IntegerNumberSystem): with ++ \spad{primes(a,b)} returns a list of all primes p with ++ \spad{a <= p <= b} == add - smallPrimes: List I := [2::I,3::I,5::I,7::I,11::I,13::I,17::I,19::I,_ - 23::I,29::I,31::I,37::I,41::I,43::I,47::I,_ - 53::I,59::I,61::I,67::I,71::I,73::I,79::I,_ - 83::I,89::I,97::I,101::I,103::I,107::I,109::I,_ - 113::I,127::I,131::I,137::I,139::I,149::I,151::I,_ - 157::I,163::I,167::I,173::I,179::I,181::I,191::I,_ - 193::I,197::I,199::I,211::I,223::I,227::I,229::I,_ - 233::I,239::I,241::I,251::I,257::I,263::I,269::I,_ - 271::I,277::I,281::I,283::I,293::I,307::I,311::I,_ - 313::I] +@ +\subsection{smallPrimes} +This is a table of all of the primes in [2..10000]. It is used by the +prime? function to check for primality. It is used by the primes function +to generate arrays of primes in a given range. Changing the range included +in this table implies changing the value of the nextSmallPrime variable. +There is a constant in the function squareFree from IntegerFactorizationPackage +that is the square of the upper bound of the table range, in this case +10000000. +<>= + smallPrimes: List I := + [2::I, 3::I, 5::I, 7::I, 11::I, 13::I, 17::I, 19::I,_ + 23::I, 29::I, 31::I, 37::I, 41::I, 43::I, 47::I, 53::I,_ + 59::I, 61::I, 67::I, 71::I, 73::I, 79::I, 83::I, 89::I,_ + 97::I, 101::I, 103::I, 107::I, 109::I, 113::I, 127::I,_ + 131::I, 137::I, 139::I, 149::I, 151::I, 157::I, 163::I,_ + 167::I, 173::I, 179::I, 181::I, 191::I, 193::I, 197::I,_ + 199::I, 211::I, 223::I, 227::I, 229::I, 233::I, 239::I,_ + 241::I, 251::I, 257::I, 263::I, 269::I, 271::I, 277::I,_ + 281::I, 283::I, 293::I, 307::I, 311::I, 313::I, 317::I,_ + 331::I, 337::I, 347::I, 349::I, 353::I, 359::I, 367::I,_ + 373::I, 379::I, 383::I, 389::I, 397::I, 401::I, 409::I,_ + 419::I, 421::I, 431::I, 433::I, 439::I, 443::I, 449::I,_ + 457::I, 461::I, 463::I, 467::I, 479::I, 487::I, 491::I,_ + 499::I, 503::I, 509::I, 521::I, 523::I, 541::I, 547::I,_ + 557::I, 563::I, 569::I, 571::I, 577::I, 587::I, 593::I,_ + 599::I, 601::I, 607::I, 613::I, 617::I, 619::I, 631::I,_ + 641::I, 643::I, 647::I, 653::I, 659::I, 661::I, 673::I,_ + 677::I, 683::I, 691::I, 701::I, 709::I, 719::I, 727::I,_ + 733::I, 739::I, 743::I, 751::I, 757::I, 761::I, 769::I,_ + 773::I, 787::I, 797::I, 809::I, 811::I, 821::I, 823::I,_ + 827::I, 829::I, 839::I, 853::I, 857::I, 859::I, 863::I,_ + 877::I, 881::I, 883::I, 887::I, 907::I, 911::I, 919::I,_ + 929::I, 937::I, 941::I, 947::I, 953::I, 967::I, 971::I,_ + 977::I, 983::I, 991::I, 997::I, 1009::I, 1013::I,_ + 1019::I, 1021::I, 1031::I, 1033::I, 1039::I, 1049::I,_ + 1051::I, 1061::I, 1063::I, 1069::I, 1087::I, 1091::I,_ + 1093::I, 1097::I, 1103::I, 1109::I, 1117::I, 1123::I,_ + 1129::I, 1151::I, 1153::I, 1163::I, 1171::I, 1181::I,_ + 1187::I, 1193::I, 1201::I, 1213::I, 1217::I, 1223::I,_ + 1229::I, 1231::I, 1237::I, 1249::I, 1259::I, 1277::I,_ + 1279::I, 1283::I, 1289::I, 1291::I, 1297::I, 1301::I,_ + 1303::I, 1307::I, 1319::I, 1321::I, 1327::I, 1361::I,_ + 1367::I, 1373::I, 1381::I, 1399::I, 1409::I, 1423::I,_ + 1427::I, 1429::I, 1433::I, 1439::I, 1447::I, 1451::I,_ + 1453::I, 1459::I, 1471::I, 1481::I, 1483::I, 1487::I,_ + 1489::I, 1493::I, 1499::I, 1511::I, 1523::I, 1531::I,_ + 1543::I, 1549::I, 1553::I, 1559::I, 1567::I, 1571::I,_ + 1579::I, 1583::I, 1597::I, 1601::I, 1607::I, 1609::I,_ + 1613::I, 1619::I, 1621::I, 1627::I, 1637::I, 1657::I,_ + 1663::I, 1667::I, 1669::I, 1693::I, 1697::I, 1699::I,_ + 1709::I, 1721::I, 1723::I, 1733::I, 1741::I, 1747::I,_ + 1753::I, 1759::I, 1777::I, 1783::I, 1787::I, 1789::I,_ + 1801::I, 1811::I, 1823::I, 1831::I, 1847::I, 1861::I,_ + 1867::I, 1871::I, 1873::I, 1877::I, 1879::I, 1889::I,_ + 1901::I, 1907::I, 1913::I, 1931::I, 1933::I, 1949::I,_ + 1951::I, 1973::I, 1979::I, 1987::I, 1993::I, 1997::I,_ + 1999::I, 2003::I, 2011::I, 2017::I, 2027::I, 2029::I,_ + 2039::I, 2053::I, 2063::I, 2069::I, 2081::I, 2083::I,_ + 2087::I, 2089::I, 2099::I, 2111::I, 2113::I, 2129::I,_ + 2131::I, 2137::I, 2141::I, 2143::I, 2153::I, 2161::I,_ + 2179::I, 2203::I, 2207::I, 2213::I, 2221::I, 2237::I,_ + 2239::I, 2243::I, 2251::I, 2267::I, 2269::I, 2273::I,_ + 2281::I, 2287::I, 2293::I, 2297::I, 2309::I, 2311::I,_ + 2333::I, 2339::I, 2341::I, 2347::I, 2351::I, 2357::I,_ + 2371::I, 2377::I, 2381::I, 2383::I, 2389::I, 2393::I,_ + 2399::I, 2411::I, 2417::I, 2423::I, 2437::I, 2441::I,_ + 2447::I, 2459::I, 2467::I, 2473::I, 2477::I, 2503::I,_ + 2521::I, 2531::I, 2539::I, 2543::I, 2549::I, 2551::I,_ + 2557::I, 2579::I, 2591::I, 2593::I, 2609::I, 2617::I,_ + 2621::I, 2633::I, 2647::I, 2657::I, 2659::I, 2663::I,_ + 2671::I, 2677::I, 2683::I, 2687::I, 2689::I, 2693::I,_ + 2699::I, 2707::I, 2711::I, 2713::I, 2719::I, 2729::I,_ + 2731::I, 2741::I, 2749::I, 2753::I, 2767::I, 2777::I,_ + 2789::I, 2791::I, 2797::I, 2801::I, 2803::I, 2819::I,_ + 2833::I, 2837::I, 2843::I, 2851::I, 2857::I, 2861::I,_ + 2879::I, 2887::I, 2897::I, 2903::I, 2909::I, 2917::I,_ + 2927::I, 2939::I, 2953::I, 2957::I, 2963::I, 2969::I,_ + 2971::I, 2999::I, 3001::I, 3011::I, 3019::I, 3023::I,_ + 3037::I, 3041::I, 3049::I, 3061::I, 3067::I, 3079::I,_ + 3083::I, 3089::I, 3109::I, 3119::I, 3121::I, 3137::I,_ + 3163::I, 3167::I, 3169::I, 3181::I, 3187::I, 3191::I,_ + 3203::I, 3209::I, 3217::I, 3221::I, 3229::I, 3251::I,_ + 3253::I, 3257::I, 3259::I, 3271::I, 3299::I, 3301::I,_ + 3307::I, 3313::I, 3319::I, 3323::I, 3329::I, 3331::I,_ + 3343::I, 3347::I, 3359::I, 3361::I, 3371::I, 3373::I,_ + 3389::I, 3391::I, 3407::I, 3413::I, 3433::I, 3449::I,_ + 3457::I, 3461::I, 3463::I, 3467::I, 3469::I, 3491::I,_ + 3499::I, 3511::I, 3517::I, 3527::I, 3529::I, 3533::I,_ + 3539::I, 3541::I, 3547::I, 3557::I, 3559::I, 3571::I,_ + 3581::I, 3583::I, 3593::I, 3607::I, 3613::I, 3617::I,_ + 3623::I, 3631::I, 3637::I, 3643::I, 3659::I, 3671::I,_ + 3673::I, 3677::I, 3691::I, 3697::I, 3701::I, 3709::I,_ + 3719::I, 3727::I, 3733::I, 3739::I, 3761::I, 3767::I,_ + 3769::I, 3779::I, 3793::I, 3797::I, 3803::I, 3821::I,_ + 3823::I, 3833::I, 3847::I, 3851::I, 3853::I, 3863::I,_ + 3877::I, 3881::I, 3889::I, 3907::I, 3911::I, 3917::I,_ + 3919::I, 3923::I, 3929::I, 3931::I, 3943::I, 3947::I,_ + 3967::I, 3989::I, 4001::I, 4003::I, 4007::I, 4013::I,_ + 4019::I, 4021::I, 4027::I, 4049::I, 4051::I, 4057::I,_ + 4073::I, 4079::I, 4091::I, 4093::I, 4099::I, 4111::I,_ + 4127::I, 4129::I, 4133::I, 4139::I, 4153::I, 4157::I,_ + 4159::I, 4177::I, 4201::I, 4211::I, 4217::I, 4219::I,_ + 4229::I, 4231::I, 4241::I, 4243::I, 4253::I, 4259::I,_ + 4261::I, 4271::I, 4273::I, 4283::I, 4289::I, 4297::I,_ + 4327::I, 4337::I, 4339::I, 4349::I, 4357::I, 4363::I,_ + 4373::I, 4391::I, 4397::I, 4409::I, 4421::I, 4423::I,_ + 4441::I, 4447::I, 4451::I, 4457::I, 4463::I, 4481::I,_ + 4483::I, 4493::I, 4507::I, 4513::I, 4517::I, 4519::I,_ + 4523::I, 4547::I, 4549::I, 4561::I, 4567::I, 4583::I,_ + 4591::I, 4597::I, 4603::I, 4621::I, 4637::I, 4639::I,_ + 4643::I, 4649::I, 4651::I, 4657::I, 4663::I, 4673::I,_ + 4679::I, 4691::I, 4703::I, 4721::I, 4723::I, 4729::I,_ + 4733::I, 4751::I, 4759::I, 4783::I, 4787::I, 4789::I,_ + 4793::I, 4799::I, 4801::I, 4813::I, 4817::I, 4831::I,_ + 4861::I, 4871::I, 4877::I, 4889::I, 4903::I, 4909::I,_ + 4919::I, 4931::I, 4933::I, 4937::I, 4943::I, 4951::I,_ + 4957::I, 4967::I, 4969::I, 4973::I, 4987::I, 4993::I,_ + 4999::I, 5003::I, 5009::I, 5011::I, 5021::I, 5023::I,_ + 5039::I, 5051::I, 5059::I, 5077::I, 5081::I, 5087::I,_ + 5099::I, 5101::I, 5107::I, 5113::I, 5119::I, 5147::I,_ + 5153::I, 5167::I, 5171::I, 5179::I, 5189::I, 5197::I,_ + 5209::I, 5227::I, 5231::I, 5233::I, 5237::I, 5261::I,_ + 5273::I, 5279::I, 5281::I, 5297::I, 5303::I, 5309::I,_ + 5323::I, 5333::I, 5347::I, 5351::I, 5381::I, 5387::I,_ + 5393::I, 5399::I, 5407::I, 5413::I, 5417::I, 5419::I,_ + 5431::I, 5437::I, 5441::I, 5443::I, 5449::I, 5471::I,_ + 5477::I, 5479::I, 5483::I, 5501::I, 5503::I, 5507::I,_ + 5519::I, 5521::I, 5527::I, 5531::I, 5557::I, 5563::I,_ + 5569::I, 5573::I, 5581::I, 5591::I, 5623::I, 5639::I,_ + 5641::I, 5647::I, 5651::I, 5653::I, 5657::I, 5659::I,_ + 5669::I, 5683::I, 5689::I, 5693::I, 5701::I, 5711::I,_ + 5717::I, 5737::I, 5741::I, 5743::I, 5749::I, 5779::I,_ + 5783::I, 5791::I, 5801::I, 5807::I, 5813::I, 5821::I,_ + 5827::I, 5839::I, 5843::I, 5849::I, 5851::I, 5857::I,_ + 5861::I, 5867::I, 5869::I, 5879::I, 5881::I, 5897::I,_ + 5903::I, 5923::I, 5927::I, 5939::I, 5953::I, 5981::I,_ + 5987::I, 6007::I, 6011::I, 6029::I, 6037::I, 6043::I,_ + 6047::I, 6053::I, 6067::I, 6073::I, 6079::I, 6089::I,_ + 6091::I, 6101::I, 6113::I, 6121::I, 6131::I, 6133::I,_ + 6143::I, 6151::I, 6163::I, 6173::I, 6197::I, 6199::I,_ + 6203::I, 6211::I, 6217::I, 6221::I, 6229::I, 6247::I,_ + 6257::I, 6263::I, 6269::I, 6271::I, 6277::I, 6287::I,_ + 6299::I, 6301::I, 6311::I, 6317::I, 6323::I, 6329::I,_ + 6337::I, 6343::I, 6353::I, 6359::I, 6361::I, 6367::I,_ + 6373::I, 6379::I, 6389::I, 6397::I, 6421::I, 6427::I,_ + 6449::I, 6451::I, 6469::I, 6473::I, 6481::I, 6491::I,_ + 6521::I, 6529::I, 6547::I, 6551::I, 6553::I, 6563::I,_ + 6569::I, 6571::I, 6577::I, 6581::I, 6599::I, 6607::I,_ + 6619::I, 6637::I, 6653::I, 6659::I, 6661::I, 6673::I,_ + 6679::I, 6689::I, 6691::I, 6701::I, 6703::I, 6709::I,_ + 6719::I, 6733::I, 6737::I, 6761::I, 6763::I, 6779::I,_ + 6781::I, 6791::I, 6793::I, 6803::I, 6823::I, 6827::I,_ + 6829::I, 6833::I, 6841::I, 6857::I, 6863::I, 6869::I,_ + 6871::I, 6883::I, 6899::I, 6907::I, 6911::I, 6917::I,_ + 6947::I, 6949::I, 6959::I, 6961::I, 6967::I, 6971::I,_ + 6977::I, 6983::I, 6991::I, 6997::I, 7001::I, 7013::I,_ + 7019::I, 7027::I, 7039::I, 7043::I, 7057::I, 7069::I,_ + 7079::I, 7103::I, 7109::I, 7121::I, 7127::I, 7129::I,_ + 7151::I, 7159::I, 7177::I, 7187::I, 7193::I, 7207::I,_ + 7211::I, 7213::I, 7219::I, 7229::I, 7237::I, 7243::I,_ + 7247::I, 7253::I, 7283::I, 7297::I, 7307::I, 7309::I,_ + 7321::I, 7331::I, 7333::I, 7349::I, 7351::I, 7369::I,_ + 7393::I, 7411::I, 7417::I, 7433::I, 7451::I, 7457::I,_ + 7459::I, 7477::I, 7481::I, 7487::I, 7489::I, 7499::I,_ + 7507::I, 7517::I, 7523::I, 7529::I, 7537::I, 7541::I,_ + 7547::I, 7549::I, 7559::I, 7561::I, 7573::I, 7577::I,_ + 7583::I, 7589::I, 7591::I, 7603::I, 7607::I, 7621::I,_ + 7639::I, 7643::I, 7649::I, 7669::I, 7673::I, 7681::I,_ + 7687::I, 7691::I, 7699::I, 7703::I, 7717::I, 7723::I,_ + 7727::I, 7741::I, 7753::I, 7757::I, 7759::I, 7789::I,_ + 7793::I, 7817::I, 7823::I, 7829::I, 7841::I, 7853::I,_ + 7867::I, 7873::I, 7877::I, 7879::I, 7883::I, 7901::I,_ + 7907::I, 7919::I, 7927::I, 7933::I, 7937::I, 7949::I,_ + 7951::I, 7963::I, 7993::I, 8009::I, 8011::I, 8017::I,_ + 8039::I, 8053::I, 8059::I, 8069::I, 8081::I, 8087::I,_ + 8089::I, 8093::I, 8101::I, 8111::I, 8117::I, 8123::I,_ + 8147::I, 8161::I, 8167::I, 8171::I, 8179::I, 8191::I,_ + 8209::I, 8219::I, 8221::I, 8231::I, 8233::I, 8237::I,_ + 8243::I, 8263::I, 8269::I, 8273::I, 8287::I, 8291::I,_ + 8293::I, 8297::I, 8311::I, 8317::I, 8329::I, 8353::I,_ + 8363::I, 8369::I, 8377::I, 8387::I, 8389::I, 8419::I,_ + 8423::I, 8429::I, 8431::I, 8443::I, 8447::I, 8461::I,_ + 8467::I, 8501::I, 8513::I, 8521::I, 8527::I, 8537::I,_ + 8539::I, 8543::I, 8563::I, 8573::I, 8581::I, 8597::I,_ + 8599::I, 8609::I, 8623::I, 8627::I, 8629::I, 8641::I,_ + 8647::I, 8663::I, 8669::I, 8677::I, 8681::I, 8689::I,_ + 8693::I, 8699::I, 8707::I, 8713::I, 8719::I, 8731::I,_ + 8737::I, 8741::I, 8747::I, 8753::I, 8761::I, 8779::I,_ + 8783::I, 8803::I, 8807::I, 8819::I, 8821::I, 8831::I,_ + 8837::I, 8839::I, 8849::I, 8861::I, 8863::I, 8867::I,_ + 8887::I, 8893::I, 8923::I, 8929::I, 8933::I, 8941::I,_ + 8951::I, 8963::I, 8969::I, 8971::I, 8999::I, 9001::I,_ + 9007::I, 9011::I, 9013::I, 9029::I, 9041::I, 9043::I,_ + 9049::I, 9059::I, 9067::I, 9091::I, 9103::I, 9109::I,_ + 9127::I, 9133::I, 9137::I, 9151::I, 9157::I, 9161::I,_ + 9173::I, 9181::I, 9187::I, 9199::I, 9203::I, 9209::I,_ + 9221::I, 9227::I, 9239::I, 9241::I, 9257::I, 9277::I,_ + 9281::I, 9283::I, 9293::I, 9311::I, 9319::I, 9323::I,_ + 9337::I, 9341::I, 9343::I, 9349::I, 9371::I, 9377::I,_ + 9391::I, 9397::I, 9403::I, 9413::I, 9419::I, 9421::I,_ + 9431::I, 9433::I, 9437::I, 9439::I, 9461::I, 9463::I,_ + 9467::I, 9473::I, 9479::I, 9491::I, 9497::I, 9511::I,_ + 9521::I, 9533::I, 9539::I, 9547::I, 9551::I, 9587::I,_ + 9601::I, 9613::I, 9619::I, 9623::I, 9629::I, 9631::I,_ + 9643::I, 9649::I, 9661::I, 9677::I, 9679::I, 9689::I,_ + 9697::I, 9719::I, 9721::I, 9733::I, 9739::I, 9743::I,_ + 9749::I, 9767::I, 9769::I, 9781::I, 9787::I, 9791::I,_ + 9803::I, 9811::I, 9817::I, 9829::I, 9833::I, 9839::I,_ + 9851::I, 9857::I, 9859::I, 9871::I, 9883::I, 9887::I,_ + 9901::I, 9907::I, 9923::I, 9929::I, 9931::I, 9941::I,_ + 9949::I, 9967::I, 9973::I] productSmallPrimes := */smallPrimes - nextSmallPrime := 317::I + nextSmallPrime := 10007::I nextSmallPrimeSquared := nextSmallPrime**2 two := 2::I tenPowerTwenty:=(10::I)**20 @@ -81,8 +284,9 @@ IntegerPrimesPackage(I:IntegerNumberSystem): with 14386156093::I, 15579919981::I, 18459366157::I, 19887974881::I, 21276028621::I ]::(List I) PomeranceLimit:=27716349961::I -- replaces (25*10**9) due to Pinch - PinchList:= [3215031751::I, 118670087467::I, 128282461501::I, 354864744877::I, - 546348519181::I, 602248359169::I, 669094855201::I ] + PinchList:= _ + [3215031751::I, 118670087467::I, 128282461501::I, 354864744877::I, + 546348519181::I, 602248359169::I, 669094855201::I ] PinchLimit:= (10**12)::I PinchList2:= [2152302898747::I, 3474749660383::I] PinchLimit2:= (10**13)::I @@ -92,6 +296,9 @@ IntegerPrimesPackage(I:IntegerNumberSystem): with count2Order:Vector NonNegativeInteger := new(1,0) -- used to check whether we observe an element of maximal two-order +@ +\subsection{primes} +<>= primes(m, n) == -- computes primes from m to n inclusive using prime? l:List(I) := @@ -107,17 +314,18 @@ IntegerPrimesPackage(I:IntegerNumberSystem): with rabinProvesCompositeSmall : (I,I,I,I,NonNegativeInteger) -> Boolean +@ +\subsection{rabinProvesCompositeSmall} +<>= rabinProvesCompositeSmall(p,n,nm1,q,k) == -- probability n prime is > 3/4 for each iteration -- for most n this probability is much greater than 3/4 t := powmod(p, q, n) -- neither of these cases tells us anything --- if not (one? t or t = nm1) then if not ((t = 1) or t = nm1) then for j in 1..k-1 repeat oldt := t t := mulmod(t, t, n) --- one? t => return true (t = 1) => return true -- we have squared someting not -1 and got 1 t = nm1 => @@ -125,18 +333,19 @@ IntegerPrimesPackage(I:IntegerNumberSystem): with not (t = nm1) => return true false +@ +\subsection{rabinProvesComposite} +<>= rabinProvesComposite(p,n,nm1,q,k) == -- probability n prime is > 3/4 for each iteration -- for most n this probability is much greater than 3/4 t := powmod(p, q, n) -- neither of these cases tells us anything if t=nm1 then count2Order(1):=count2Order(1)+1 --- if not (one? t or t = nm1) then if not ((t = 1) or t = nm1) then for j in 1..k-1 repeat oldt := t t := mulmod(t, t, n) --- one? t => return true (t = 1) => return true -- we have squared someting not -1 and got 1 t = nm1 => @@ -147,10 +356,12 @@ IntegerPrimesPackage(I:IntegerNumberSystem): with # rootsMinus1 > 2 => true -- Z/nZ can't be a field false +@ +\subsection{prime?} +<>= prime? n == n < two => false n < nextSmallPrime => member?(n, smallPrimes) --- not one? gcd(n, productSmallPrimes) => false not (gcd(n, productSmallPrimes) = 1) => false n < nextSmallPrimeSquared => true @@ -202,6 +413,9 @@ IntegerPrimesPackage(I:IntegerNumberSystem): with rabinProvesComposite(currPrime,n,nm1,q,k) => return false true +@ +\subsection{nextPrime} +<>= nextPrime n == -- computes the first prime after n n < two => two @@ -209,6 +423,9 @@ IntegerPrimesPackage(I:IntegerNumberSystem): with while not prime? n repeat n := n + two n +@ +\subsection{prevPrime} +<>= prevPrime n == -- computes the first prime before n n < 3::I => error "no primes less than 2" @@ -270,13 +487,21 @@ IntegerRoots(I:IntegerNumberSystem): Exports == Implementation where 52::I,64::I,73::I,81::I,97::I,100::I,112::I,121::I] two := 2::I - +@ +\subsection{perfectSquare?} +<>= perfectSquare? a == (perfectSqrt a) case I + +@ +\subsection{perfectNthPower?} +<>= perfectNthPower?(b, n) == perfectNthRoot(b, n) case I +@ +\subsection{perfectNthRoot} +<>= perfectNthRoot n == -- complexity (log log n)**2 (log n)**2 m:NNI --- one? n or zero? n or n = -1 => [n, 1] (n = 1) or zero? n or n = -1 => [n, 1] e:NNI := 1 p:NNI := 2 @@ -287,16 +512,17 @@ IntegerRoots(I:IntegerNumberSystem): Exports == Implementation where p := convert(nextPrime(p::I))@Integer :: NNI [n, e] +@ +\subsection{approxNthRoot} +<>= approxNthRoot(a, n) == -- complexity (log log n) (log n)**2 zero? n => error "invalid arguments" --- one? n => a (n = 1) => a n=2 => approxSqrt a negative? a => odd? n => - approxNthRoot(-a, n) 0 zero? a => 0 --- one? a => 1 (a = 1) => 1 -- quick check for case of large n ((3*n) quo 2)::I >= (l := length a) => two @@ -311,15 +537,24 @@ IntegerRoots(I:IntegerNumberSystem): Exports == Implementation where z := x-y x +@ +\subsection{perfectNthRoot} +<>= perfectNthRoot(b, n) == (r := approxNthRoot(b, n)) ** n = b => r "failed" +@ +\subsection{perfectSqrt} +<>= perfectSqrt a == a < 0 or not member?(a rem (144::I), resMod144) => "failed" (s := approxSqrt a) * s = a => s "failed" +@ +\subsection{approxSqrt} +<>= approxSqrt a == a < 1 => 0 if (n := length a) > (100::I) then @@ -369,9 +604,12 @@ IntegerFactorizationPackage(I): Exports == Implementation where Implementation ==> add import IntegerRoots(I) - + BasicSieve: (I, I) -> FF +@ +\subsection{squareFree} +<>= squareFree(n:I):FF == u:I if n<0 then (m := -n; u := -1) @@ -381,18 +619,119 @@ IntegerFactorizationPackage(I): Exports == Implementation where rec.xpnt := 2 * rec.xpnt makeFR(u * unit sv, l) -- avoid using basic sieve when the lim is too big - lim := 1 + approxNthRoot(m,3) - lim > (100000::I) => makeFR(u, factorList factor m) + -- we know the sieve constants up to sqrt(100000000) + lim := 1 + approxSqrt(m) + lim > (100000000::I) => makeFR(u, factorList factor m) x := BasicSieve(m, lim) y := --- one?(m:= unit x) => factorList x ((m:= unit x) = 1) => factorList x (v := perfectSqrt m) case I => concat_!(factorList x, ["sqfr",v,2]$FFE) concat_!(factorList x, ["sqfr",m,1]$FFE) makeFR(u, y) - -- Pfun(y: I,n: I): I == (y**2 + 5) rem n +@ +\subsection{PollardSmallFactor} +This is Brent's\cite{1} optimization of Pollard's\cite{2} rho factoring. +Brent's algorithm is about 24 percent faster than Pollard's. Pollard;s +algorithm has complexity $O(p^{1/2})$ where $p$ is the smallest prime +factor of the composite number $N$. + +Pollard's idea is based on the observation that two numbers $x$ and $y$ +are congruent modulo $p$ with probability 0.5 after $1.177*\sqrt{p}$ numbers +have been randomly chosen. If we try to factor $n$ and $p$ is a factor of +$n$, then +$$1 < gcd(\vert x-y\vert,n) \le n$$ since $p$ divides both $\vert x-y\vert$ +and $n$. + +Given a function $f$ which generates a pseudo-random sequence of numbers +we allow $x$ to walk the sequence in order and $y$ to walk the sequence +at twice the rate. At each cycle we compute $gcd(\vert x-y\vert,n)$. +If this GCD ever equals $n$ then $x=y$ which means that we have walked +"all the way around the pseudo-random cycle" and we terminate with failure. + +This algorithm returns failure on all primes but also fails on some +composite numbers. + +Quoting Brent's back-tracking idea: +\begin{quote} +The best-known algorithm for finding GCDs is the Euclidean algorithm +which takes $O(\log N)$ times as long as one multiplication mod $N$. Pollard +showed that most of the GCD computations in Floyd's algorithm could be +dispensed with. ... The idea is simple: if $P_F$ computes $GCD(z_1,N)$, +$GCD(z_2,N)$,$\ldots$, then we compute +$$q_i=\prod_{j=1}^i{z_j}(\textrm{mod }N)$$ +and only compute $GCD(q_i,N)$ when $i$ is a multiple of $m$, where +$\log N < < m < < N^{1/4}$. Since $q_{i+1}=q_i \times z_{i+1}(\textrm{mod }N)$, +the work required for each GCD computation in algorithm $P_F$ is effectively +reduced to that for a multiplication mod $N$ in the modified algorithm. +The probability of the algorithm failing because $q_i=0$ increases, so it +is best not to choose $m$ too large. This problem can be minimized by +backtracking to the state after the previous GCD computation and setting +$m=1$. +\end{quote} +Brent incorporates back-tracking, omits the random choice of u, and +makes some minor modifications. His algorithm (p192-183) reads: + +\noindent +$y:=x_0; r:=1; q:=1;$ + +\noindent +\hbox{\hskip 0.5cm}{\bf repeat} $x:=y;$ + +\noindent +\hbox{\hskip 1.0cm}{\bf for} $i:=1$ {\bf to} $r$ {\bf do} $y:=f(y); k:=0;$ + +\noindent +\hbox{\hskip 1.0cm}{\bf repeat} $ys:=y;$ + +\noindent +\hbox{\hskip 1.5cm}{\bf for} $i:=1$ {\bf to} $min(m,r-k)$ {\bf do} + +\noindent +\hbox{\hskip 2.0cm}{\bf begin} $y:=f(y); q:=q*\vert x-y\vert mod N$ + +\noindent +\hbox{\hskip 2.0cm}{\bf end}; + +\noindent +\hbox{\hskip 1.5cm}$G:=GCD(q,N); k:=k+m$ + +\noindent +\hbox{\hskip 1.0cm}{\bf until} $(k \ge r)$ {\bf or} $(G > 1); r:=2*r$ + +\noindent +\hbox{\hskip 0.5cm}{\bf until} $G > 1$; + +\noindent +\hbox{\hskip 0.5cm}{\bf if} $G=N$ {\bf then} + +\noindent +\hbox{\hskip 1.0cm}{\bf repeat} $ys:=f(ys); G:=GCD(\vert y-yx\vert,N)$ + +\noindent +\hbox{\hskip 1.0cm}{\bf until} $G > 1$; + +\noindent +\hbox{\hskip 0.5cm}{\bf if} $G=N$ {\bf then} failure {\bf else} success + +Here we use the function +$$(y*y+5::I)~{\textrm rem}~ n$$ +as our pseudo-random sequence with a random starting value for y. + +On possible optimization to explore is to keep a hash table for the +computed values of the function $y_{i+1}:=f(y_i)$ since we effectively +walk the sequence several times. And we walk the sequence in a loop +many times. But because we are generating a very large number of +numbers the array can be a simple array of fixed size that captures +the last n values. So if we make a fixed array F of, say $2^q$ +elements we can store $f(y_i)$ in F[$y_i$ mod $2^q$]. + +One property that this algorithm assumes is that the function used +to generate the numbers has a long, hopefully complete, period. It +is not clear that the recommended function has that property. + +<>= PollardSmallFactor(n:I):Union(I,"failed") == -- Use the Brent variation x0 := random()$I @@ -405,7 +744,6 @@ IntegerFactorizationPackage(I): Exports == Implementation where x := y for i in 1..convert(r)@Integer repeat y := (y*y+5::I) rem n - q := (q*abs(x-y)) rem n k:I := 0 until (k>=r) or (G>1) repeat ys := y @@ -422,22 +760,55 @@ IntegerFactorizationPackage(I): Exports == Implementation where G=n => "failed" G - BasicSieve(r, lim) == - l:List(I) := - [1::I,2::I,2::I,4::I,2::I,4::I,2::I,4::I,6::I,2::I,6::I] - concat_!(l, rest(l, 3)) - d := 2::I - n := r +@ +\subsection{BasicSieve} +We create a list of prime numbers up to the limit given. The prior code +used a circular list but tests of that list show that on average more +than 50% of those numbers are not prime. Now we call primes to generate +the required prime numbers. Overall this is a small percentage of the +time needed to factor. + +This loop uses three pieces of information +\begin{enumerate} +\item n which is the number we are testing +\item d which is the current prime to test +\item lim which is the upper limit of the primes to test +\end{enumerate} + +We loop d over the list of primes. If the remaining number n is +smaller than the square of d then n must be prime and if it is +not one, we add it to the list of primes. If the remaining number +is larger than the square of d we remove all factors of d, reducing +n each time. Then we add a record of the new factor and its multiplicity, m. +We continue the loop until we run out of primes. + +Annoyingly enough, primes does not return an ordered list so we fix this. + +The sieve works up to a given limit, reducing out the factors that it +finds. If it can find all of the factors than it returns a factored +result where the first element is the unit 1. If there is still a +part of the number unfactored it returns the number and a list of +the factors found and their multiplicity. + +Basically we just loop thru the prime factors checking to see if +they are a component of the number, n. If so, we remove the factor from +the number n (possibly m times) and continue thru the list of primes. +<>= + BasicSieve(n, lim) == + p:=primes(1::I,lim::I)$IntegerPrimesPackage(I) + l:List(I) := append([first p],reverse rest p) ls := empty()$List(FFE) - for s in l repeat - d > lim => return makeFR(n, ls) + for d in l repeat if n1 then ls := concat_!(ls, ["prime",n,1]$FFE) return makeFR(1, ls) for m in 0.. while zero?(n rem d) repeat n := n quo d if m>0 then ls := concat_!(ls, ["prime",d,convert m]$FFE) - d := d+s + makeFR(n,ls) +@ +\subsection{BasicMethod} +<>= BasicMethod n == u:I if n<0 then (m := -n; u := -1) @@ -445,6 +816,38 @@ IntegerFactorizationPackage(I): Exports == Implementation where x := BasicSieve(m, 1 + approxSqrt m) makeFR(u, factorList x) +@ +\subsection{factor} +The factor function is many orders of magnitude slower than the results +of other systems. A posting on sci.math.symbolic showed that NTL could +factor the final value (t6) in about 11 seconds. Axiom takes about 8 hours. +\begin{verbatim} +a1:=101 +a2:=109 +t1:=a1*a2 +factor t1 + +a3:=21525175387 +t2:=t1*a3 +factor t2 + +a4:=218301576858349 +t3:=t2*a4 +factor t3 + +a5:=13731482973783137 +t4:=t3*a5 +factor t4 + +a6:=23326138687706820109 +t5:=t4*a6 +factor t5 + +a7:=4328240801173188438252813716944518369161 +t6:=t5*a7 +factor t6 +\end{verbatim} +<>= factor m == u:I zero? m => 0 @@ -452,7 +855,6 @@ IntegerFactorizationPackage(I): Exports == Implementation where else (n := m; u := 1) b := BasicSieve(n, 10000::I) flb := factorList b --- one?(n := unit b) => makeFR(u, flb) ((n := unit b) = 1) => makeFR(u, flb) a:LMI := dictionary() -- numbers yet to be factored b:LMI := dictionary() -- prime factors found @@ -465,7 +867,7 @@ IntegerFactorizationPackage(I): Exports == Implementation where (s := perfectNthRoot n).exponent > 1 => insert_!(s.base, a, c * s.exponent) -- test for a difference of square - x:=approxSqrt n; + x:=approxSqrt n if (x**2 insert_!(x+y,a,c) @@ -521,14 +923,16 @@ IntegerFactorizationPackage(I): Exports == Implementation where --SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. @ <<*>>= -<> - <> <> <> @ \eject \begin{thebibliography}{99} -\bibitem{1} nothing +\bibitem{1} Brent, Richard, ``An Improved Monte Carlo Factorization +Algorithm'', BIT 20, 1980, pp176-184, +http://web.comlab.ox.ac.uk/oucl/work/richard.brent/pd/rpb051i.pdf +\bibitem{2} Pollard, J.M., ``A Monte Carlo method for factorization'' +BIT Numerical Mathematics 15(3), 1975, pp331-334 \end{thebibliography} \end{document}