Introduction

Every positive integer can be written as a product of primes, and finding that product is called integer factoring. It sounds elementary — yet no one knows a fast general algorithm, and the hardest instances safeguard modern encryption.

For most of history the best strategy was trial division: divide by 2, then 3, then 5, and so on. That works for small numbers, but for a 100-digit number it would take longer than the age of the universe.

In 1981, Carl Pomerance invented the Quadratic Sieve (QS) — an algorithm that runs in sub-exponential time, far faster than trial division, yet (so far as anyone knows) still not polynomial. It held the world speed record for factoring large integers through the 1980s and into the 1990s, and it still reigns for numbers up to roughly 100 digits.

The key insight is elegant: instead of dividing, hunt for pairs of numbers whose squares are congruent modulo NN — and build a system of equations over F2\mathbb{F}_2 that, when solved, hands you a factor almost for free. See how factoring relates to public-key cryptography, and how P vs NP looms in the background.

Build the Relation Matrix

The demo sieves values xx near ⌊N⌋\lfloor\sqrt{N}\rfloor, tests whether x2 mod Nx^2 \bmod N factors completely over a small factor base of primes, and records each successful factorization as a row in a binary matrix. When enough smooth relations are found it runs Gaussian elimination over F2\mathbb{F}_2 to find a dependency — and that dependency yields a non-trivial factor of NN.

<p class="hint">Pick a number <b>N</b> to factor. The sieve tests values of <em>x</em>
near √N, checks whether <em>x</em>² mod N factors completely over the factor base,
and collects those <b>smooth relations</b> into a binary matrix. Once enough rows
accumulate, Gaussian elimination over 𝔽₂ finds a dependency that yields a factor.</p>
<div class="controls">
  <label>N = <select id="nsel">
    <option value="77">77 = ?</option>
    <option value="143">143 = ?</option>
    <option value="323" selected>323 = ?</option>
    <option value="667">667 = ?</option>
    <option value="1147">1147 = ?</option>
  </select></label>
  <button id="sieve" type="button">Run Sieve</button>
  <button id="solve" type="button" disabled>Solve Matrix</button>
  <button id="reset" type="button" class="ghost">Reset</button>
