Skip to content

Commit

Permalink
chore: reorganise Miller-Rabin code into seperate files
Browse files Browse the repository at this point in the history
  • Loading branch information
h5law committed May 18, 2024
1 parent e8d8378 commit ab40861
Show file tree
Hide file tree
Showing 3 changed files with 108 additions and 100 deletions.
30 changes: 30 additions & 0 deletions generators.go
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
package primality

// SieveOfEratosthenes is an implementation of the ancient sieve of Eratosthenes
// to generate a list of primes up to the provided integer.
func SieveOfEratosthenes(n uint64) []byte {
if n < 2 {
return nil // 1 is not prime
}
// Initialise a boolean slice with all false values
b := make([]bool, n+1)
b[0], b[1] = true, true
// Set all multiples of numbers up to sqrt(n) to true
for i := uint64(2); i*i <= n; i++ {
if b[i] {
continue
}
for j := i * i; j <= n; j += i {
b[j] = true
}
}
// Initialise a slice to store the primes we generate
primes := make([]byte, n+1)
// Primes are remaining false values indexes
for i, v := range b {
if !v {
primes[i] = 0x1
}
}
return primes
}
100 changes: 0 additions & 100 deletions miller-rabin.go
Original file line number Diff line number Diff line change
Expand Up @@ -3,113 +3,13 @@ package primality
import (
"math/big"
"math/rand"
"slices"
"time"
)

// SieveOfEratosthenes is an implementation of the ancient sieve of Eratosthenes
// to generate a list of primes up to the provided integer.
func SieveOfEratosthenes(n uint64) []uint64 {
if n < 2 {
return nil // 1 is not prime
}
// Initialise a slice to store the primes we generate
primes := make([]uint64, 0, n)
// Initialise a boolean slice with all false values
b := make([]bool, n+1)
b[0], b[1] = true, true
// Set all multiples of numbers up to sqrt(n) to true
for i := uint64(2); i*i <= n; i++ {
if b[i] {
continue
}
for j := i * i; j <= n; j += i {
b[j] = true
}
}
// Primes are remaining false values indexes
for i, v := range b {
if !v {
primes = append(primes, uint64(i))
}
}
return primes
}

var zero = big.NewInt(0)
var one = big.NewInt(1)
var two = big.NewInt(2)

// Return new slice containing only the elements in a that are not in b.
//
// d = {x : x ∈ a and x ∉ b}
func difference(a, b []uint64) []uint64 {
diff := make([]uint64, 0, len(a))
for _, v := range a {
if found := slices.Contains(b, v); !found {
diff = append(diff, v)
}
}
return diff
}

// Generate a bitmask for primes from lo->hi, the bitmask will index the primes
// at their relative bit indexes for quick comparison and evaluation of primes.
func primemask(lo, hi uint64) *big.Int {
low := make([]uint64, lo)
high := make([]uint64, hi)
// Get primes 2->hi
high = SieveOfEratosthenes(hi)
// Only get lower limit primes if lo > 2
if lo > 2 {
low = SieveOfEratosthenes(lo)
}
// Use only difference between high and low slices
diff := difference(high, low)
// Create the mask by performing mask |= 1<<prime
mask := new(big.Int)
one := big.NewInt(1)
for i := 0; i < len(diff); i++ {
prime := new(big.Int)
// Left shift
prime = prime.Lsh(one, uint(diff[i]))
mask = mask.Or(mask, prime)
}
return mask
}

// bigPowMod uses the binary exponentiation of the exponent provided to
// compute the base raised to the exponent mod the provided modulus.
func bigPowMod(b, e, m *big.Int) *big.Int {
// Special case
if m.Cmp(one) == 0 {
return zero
}
// Copy arguments
base := new(big.Int)
base = base.Set(b)
exp := new(big.Int)
exp = exp.Set(e)
mod := new(big.Int)
mod = mod.Set(m)
// Initialise remainder
r := big.NewInt(int64(1))
// Use the binary exponentiation of e to successively calculate b^e (mod m)
base = base.Mod(base, mod)
em := new(big.Int)
for exp.Cmp(zero) > 0 {
em = em.Mod(exp, two)
if em.Cmp(one) == 0 {
r = r.Mul(r, base).Mod(r, mod)
}
// Divide exponent by 2
exp = exp.Rsh(exp, uint(1))
// Repeatedly square b to achieve b=b^(2^i) (mod m)
base = base.Mul(base, base).Mod(base, mod)
}
return r
}

// MillerRabin is an implementation of the Miller-Rabin primality test, it is
// a probabilistic method - and as such using 25 repetitions/rounds is
// recommended as it ensures a higher probability of accuracy of the result.
Expand Down
78 changes: 78 additions & 0 deletions utils.go
Original file line number Diff line number Diff line change
@@ -0,0 +1,78 @@
package primality

import (
"bytes"
"math/big"
)

// Return new slice containing only the elements in a that are not in b.
//
// d = {x : x ∈ a and x ∉ b}
func difference(a, b []byte) []byte {
diff := make([]byte, 0, len(a))
for _, v := range a {
if found := bytes.Contains(b, []byte{v}); !found {
diff = append(diff, v)
}
}
return diff
}

// Generate a bitmask for primes from lo->hi, the bitmask will index the primes
// at their relative bit indexes for quick comparison and evaluation of primes.
func primemask(lo, hi uint64) *big.Int {
low := make([]byte, lo)
high := make([]byte, hi)
// Get primes 2->hi
high = SieveOfEratosthenes(hi)
// Only get lower limit primes if lo > 2
if lo > 2 {
low = SieveOfEratosthenes(lo)
}
// Use only difference between high and low slices
diff := difference(high, low)
// Create the mask by performing mask |= 1<<prime
mask := new(big.Int)
one := big.NewInt(1)
for i := 0; i < len(diff); i++ {
prime := new(big.Int)
// Left shift
if diff[i] == 0x1 {
prime = prime.Lsh(one, uint(i))
mask = mask.Or(mask, prime)
}
}
return mask
}

// bigPowMod uses the binary exponentiation of the exponent provided to
// compute the base raised to the exponent mod the provided modulus.
func bigPowMod(b, e, m *big.Int) *big.Int {
// Special case
if m.Cmp(one) == 0 {
return zero
}
// Copy arguments
base := new(big.Int)
base = base.Set(b)
exp := new(big.Int)
exp = exp.Set(e)
mod := new(big.Int)
mod = mod.Set(m)
// Initialise remainder
r := big.NewInt(int64(1))
// Use the binary exponentiation of e to successively calculate b^e (mod m)
base = base.Mod(base, mod)
em := new(big.Int)
for exp.Cmp(zero) > 0 {
em = em.Mod(exp, two)
if em.Cmp(one) == 0 {
r = r.Mul(r, base).Mod(r, mod)
}
// Divide exponent by 2
exp = exp.Rsh(exp, uint(1))
// Repeatedly square b to achieve b=b^(2^i) (mod m)
base = base.Mul(base, base).Mod(base, mod)
}
return r
}

0 comments on commit ab40861

Please sign in to comment.