34 people like it.

Miller–Rabin primality test

Miller–Rabin primality test is an algorithm which determines whether a given number is probable prime. For more information go to http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test

 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: 

open System.Numerics

///This implementation is based on the Miller-Rabin Haskell implementation 
///from http://www.haskell.org/haskellwiki/Testing_primality
let pow' mul sq x' n' = 
    let rec f x n y = 
        if n = 1I then
            mul x y
        else
            let (q,r) = BigInteger.DivRem(n, 2I)
            let x2 = sq x
            if r = 0I then
                f x2 q y
            else
                f x2 q (mul x y)
    f x' n' 1I
        
let mulMod (a :bigint) b c = (b * c) % a
let squareMod (a :bigint) b = (b * b) % a
let powMod m = pow' (mulMod m) (squareMod m)
let iterate f = Seq.unfold(fun x -> let fx = f x in Some(x,fx))

///See: http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
let millerRabinPrimality n a =
    let find2km n = 
        let rec f k m = 
            let (q,r) = BigInteger.DivRem(m, 2I)
            if r = 1I then
                (k,m)
            else
                f (k+1I) q
        f 0I n
    let n' = n - 1I
    let iter = Seq.tryPick(fun x -> if x = 1I then Some(false) elif x = n' then Some(true) else None)
    let (k,m) = find2km n'
    let b0 = powMod n a m

    match (a,n) with
        | _ when a <= 1I && a >= n' -> 
            failwith (sprintf "millerRabinPrimality: a out of range (%A for %A)" a n)
        | _ when b0 = 1I || b0 = n' -> true
        | _  -> b0 
                 |> iterate (squareMod n) 
                 |> Seq.take(int k)
                 |> Seq.skip 1 
                 |> iter 
                 |> Option.exists id 

///For Miller-Rabin the witnesses need to be selected at random from the interval [2, n - 2]. 
///More witnesses => better accuracy of the test.
///Also, remember that if Miller-Rabin returns true, then the number is _probable_ prime. 
///If it returns false the number is composite.
let isPrimeW witnesses = function
    | n when n < 2I -> false
    | n when n = 2I -> true
    | n when n = 3I -> true
    | n when n % 2I = 0I -> false
    | n             -> witnesses |> Seq.forall(millerRabinPrimality n)

///Example
let isPrime = isPrimeW [2I;3I]
namespace System
namespace System.Numerics
val pow' : mul:('a -> BigInteger -> BigInteger) -> sq:('a -> 'a) -> x':'a -> n':BigInteger -> BigInteger

Full name: Script.pow'


This implementation is based on the Miller-Rabin Haskell implementation
from http://www.haskell.org/haskellwiki/Testing_primality
val mul : ('a -> BigInteger -> BigInteger)
val sq : ('a -> 'a)
val x' : 'a
val n' : BigInteger
val f : ('a -> BigInteger -> BigInteger -> BigInteger)
val x : 'a
val n : BigInteger
val y : BigInteger
val q : BigInteger
val r : BigInteger
Multiple items
type BigInteger =
  struct
    new : value:int -> BigInteger + 7 overloads
    member CompareTo : other:int64 -> int + 3 overloads
    member Equals : obj:obj -> bool + 3 overloads
    member GetHashCode : unit -> int
    member IsEven : bool
    member IsOne : bool
    member IsPowerOfTwo : bool
    member IsZero : bool
    member Sign : int
    member ToByteArray : unit -> byte[]
    ...
  end

Full name: System.Numerics.BigInteger

--------------------
BigInteger()
BigInteger(value: int) : unit
BigInteger(value: uint32) : unit
BigInteger(value: int64) : unit
BigInteger(value: uint64) : unit
BigInteger(value: float32) : unit
BigInteger(value: float) : unit
BigInteger(value: decimal) : unit
BigInteger(value: byte []) : unit
BigInteger.DivRem(dividend: BigInteger, divisor: BigInteger, remainder: byref<BigInteger>) : BigInteger
val x2 : 'a
val mulMod : a:bigint -> b:BigInteger -> c:BigInteger -> BigInteger

Full name: Script.mulMod
val a : bigint
type bigint = BigInteger

Full name: Microsoft.FSharp.Core.bigint
val b : BigInteger
val c : BigInteger
val squareMod : a:bigint -> b:BigInteger -> BigInteger

Full name: Script.squareMod
val powMod : m:bigint -> (BigInteger -> BigInteger -> BigInteger)

Full name: Script.powMod
val m : bigint
val iterate : f:('a -> 'a) -> ('a -> seq<'a>)

Full name: Script.iterate
val f : ('a -> 'a)
module Seq

from Microsoft.FSharp.Collections
val unfold : generator:('State -> ('T * 'State) option) -> state:'State -> seq<'T>

Full name: Microsoft.FSharp.Collections.Seq.unfold
val fx : 'a
union case Option.Some: Value: 'T -> Option<'T>
val millerRabinPrimality : n:BigInteger -> a:BigInteger -> bool

Full name: Script.millerRabinPrimality


See: http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
val a : BigInteger
val find2km : (BigInteger -> BigInteger * BigInteger)
val f : (BigInteger -> BigInteger -> BigInteger * BigInteger)
val k : BigInteger
val m : BigInteger
val iter : (seq<BigInteger> -> bool option)
val tryPick : chooser:('T -> 'U option) -> source:seq<'T> -> 'U option

Full name: Microsoft.FSharp.Collections.Seq.tryPick
val x : BigInteger
union case Option.None: Option<'T>
val b0 : BigInteger
val failwith : message:string -> 'T

Full name: Microsoft.FSharp.Core.Operators.failwith
val sprintf : format:Printf.StringFormat<'T> -> 'T

Full name: Microsoft.FSharp.Core.ExtraTopLevelOperators.sprintf
val take : count:int -> source:seq<'T> -> seq<'T>

Full name: Microsoft.FSharp.Collections.Seq.take
Multiple items
val int : value:'T -> int (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.int

--------------------
type int = int32

Full name: Microsoft.FSharp.Core.int

--------------------
type int<'Measure> = int

Full name: Microsoft.FSharp.Core.int<_>
val skip : count:int -> source:seq<'T> -> seq<'T>

Full name: Microsoft.FSharp.Collections.Seq.skip
module Option

from Microsoft.FSharp.Core
val exists : predicate:('T -> bool) -> option:'T option -> bool

Full name: Microsoft.FSharp.Core.Option.exists
val id : x:'T -> 'T

Full name: Microsoft.FSharp.Core.Operators.id
val isPrimeW : witnesses:seq<BigInteger> -> _arg1:BigInteger -> bool

Full name: Script.isPrimeW


For Miller-Rabin the witnesses need to be selected at random from the interval [2, n - 2].
More witnesses => better accuracy of the test.
Also, remember that if Miller-Rabin returns true, then the number is _probable_ prime.
If it returns false the number is composite.
val witnesses : seq<BigInteger>
val forall : predicate:('T -> bool) -> source:seq<'T> -> bool

Full name: Microsoft.FSharp.Collections.Seq.forall
val isPrime : (BigInteger -> bool)

Full name: Script.isPrime


Example
Next Version Raw view Test code New version

More information

Link:http://fssnip.net/2w
Posted:13 years ago
Author:Cesar Mendoza
Tags: prime testing , primes

This web site is created using F# and Suave web server. It is hosted on Azure and the source code is on GitHub. Contributions are welcome!

The first version of fssnip.net has been created by @tomaspetricek back in 2010. This web site is a new, open-source and contribution-friendly version.