</div>
<div id="fb-row" class="fb-row"></div>
<div id="matrix-wrap" class="matrix-wrap"></div>
<div class="status" id="status">Choose N and click Run Sieve.</div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; font-size: 14px; }
.hint { font-size: .88rem; color: #444; margin: 0 0 .7rem; line-height: 1.45; }
.controls { display: flex; gap: .5rem; align-items: center; flex-wrap: wrap; margin-bottom: .6rem; }
label { font-size: .9rem; }
select { font-size: .9rem; padding: .2rem .4rem; border-radius: 6px; border: 1px solid #adb1b8; }
button { font: 600 13px system-ui; padding: .4rem .8rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .45; cursor: default; }
.fb-row { font-size: .8rem; color: #555; margin-bottom: .4rem; min-height: 1.2em; }
.matrix-wrap { overflow-x: auto; margin-bottom: .5rem; }
table { border-collapse: collapse; font-size: .82rem; font-family: ui-monospace, monospace; }
th { background: #1d3557; color: #fff; padding: 3px 7px; text-align: center; white-space: nowrap; }
td { padding: 3px 7px; text-align: center; border-bottom: 1px solid #e0e4e8; }
tr:nth-child(even) td { background: #f5f7f9; }
td.one { color: #e63946; font-weight: 700; }
td.dep { background: #fff3cd !important; }
.status { font-size: .95rem; font-weight: 600; margin-top: .3rem; min-height: 1.4em; }
.status.ok { color: #0a7d33; }
.status.bad { color: #c92f3c; }
.status.info { color: #1d3557; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
/* Quadratic Sieve demo: collect smooth relations, build the binary exponent
matrix (mod 2), then use Gaussian elimination over F_2 to find a linear
dependency. That dependency gives x^2 ≡ y^2 (mod N) => factor via gcd. */

function gcd(a, b) { while (b) { [a, b] = [b, a % b]; } return a; }

// Small primes up to 30 as factor base
const PRIMES = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29];

// Trial-divide v completely over factor base; return exponent vector or null
function smoothDecompose(v, primes) {
  if (v <= 0) return null;
  const exps = new Array(primes.length).fill(0);
  let rem = v;
  for (let i = 0; i < primes.length; i++) {
    while (rem % primes[i] === 0) { rem = rem / primes[i]; exps[i]++; }
  }
  return rem === 1 ? exps : null;
}

// Gaussian elimination over F_2; returns list of dependency rows (indices)
function f2Elimination(matrix) {
  const rows = matrix.map(r => r.slice());
  const n = rows.length, m = rows[0].length;
  const pivot = new Array(m).fill(-1);
  const rowOf = []; // rowOf[col] = pivot row index
  for (let col = 0; col < m; col++) {
    let found = -1;
    for (let row = 0; row < n; row++) {
      if (rows[row][col] === 1 && pivot.indexOf(row) === -1 && rowOf.indexOf(row) === -1) {
        found = row; break;
      }
    }
    if (found === -1) continue;
    rowOf[col] = found;
    pivot[found] = col;
    for (let row = 0; row < n; row++) {
      if (row !== found && rows[row][col] === 1) {
        for (let c = 0; c < m; c++) rows[row][c] ^= rows[found][c];
      }
    }
  }
  // Rows not chosen as pivots are linear dependencies
  const deps = [];
  for (let row = 0; row < n; row++) {
    if (pivot.indexOf(row) === -1) deps.push(row);
  }
  return deps;
}

let N = 323, relations = [], factorBase = [], solveEnabled = false;

function reset() {
  relations = [];
  solveEnabled = false;
  document.getElementById('solve').disabled = true;
  document.getElementById('matrix-wrap').innerHTML = '';
  document.getElementById('fb-row').textContent = '';
  setStatus(`Choose N and click Run Sieve.`, '');
}

function setStatus(t, cls) {
  const s = document.getElementById('status');
  s.textContent = t; s.className = 'status ' + cls;
}

function renderMatrix() {
  if (relations.length === 0) { document.getElementById('matrix-wrap').innerHTML = ''; return; }
  const fb = factorBase;
  let html = '<table><thead><tr><th>x</th><th>x²mod N</th>';
  fb.forEach(p => { html += '<th>' + p + '</th>'; });
  html += '</tr></thead><tbody>';
  relations.forEach(rel => {
    html += '<tr><td>' + rel.x + '</td><td>' + rel.val + '</td>';
    rel.exps.forEach(e => {
      const bit = e % 2;
      html += '<td class="' + (bit ? 'one' : '') + '">' + bit + '</td>';
    });
    html += '</tr>';
  });
  html += '</tbody></table>';
  document.getElementById('matrix-wrap').innerHTML = html;
}

function runSieve() {
  N = parseInt(document.getElementById('nsel').value);
  factorBase = PRIMES.filter(p => p < Math.ceil(Math.sqrt(N)) + 5);
  document.getElementById('fb-row').textContent =
    fmt(`Factor base: {%factorBase%} (primes ≤ %v%)`, { factorBase: factorBase.join(', '), v: factorBase[factorBase.length-1] });
  relations = [];

  const root = Math.floor(Math.sqrt(N));
  const needed = factorBase.length + 2; // a few extra for dependencies
  let found = 0, tried = 0;

  for (let x = root; x <= root + 200 && found < needed; x++) {
    const val = (x * x) % N;
    tried++;
    const exps = smoothDecompose(val, factorBase);
    if (exps) {
      relations.push({ x, val, exps });
      found++;
    }
  }
  // Also try negative side: x from root down
  for (let x = root - 1; x >= 2 && relations.length < needed + 2; x--) {
    const val = (x * x) % N;
    tried++;
    const exps = smoothDecompose(val, factorBase);
    if (exps) {
      if (!relations.find(r => r.x === x)) relations.push({ x, val, exps });
    }
  }
  relations.sort((a, b) => a.x - b.x);

  renderMatrix();
  if (relations.length >= factorBase.length + 1) {
    document.getElementById('solve').disabled = false;
    setStatus(
      fmt(`Collected %length% smooth relations (need ≥ %length2%%v% ). Ready to solve!`, { length: relations.length, length2: factorBase.length, v: 1 }), 'info'
    );
  } else {
    setStatus(
      fmt(`Only %length% smooth relations found — not enough. Try a larger N.`, { length: relations.length }), 'bad'
    );
  }
}

function solveMatrix() {
  if (relations.length === 0) return;
  const fb = factorBase;
  // Build exponent matrix mod 2
  const matrix = relations.map(r => r.exps.map(e => e % 2));
  const deps = f2Elimination(matrix);

  if (deps.length === 0) {
    setStatus(`No dependency found in this set of relations — run sieve again for more.`, 'bad');
    return;
  }

  // Highlight dependency rows
  const depIdx = deps[0];
  const rows = document.querySelectorAll('#matrix-wrap tbody tr');
  relations.forEach((_, i) => {
    if (rows[i]) rows[i].querySelectorAll('td').forEach(td => td.classList.remove('dep'));
  });

  // For each dependency, attempt to extract a factor
  let factored = false;
  for (const depRow of deps) {
    /* The dependency is a zero row in the eliminated matrix, meaning the
XOR of some subset of original rows is zero — find that subset. */
    // We just use the single relation row directly for the toy demo.
    const rel = relations[depRow];
    const x = rel.x;
    /* y^2 = product of x^2 mod N for subset; in our toy, subset is just {rel}
so y = x, and x^2 mod N = rel.val. We need gcd(x - sqrt(rel.val), N).
For a real factor we need a nontrivial congruence; try gcd(x-1, N) as proxy.
More correctly: combine all rows in subset. */
    // Build subset from the dependency (rows that contribute to zero row):
    // For simplicity, iterate over pairs of relations and check x1^2*x2^2 ≡ square.
    for (let i = 0; i < relations.length && !factored; i++) {
      for (let j = i + 1; j < relations.length && !factored; j++) {
        const ri = relations[i], rj = relations[j];
        // Combined exponents
        const combined = ri.exps.map((e, k) => e + rj.exps[k]);
        const allEven = combined.every(e => e % 2 === 0);
        if (allEven) {
          // x = ri.x * rj.x mod N; y = product of p^(exp/2)
          const xVal = (ri.x * rj.x) % N;
          let yVal = 1;
          combined.forEach((e, k) => { for (let t = 0; t < e / 2; t++) yVal *= fb[k]; });
          yVal = yVal % N;
          const diff = Math.abs(xVal - yVal);
          const f = gcd(diff, N);
          if (f > 1 && f < N) {
            if (rows[i]) rows[i].querySelectorAll('td').forEach(td => td.classList.add('dep'));
            if (rows[j]) rows[j].querySelectorAll('td').forEach(td => td.classList.add('dep'));
            setStatus(
              fmt(`Dependency: rows for x= x=%x% and x= x=%x2% → %xVal%² ≡ %yVal%² (mod %N%) → gcd(%diff%, %N%) = %f% ✓  So %N% = %f% × %v%`, { x: ri.x, x2: rj.x, xVal: xVal, yVal: yVal, N: N, diff: diff, f: f, v: N / f }), 'ok'
            );
            factored = true;
          }
        }
      }
    }
    if (factored) break;
  }
  if (!factored) {
    // fallback: show what we found via gcd of differences
    // Try single-row perfect square
    let fallbackMsg = 'Matrix solved but all gcd values were trivial (1 or N). ';
    for (const rel of relations) {
      const sqRoot = Math.round(Math.sqrt(rel.val));
      if (sqRoot * sqRoot === rel.val) {
        const f = gcd(Math.abs(rel.x - sqRoot), N);
        if (f > 1 && f < N) {
          setStatus(
            rel.x + '² ≡ ' + sqRoot + '² (mod ' + N + ')' +
            ' → gcd(' + Math.abs(rel.x - sqRoot) + ', ' + N + ') = ' + f +
            ' ✓  ' + N + ' = ' + f + ' × ' + (N / f), 'ok'
          );
          factored = true;
          break;
        }
      }
    }
    if (!factored) setStatus(fmt(`%fallbackMsg%Try running the sieve again for more relations.`, { fallbackMsg: fallbackMsg }), 'bad');
  }
}

document.getElementById('sieve').addEventListener('click', runSieve);
document.getElementById('solve').addEventListener('click', solveMatrix);
document.getElementById('reset').addEventListener('click', reset);
document.getElementById('nsel').addEventListener('change', reset);
reset();
Result
Console

Notice: collecting smooth relations is the sieve (fast), while solving the matrix is linear algebra (also fast). The bottleneck is finding enough smooth numbers — and the sub-exponential magic is that smooth numbers appear far more often than you might expect, thanks to the structure of x2 mod Nx^2 \bmod N.

The Real Complexity

The Quadratic Sieve is not just "better than trial division" — its complexity falls into a mathematically precise middle ground.

  • Trial division needs O(N)O(\sqrt{N}) steps — exponential in the number of digits.
  • The Quadratic Sieve runs in LN ⁣[12, 1]=e(1+o(1))ln⁡Nln⁡ln⁡NL_N\!\left[\tfrac{1}{2},\,1\right] = e^{(1+o(1))\sqrt{\ln N \ln \ln N}} steps — sub-exponential: faster than any pure exponential, slower than any polynomial.
  • The General Number Field Sieve (GNFS), invented in 1993, improved the exponent to LN ⁣[13, 1.923]L_N\!\left[\tfrac{1}{3},\,1.923\right] and dethroned QS for large numbers. QS still wins below roughly 1010010^{100}.
  • Status: open. No polynomial-time classical algorithm for factoring is known. If one existed, it would break RSA and imply P=NP\text{P} = \text{NP} is at least consistent — though factoring is not known to be NP-complete. Shor's algorithm (1994) factors in polynomial time on a quantum computer, but no machine large enough to threaten real keys exists yet.

The practical punchline: RSA-768 (232 digits) fell in 2009 after two years of distributed computation using the GNFS — the same family of ideas that the Quadratic Sieve pioneered.

Where It Matters

The Quadratic Sieve is not a curiosity — it is the reason modern cryptographic key sizes are what they are:

  • Breaking RSA keys: QS factored RSA-129 (429 bits) in 1994 after 8 months of distributed work, publicly embarrassing a 1977 bet that it would take 40 quadrillion years.
  • Inspiring the GNFS: the relation-matrix idea from QS was directly generalized into the Number Field Sieve, which is still the best known algorithm for large factors and the tool that cracked RSA-768.
  • Setting key-size floors: every recommendation to use 2048-bit or 4096-bit RSA keys traces back to complexity analyses built on QS and GNFS running times.
  • Smooth-number theory: the smooth-number estimates QS relies on feed into random factoring algorithms, primality tests, and discrete-log computations — the engine under discrete logarithms and elliptic-curve cryptography.
  • Educational value: QS is the cleanest path from "factoring is hard" to "here is why" — it makes the sub-exponential gap tangible and shows how linear algebra over F2\mathbb{F}_2 can crack a number-theory problem.

Understand the Quadratic Sieve and you hold the key to why factoring is believed hard — and what it would mean for cryptography if it were not.

Conclusion

The Quadratic Sieve is a beautiful idea at the boundary of possibility. It does not factor by dividing — it factors by finding coincidences, turning a number-theory problem into a linear-algebra problem that a computer can solve in sub-exponential time.

Its legacy lives in every TLS handshake: the 2048-bit RSA minimum exists because QS and its descendant, the GNFS, are the fastest tools we have — and they are still far slower than polynomial. If someone found a polynomial-time factoring algorithm tomorrow, the world's encrypted traffic would crumble overnight.

Until that day, the Quadratic Sieve stands as proof that mathematics can shave an exponential down to sub-exponential — and a reminder that shaving it all the way to polynomial remains the open question behind P vs NP.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/quadratic-sieve/Content licensed under CC BY-NC 4.0.