diff --git a/generators.go b/generators.go
new file mode 100644
index 0000000..e3d2777
--- /dev/null
+++ b/generators.go
@@ -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
+}
diff --git a/miller-rabin.go b/miller-rabin.go
index 8264c15..738dcc1 100644
--- a/miller-rabin.go
+++ b/miller-rabin.go
@@ -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.
diff --git a/utils.go b/utils.go
new file mode 100644
index 0000000..42840ef
--- /dev/null
+++ b/utils.go
@@ -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
+}