Quadratic Sieve

qsieve.js

// Writing this in JavaScript was a mistake.
const B = 10000;

const N = 9329629409n * 2315959301n;

console.time('exec');

const bigintSqrt = (v, ceil = false) => {
  if (v < 2n) return v;
  let iter = v >> 1n;
  while (true) {
    const next = (iter + (v / iter)) >> 1n;
    if (iter === next || next === iter + 1n) break;
    iter = next;
  }
  if (ceil && iter * iter < v) return iter + 1n;
  return iter;
};

const primeBasisB = [];
const composite = new Uint8Array(B + 6 >> 3);
const isComposite = i => (composite[i - 2 >> 3] >> ((i - 2) & 7)) & 1;
for (let i = 2; i <= B; ++i) {
  if (!isComposite(i)) {
    primeBasisB.push(BigInt(i));
    for (let j = i * 2; j < B; j += i) {
      composite[j - 2 >> 3] |= 1 << ((j - 2) & 7);
    }
  }
}

const gcd = (a, b) => b === 0n ? a : gcd(b, a % b);

const primeCountB = primeBasisB.length;
const aBase = bigintSqrt(N);
const aOffsets = new BigUint64Array(primeCountB + 1);

const rowBytes = primeCountB + 7 >> 3;
const basisCoeff = new Uint8Array(rowBytes);
const mat = new Uint8Array((primeCountB + 1) * rowBytes);
console.timeLog('exec', 'init');
for (let aOffset = 1n, i = 0; i <= primeCountB; aOffset += 1n) {
  let a = aOffset + aBase;
  let b = a * a - N;

  for (let k = 0; k < primeCountB && b > 1n; ++k) {
    const basisVal = primeBasisB[k];

    while (true) {
      const div = b / basisVal;
      if (div * basisVal !== b) break;
      basisCoeff[k >> 3] ^= 1 << (k & 7);
      b = div;
    }
  }

  if (aOffset % BigInt(B) === 0n) {
    console.timeLog('exec', 'scan', aOffset)
  }

  if (b === 1n) {
    if (i % Math.floor(primeCountB / 10) === 0) {
      console.timeLog('exec', 'base', i)
    }
    mat.set(basisCoeff, rowBytes * i);
    aOffsets[i++] = aOffset;
  }
  basisCoeff.fill(0);
}

console.timeLog('exec', 'bases');

const workingMat = mat.slice();
// augMat@mat == workingMat
// any rows of zeros in workingMat after row reduction correspond to rows
// in augMat where each value in the row acts like a coefficient to a row
// in mat. Basically the row in augMat is in the left-nullspace of mat.
const augRowBytes = primeCountB + 8 >> 3;
const augMat = new Uint8Array((primeCountB + 1) * augRowBytes);
for (let i = 0; i <= primeCountB; ++i) {
  augMat[i * augRowBytes + (i >> 3)] |= 1 << (i & 7);
}

const rowXOR = (i1, i2) => {
  const r1 = i1 * rowBytes, r2 = i2 * rowBytes;
  const ar1 = i1 * augRowBytes, ar2 = i2 * augRowBytes;
  for (let i = 0; i < rowBytes; ++i) {
    workingMat[r1 + i] ^= workingMat[r2 + i];
    augMat[ar1 + i] ^= augMat[ar2 + i];
  }
}

let curRow = 0;
for (let pivot = 0; pivot < primeCountB; ++pivot) {
  let leadingRow = curRow;
  for (; leadingRow <= primeCountB; ++leadingRow) {
    if ((workingMat[leadingRow * rowBytes + (pivot >> 3)] >> (pivot & 7)) & 1) {
      break;
    }
  }
  if (leadingRow > primeCountB) {
    continue;
  }
  if (leadingRow !== curRow) rowXOR(curRow, leadingRow);
  for (let row = curRow + 1; row <= primeCountB; ++row) {
    if ((workingMat[row * rowBytes + (pivot >> 3)] >> (pivot & 7)) & 1) {
      rowXOR(row, curRow);
    }
  }
  ++curRow;
}

console.timeLog('exec', 'rref');

for (let i = curRow; i <= primeCountB; ++i) {
  let x = 1n;
  let ySq = 1n;
  for (let j = 0; j <= primeCountB; ++j) {
    if ((augMat[i * augRowBytes + (j >> 3)] >> (j & 7)) & 1) {
      const a = aBase + aOffsets[j];
      x = (x * a) % N;
      ySq = (ySq * (a * a - N));
    }
  }
  const y = bigintSqrt(ySq) % N;
  const factor = gcd(x + y, N);
  if (factor !== 1n && factor !== N) {
    console.log('Found factor', factor, '|', N);
    console.log('rest =', N / factor);
    break;
  }
}

Apparently this is slow not because I wrote this in JavaScript, but rather because I'm picking my as rather stupidly. The reason they call it a "quadratic sieve" is because you're supposed to pick both your prime factors by solving some strange modular quadratic congruence. (It's probably also JavaScript though.)