Skip to content

Instantly share code, notes, and snippets.

@lancejpollard
Created January 12, 2022 15:35
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save lancejpollard/8f0788feea5e74b9974c7d4a9a8039a6 to your computer and use it in GitHub Desktop.
Save lancejpollard/8f0788feea5e74b9974c7d4a9a8039a6 to your computer and use it in GitHub Desktop.
AKS + Miller Prime Test Combination in JavaScript
// function binomialCoeff(n, k) {
// if ((k === 0) || (k === n)) {
// return 1n
// } else {
// return binomialCoeff(n - 1n, k - 1n) + binomialCoeff(n - 1n, k)
// }
// }
// function isDivisible(value) {
// return (value % this == 0n)
// }
// function isPrime(p) {
// if ((p === 0n) || (p === 1n)) {
// return true
// } else {
// let coeff = Array.from(Array(p + 1n), (x, index) => binomialCoeff(p, index))
// coeff.shift()
// coeff.splice(-1, 1)
// return coeff.every(isDivisible, p)
// }
// }
// array used to store coefficients .
let c = [];
// function to calculate the coefficients
// of (x - 1)^n - (x^n - 1) with the help
// of Pascal's triangle .
function coef(n)
{
c[0] = 1n;
for (let i = 0n; i < n; c[0] = -c[0], i++) {
c[1n + i] = 1n;
for (let j = i; j > 0n; j--)
c[j] = c[j - 1n] - c[j];
}
}
// function to check whether
// the number is prime or not
function isPrimeAKS(n)
{
// Calculating all the coefficients by
// the function coef and storing all
// the coefficients in c array .
coef(n);
// subtracting c[n] and adding c[0] by 1
// as ( x - 1 )^n - ( x^n - 1), here we
// are subtracting c[n] by 1 and adding
// 1 in expression.
c[0]++;
c[n]--;
// checking all the coefficients whether
// they are divisible by n or not.
// if n is not prime, then loop breaks
// and (i > 0).
let i = n;
while ((i--) > 0n && c[i] % n == 0n)
;
// Return true if all coefficients are
// divisible by n.
return i < 0;
}
function isPrime(n) {
if (isPrimeMiller(n)) {
return isPrimeAKS(n)
}
return false
}
module.exports = isPrime
// Javascript program Miller-Rabin primality test
// based on JavaScript code found at https://www.geeksforgeeks.org/primality-test-set-3-miller-rabin/
// Utility function to do
// modular exponentiation.
// It returns (x^y) % p
function power(x, y, p)
{
// Initialize result
// (JML- all literal integers converted to use n suffix denoting BigInt)
let res = 1n;
// Update x if it is more than or
// equal to p
x = x % p;
while (y > 0n)
{
// If y is odd, multiply
// x with result
if (y & 1n)
res = (res*x) % p;
// y must be even now
y = y/2n; // (JML- original code used a shift operator, but division is clearer)
x = (x*x) % p;
}
return res;
}
// This function is called
// for all k trials. It returns
// false if n is composite and
// returns false if n is
// probably prime. d is an odd
// number such that d*2<sup>r</sup> = n-1
// for some r >= 1
function millerTest(d, n)
{
// (JML- all literal integers converted to use n suffix denoting BigInt)
// Pick a random number in [2..n-2]
// Corner cases make sure that n > 4
/*
JML- I can't mix the Number returned by Math.random with
operations involving BigInt. The workaround is to create a random integer
with precision 6 and convert it to a BigInt.
*/
const r = BigInt(Math.floor(Math.random() * 100_000))
// JML- now I have to divide by the multiplier used above (BigInt version)
const y = r*(n-2n)/100_000n
let a = 2n + y % (n - 4n);
// Compute a^d % n
let x = power(a, d, n);
if (x == 1n || x == n-1n)
return true;
// Keep squaring x while one
// of the following doesn't
// happen
// (i) d does not reach n-1
// (ii) (x^2) % n is not 1
// (iii) (x^2) % n is not n-1
while (d != n-1n)
{
x = (x * x) % n;
d *= 2n;
if (x == 1n)
return false;
if (x == n-1n)
return true;
}
// Return composite
return false;
}
// It returns false if n is
// composite and returns true if n
// is probably prime. k is an
// input parameter that determines
// accuracy level. Higher value of
// k indicates more accuracy.
function isPrimeMiller( n, k=40)
{
// (JML- all literal integers converted to use n suffix denoting BigInt)
// Corner cases
if (n <= 1n || n == 4n) return false;
if (n <= 3n) return true;
// Find r such that n =
// 2^d * r + 1 for some r >= 1
let d = n - 1n;
while (d % 2n == 0n)
d /= 2n;
// Iterate given nber of 'k' times
for (let i = 0; i < k; i++)
if (!millerTest(d, n))
return false;
return true;
}
module.exports = isPrime
// // small
// logBetween(0n, 10n)
// // normal but bigger than 65536
// logBetween(0n, 100000n)
// // pretty big but still within JS 2 ** 64
// logBetween(0n, 100000100000100000100000n)
// // big, outside of 2 ** 64
// logBetween(0n, 100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000n)
// // huge, way outside 2 ** 64
// logBetween(0n, 100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000n)
// // between normal and pretty big
// logBetween(100000n, 100000100000100000100000n)
// // between pretty big and big
// logBetween(100000100000100000100000n, 100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000100000n)
// function logBetween(a, b) {
// console.log(randomBigIntBetween(a, b), `${a} <=> ${b}`)
// }
module.exports = randomBigIntBetween
function randomBigIntBetween(minInclusive, maxExclusive) {
var maxInclusive = (maxExclusive - minInclusive) - BigInt(1)
var x = BigInt(1)
var y = BigInt(0)
while(true) {
x = x * BigInt(2)
var randomBit = BigInt(Math.random()<0.5 ? 1 : 0)
y = y * BigInt(2) + randomBit
if(x > maxInclusive) {
if (y <= maxInclusive) { return y + minInclusive }
// Rejection
x = x - maxInclusive - BigInt(1)
y = y - maxInclusive - BigInt(1)
}
}
}
const fs = require('fs')
const isPrime = require('./findPrime')
const rbigint = require('./rbigint')
const stream = fs.createWriteStream('tmp/prime.large.csv', { flags: 'a+' })
const wait = () => new Promise((res) => stream.once('drain', res))
start()
async function start() {
stream.write(`prime,\n`)
let i = 0
while (true) {
let n = rbigint(
1000000000100000000010000000001000000000100000000010000000001000000000n,
10000000001000000000100000000010000000001000000000100000000010000000001000000000n
)
if (isPrime(n) && n % 4n === 3n) {
console.log(String(n))
const canContinue = stream.write(`${n},\n`)
if (!canContinue) {
await wait();
}
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment