Page 1 of 2 12 LastLast
Results 1 to 10 of 16

Thread: Miller-Rabin primality test

  1. #1
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152

    Miller-Rabin primality test

    I tried to implement the Miller-Rabin primality test, in Racket (Scheme).

    http://en.wikipedia.org/wiki/Miller%...primality_test

    Why did I use Racket? Because, I like it.

    The test is probabilistic, not deterministic. It doesn't absolutely determine if an integer is prime. You set the probability, P, and if the test says that a particular integer is prime, then, the likelihood is at least P, that it is actually prime.

    My experience with my implementation of it, is that it sometimes makes errors. So, maybe I did something wrong.

    There are 109 primes less than 600, and the test said that 24 of them are not prime. However, it seems that as the candidate integer gets bigger, the procedure becomes more accurate.

    Below is shown the Racket code, and some tests I did on integers in the Racket REPL (read-eval-print loop).

    For each test, the output is either #t (true), or #f (false), and the time in seconds to do the calculation, is also indicated.

    For the tests, I used primes which I got from,

    http://en.wikipedia.org/wiki/List_of_prime_numbers .

    I did 8 tests, on primes from the web page.

    The first 6 are Carol primes, and the next 2 are Bell number primes.

    The procedure correctly indicated that all 8 are prime.

    Then, to make sure that the procedure wasn't just answering true for every large integer, I did the 8 tests again.

    This time, I added the value 2, to each of the original 8 primes, and tested those integers.

    For the second set of 8 tests, the procedure indicated that the first 7 were not prime, as I hoped it would.

    However, for the 8th test in the second set, the procedure indicated that the integer is prime.

    At that point I thought that the procedure was no good.

    But, I decided to try one more thing.

    I went to,

    http://www.wolframalpha.com/ .

    I executed the command,

    "PrimeQ[359334085968622831041960188598043661065388726959079839]".

    "Super-Mathematica" indicated, that the number is indeed prime.

    So, I felt better.

    In all of the tests, I set the probability, P, equal to, 0.999.



    ; code ---------------------------------------------------------------------------------------------
    
    #lang racket
    
    (define ran-max 4294967087)
    
    (define (square x) (* x x))
    
    (define (expmod base exp m)
      (cond ((= exp 0) 1)
            ((even? exp) (remainder (square (expmod base (/ exp 2) m)) m))
            (else (remainder (* base (expmod base (- exp 1 ) m)) m))))
    
    (define (power-of-two-factor x)
      (define (power-of-two-internal x n)
        (if (= (modulo x 2) 1)
            n
            (power-of-two-internal (/ x 2) (+ n 1))))
      (power-of-two-internal x 0))
    
    (define (factor-even-number x)
      (let ((po2 (power-of-two-factor x)))
        (values po2 (/ x (expt 2 po2)))))
    
    (define (big-random n)
      (let ((r (round (/ (* n (random ran-max)) ran-max))))
        (if (= r 1) (big-random n) r)))
    
    (define (num-repetitions p)
      (define 1-p (- 1 p))
      (ceiling (- (/ (log 1-p) (log 4)))))
    
    (define (is-prime? n prob)
      (cond
        ((< n 5) #f)
        ((even? n) #f)
        (else
         (define-values (s d) (factor-even-number (- n 1)))
         (define k (num-repetitions prob))    
         (define (outer-loop k)
           (define a (big-random (- n 2)))
           (define x (expmod a d n))
           (cond
             ((= k 0) #t)
             ((or (= x 1) (= x (- n 1))) (outer-loop (- k 1)))
             (else (inner-loop (- s 1) x))))    
         (define (inner-loop i y)
           (set! y (modulo (* y y) n))
           (cond
             ((or (= y 1) (= i 0)) #f)
             ((= y (- n 1)) (outer-loop (- k 1)))
             (else (inner-loop (- i 1) y))))     
         (outer-loop k))))
    
    (define (time-is-prime? n p)
      (let ((t1 (current-milliseconds)))
        (let ((prime? (is-prime? n p)))
          (let ((t2 (current-milliseconds)))
            (let ((tt (/ (- t2 t1) 1000.0)))
              (values prime? tt))))))
    
    ; REPL interactions --------------------------------------------------------------------------------
    
    Welcome to DrRacket, version 5.1 [3m].
    Language: racket.
    > (time-is-prime? 68718952447 0.999)
    #t
    0.141
    > (time-is-prime? 274876858367 0.999)
    #t
    0.122
    > (time-is-prime? 4398042316799 0.999)
    #t
    0.108
    > (time-is-prime? 1125899839733759 0.999)
    #t
    0.113
    > (time-is-prime? 18014398241046527 0.999)
    #t
    0.097
    > (time-is-prime? 1298074214633706835075030044377087 0.999)
    #t
    0.274
    > (time-is-prime? 35742549198872617291353508656626642567 0.999)
    #t
    0.17
    > (time-is-prime? 359334085968622831041960188598043661065388726959079837 0.999)
    #t
    1.488
    > (time-is-prime? 68718952449 0.999)
    #f
    0
    > (time-is-prime? 274876858369 0.999)
    #f
    0.001
    > (time-is-prime? 4398042316801 0.999)
    #f
    0.123
    > (time-is-prime? 1125899839733761 0.999)
    #f
    0.001
    > (time-is-prime? 18014398241046529 0.999)
    #f
    0.001
    > (time-is-prime? 1298074214633706835075030044377089 0.999)
    #f
    0.002
    > (time-is-prime? 35742549198872617291353508656626642569 0.999)
    #f
    0.118
    > (time-is-prime? 359334085968622831041960188598043661065388726959079839 0.999)
    #t
    3.527
    >
    

    Last edited by danbaron; 18-04-2011 at 02:37.
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  2. #2
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    I tried a 100 digit prime, which I got from,

    http://primes.utm.edu/lists/small/small.html#80 .

    Then, I added 2 to it, and tried that.

    It seems to work OK.

    Welcome to DrRacket, version 5.1 [3m].
    Language: racket.
    > (time-is-prime? 2074722246773485207821695222107608587480996474721117292752992589912196684750549658310084416732550077 0.999)
    #t
    5.463
    > (time-is-prime? 2074722246773485207821695222107608587480996474721117292752992589912196684750549658310084416732550079 0.999)
    #f
    0.086
    >
    
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  3. #3
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    I tried a 200 digit prime, which I got from,

    http://primes.utm.edu/lists/small/small2.html .

    Then, I added 2 to it, and tried that.

    It still seems to work OK.

    Welcome to DrRacket, version 5.1 [3m].
    Language: racket.
    > (time-is-prime? 58021664585639791181184025950440248398226136069516938232493687505822471836536824298822733710342250697739996825938232641940670857624514103125986134050997697160127301547995788468137887651823707102007839 0.999)
    #t
    3.899
    > (time-is-prime? 58021664585639791181184025950440248398226136069516938232493687505822471836536824298822733710342250697739996825938232641940670857624514103125986134050997697160127301547995788468137887651823707102007841 0.999)
    #f
    0.165
    >
    

    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  4. #4
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    I tried a 1050 digit prime, which I got from,

    http://primes.utm.edu/primes/lists/all.txt .

    It is #5706 in the list, 11^1008 + 998672782.

    Then, since the last digit of the number is 3, I tried the number plus 4.

    It still seems to work OK.

    Welcome to DrRacket, version 5.1 [3m].
    Language: racket.
    > (define a (+ (expt 11 1008) 998672782))
    > a
    529452056448793684463171701371965165599024891268044535141171238938714197123500025362132544081526670388715712690296898616337461300706986447973191761836855880775521428076605793975641811081746047876506709541924966850486767835110575769384906730749667777016471469177953052696916660058190262098621182600138065372500625688141864409509834051913316703890710975152285824832755202082018739474065720249911588366808071879632022058226980343003483610162880483389826691629327795587270718562153925742749978926314647876232385694302214489674918510983781514194149348746847611591937897776475529798214185181022673653771123333515766636914861846761384172063717773371541993327403277856399026738370062644554264826051111760730911233180455686309777627985916941132711706748812406408935146311575975928657195294322425757348388741593204454372595463565373254895887022065283560617002190260565624775467971401594792170355772678679741987491858658211032034438449802653787330079926121284263030291537500873819738703978995623111675987883052680605464876084482213424458493450872074280702291663
    > (define a+4 (+ a 4))
    > (time-is-prime? a 0.999)
    #t
    15.988
    > (time-is-prime? a+4 0.999)
    #f
    2.402
    >
    
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  5. #5
    Member Johannes's Avatar
    Join Date
    Nov 2010
    Location
    Wuustwezel, Belgium
    Age
    56
    Posts
    95
    Rep Power
    25
    Dan,

    My program (below) agrees with all of your checks. What is striking is that my script is much faster for the first examples, performs about equally for the 200-digit numbers, and is painfullly slower for the 1050-digit numbers. At a guess Racket is using FFT multiplication for its big integers.

    One question: when you calculate the number of loops, the -log(1-p)/log(4) comes out as 4.98. Is that result rounded or truncated? Because if it is truncated your program is incorrect with respect to the probability.
    ' > BigInt_PrimeMillerRabin
     
    Uses "Console"
    Module "BigInt"
    Alias String As BigInt
     
    Randomize
     
    TimeMR("68718952447")
    TimeMR("274876858367")
    TimeMR("4398042316799")
    TimeMR("1125899839733759")
    TimeMR("18014398241046527")
    TimeMR("1298074214633706835075030044377087")
    TimeMR("35742549198872617291353508656626642567")
    TimeMR("359334085968622831041960188598043661065388726959079837")
     
    TimeMR("68718952449")
    TimeMR("274876858369")
    TimeMR("4398042316801")
    TimeMR("1125899839733761")
    TimeMR("18014398241046529")
    TimeMR("1298074214633706835075030044377089")
    TimeMR("35742549198872617291353508656626642569")
    TimeMR("359334085968622831041960188598043661065388726959079839")
     
    TimeMR("2074722246773485207821695222107608587480996474721117292752992589912196684750549658310084416732550077")
    TimeMR("2074722246773485207821695222107608587480996474721117292752992589912196684750549658310084416732550079")
     
    TimeMR("58021664585639791181184025950440248398226136069516938232493687505822471836536824298822733710342250697739996825938232641940670857624514103125986134050997697160127301547995788468137887651823707102007839")
    TimeMR("58021664585639791181184025950440248398226136069516938232493687505822471836536824298822733710342250697739996825938232641940670857624514103125986134050997697160127301547995788468137887651823707102007841")
     
    TimeMR("529452056448793684463171701371965165599024891268044535141171238938714197123500025362132544081526670388715712690296898616337461300706986447973191761836855880775521428076605793975641811081746047876506709541924966850486767835110575769384906730749667777016471469177953052696916660058190262098621182600138065372500625688141864409509834051913316703890710975152285824832755202082018739474065720249911588366808071879632022058226980343003483610162880483389826691629327795587270718562153925742749978926314647876232385694302214489674918510983781514194149348746847611591937897776475529798214185181022673653771123333515766636914861846761384172063717773371541993327403277856399026738370062644554264826051111760730911233180455686309777627985916941132711706748812406408935146311575975928657195294322425757348388741593204454372595463565373254895887022065283560617002190260565624775467971401594792170355772678679741987491858658211032034438449802653787330079926121284263030291537500873819738703978995623111675987883052680605464876084482213424458493450872074280702291663")
    TimeMR("529452056448793684463171701371965165599024891268044535141171238938714197123500025362132544081526670388715712690296898616337461300706986447973191761836855880775521428076605793975641811081746047876506709541924966850486767835110575769384906730749667777016471469177953052696916660058190262098621182600138065372500625688141864409509834051913316703890710975152285824832755202082018739474065720249911588366808071879632022058226980343003483610162880483389826691629327795587270718562153925742749978926314647876232385694302214489674918510983781514194149348746847611591937897776475529798214185181022673653771123333515766636914861846761384172063717773371541993327403277856399026738370062644554264826051111760730911233180455686309777627985916941132711706748812406408935146311575975928657195294322425757348388741593204454372595463565373254895887022065283560617002190260565624775467971401594792170355772678679741987491858658211032034438449802653787330079926121284263030291537500873819738703978995623111675987883052680605464876084482213424458493450872074280702291667")
     
    WaitKey
     
    Sub TimeMR(n As BigInt,Optional p As Ext)
    Local f As Boolean
    Local t As Quad
    n = BigInt_FromString(n)
    If p=0 Then p=.999
    HiResTimer_Init
    HiResTimer_Get
    f = MillerRabin(n,p)
    t=HiResTimer_Delta
    PrintL BigInt_ToString(n)
    If f Then
      Print  "PRIME "
    Else
      Print  "      "
    EndIf
    PrintL Format$(t/1000000,"0.000")
    End Sub
     
    Function MillerRabin(n As BigInt,p As Ext) As Boolean
    Local a,d,m,n1,n2,t,w,x As BigInt
    Local i,j,s As DWord
    Local b As Byte
    Local q As Ext = 1
    n1 = BigInt_Dec(n)
    n2 = BigInt_Dec(n1)
    d = n1
    While BigInt_IfEven(d)
      s += 1
      d = BigInt_Div(d,2)
    Wend
    w = BigInt_FromInteger(1)
    t = BigInt_FromInteger(2)
    m = BigInt_FromInteger(1000000000)
    p = 1-p
    Do
      q /= 4
      a = BigInt_Div(BigInt_Mul(n2,BigInt_FromInteger(Rnd(1,1000000000))),m)
      a = BigInt_Max(a,t)
      x = w
      For i=2 To Len(d)
        b = Asc(Mid$(d,i,1))
        For j=1 To 8
          If (b And 1)=1 Then x = BigInt_Mod(BigInt_Mul(x,a),n)
          a = BigInt_Mod(BigInt_Mul(a,a),n)
          SHIFT RIGHT b,1
        Next j
      Next i
      If BigInt_IfEQ(x,w) Then Iterate Do
      If BigInt_IfEQ(x,n1) Then Iterate Do
      If s>1 Then
        For i=2 To s
          x = BigInt_Mod(BigInt_Mul(x,x),n)
          If BigInt_IfEQ(x,w) Then Return FALSE
          If BigInt_IfEQ(x,n1) Then Iterate Do
        Next i
      EndIf
      Return FALSE
    Loop Until q<=p
    Return TRUE
    End Function
    
    Last edited by Johannes; 18-04-2011 at 16:50.
    Boole and Turing, help me!

    Primary programming: 200 MHz ARM StrongARM, RISC OS 4.02, BASIC V, ARM assembler.
    Secondary programming: 3.16 GHz Intel Core 2 Duo E8500, Vista Home Premium SP2, thinBasic, x86 assembler.

  6. #6
    Member Johannes's Avatar
    Join Date
    Nov 2010
    Location
    Wuustwezel, Belgium
    Age
    56
    Posts
    95
    Rep Power
    25
    Dan, remember this one?
    Quote Originally Posted by danbaron View Post
    Factor the following composite integer.

    21130819027319923178509622632727980622577522491335986162202750512019525136688895741633235387039306968582318693021114825649130164135868576427684517893
    Factoring is still "impossible" but with the Miller-Rabin test it takes about half a second on my system to determine that that number is, in fact, composite.
    Boole and Turing, help me!

    Primary programming: 200 MHz ARM StrongARM, RISC OS 4.02, BASIC V, ARM assembler.
    Secondary programming: 3.16 GHz Intel Core 2 Duo E8500, Vista Home Premium SP2, thinBasic, x86 assembler.

  7. #7
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    (I don't think anyone will call you lazy, Johannes.)

    "ceiling", so, hopefully, 4.98 --> 5.

    (But, it should do at least one repetition no matter what the value of p is. As it is currently written, if you enter p as 0, then it will do no repetitions.)

    (define (num-repetitions p)
      (define 1-p (- 1 p))
      (ceiling (- (/ (log 1-p) (log 4)))))
    

    I remember it, so, I'll do it (or, more accurately, Racket will do it).

      Welcome to DrRacket, version 5.1 [3m].
      Language: racket.
      > (time-is-prime? 21130819027319923178509622632727980622577522491335986162202750512019525136688895741633235387039306968582318693021114825649130164135868576427684517893 0.999)
      #f
      0.133
    >
    

    I ran two more numbers from the list,

    http://primes.utm.edu/primes/lists/all.txt .

    #5680, 546! - 1, 1260 digits

    #5647, 83 * 2^5318 - 1, 1603 digits

    Welcome to DrRacket, version 5.1 [3m].
    Language: racket.
    > (define a (- (fact 546) 1))
    > (time-is-prime? a 0.999)
    #t
    29.135
    > a
    141302009261418325453762621393526286698658617796079910071353871488907640504160743565501241236107733762080837422133568434484213467457560621538131284091550631968816592748754460044234605778793037038634578162615257746896342308546754776364849349150715175577496877924111556740124186753117294177162303783405503829723023736326990734833969973121886593758921817486418079092325492863765948445727040946819291138129804200523981866521636537165676001267152151758174473660767772438139761519901508251733325317410707450984627544775891183776230678501637671689726049915871033588824272610663202323719389931283164345628668801970555064368945390634202638535930076803929062184736870310291943624402875877497305289749392357778468456405062145336924250773407066015135454416984550131804906976797622941301057536164271657477854025104530247400241587703885280095204104624235632766825184645847127845152855411530623429061372120607469651273815795437429759478958062539197279629575104299984571584633649528868101236905793259525917436877798780909275667980803203926354424450148098596386681966959342498581232330994696677186100809764161773397785695845695608487026255462399999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999
    > (define b (+ a 2))
    > (time-is-prime? b 0.999)
    #f
    3.927
    > (define c (- (* 83 (expt 2 5318)) 1))
    > (time-is-prime? c 0.999)
    #t
    66.261
    > c
    6260298383829657225286112584739828350961566071813614817090655661726702203000727443155352137677415727639270005152424260690857553900368341450346970741308511417934442673066866644403367228220364439878425897186245253177691894131530314490537670719045204854370814516293642618537389174767012520242616529405673247696080934123952725656223766315340714614855806548932055139425179151479881206038701763802779138916668883306784294260180704008909313195142729581852896420768539776228014321360869785451249902588785914955877269128146778963486558120668354631060167501141743903122509454201202547688715196727040523077806030531988397663538853730963809344267934674614590712227267805424158243982282344265009362357321981931239055386336840472883326309867734816135908681990888503697042094651571968230153780164618679215781152726244096578224696053246041697606541689435068727931476353894324623910695079566421040585311985103805812421650786510353427759898861776000121326090856058979212077591696177260616037165555751507285878767478858118625647588877281294325434409067338807313030610213165442239480105674376617197846030146545262330219853080619480300992036754949135625216566071983425953282414885683202259216717086223201216107231214643841453821800479705443657770734535275127870398933661350313638594922371398078999763881759388037191499672886496845677727495900732005269358724467804503279937265707718167486322051615118497652040645286382390708500068385222336755432189454238226603697941580242718059483611416057744941499700978593848565252669806462644516538399273148998665777617816962249185512503379530702903585164261671957932841026458847073533951
    > (define d (+ c 2))
    > (time-is-prime? d 0.999)
    #f
    4.931
    >
    

    When I tried numbers bigger than #5647, something went wrong, because Racket calculated forever.

    We probably agree that this way is much much better than the deterministic way.

    I agree, Racket probably uses FFT to multiply integers. It is amazingly fast.



    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  8. #8
    Member Johannes's Avatar
    Join Date
    Nov 2010
    Location
    Wuustwezel, Belgium
    Age
    56
    Posts
    95
    Rep Power
    25
    Quote Originally Posted by danbaron View Post
    "ceiling", so, hopefully, 4.98 --> 5.
    Oh snap!

    Should have seen that.

    I was wondering what you were doing there and then I thought to check Wikipedia to see what the reliability was. Amazing that with only 5 random tries you can get less than a one-in-a-thousand error. Ten tries and it's less than one in a million.

    I did the repetitive division for the probability because it is probably (ahem) faster than the two LOGs. But you get the prize for cryptic code on that one.
    Boole and Turing, help me!

    Primary programming: 200 MHz ARM StrongARM, RISC OS 4.02, BASIC V, ARM assembler.
    Secondary programming: 3.16 GHz Intel Core 2 Duo E8500, Vista Home Premium SP2, thinBasic, x86 assembler.

  9. #9
    Member Johannes's Avatar
    Join Date
    Nov 2010
    Location
    Wuustwezel, Belgium
    Age
    56
    Posts
    95
    Rep Power
    25
    Quote Originally Posted by danbaron View Post
    When I tried numbers bigger than #5647, something went wrong, because Racket calculated forever.
    I'm running #5551 ((2^10211-1)/306772303457009724362047724636324707614338377) with 3,030 digts now. The calculation of the number itself took almost no time at all. The Miller-Rabin test is... running.
    Boole and Turing, help me!

    Primary programming: 200 MHz ARM StrongARM, RISC OS 4.02, BASIC V, ARM assembler.
    Secondary programming: 3.16 GHz Intel Core 2 Duo E8500, Vista Home Premium SP2, thinBasic, x86 assembler.

  10. #10
    Member Johannes's Avatar
    Join Date
    Nov 2010
    Location
    Wuustwezel, Belgium
    Age
    56
    Posts
    95
    Rep Power
    25
    Quote Originally Posted by Johannes View Post
    The Miller-Rabin test is... running.
    Two thousand seven hundred and nineteen seconds. And the number is prime with a certainty of at least 99.9%.

    Dan, most likely Racket is running into memory problems with the larger numbers. I think I remember you having to define the amount of memory you want to work with.
    Boole and Turing, help me!

    Primary programming: 200 MHz ARM StrongARM, RISC OS 4.02, BASIC V, ARM assembler.
    Secondary programming: 3.16 GHz Intel Core 2 Duo E8500, Vista Home Premium SP2, thinBasic, x86 assembler.

Page 1 of 2 12 LastLast

Members who have read this thread: 0

There are no members to list at the moment.

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •