[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Axiom-developer] 20080312.01.tpd.patch (BasicSieve, primes, intfact doc
|
From: |
daly |
|
Subject: |
[Axiom-developer] 20080312.01.tpd.patch (BasicSieve, primes, intfact documentation) |
|
Date: |
Wed, 12 Mar 2008 19:35:10 -0600 |
This patch updates the documentation for intfact.spad.
It expands the list of pre-computed primes to cover the range [2..10000].
It documents the PollardSmallFactor function and fixes it to conform
to Brent's published algorithm (an improvement on Pollard's).
It documents and optimizes BasicSieve. The original function used a
circular list to generate the prime sieve but testing showed that the
method generates primes less than half the time so most of the work
was wasted. This now uses the primes function which has been expanded
to improve the speed of small prime generation.
The factor function, which led to these changes, is also documented with
the original problem from sci.math.symbolic. It appears that Brent's
algorithm needs to be replaced for factoring larger primes. Axiom takes
about 10 hours to factor the given number. The algorithm under development
can do the factoring in about 7 seconds. Due to the complexity of the
code it is not included in this patch.
Tim
=========================================================================
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.
<<package PRIMES IntegerPrimesPackage>>=
)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.
+<<package PRIMES IntegerPrimesPackage>>=
+ 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}
+<<package PRIMES IntegerPrimesPackage>>=
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}
+<<package PRIMES IntegerPrimesPackage>>=
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}
+<<package PRIMES IntegerPrimesPackage>>=
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?}
+<<package PRIMES IntegerPrimesPackage>>=
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}
+<<package PRIMES IntegerPrimesPackage>>=
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}
+<<package PRIMES IntegerPrimesPackage>>=
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?}
+<<package IROOT IntegerRoots>>=
perfectSquare? a == (perfectSqrt a) case I
+
+@
+\subsection{perfectNthPower?}
+<<package IROOT IntegerRoots>>=
perfectNthPower?(b, n) == perfectNthRoot(b, n) case I
+@
+\subsection{perfectNthRoot}
+<<package IROOT IntegerRoots>>=
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}
+<<package IROOT IntegerRoots>>=
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}
+<<package IROOT IntegerRoots>>=
perfectNthRoot(b, n) ==
(r := approxNthRoot(b, n)) ** n = b => r
"failed"
+@
+\subsection{perfectSqrt}
+<<package IROOT IntegerRoots>>=
perfectSqrt a ==
a < 0 or not member?(a rem (144::I), resMod144) => "failed"
(s := approxSqrt a) * s = a => s
"failed"
+@
+\subsection{approxSqrt}
+<<package IROOT IntegerRoots>>=
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}
+<<package INTFACT IntegerFactorizationPackage>>=
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.
+
+<<package INTFACT IntegerFactorizationPackage>>=
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.
+<<package INTFACT IntegerFactorizationPackage>>=
+ 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 n<d*d then
if n>1 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}
+<<package INTFACT IntegerFactorizationPackage>>=
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}
+<<package INTFACT IntegerFactorizationPackage>>=
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<n) then x:=x+1
(y:=perfectSqrt (x**2-n)) case I =>
insert_!(x+y,a,c)
@@ -521,14 +923,16 @@ IntegerFactorizationPackage(I): Exports == Implementation
where
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
-<<license>>
-
<<package PRIMES IntegerPrimesPackage>>
<<package IROOT IntegerRoots>>
<<package INTFACT IntegerFactorizationPackage>>
@
\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}
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Axiom-developer] 20080312.01.tpd.patch (BasicSieve, primes, intfact documentation),
daly <=