Primality Testing in Polynomial Time (Ⅱ)

It turns out the conditions of the AKS theorem are stronger than they appear; they themselves imply that \(n\) is prime. To show this, we need the following theorem, which we'll state without proof: We also need a couple of lemmas:

Then, we can prove the strong version of the AKS theorem:

The only remaining loose end is to show that there exists an \(r\) with \(o_r(n) \gt \lg^2 n\) and that it's small enough (i.e., polylog in \(n\)). The existence of \(r\) is easy to see; we can simply pick the smallest \(r\) that is co-prime to \(n\) and greater than \(n^{\lg^2 n}\). But that's obviously too big. We can do better:

Also, it turns out that about half the time, we can do better. We'll state this theorem without proof:

Without further ado, here is a simple version of the AKS algorithm, given \(n \ge 2\):

1. Starting from \(\lceil \lg n \rceil^2 + 2\), find an \(r\) such that \(\gcd(r, n) = 1\) and \(o_r(n) \gt \lceil \lg n \rceil^2\).
2. Compute \(M = \lfloor \sqrt{r - 1} \rfloor \lceil \lg n \rceil + 1\).
3. Search for a prime factor of \(n\) less than \(M\). If one is found, return “composite”. If none are found and \(M \ge \lfloor \sqrt{n} \rfloor\), return “prime”.
4. For each \(1 \le a \lt M\), compute \((X + a)^n\), reducing coefficients mod \(n\) and powers mod \(r\). If the result is not equal to \(X^{n\text{ mod }r} + a\), return “composite”.
5. Otherwise, return “prime”.

Below is step 1 implemented in Javascript; however, here we bound \(r\) explicitly to be able to detect bugs quickly.^[5]

// Returns an upper bound for r such that o_r(n) > ceil(lg(n))^2 that
// is polylog in n.
function calculateAKSModulusUpperBound(n) {
  n = SNat.cast(n);
  var ceilLgN = new SNat(n.ceilLg());
  var rUpperBound = ceilLgN.pow(5).max(3);
  var nMod8 = n.mod(8);
  if (nMod8.eq(3) || nMod8.eq(5)) {
    rUpperBound = rUpperBound.min(ceilLgN.pow(2).times(8));
  }
  return rUpperBound;
}

// Returns the least r such that o_r(n) > ceil(lg(n))^2 >= ceil(lg(n)^2).
function calculateAKSModulus(n, multiplicativeOrderCalculator) {
  n = SNat.cast(n);
  multiplicativeOrderCalculator =
    multiplicativeOrderCalculator || calculateMultiplicativeOrderCRT;

  var ceilLgN = new SNat(n.ceilLg());
  var ceilLgNSq = ceilLgN.pow(2);
  var rLowerBound = ceilLgNSq.plus(2);
  var rUpperBound = calculateAKSModulusUpperBound(n);

  for (var r = rLowerBound; r.le(rUpperBound); r = r.plus(1)) {
    if (n.gcd(r).ne(1)) {
      continue;
    }
    var o = multiplicativeOrderCalculator(n, r);
    if (o.gt(ceilLgNSq)) {
      return r;
    }
  }

  throw new Error('Could not find AKS modulus');
}

Here is step 2 implemented in Javascript:

// Returns floor(sqrt(r-1)) * ceil(lg(n)) + 1 > floor(sqrt(Phi(r))) * lg(n).
function calculateAKSUpperBoundSimple(n, r) {
  n = SNat.cast(n);
  r = SNat.cast(r);

  // Use r - 1 instead of calculating Phi(r).
  return r.minus(1).floorRoot(2).times(n.ceilLg()).plus(1);
}

Here is part of step 3 implemented in Javascript, along with the comments for the functions used in trial division:

// Given a number n, a generator function getNextDivisor, and a
// processing function processPrimeFactor, factors n using the
// divisors returned by genNextDivisor and passes each prime factor
// with its multiplicity to processPrimeFactor.
//
// getNextDivisor is passed the current unfactorized part of n and it
// should return the next divisor to try, or null if there are no more
// divisors to generate (although processPrimeFactor may still be
// called).  processPrimeFactor is called with each non-trivial prime
// factor and its multiplicity.  If it returns a false value, it won't
// be called anymore.
function trialDivide(n, getNextDivisor, processPrimeFactor) {
  ...
}

// Returns a generator that generates primes up to 7, then odd numbers
// up to floor(sqrt(n)), using a mod-30 wheel to eliminate odd numbers
// that are known composite (roughly half).
function makeMod30WheelDivisorGenerator() {
  ...
}

// Returns the first factor of n < m from generator, or null if there
// is no such factor.
function getFirstFactorBelow(n, M, generator) {
  n = SNat.cast(n);
  M = SNat.cast(M);
  generator = generator || makeMod30WheelDivisorGenerator();

  var boundedGenerator = function(n) {
    var d = generator(n);
    return (d && d.lt(M)) ? d : null;
  };
  var factor = null;
  trialDivide(n, boundedGenerator, function(p, k) {
    if (p.lt(M.min(n))) {
      factor = p;
    }
    return false;
  });
  return factor;
}

