Skip to content

Instantly share code, notes, and snippets.

@primo-ppcg
Last active July 23, 2023 19:28
Show Gist options
  • Save primo-ppcg/31a31ecede0296457b6c57fe1f0c8020 to your computer and use it in GitHub Desktop.
Save primo-ppcg/31a31ecede0296457b6c57fe1f0c8020 to your computer and use it in GitHub Desktop.
(use spork/test)
(use spork/misc)
(use spork/math)
(def- small-primes
"primes less than 212"
[2 3 5 7 11 13 17 19 23 29 31 37
41 43 47 53 59 61 67 71 73 79 83 89
97 101 103 107 109 113 127 131 137 139 149 151
157 163 167 173 179 181 191 193 197 199 211])
(def- soe-indices
"pre-calced sieve of eratosthenes for 2 3 5 7"
[1 11 13 17 19 23 29 31 37 41 43 47
53 59 61 67 71 73 79 83 89 97 101 103
107 109 113 121 127 131 137 139 143 149 151 157
163 167 169 173 179 181 187 191 193 197 199 209])
(def- soe-offsets
"distances between sieve values"
[10 2 4 2 4 6 2 6 4 2 4 6
6 2 6 4 2 6 4 6 8 4 2 4
2 4 8 6 4 6 2 4 6 2 6 6
4 2 4 6 2 6 4 2 4 2 10 2])
(def- soe-table
"tabulated offsets mod 105"
[0 10 2 0 4 0 0 0 8 0 0 2 0 4 0
0 6 2 0 4 0 0 4 6 0 0 6 0 0 2
0 6 2 0 4 0 0 4 6 0 0 2 0 4 2
0 6 6 0 0 0 0 6 6 0 0 0 0 4 2
0 6 2 0 4 0 0 4 6 0 0 2 0 6 2
0 6 0 0 4 0 0 4 6 0 0 2 0 4 8
0 0 2 0 10 0 0 4 0 0 0 2 0 4 2])
(def- rm-bases-32
"bases for rabin-miller determinism to 2^32"
[2 7 61])
(def- rm-bases-64
"bases for rabin-miller determinism to 2^64"
[2 325 9375 28178 450775 9780504 1795265022])
(def- prime-prods
#(* ;(take 13 small-primes))
304250263527210)
(defn rabin-miller-prp?
[n & bases]
(def ps
(if (empty? bases)
(if (int? n) rm-bases-32 rm-bases-64)
bases))
(var d (- n 1))
(var s 0)
(while (zero? (mod d 2))
(++ s)
(set d (div d 2)))
(label result
(each p ps
(var x (powmod p d n))
(when (compare< 1 x (- n 1))
(for r 1 s
(set x (mulmod x x n))
(if-not (compare< 1 x (- n 1)) (break)))
(if-not (compare= x (- n 1))
(return result false))))
true))
(defn prime?
[n]
(cond
(compare<= n 211) (do
(def m (binary-search n small-primes compare<))
(compare= n (in small-primes m)))
(= 0 (jacobi n prime-prods)) false
(rabin-miller-prp? n)))
(defn next-prime
"Returns the next prime number strictly larger than `n`."
[n]
(label result
(def x (+ n 1))
(if (compare<= x 2) (return result 2))
(when (compare<= x 211)
(def m (binary-search x small-primes compare<))
(return result (in small-primes m)))
(def j (mod x 210))
(def m (binary-search j soe-indices compare<))
(var i (+ x (- (in soe-indices m) j)))
(def offs [;(slice soe-offsets m) ;(slice soe-offsets 0 m)])
(forever
(each o offs
(if (prime? i)
(return result i))
(+= i o)))))
(defn primes
"A boundless prime number generator."
[]
(coro
(each n small-primes (yield n))
(def sieve @{221 13 253 11})
(def pg (drop 6 (primes)))
(var p (resume pg))
(var q (* p p))
(var n 211)
(forever
(each offset soe-offsets
(+= n offset)
(def stp (get sieve n))
(cond
stp (do
(put sieve n nil)
(var nxt (div n stp))
(+= nxt (in soe-table (mod nxt 105)))
(while (get sieve (* nxt stp))
(+= nxt (in soe-table (mod nxt 105))))
(put sieve (* nxt stp) stp))
(< n q) (yield n)
(do
(put sieve (+ q (* p (in soe-table (mod p 105)))) p)
(set p (resume pg))
(set q (* p p))))))))
(do
"tests"
(defn trial-division
[n]
(cond
(< n 2) false
(= n 2) true
(= 0 (mod n 2)) false
(label result
(var p 3)
(while (<= (* p p) n)
(if (= 0 (mod n p)) (return result false))
(+= p 2))
true)))
(def big-primes
[100123456789
100529784361
101103163367
101107157131
101111111111
101234567897
101601701401
101740496633
103723971119
105840677923
107928278317
108280301737
109297270343
109661199601
111160669681
111161191111
111417192023
111696081881
112346958007
113013596393
113131311401
113391385603
115829122963
118116011609
118233435679
119218851371
121392937879
123462985007
123511311277
124567987631
125411328001
126704222713
127397154761
129866728583
130147795189
131681894401
131707310437
134141142143
136363636361
137168442221
137438691329
137438953481
139149151157
141516182021
141592653589
142112242123
142356474523
142915724779
143014298809
145679876431
150614187107
152164007921
157119131933
158033851523
162536496481
162849681287
166400805323
166425493681
166425813649
166661866681
168409140869
171190210231
171727482881
172573565537
181030208131
183208285259
186125268239
188886809101
189686906191
189795640003
191198888863
191373251117
194036151289
196719631961
198765432101
200560490131
201820192019
205235235391
220123456789
220999999999
232911191713
235375676171
239651440411
248857367251
251117233141
252097800623
255255255251
255339610381
255339670841
257778788987
260389232731
262364233421
274860381259
275311670611
277535577223
285311670733
293826343073
298999999999
308457624821
313191181151
313473008141
315569251763
316234143227
323866659173
328628315459
333227777777
333332333333
334473276937
337016056721
342054335261
344980016453
351725765537
352573616771
353373727757
355323522737
355711131719
357315397211
358729059641
381654729067
393142151459
411379717319
415926535897
433123191173
445317119867
454280348267
454935396569
463479467461
483148266971
507995698619
512463676147
526858348381
539423223413
542939080319
547716131821
548535559133
553559562581
555553332211
555555555551
555555577777
576529484441
581485876661
582017570143
597325496783
608888888809
608981813029
619010109889
619737131179
635365536563
655372571753
655373525717
657835997711
661991890601
662626111933
665545362821
666666555551
674960891221
676703534071
689101181569
698908090691
700123456789
701234567897
737797911337
744509905409
753352617167
753676171253
757777777777
762394144181
762939453127
765430123567
774004377953
777775777777
816425361649
824633702441
825353008489
837973716761
902659997773
908209935089
921023456789
964991203151
976543312183
977973373171
978058181203
981563119939
999998999999
999999000001
1000000000039
1000000000063
1000000001159
1000000002217
1000000999999
1000008000001
1000123457689
1000123465987
1012345678901
1012346665879
1024688864201
1050100100501
1078834318169
1099511628401
1100011100011
1110110110111
1110111110111
1111111999999
1111118111111
1111120211111
1116295198451
1118810161681
1123529253211
1131313515313
1133328809969
1146890986411
1161468891953
1161880189181
1169769749219
1196980996691
1218131412181
1226280710981
1231507051321
1252099518401
1312961745803
1344409044431
1346140268867
1349623494409
1398023584459
1428571758241
1431535622177
1468986414689
1481754620933
1486287894323
1497631652873
1534139560403
1534696030979
1537267083749
1597233891353
1606081806061
1616110116161
1618033308161
1618033988749
1625800359439
1666666666661
1686910196861
1719004000057
1732565537257
1818116091619
1835211125381
1886881868881
1888081808881
1888691180161
1906101090199
1915318120519
1918960000861
1919110119191
1988006860681
2004006004001
2090609914951
2328316694653
2335557777577
2351123294153
2651541199513
2705863197041
2748779069441
2760889966651
2822916096001
2999999899999
3059220303001
3111111111113
3192202523137
3238428376721
3273534353723
3278744415797
3484606209571
3525357323233
3531577135439
3571113171931
3651980735243
3664533202453
3763863863761
3889108085627
4016465016163
4398042316799
4615392897979
4643668177711
4666666666667
4913685912167
5237732522753
5298696301807
5335039446259
5352640789913
5510987654321
5555555555551
5599297703999
5783272917893
5942636062289
6056529316217
6123456789101
6123584726269
6435409832617
6440261849537
6660000000001
6746328388801
6987191424553
7073313533333
7142857142857
7177111117717
7184500054817
7284717174827
7625597485003
7666466646667
7762868682677
7777774777777
7777775552353
7929188600987
8212668876301
8284590452353
8313667333393
8401414020133
9095665192937
9203225223029
9740985827329
9876543244501
9977770001777
9999919199999
9999991111111
9999999998987
10000000000283
10006181886601
10011111122113
10012345678901
10102323454577
10110111011111
10861196119801
10987654321001
11091501631241
11108452651921
11111333333333
11236129794797
11376747235397
11410337850553
12326713558771
12345678976541
12345678987431
12456898765321
12702232868389
12966235728587
13172937415361
13301522971817
13334390952317
13791315212531
14014014014159
14444999999999
15154262241479
15261281789861
15275163131153
15307263442931
15423094826093
16000611609061
18197520028579
18285670562881
18640889198609
18691113008663
18826507658281
19489355557003
19981091998181
20215043836789
20325153628489
20479999999999
20991129234731
21732382677641
22333555577777
23513892331597
24738041398529
27064032706411
27292996856087
27985032798461
28116440335967
28722900390631
28862180229031
28862180229121
28999134198739
29065965967667
29202114663949
29313741434753
29475832947601
29996224275833
32267562876001
32353461605953
33670369817243
34127252370851
34445333343437
34524689549219
35557777777987
35711101131151
35795175733399
37383940414243
38332212254719
39153010938487
41434547495153
43308631298509
45678910111213
47101299815623
50000010000001
52529393474711
55555544444441
57757791939599
58550453144611
58575551494639
59616771737983
60000000000013
60111111111109
60419650387267
60818091990661
62481801147341
63604045061911
67280421310721
68289980006401
69010616901091
69747782858687
73275315729173
73716967656361
73727170696867
74596893730427
75022592087629
78456580281239
79012338765433
81744303091421
83695120256591
86812553978993
87255058923913
88343532313733
88888888888889
88929267773197
89637484042681
89888786858483
90442568793803
92200000000001
93917775575551
97012345678979
97797371371331
99999999944441
100000000105583
100033000330001
100110101011001
100352695899791
101012345678789
101012345678989
101103107109113
101103107109307
101741582076581
105010123456789
105120615024773
106111115118119
106800081666611
109139149179199
111111151111111
116315256993601
116666666666611
119315717513911
119755797430727
122243323466323
122334444555553
123123454321321
123456136101521
123456789101213
123456797654321
123467898764321
123571113171923
125216512343729
131175323571113
134567897654321
135791113151719
139239739439839
140737488356207
141616531685159
144121100816449
145661743208479
146263628637547
147625901476283
151515151515151
166986101689661
170669145704411
176860696068671
181395559296673
182697604666543
189811909119191
190296921999337
191373383373191
192549818945291
197198199200201
199689672115897
231917131175321
233444000011111
253931039382791
264575131106459
274988895429719
275604952420573
277777788888989
278801239404553
281474993487871
291598227841757
302875106592241
303160419086407
307031012708191
309427398372577
310950284569927
311111111111113
313353373353313
318532111235813
320255973501901
324578696875423
333333373333333
333555577577777
344555666677777
347811194367163
348413484334847
355693655479801
399899999999999
409881449991499
440653748357699
468395662504823
484511389338941
485398038695407
493935332521159
518463576809201
535006138814359
555554444332213
578415690713087
597655503030737
666666666666631
688846502588399
691096906169011
727777887889889
733709701389359
736353433323133
737373737373737
766788787887667
777737777777777
777773333999551
777777355553753
777777715963511
789012345678901
799333555511111
809709509409109
810000006707603
870530414842019
888888877777777
889292677731979
914148152112161
923291713111753
930348707843039
953467954114363
973369606963379
977555333311111
979853562951413
980835832582657
987654321346789
988666444411111
990377764891511
999888866553221
999998727899999
1000000000100011
1003229774283941
1011001110001111
1084051191974761
1086868110868699
1090109110921093
1106098069699111
1111235916285193
1181118711931201
1208152477719361
1211221211212211
1235711131175321
1238926361552897
1255571292290713
1301477951890177
1311753223571113
1311870831664661
1333333333333333
1379131521253133
1391098979592919
1423214346574567
1510553619999637
1593350922240001
1609161918110111
1616668118999191
1657688103077659
1680588011350901
1693182318746371
1737476797107127
1867357583296451
1886999912279359
1889080110806881
1989530586646177
1990474529917009
2035802523820057
2067078579454279
2222226095589337
2231588810593399
2345678911111111
2355457523880889
2357275332573527
2357353373727757
2380072466365871
2468642135797531
2573253757232573
2744337540964913
2828074326968519
2856646544865959
2956667885147129
3129313031303131
3325997869054417
3391382115599173
3637383940414243
3733799929399999
3931520917431241
3940518300860411
4429978144299823
4444280714420857
4547495153555759
4564564564564561
4847464544434241
5111111111111119
5944066965503999
5953474341373129
5999999999899999
6171054912832631
6241156164232441
6435616333396997
6664666366626661
6735249118018991
6988699669998001
7005264275346131
7190597297273099
7523725352733257
7753757725325377
7897466719774591
7897897897897897
8008808808808813
8343656381177203
8690333381690951
8778405442862239
9007199254740847
9007199254740881])
(def bigger-primes
["9293787934331213"
"9436835835813811"
"9779737137133111"
"9949370777987917"
"9999999900000001"
"10000000002065383"
"10106099866009181"
"10112269203786181"
"10234567876543201"
"10269797835402631"
"10939058860032031"
"10968523783746901"
"11018881818881011"
"11120933330250889"
"11333555557777777"
"11417969834487023"
"12223242526272829"
"12348516587712457"
"12962962592588887"
"13000000000000019"
"13107181911273173"
"13173779397737131"
"13176245766932173"
"13311464115101051"
"13457986268975431"
"13661181333262459"
"13666666666666613"
"13666666666666631"
"15125111011152151"
"15795210542132999"
"17000000000000071"
"19019684767739993"
"19149376019158129"
"19611091991096119"
"22222223333355757"
"22263490219831439"
"22439962446379651"
"23456789101112123"
"27714215764134587"
"30000006160000003"
"31896090547101491"
"32494189635056477"
"33334444777777777"
"33832795028841971"
"36028797018963913"
"37033804397792473"
"37124508045065437"
"43142746595714191"
"44444446666688899"
"47784853865664217"
"55350776431903243"
"55453628211510631"
"55555555555555999"
"59604644783353249"
"59649589127497217"
"66555444443333333"
"66606660666066601"
"69668002914515347"
"70000006660000007"
"71117012721071117"
"71828182828182817"
"74747474747474747"
"77777777977777777"
"81615141210198641"
"91521253335394951"
"94350482728405349"
"97654321012345679"
"98764321012345789"
"98765321012456789"
"98765432101456789"
"98823199699242779"
"99194853094755497"
"99999988898898889"
"99999999999899999"
"99999999999999997"
"100109100129100151"
"102598800232111471"
"111101234567891111"
"113127131137139149"
"115578717622022981"
"117116115114112111"
"124277066894534299"
"137153163127255511"
"139717391739171397"
"147578905723456789"
"153722867280912929"
"159067808851610657"
"160409920439929091"
"161111111111111111"
"165678739293946997"
"166667166667000003"
"169177178183185187"
"180811181061181081"
"191128877173556587"
"212345678987654321"
"220724578848317699"
"222222222222222221"
"222333555557777777"
"255455445544554553"
"261798184036870849"
"277233211975331753"
"280743771536011781"
"300000224101777931"
"302301299298295291"
"311172329414761211"
"316227766016837933"
"333333331241623291"
"352521223451364323"
"356112891082412207"
"359681422627177289"
"369121518212427307"
"421538917598915629"
"443231282319131071"
"444444666868899899"
"449317973725128511"
"496481100121144169"
"581321345589144233"
"599999999999899999"
"611091996016611019"
"618681508598750923"
"618819619619660099"
"666666555557777777"
"707274757677788081"
"777777722155555333"
"789795449254776509"
"909090909090909091"
"953947941937929919"
"973006384792642181"
"974383618913296759"
"986143414027351997"
"989466010702279111"
"1000000000000000003"
"1000000000000000009"
"1010203040506070809"
"1010999888887777777"
"1011112131415161719"
"1018264464644628101"
"1023456987896543201"
"1023458697968543201"
"1023489657569843201"
"1023975468645793201"
"1023985674765893201"
"1029001476404750269"
"1034578692968754301"
"1034725698965274301"
"1056208858880824379"
"1109011810160101609"
"1111001101011001111"
"1111111111111111111"
"1111111111111111417"
"1111111112345679001"
"1111918171614101513"
"1111918171614121013"
"1119416189101109149"
"1122334455667788991"
"1147797409030815779"
"1182932213481679513"
"1234567654321234567"
"1234567894987654321"
"1241111124211111421"
"1294584340434854921"
"1341351362137138139"
"1425172824437699411"
"1524157875019052101"
"1618964990108856391"
"1640293473202755797"
"1712329866165608783"
"1732050807568877293"
"1799999999999999999"
"1801241230056600523"
"1971179391939711791"
"1981856124216581891"
"2015121110987654321"
"2223243435546756677"
"2305843009213693951"
"2305843011361841323"
"2327074306453592351"
"2445566867868888899"
"2999998999999999999"
"3203000719597029781"
"3294760066632959081"
"3319393725273939133"
"3331113965338635107"
"3399714628553118049"
"3441173392933731913"
"3842148274728412483"
"3868168229228618683"
"3874204890000000001"
"4968487560233624827"
"5449909581264135103"
"5555555555555555533"
"6082394749206781697"
"6284968128335503367"
"6351332252927568989"
"6568408355712890627"
"6787988999657777797"
"7393913311133193937"
"7532753275327532753"
"7654321012345678987"
"8093914354023690019"
"8223335555577777779"
"8230452606740808761"
"8512677386048191063"
"8888888897888888899"
"9181531581341931811"])
(def big-pseudoprimes
[1373653
1530787
1987021
2284453
3116107
5173601
6787327
11541307
13694761
15978007
16070429
16879501
25326001
27509653
27664033
28527049
54029741
61832377
66096253
74927161
80375707
101649241
161304001
960946321
1157839381
3215031751
3697278427
5764643587
6770862367
14386156093
15579919981
18459366157
19887974881
21276028621
27716349961
29118033181
37131467521
41752650241
42550716781
43536545821
118670087467
307768373641
315962312077
354864744877
457453568161
528929554561
546348519181
602248359169
1362242655901
1871186716981
2152302898747
2273312197621
2366338900801
3343433905957
3461715915661
3474749660383
3477707481751
4341937413061
4777422165601
5537838510751])
(def jacobi-tests
[0 1 1
0 585213579152349 0
8 6720259138994321 1
4 3529872439860193 1
1 8865951665456613 1
2 3751179396445577 1
3 8857077568716611 1
5 6724496004486515 0
6 8046472257518507 -1
7 1018578125212317 -1
2883827 176912159547505 0
28247691 946201408616763 0
19708837 522602379673365 0
13275491 879901654283069 0
60217093 5240865403156329 0
47257397 4366099047223353 0
43746771 3450465360892371 0
28218249 3334597467714585 0
11487823 1239188652493365 0
12927767599 8672338641 1
60316482065 10376205733 1
12718160572 40414256491 -1
120812815652 231360065535 0
117917744387 386108597301 -1
1972212248293 500091801737 -1
6182825798573 657084002447 -1
7724453058855 4314999370589 1
7587489170907 8474042507595 0
9696584664327 8446084420521 0
28511645797572 26659633033591 1
51420282718547 82865707693791 1
21704899141827 30086928426459 0
488044769925626 1 1
260979839748477 5 -1
694655438452973 17 1
207516837635459 5084893 0
396249202456023 14136801 0
121253238725769 94717919 0
357644815400260 27389419382869 0
114127918676057 223507053233071 1
533840837447381 410811659094573 1
444319204437207 407806165804965 0
859875223300786 587409117327703 0
458664366106310 542012645716555 0
424471304525198 2351874079908277 1
557014689515323 1335410288492449 0
153451963855887 1791234085341119 0
7601489827875349 3 1
2448829697631485 9 1
6402808835167554 7 -1
8658574414874615 15 0
3057215490220020 11 -1
3667924217114464 13 -1
1636349321668428 67995153 0
3468280965012603 133198631 0
3447863313944205 127231715 0
4774734459975410 107930263 0
4508174835245718 114427053 0
4341041231878004 107549279 0
5914682029309478 4561197873371 1
2592183347075284 6758805834169 1
8119957211750944 81950119387153 -1
3928288016132477 47031755147963 -1
1567497221445441 676271573590573 1
3849432154013499 124197603917203 1
8989478004945231 813421418837009 1
2352560480455240 546175802817609 0
1397775408790930 301656063758065 0
7300132779063743 923821524336799 -1
2172713837802427 736479480290711 -1
1168603720789234 7699546778746913 1
2323469328644027 2276564181584985 1
4560278553513941 3783685953479921 1
6596617150158661 5927630295508915 1
6385868191676153 4426089584466073 1
8723741114301869 1559345136119577 1
4622568476421908 4170463869060991 1
7126749286305196 6439298346766197 1
6307994731864734 3488208932992171 1
2292810700693211 8152126868198791 1
5615410673045322 4869530827612931 1
7452335230932670 2019154359110699 1
5130615933238582 2882556393517377 1
3525582389132045 3515179569082871 1
3845922998562750 1825997838443061 0
3933618783354591 5378276440500211 0
3851405698358413 4880238089886453 0
5878988890756469 2861006357180511 0
8874876052301087 6136905421930341 0
1435091917964529 2898547151934585 0
6680479536006686 1853547664257451 -1
3833863307507711 1780254117051967 -1
1437351394353646 1828705433408745 -1
1922007955867708 1255085380091051 -1
6837057651372637 5930513687511503 -1
8855677375508641 7328760412908381 -1
3387748941432134 8519439313817431 -1
3550609424633387 1840750301656275 -1])
(def pg (take 10000 (primes)))
(timeit-loop [p :in pg] (assert (prime? p) (string "prime " p)))
(assert (= 104729 (last pg)) "10000th prime")
(assert (all = (map next-prime pg) (drop 1 pg)) "next-prime")
(for n 1 100000
(assert (= (trial-division n) (prime? n)) n))
(assert (all (fn [[a m o]] (= (jacobi a m) o)) (partition 3 jacobi-tests)) "jacobi")
(timeit-loop [p :in big-primes] (assert (prime? p) (string "big primes " p)))
(timeit-loop [p :in (map int/s64 big-primes)] (assert (prime? p) (string "big primes s64 " p)))
(timeit-loop [p :in (map int/s64 bigger-primes)] (assert (prime? p) (string "bigger primes s64 " p)))
(timeit-loop [p :in (map int/u64 bigger-primes)] (assert (prime? p) (string "bigger primes u64 " p)))
(assert (all |(rabin-miller-prp? $ 2 3) big-pseudoprimes) "big pseudoprimes base 2 3")
(assert (not (some rabin-miller-prp? big-pseudoprimes)) "big pseudoprimes rabin-miller"))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment