Find perfect numbers

perfect-numbers.js

/**
 * 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

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:

• draw a random number
• compute the Jacobi Symbol
• compute the power of a large number to a large number

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:

19th number is found in ~35s.