Below is a function that ties steps 1 to 3 together; it is useful for testing purposes to separate it from the other steps. (Actually, we use a different function to compute \(M\) which computes \(φ(r)\) instead of using \(r - 1\) so that we always have the tightest bound possible for \(M\).)

// The getAKSParameters* functions below return a parameters object
// with the following fields:
//
//   n: the number the parameters are for.
//
//   factor: A prime factor of n.  If present, the fields below may
//           not be present.
//
//   isPrime: if set, n is prime.  If present, the fields below may
//            not be present.
//
//   r: the AKS modulus for n.
//
//   M: the AKS upper bound for n.

function getAKSParametersSimple(n) {
  n = SNat.cast(n);

  var r = calculateAKSModulus(n);
  var M = calculateAKSUpperBound(n, r);
  var parameters = {
    n: n,
    r: r,
    M: M
  };

  var factor = getFirstFactorBelow(n, M);
  if (factor) {
    parameters.factor = factor;
  } else if (M.gt(n.floorRoot(2))) {
    parameters.isPrime = true;
  }

  return parameters;
}

Finally, here is step 4 implemented in Javascript:

// Returns whether (X + a)^n = X^n + a mod (X^r - 1, n).
function isAKSWitness(n, r, a) {
  n = SNat.cast(n);
  r = SNat.cast(r);
  a = SNat.cast(a);

  function reduceAKS(p) {
    return p.modPow(r).mod(n);
  }

  function prodAKS(x, y) {
    return reduceAKS(x.times(y));
  };

  var one = new SPoly(new SNat(1));
  var xn = one.shiftLeft(n.mod(r));
  var ap = new SPoly(a);
  var lhs = one.shiftLeft(1).plus(ap).pow(n, prodAKS);
  var rhs = reduceAKS(one.shiftLeft(n).plus(ap));
  return lhs.ne(rhs);
}

// Returns the first a < M that is an AKS witness for n, or null if
// there isn't one.
function getFirstAKSWitness(n, r, M) {
  for (var a = new SNat(1); a.lt(M); a = a.plus(1)) {
    if (isAKSWitness(n, r, a)) {
      return a;
    }
  }
  return null;
}

Here's the code that ties it all together:

// Returns whether n is prime or not using the AKS primality test.
function isPrimeByAKS(n) {
  n = SNat.cast(n);

  var parameters = getAKSParameters(n);
  if (parameters.factor) {
    return false;
  }
  if (parameters.isPrime) {
    return true;
  }
  return (getFirstAKSWitness(n, parameters.r, parameters.M) == null);
}

Here is a slightly more complicated version of the AKS algorithm. Again given \(n \ge 2\):

1. Search for a prime factor of \(n\) less than \(\lceil \lg n \rceil^2 + 2\). If one is found, return “composite”.
2. For each \(r\) from \(\lceil \lg n \rceil^2 + 2\):
   1. If \(r \gt \lfloor \sqrt{n} \rfloor\), return “prime”.
   2. If \(r\) divides \(n\), return “composite”.
   3. Otherwise, factorize \(r\).
   4. Compute \(o_r(n)\) using \(r\)'s prime factors. If it is less than or equal to \(\lceil \lg n \rceil^2\), jump back to the top of the loop with the next \(r\).
   5. Otherwise, compute \(φ(r)\) using \(r\)'s prime factors.
   6. Compute \(M = \lfloor \sqrt{φ(r)} \rfloor \lceil \lg n \rceil + 1\), and break out of the loop.
3. For each \(1 \le a \lt M\), compute \((X + a)^n\), reducing coefficients mod \(n\) and powers mod \(r\). If the result is not equal to \(X^{n\text{ mod }r} + a\), return “composite”.
4. Otherwise, return “prime”.

Here are steps 1 and 2 of the above algorithm, implemented in Javascript:

function getAKSParameters(n, factorizer) {
  n = SNat.cast(n);
  factorizer = factorizer || defaultFactorizer;

  var ceilLgN = new SNat(n.ceilLg());
  var ceilLgNSq = ceilLgN.pow(2);
  var floorSqrtN = n.floorRoot(2);

  var rLowerBound = ceilLgNSq.plus(2);
  var rUpperBound = calculateAKSModulusUpperBound(n).min(floorSqrtN);

  var parameters = {
    n: n
  };

  var factor = getFirstFactorBelow(n, rLowerBound);
  if (factor) {
    parameters.factor = factor;
    return parameters;
  }

  for (var r = rLowerBound; r.le(rUpperBound); r = r.plus(1)) {
    if (n.mod(r).isZero()) {
      parameters.factor = d;
      return parameters;
    }

    var rFactors = getFactors(r, factorizer);
    var o = calculateMultiplicativeOrderCRTFactors(n, rFactors, factorizer);
    if (o.gt(ceilLgNSq)) {
      parameters.r = r;
      parameters.M = calculateAKSUpperBoundFactors(n, rFactors);
      return parameters;
    }
  }

  if (rUpperBound.eq(floorSqrtN)) {
    parameters.isPrime = true;
    return parameters;
  }

  throw new Error('Could not find AKS modulus');
}