Last active
November 21, 2018 08:46
-
-
Save GMartigny/632e643114c5d6d6ba3c806695b801c8 to your computer and use it in GitHub Desktop.
Find perfect numbers
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/** | |
* Do the work | |
* @param {Object} parameters - Selected value on inputs | |
*/ | |
function work (parameters) { | |
const limit = +parameters.limit; | |
/** | |
* Compute that weird jacobi thing | |
* from: http://2000clicks.com/mathhelp/NumberTh27JacobiSymbolAlgorithm.aspx | |
*/ | |
const jacobi = (a, b) => { | |
if (b % 2n == 0) | |
return 0; | |
let j = 1n; | |
while (a) { | |
while (a % 2n == 0) { | |
a >>= 1n; | |
if (b % 8n == 3 || b % 8n == 5) | |
j *= -1n; | |
} | |
[a, b] = [b, a]; | |
if (a % 4n == 3 && b % 4n == 3) | |
j *= -1n; | |
a %= b; | |
} | |
if (b == 1) | |
return j; | |
return 0; | |
}; | |
/** | |
* Compute a power between large number without braking the memory | |
* from: https://www.wikiwand.com/en/Modular_exponentiation#/Left-to-right_binary_method | |
*/ | |
const modularPow = (base, exp, mod) => { | |
let res = 1n; | |
while (exp) { | |
res **= 2n; | |
if (exp % 2n == 1) | |
res = (res * base) % mod; | |
exp >>= 1n; | |
} | |
return res; | |
}; | |
const r = Math.random; | |
/** | |
* Return random number up to a BigInt upper bound (this is not correct randomness, but good enough) | |
*/ | |
const random = (to) => { | |
const str = to.toString(); | |
return BigInt(r().toString().substr(2, str.length - 1)) * BigInt(str[0]); | |
}; | |
const primes = []; | |
/** | |
* Check if integer is a prime number | |
*/ | |
const isPrime = (n) => { | |
if (n < 4) return n > 1; | |
// Check against already known primes | |
if (primes.some(p => (n % p) == 0)) return false; | |
// from: https://www.wikiwand.com/en/Solovay%E2%80%93Strassen_primality_test#/Algorithm_and_running_time | |
let i = 15; // 0.003% chance of false positive | |
while (i--) { | |
const a = random(n - 4n) + 2n; | |
const x = jacobi(a, n); | |
if (x == 0 || modularPow(a, (n - 1n) / 2n, n) != (x + n) % n) | |
return false; | |
} | |
return true; | |
}; | |
/** | |
* Return a generator giving only prime numbers | |
* from: https://www.wikiwand.com/en/Sieve_of_Eratosthenes | |
*/ | |
function* primeGenerator (l) { | |
primes.push(2n); | |
yield 2n; | |
let possibilities = (new Array(l / 2)).fill().map((_, i) => 3n + (2n * BigInt(i))); | |
while (possibilities.length) { | |
const prime = possibilities[0]; | |
possibilities = possibilities.filter(x => x % prime); | |
primes.push(prime); | |
progression(1 - (possibilities.length / l)); | |
yield prime; | |
} | |
} | |
/** | |
* Compute the number of form b(3) => b11100, b(n) => b1(n time)0(n-1 time) | |
*/ | |
const buildPerfect = (n) => ((2n ** n) - 1n) * (2n ** (n - 1n)); | |
const generator = primeGenerator(limit); | |
let prime = generator.next(); | |
const results = []; | |
const exp = []; | |
while (!prime.done) { | |
if (isPrime((2n ** prime.value) - 1n)) { | |
exp.push(prime.value); | |
results.push(buildPerfect(prime.value)); | |
} | |
prime = generator.next(); | |
} | |
send(`Found ${results.length} perfect numbers.`); | |
send("Exponents: " + exp.join(", ")); | |
return results.map((n, i) => `${i + 1}. ${n}`); | |
} | |
work({limit: 2300}); // Will go up to 2300, giving the 17nth result |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
This script output Perfect Numbers after having checked all prime up to a limit passed in parameter.
Perfect Numbers tend to grow fast, the 8th is 19 digits long and the 13th is 314 digits long. Javascript can't natively handle such large number without going overflow. Hopefully, using BigInt allow you to manipulate arbitrary large number. Only problem, they are not compatible with native number. All calculations need to be done using BigInt.
This is based on the assertion that in binary, perfect numbers always have the same form:
1(n time)0(n-1 time)
.Thanks to OEIS, I found out that
n
follow the Mersenne exponents list. A Mersenne exponents is a prime p such that 2^p - 1 is prime.First task is to generate prime numbers, I use the Sieve of Eratosthenes for its simplicity and speed.
Then check if this prime is still prime after
2^p-1
. This check can be done with the Solovay–Strassen primality test. This method, however, has 3 caveats:The first one is tricky without floating numbers, I came out with a "good enough" solution.
The jacobi symbol can be solved with any algorithm found online.
Finally, we can use Modular exponentiation to compute powers without overflowing the memory. This last point is the main performance bottleneck.
Results using my MathWorker on a not so performant computer: