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));