below is Pari/GP code which prints out a table for the optimal amount of trial division to do before calling ispseudoprime on large primes, automating some of previously discussed. the table is specific to the computer on which the code is run, specific to the relative speeds of trial division and modular exponentiation.
ispseudoprime internally does trial division only up to 101, so parithreshold is 102. to confirm this, we use Mersenne prime 2^19937-1. the first call to ispseudoprome returns immediately.
? #
timer = 1 (on)
? y=2^19937-1;
? ispseudoprime(101*y)
0
? ispseudoprime(103*y)
cpu time = 1,313 ms, real time = 1,316 ms.
0
we construct a composite which will pass ispseudoprime's internal trial division by multiplying large random primes (primes approximately 2^1000 = growthrate). (unfixed bug: asking for a random prime less than 2^1000 might, with miniscule probability, generate a prime less than 102.)
we do some fine adjustment of growthrate so that the composites are approximately 1000n+buffer bits (buffer = 5) to make the table pretty. calling randomprime with argument 2^(x+1) will, on average, provide a prime of size 2^x.
the variable "dummy" is to make sure the calls to mod and ispseudoprime do not get optimized away. this might be more cautious than necessary.
after its small amount of trial division, ispseudoprime does the BPSW primality test, which starts with a Miller-Rabin test with base 2. when given a composite, with high probability this first test fails, and we assume that this is always the case. if we are unlucky and hit a composite which has to do the rest of BPSW, then that line's entry in the second column will be too large by a factor of 3 or 4.
the output is probabilistic; do multiple runs to assure yourself of results.
the columns are (1) width in bits of the composite being tested, (2) milliseconds for a failing ispseudoprime test, (3) milliseconds for one trial division, (4) how far (with forprime) you should do trial division for numbers of that size. for example, if you are testing primality of 10000-bit numbers on the same computer as tested below, you should define a function like this:
mypseudoprime(n) = forprime(p=100, 301720, if(n%p==0,return(0))); ispseudoprime(n)
future work: fit output to a curve; interpolate and extrapolate.
the code:
? print("(bits) (miller-rabin ms) (trial division ms) (forprime limit)"); buffer=5; growthrate=1000; trialdivisioniterations=500000; parithreshold=102; endprime= prime(primepi(parithreshold)+trialdivisioniterations); millerrabiniterations=10; k=randomprime(2^(growthrate+buffer+1)); for(i=2,+oo,desired=growthrate*i+buffer; thisgrow=desired-round(log(k)/log(2)); k*=randomprime(2^(thisgrow+1)); dummy=0; gettime; for(j=1,millerrabiniterations,dummy+=ispseudoprime(k)); timemillerrabin=gettime/millerrabiniterations; forprime(p=parithreshold,endprime,if(k%p==0,dummy+=1)); timetrialdivision=gettime/trialdivisioniterations; printf("%.1f %.1f (%.2e) %.1f\n", log(k)/log(2), timemillerrabin, timetrialdivision, timemillerrabin/timetrialdivision))
example output:
(bits) (miller-rabin ms) (trial division ms) (forprime limit)
2004.2 4.1 (3.12 e-4) 13141.0
3005.2 13.2 (3.84 e-4) 34375.0
4004.5 25.1 (4.50 e-4) 55777.8
5005.4 44.7 (5.20 e-4) 85961.5
6005.5 73.2 (5.92 e-4) 123648.6
7005.7 110.6 (6.60 e-4) 167575.8
8005.7 149.9 (7.26 e-4) 206473.8
8995.9 199.4 (8.00 e-4) 249250.0
10005.2 263.1 (8.72 e-4) 301720.2
11003.1 335.0 (9.40 e-4) 356383.0
12000.7 420.4 (1.01 e-3) 417892.6
13005.2 519.3 (1.07 e-3) 484421.6
14003.4 612.3 (1.14 e-3) 535227.3
15005.9 724.0 (1.22 e-3) 595394.7
16005.6 861.0 (1.29 e-3) 668478.3
17004.4 1011.3 (1.35 e-3) 746898.1
18005.0 1173.8 (1.43 e-3) 818549.5
19003.5 1351.4 (1.49 e-3) 904551.5
20006.1 1556.2 (1.58 e-3) 986185.0
21005.0 1732.7 (1.63 e-3) 1064312.0
22005.3 1909.5 (1.70 e-3) 1123235.3
23005.8 2152.0 (1.77 e-3) 1213077.8
24002.6 2414.8 (1.84 e-3) 1310966.3
25004.9 2666.9 (1.91 e-3) 1399213.0
26005.4 2958.1 (1.98 e-3) 1495500.5
27006.0 3245.6 (2.04 e-3) 1589422.1
28006.0 3560.3 (2.12 e-3) 1677804.0
29005.8 3913.1 (2.19 e-3) 1790073.2
30004.0 4247.4 (2.25 e-3) 1886056.8
31004.0 4618.4 (2.32 e-3) 1987263.3
32004.6 5003.2 (2.39 e-3) 2091638.8
33000.0 5349.2 (2.46 e-3) 2170941.6
34003.2 5781.2 (2.55 e-3) 2265360.5
35004.9 6213.4 (2.60 e-3) 2386098.3
36003.3 6713.6 (2.67 e-3) 2512574.9
37002.1 7305.4 (2.75 e-3) 2656509.1
38005.1 8087.1 (2.81 e-3) 2875924.6
39003.6 8526.0 (2.88 e-3) 2960416.7
39997.8 8835.2 (3.00 e-3) 2943104.6
41005.1 9398.1 (3.02 e-3) 3114015.9
42005.8 9832.8 (3.09 e-3) 3184196.9
43003.1 10339.4 (3.16 e-3) 3267825.5
44005.6 11007.1 (3.23 e-3) 3409882.3
45003.5 11644.2 (3.29 e-3) 3537120.3
46004.8 12352.8 (3.37 e-3) 3667696.0
47005.6 13060.3 (3.44 e-3) 3798807.4
48005.5 13850.3 (3.50 e-3) 3952711.2
49003.3 14825.6 (3.58 e-3) 4145861.3
50005.9 15435.8 (3.65 e-3) 4231304.8
51005.6 16128.1 (3.71 e-3) 4347196.8
52004.8 16998.9 (3.78 e-3) 4497063.5
53005.0 17837.7 (3.85 e-3) 4637987.5
54004.5 18631.3 (3.92 e-3) 4752882.7
55004.8 19613.3 (4.00 e-3) 4903325.0
56004.1 20424.9 (4.06 e-3) 5025812.0
57005.9 21383.3 (4.13 e-3) 5177554.5
58004.2 22470.3 (4.19 e-3) 5357725.3
59003.3 23287.9 (4.27 e-3) 5458954.5
60006.3 24412.8 (4.34 e-3) 5627662.5
61005.9 25196.2 (4.40 e-3) 5726409.1
62003.8 25857.4 (4.48 e-3) 5766592.3
63004.8 26956.6 (4.54 e-3) 5932350.4
64005.4 27839.4 (4.61 e-3) 6038915.4
65005.5 29078.7 (4.70 e-3) 6192227.4
66005.4 29854.0 (4.76 e-3) 6266582.7
67003.5 30966.1 (4.82 e-3) 6424502.1
68006.4 32222.3 (4.90 e-3) 6570615.8
69005.0 33446.1 (4.97 e-3) 6732306.8
70005.6 34635.7 (5.05 e-3) 6861271.8
71005.0 35762.7 (5.10 e-3) 7012294.1
72005.1 36893.8 (5.17 e-3) 7130614.6
73003.9 38375.0 (5.23 e-3) 7334671.3
74005.5 39864.7 (5.32 e-3) 7496182.8
^C
*** at top-level: ...errabiniterations,dummy+=ispseudoprime(k));ti
*** ^--------------------
*** ispseudoprime: user interrupt after 2h, 29min, 5,695 ms
*** Break loop: <Return> to continue; 'break' to go back to GP prompt
break> break
No comments :
Post a Comment