Finding the 4 squares that sums up to a given number

Question

I wrote this algorithm in JS with the main goal that it will give an answer for 10^750-10^1000 in 2-3 seconds but it solves 10^150 in 2-3 seconds so I am a little far from my goal, the problem is that I sieved through the algorithm couple of times, but nothing I do seems to speed the algorithm significantly.

I would appreciate if any one could enhance my code or give a new one in JS to reach my goal of solving 10^750-10^1000 in under 3 seconds

I implemented Michael O. Rabin and Jeffrey Shallit algorithm for finding the 4 squares summing up to given number n and Miller-Rabin for primality test, I think that these are the very best and fast, so why is my algorithm is not fast as "theoretically could be" ?!

Edit : I uploaded first to stack overflow but was advised to come to this site.

Edit : after a little work i managed to make it run in 0.4 second for numbers like 10^300 but i still get 10.6second for 10^750 which is my goal (but now the code is much faster), i also fragmented the code and debugged it to see what process is time consuming and its in the fourSquares where i try to find a prime that is -1 of a multiple of my number n.

 

let primes = [2n, 3n, 5n, 7n];

function power(x, y, p) {
    let res = 1n;
    x %= p;
    while (y > 0n) {
        if (y % 2n == 1n)
            res = (res * x) % p;
        x = (x * x) % p;
        y /= 2n;
    }
    return res;
}

function millerTest(a, d, n) {
    let x = power(a, d, n);
    if (x == 1n || x == n - 1n)
        return true;
    while (d != n - 1n) {
        x = (x * x) % n;
        d *= 2n;
        if (x == 1n)
            return false;
        if (x == n - 1n)
            return true;
    }
    return false;
}

function isPrime(n) {
    if (n <= 1n || n == 4n) return false;
    for (let i = 0; i < primes.length; i++) if (n == primes[i]) return true;
    let d = n - 1n;
    while (d % 2n == 0)
        d /= 2n;
    for (let i = 0; i < primes.length; i++)
        if (!millerTest(primes[i], d, n))
            return false;
    return true;
}

function abs(a) {
    if (a < 0n) return -a;
    return a;
}

function round(a, b) {
    if (b == 0n) return undefined;
    if (a % b == 0n) return a / b;
    let d = a / b;
    let x = abs(a - (d - 1n) * b);
    let y = abs(a - d * b);
    let z = abs(a - (d + 1n) * b);
    if (x <= y && x <= z) return d - 1n;
    if (z <= y && z <= x) return d + 1n;
    return d;
}

function complexIsZero(c) {
    if (c[0] == 0n && c[1] == 0n) return true;
    return false;
}

function complexNorm(c) {
    return c[0] * c[0] + c[1] * c[1];
}

function complexGCD(c1, c2) {
    if (complexIsZero(c1)) return c2;
    if (complexIsZero(c2)) return c1;
    let u = c1[0];
    let v = c1[1];
    let x = c2[0];
    let y = c2[1];
    let a = round(u * x + v * y, x * x + y * y);
    let b = round(v * x - u * y, x * x + y * y);
    while (!complexIsZero(c2)) {
        c1 = c2;
        c2 = [u - a * x + b * y, v - b * x - a * y];
        u = c1[0];
        v = c1[1];
        x = c2[0];
        y = c2[1];
        a = round(u * x + v * y, x * x + y * y);
        b = round(v * x - u * y, x * x + y * y);
    }
    return c1;
}

function twoSquares(p) {
    if (isPrime(p) && p % 4n == 1n) {
        for (let i = 0; i < primes.length; i++)
            if (power(primes[i], (p - 1n) / 2n, p) == p - 1n) {
                let x = power(primes[i], (p - 1n) / 4n, p);
                return complexGCD([x, 1n], [p, 0n]);
            }
    }
    return [0n, 0n];
}


//Quaternion Class
function QIsZero(q1) {
    for (let i = 0; i < q1.length; i++) q1[i] = BigInt(q1[i]);
    for (let i = 0; i < q1.length; i++) if (q1[i] != 0n) return false;
    return true;
}

function Qadd(q1, q2) {
    for (let i = 0; i < q1.length; i++) q1[i] = BigInt(q1[i]);
    for (let i = 0; i < q2.length; i++) q2[i] = BigInt(q2[i]);
    let a1 = q1[0];
    let b1 = q1[1];
    let c1 = q1[2];
    let d1 = q1[3];
    let a2 = q2[0];
    let b2 = q2[1];
    let c2 = q2[2];
    let d2 = q2[3];
    return [a1 + a2, b1 + b2, c1 + c2, d1 + d2];
}

function Qsub(q1, q2) {
    for (let i = 0; i < q1.length; i++) q1[i] = BigInt(q1[i]);
    for (let i = 0; i < q2.length; i++) q2[i] = BigInt(q2[i]);
    let a1 = q1[0];
    let b1 = q1[1];
    let c1 = q1[2];
    let d1 = q1[3];
    let a2 = q2[0];
    let b2 = q2[1];
    let c2 = q2[2];
    let d2 = q2[3];
    return [a1 - a2, b1 - b2, c1 - c2, d1 - d2];
}

function Qconj(q1) {
    for (let i = 0; i < q1.length; i++) q1[i] = BigInt(q1[i]);
    let a1 = q1[0];
    let b1 = q1[1];
    let c1 = q1[2];
    let d1 = q1[3];
    return [a1, -b1, -c1, -d1];
}

