[ create a new paste ] login | about

Project: programmingpraxis
Link: http://programmingpraxis.codepad.org/mbEWx4S1    [ raw code | output | fork ]

Scheme, pasted on Nov 11:
; 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)


Output:
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


Create a new paste based on this one


Comments: