Project:
 ```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 ``` ```; frobenius pseudoprimality test (define rand #f) (define randint #f) (let ((two31 #x80000000) (a (make-vector 56 -1)) (fptr #f)) (define (mod-diff x y) (modulo (- x y) two31)) ; generic version ; (define (mod-diff x y) (logand (- x y) #x7FFFFFFF)) ; fast version (define (flip-cycle) (do ((ii 1 (+ ii 1)) (jj 32 (+ jj 1))) ((< 55 jj)) (vector-set! a ii (mod-diff (vector-ref a ii) (vector-ref a jj)))) (do ((ii 25 (+ ii 1)) (jj 1 (+ jj 1))) ((< 55 ii)) (vector-set! a ii (mod-diff (vector-ref a ii) (vector-ref a jj)))) (set! fptr 54) (vector-ref a 55)) (define (init-rand seed) (let* ((seed (mod-diff seed 0)) (prev seed) (next 1)) (vector-set! a 55 prev) (do ((i 21 (modulo (+ i 21) 55))) ((zero? i)) (vector-set! a i next) (set! next (mod-diff prev next)) (set! seed (+ (quotient seed 2) (if (odd? seed) #x40000000 0))) (set! next (mod-diff next seed)) (set! prev (vector-ref a i))) (flip-cycle) (flip-cycle) (flip-cycle) (flip-cycle) (flip-cycle))) (define (next-rand) (if (negative? (vector-ref a fptr)) (flip-cycle) (let ((next (vector-ref a fptr))) (set! fptr (- fptr 1)) next))) (define (unif-rand m) (let ((t (- two31 (modulo two31 m)))) (let loop ((r (next-rand))) (if (<= t r) (loop (next-rand)) (modulo r m))))) (init-rand 19380110) ; happy birthday donald e knuth (set! rand (lambda seed (cond ((null? seed) (/ (next-rand) two31)) ((eq? (car seed) 'get) (cons fptr (vector->list a))) ((eq? (car seed) 'set) (set! fptr (caadr seed)) (set! a (list->vector (cdadr seed)))) (else (/ (init-rand (modulo (numerator (inexact->exact (car seed))) two31)) two31))))) (set! randint (lambda args (cond ((null? (cdr args)) (if (< (car args) two31) (unif-rand (car args)) (floor (* (next-rand) (car args))))) ((< (car args) (cadr args)) (let ((span (- (cadr args) (car args)))) (+ (car args) (if (< span two31) (unif-rand span) (floor (* (next-rand) span)))))) (else (let ((span (- (car args) (cadr args)))) (- (car args) (if (< span two31) (unif-rand span) (floor (* (next-rand) span)))))))))) (define (expm b e m) (define (m* x y) (modulo (* x y) m)) (cond ((zero? e) 1) ((even? e) (expm (m* b b) (/ e 2) m)) (else (m* b (expm (m* b b) (/ (- e 1) 2) m))))) (define (digits n . args) (let ((b (if (null? args) 10 (car args)))) (let loop ((n n) (d '())) (if (zero? n) d (loop (quotient n b) (cons (modulo n b) d)))))) (define square? (let ((q11 (make-vector 11 #f)) (q63 (make-vector 63 #f)) (q64 (make-vector 64 #f)) (q65 (make-vector 65 #f))) (do ((k 0 (+ k 1))) ((< 5 k)) (vector-set! q11 (modulo (* k k) 11) #t)) (do ((k 0 (+ k 1))) ((< 31 k)) (vector-set! q63 (modulo (* k k) 63) #t)) (do ((k 0 (+ k 1))) ((< 31 k)) (vector-set! q64 (modulo (* k k) 64) #t)) (do ((k 0 (+ k 1))) ((< 32 k)) (vector-set! q65 (modulo (* k k) 65) #t)) (lambda (n) (if (not (integer? n)) (error 'square? "must be integer") (if (< n 1) #f (if (not (vector-ref q64 (modulo n 64))) #f (let ((r (modulo n 45045))) (if (not (vector-ref q63 (modulo r 63))) #f (if (not (vector-ref q65 (modulo r 65))) #f (if (not (vector-ref q11 (modulo r 11))) #f (let ((q (isqrt n))) (if (= (* q q) n) q #f)))))))))))) (define (jacobi a n) (if (not (and (integer? a) (integer? n) (positive? n) (odd? n))) (error 'jacobi "modulus must be positive odd integer") (let jacobi ((a a) (n n)) (cond ((= a 0) 0) ((= a 1) 1) ((= a 2) (case (modulo n 8) ((1 7) 1) ((3 5) -1))) ((even? a) (* (jacobi 2 n) (jacobi (quotient a 2) n))) ((< n a) (jacobi (modulo a n) n)) ((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (jacobi n a))) (else (jacobi n a)))))) (define (inverse x m) (let loop ((x x) (b m) (a 0) (u 1)) (if (zero? x) (if (= b 1) (modulo a m) (error 'inverse "must be coprime")) (let* ((q (quotient b x))) (loop (modulo b x) x u (- a (* u q))))))) (define (frobenius-pseudoprime? n) (let loop ((a (randint 1 n)) (b (randint 1 n))) (display "a = ") (display a) (newline) (display "b = ") (display b) (newline) (let ((d (- (* a a) (* 4 b)))) (display "d = ") (display d) (newline) (if (or (= a b) (square? d) (not (= (gcd (* 2 a b d) n) 1))) (loop (randint 1 n) (randint 1 n)) (let* ((m (/ (- n (jacobi d n)) 2)) (w (modulo (- (* a a (inverse b n)) 2) n))) (display "m = ") (display m) (newline) (display "w = ") (display w) (newline) (let loop ((u 2) (v w) (ms (digits m 2))) (display "u, v = ") (display u) (display " ") (display v) (newline) (if (pair? ms) (if (zero? (car ms)) (loop (modulo (- (* u u) 2) n) (modulo (- (* u v) w) n) (cdr ms)) (loop (modulo (- (* u v) w) n) (modulo (- (* v v) 2) n) (cdr ms))) (if (positive? (modulo (- (* u w) (* 2 v)) n)) #f (let ((x (expm b (/ (- n 1) 2) n))) (zero? (modulo (- (* u x) 2) n))))))))))) (define (strong-pseudoprime? n a) (let loop ((r 0) (s (- n 1))) (if (even? s) (loop (+ r 1) (/ s 2)) (if (= (expm a s n) 1) #t (let loop ((r r) (s s)) (cond ((zero? r) #f) ((= (expm a s n) (- n 1)) #t) (else (loop (- r 1) (* s 2))))))))) (define (prime? n) (let loop ((ps (list 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))) (if (pair? n) (if (zero? (modulo n (car ps))) (= n (car ps)) (loop (cdr ps))) (and (strong-pseudoprime? n 2) (strong-pseudoprime? n 3) (frobenius-pseudoprime? n))))) (display (prime? 3825123056546413051)) (newline) (display (prime? (- (expt 2 89) 1))) (newline) ```
 ```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 ``` ```a = 165967461024650441079709501 b = 7001028712697292951128644051 d = 27545198118968863237666585238403300920261682035092797 m = 1912561528273206525 w = 2213324777930859275 u, v = 2 2213324777930859275 u, v = 2213324777930859275 2270875488901952443 u, v = 2007894263127594961 1687991195430853439 u, v = 2239286715426570691 476815029326724428 u, v = 247086105996717166 3806344061966365995 u, v = 1462298248057697650 2583137650687404949 u, v = 2821499742648653877 2984459978153769893 u, v = 3427330740843507973 420094092277495933 u, v = 1112507905812958328 1081102756955841060 u, v = 3779295351914804736 2306888088755405988 u, v = 1175302810330265884 3368103141916367573 u, v = 2720568835640293085 2811135225279095600 u, v = 571680156765576615 3573113230820201784 u, v = 29854169404030427 1554230182923239672 u, v = 2366505738093310667 112015345582914815 u, v = 124314474894071747 3614378250012894611 u, v = 638203658611125557 2234915036850351448 u, v = 2619730335858260597 3685888905819189923 u, v = 190955847588099271 244254454003494807 u, v = 354219928337041106 360631916832940051 u, v = 656127758148161019 535746034812376895 u, v = 3133659866379517321 3281598211033578738 u, v = 3134258387534226122 2216419471685290835 u, v = 3441336101617737312 758628971702849920 u, v = 1252744110529130392 2317020009311125011 u, v = 463308875267645875 2315653402677417967 u, v = 2437472882617568381 3329765746886252669 u, v = 1195581244071628614 1885965194842555963 u, v = 815910564578299445 761600033230221829 u, v = 2554849391985244961 1684611457422367497 u, v = 1884889623548593373 1168084421426390107 u, v = 723257737991525054 2738050095009420346 u, v = 3461238375413078921 757738372030983585 u, v = 814681908936713852 2685241255179823957 u, v = 1285540921692471745 2797902866954000563 u, v = 502015006379838323 181039916727740011 u, v = 1498727703948185907 2085803332842892739 u, v = 3584496407325939422 2606324993105816375 u, v = 875388091091613033 1526995320229481309 u, v = 1737311941415014842 2635016795040862007 u, v = 2604245554657808676 450351923215198852 u, v = 2113424805407903020 3341939848305923324 u, v = 126271525932539236 1591730811051306670 u, v = 1300137571171722631 442926449426522190 u, v = 3268569301359976101 2132769228941037542 u, v = 160988969903037765 1420853032314982931 u, v = 2794284400777774734 2650427347995800224 u, v = 3050992731629812822 1519340126424836675 u, v = 3153899986792930395 963845504333573250 u, v = 1393385294753283391 2354753136787048722 u, v = 3425009176660774384 2062528201300819569 u, v = 1192798244768972750 1320419794991335810 u, v = 1596748828049567153 484224608426600773 u, v = 2607981864631459279 1452653866977478065 u, v = 1668879768778497193 2121928893053266912 u, v = 3336417290281957086 2893888326335817620 u, v = 2632615708435959484 3099090873433722937 u, v = 1780180586542791111 1189735485790610065 u, v = 565715865059161893 877496693304216728 u, v = 2416196454755507897 2262497401412406207 u, v = 283400770974670135 2616942119369586582 u, v = 3683053673758867278 525037557735457615 #f a = 914083492632061994498674630635828301 b = 1158281428926235899874054270740091141 d = 835548631502428934583168956376458821790355103078220574397467805392182037 m = 309485009821345068724781056 w = 586275862631858978966654741 u, v = 2 586275862631858978966654741 u, v = 586275862631858978966654741 222818930194943302282193037 u, v = 222818930194943302282193037 472448015783325842703261055 u, v = 608221536678616678000770884 185031263961022821664601546 u, v = 231000882874457832498928448 86679538593585339800410236 u, v = 575419614636922147489983213 303804712254463661234872018 u, v = 583377248176362532165524082 182268758263599343945433348 u, v = 243332605159384473961528403 541960925303064397859732737 u, v = 533884114922593974817165737 113708178145274743680529536 u, v = 120598983929906653827933772 203578926045793448107208835 u, v = 152670204332325388404086334 528346861251835948446382290 u, v = 524733476303755191928318227 559972895154860817422589365 u, v = 502858817618152163566686904 402866011170952487841956440 u, v = 252731397871885496699659378 353792222169477055594085925 u, v = 20239631092543810217580616 156172885859107064092707127 u, v = 108500424028539798330421624 386907417080473180698479744 u, v = 402497375382157821306660581 152987881961959729155355293 u, v = 468949715032562161239554941 439336459037533738016251372 u, v = 322138411401626962915114516 337629406124665340496059662 u, v = 35117208844729618711853894 475971093452760102271154626 u, v = 414498844117477900439509913 404610073682792944711946248 u, v = 93232056164360779098620083 614087033948038824191008708 u, v = 561746353032347998040788583 44345896842760657986148596 u, v = 345679684103324143302459707 514469134689958720950661086 u, v = 178456670777488872064520988 531725193651281466549146425 u, v = 43248382760419421436179442 319590580507104130444145725 u, v = 185644152642008105720968955 203100210385000430907929077 u, v = 42246687381995543671951745 565301615197387230734176361 u, v = 121634640024192876242090775 439214727601671891188535276 u, v = 581042888334574862714072800 207519258949365447648741813 u, v = 257882220502393960426383961 439195480949173726241016007 u, v = 72796502827656116240532658 274920237267095047553403422 u, v = 192512313759206296780953156 12730555809530004053844484 u, v = 137378001358651000999179804 595502614618067240344552746 u, v = 5306781747643928509952333 22930964063821211122926619 u, v = 31647249704115216071043950 546440196530502965453256380 u, v = 349382349085327612444907591 10880140381866697208284237 u, v = 602512363726227940019647144 355942326745655048321881830 u, v = 566866788129004968968502493 402494326786756541330285298 u, v = 421204797995896354011736629 523516856899995139072136589 u, v = 107038802935660997511809113 183727771346768775073657766 u, v = 142073608212538095579103859 460991736506446276414255270 u, v = 596282407708249792522875869 76816240437964608427428270 u, v = 147677136255624845695316238 92566507221375709768338393 u, v = 70744199768718234362305618 277188238799938103486843732 u, v = 373773981737788761451432714 614017252472516521413240768 u, v = 255677975746776137288568642 189283308294216110169629150 u, v = 93085320541245517299619774 5631787617882395529445936 u, v = 149522066443119121057487787 423817947028916381825632184 u, v = 262103703520689899985856720 442874579868298284834961681 u, v = 242111344091848421265535266 178921122160709602631880395 u, v = 339342656651553528284810034 126553337282104772244464109 u, v = 17996105943286784488675969 55495115347931865875927299 u, v = 109171394508717662628567835 608043096846806442221523438 u, v = 572132580670010213141390217 494764146356894402802938202 u, v = 15376648476182311545974234 168381012423552332181715108 u, v = 124010230354554384776727076 226132925870641626015761917 u, v = 404015955584838738938921120 403250195004505731753920744 u, v = 409049912911409888925227415 218469003585865917437524808 u, v = 83916026489839527015329338 73044332107185518473795006 u, v = 90530710286622429751519544 305081606875290446603767905 u, v = 68903624414125023878092436 339034140052122395578862370 u, v = 194421710740941861272485201 153301850664856298900165097 u, v = 390099320673176995051444373 205379459735199206280069788 u, v = 539146489093062647936046750 297941641589569752750137725 u, v = 456227527134349108981108860 180960419325054857957703440 u, v = 395083570046511465239138240 494557825876180209312514712 u, v = 403386164259153140107341915 228029856437898736492402699 u, v = 277742636415709085384891628 170657539623993910266185076 u, v = 147854539654622100061648694 538715612041615389331230398 u, v = 496000786316959896649290117 359612327589950726488998756 u, v = 269946114720766867952294779 115787439213775592435477887 u, v = 82667829751602586256008147 128026318930869743942547150 u, v = 423631779015658699100026514 112778465103833074653642045 u, v = 391483091139517355082097370 378650394108506984963553312 u, v = 369034084895607829940757934 150726000577289492954695552 u, v = 95363944458379653483098912 433611047043764054052344408 u, v = 384476859026094186848976519 527283986220003536619145495 u, v = 559385705224714492934336595 168961110990735490351097314 u, v = 233840598130087066597597741 99939359127185672407443338 u, v = 51436093599450562821421507 461155191028561838665105500 u, v = 555372695984532823653348101 595235604532038907064525717 u, v = 203004199961331080879316213 467085185982919729051489829 u, v = 386196223548689448070578677 99287870137644420437277396 u, v = 553043523520208682228768713 70776060488636309173459378 u, v = 476813087387959183305530147 517497328628547721586974348 u, v = 68271965619881530889857423 383027336777526488370938620 u, v = 35184372088832 48697543093884042116867575 u, v = 0 589302486790599960080862723 u, v = 618970019642690137449562109 32694157010831158482907370 #t ```