function Qmul(q1, q2) {
    for (let i = 0; i < q1.length; i++) q1[i] = BigInt(q1[i]);
    for (let i = 0; i < q2.length; i++) q2[i] = BigInt(q2[i]);
    let a1 = q1[0];
    let b1 = q1[1];
    let c1 = q1[2];
    let d1 = q1[3];
    let a2 = q2[0];
    let b2 = q2[1];
    let c2 = q2[2];
    let d2 = q2[3];
    let x = a1 * a2 - b1 * b2 - c1 * c2 - d1 * d2;
    let y = a1 * b2 + b1 * a2 + c1 * d2 - d1 * c2;
    let z = a1 * c2 - b1 * d2 + c1 * a2 + d1 * b2;
    let w = a1 * d2 + b1 * c2 - c1 * b2 + d1 * a2;
    return [x, y, z, w];
}

function Qnorm(q1) {
    for (let i = 0; i < q1.length; i++) q1[i] = BigInt(q1[i]);
    let a1 = q1[0];
    let b1 = q1[1];
    let c1 = q1[2];
    let d1 = q1[3];
    return a1 * a1 + b1 * b1 + c1 * c1 + d1 * d1;
}

function Qdiv(q, r) {
    let q0 = q[0];
    let q1 = q[1];
    let q2 = q[2];
    let q3 = q[3];
    let r0 = r[0];
    let r1 = r[1];
    let r2 = r[2];
    let r3 = r[3];
    let t0 = r0 * q0 + r1 * q1 + r2 * q2 + r3 * q3;
    let t1 = r0 * q1 - r1 * q0 - r2 * q3 + r3 * q2;
    let t2 = r0 * q2 + r1 * q3 - r2 * q0 - r3 * q1;
    let t3 = r0 * q3 - r1 * q2 + r2 * q1 - r3 * q0;
    let norm = Qnorm(r);
    let a = round(t0, norm);
    let b = round(t1, norm);
    let c = round(t2, norm);
    let d = round(t3, norm);
    return [a, b, c, d];
}

function gcdQuaternion(a, b) {
    if (QIsZero(a)) return b;
    if (QIsZero(b)) return a;
    let c, t;
    while (!QIsZero(b)) {
        c = Qdiv(a, b);
        t = Qsub(a, Qmul(b, c));
        a = b;
        b = t;
    }
    return a;
}

function fourSquares(n) {
    if (n < 0n) return undefined;
    if (n == 0n) return [0n, 0n, 0n, 0n];
    let d = 0n;
    while (n % 2n == 0) {
        n /= 2n;
        d += 1n;
    }
    let e = (-4n) ** (d / 4n);
    let f = d % 4n;
    let M = 30n;
    let k = 1n;
    let p = M * n * k - 1n;
    let add = 2n * M * n;
    while (!isPrime(p))
        p += add;
    let AB = twoSquares(p);
    let A = AB[0];
    let B = AB[1];
    let Q = gcdQuaternion([A, B, 1n, 0n], [n, 0n, 0n, 0n]);
    for (let i = 0; i < Q.length; i++) Q[i] *= e;
    if (f == 1n) Q = Qmul([1n, 1n, 0n, 0n], Q);
    if (f == 2n) Q = Qmul([0n, 2n, 0n, 0n], Q);
    if (f == 3n) Q = Qmul([-2n, 2n, 0n, 0n], Q);
    return Q;
}

console.log(fourSquares(10n ** 300n));

Answer 1

Avoid linear scans

    for (let i = 0; i < primes.length; i++) if (n == primes[i]) return true;

So if n is 6, we should compare it to every value in primes just to find out that it's not one of them? At minimum, maintain primes in sorted order and do a binary search. And it's quite possible that you should keep primes in a set that can answer this question directly.

There's only one even prime

function nextPrime(n) {
    n++;
    while (!isPrimeSimple(n)) n++;
    return n;
}

So if n is 1, we should check 2 (OK), 3 (OK), 4? Why 4, 6, 8, etc.? They will never be prime. Consider

function nextPrime(n) {
    if (n < 2) {
        return 2;
    }

    // if even, add one to make odd
    // if odd, add two to get the next odd
    n += 1 + (n % 2);

    while (!isPrimeSimple(n)) {
        n += 2;
    }

    return n;
}

The parentheses and curly brackets are not necessary, but I prefer them. You can leave out the initial check if you only call this with values of n at least 2.

This only increments half as many times.

You don't need to multiply

    if (n == 2 || n == 3) return true;
    for (let i = 2; i * i <= n; i++) if (n % i == 0) return false;

Well, first, there's no need to increment by one or check three separately.

    if (n % 2 === 0) {
        return n === 2;
    }

    for (let i = 3; i <= n / i; i += 2) {
        if (n % i === 0) {
            return false;
        }
    }

Three will correctly return true after this code.

We don't check any even potential factors other than two itself, as we would already have returned after checking two.

Now, in terms of multiplying versus division, you might ask why it's easier to calculate n / i than i * i. In general, it isn't. But in this specific case you are calculating n % i. And calculating n / i and n % i are often the same operation. So you might as well use the n / i that you probably already have.

It's also possible that your interpreter/compiler won't perform that optimization. If not and you can't switch to a better compiler, calculating the square root of n once would probably be cheaper than repeatedly squaring values of i.

Profiling

These are just things that jumped out at me. A better way of proceeding would be to run the code in a profiling debugger and have it tell you where you are spending most of your time. Is it in isSimplePrime? isPrime? Perhaps the BigInt constructor is slowing you down. Find out where your code is slow and then look for ways to optimize that.

Some common advice would be to generate nextPrime using a sieve, e.g. Sieve of Eratosthenes. Then you could toss isSimplePrime entirely.

Compare your Miller test to a straight up test. Or to one of the alternatives.