Primality Testing in Polynomial Time (Ⅱ)

Fred Akalin

(Note: this article isn't fully polished yet, but I thought it would be a shame to let it languish during my sabbatical. Happy new year!)

5. Strengthening the AKS theorem

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:

(Lenstra's squarefree test.) If \(a^n \equiv a \; (\mathrm{mod} \; n)\) for \(1 \le a \lt \ln^2 n\), then \(n\) is squarefree.[1]

We also need a couple of lemmas:

(Lemma 1.) For \(0 \le a \lt n\) and \(r \gt 1\), let \[ (X + a)^n \equiv X^n + a \; (\mathrm{mod} \; X^r - 1, n)\text{.} \] Then \[ (a + 1)^n = a + 1 \; (\mathrm{mod} \; n)\text{.} \]

Proof. By definition, \((X + a)^n - (X^n + a) = k(X) \cdot (X^r - 1) \; (\mathrm{mod} \; n)\). Treating both sides as a function of \(x\) and substituting in \(1\), we immediately get \((1 + a)^n - (1 + a) = 0 \; (\mathrm{mod} \; n)\). ∎

(Lemma 2.) For \(n \gt 1\), \(\lfloor \lg n \rfloor \cdot \lg n \gt \ln^2 n\).

Proof. Since \(\ln n = \frac{\lg n}{\lg e}\) and \(e \gt 2\), \(\lg n \gt \ln n\) for \(n \gt 1\).

Letting \(k = \lfloor \lg n \rfloor\), \(\ln n \lt \frac{k + 1}{\lg e}\), so if \(\frac{k + 1}{\lg e} \lt k\), that implies that \(\ln n \lt k\). Solving for \(k\), we get that \(k \gt \frac{1}{\lg e - 1}\), which is true when \(n \ge 8\).

So if \(n \ge 8\), then \(\ln n \lt \lfloor \lg n \rfloor\). Checking manually, we find that \(\ln n \lt \lfloor \lg n \rfloor\) holds also for \(n \in \{ 2, 4, 5, 6, 7 \}\), immediately implying the lemma for all \(n \gt 1\) except \(3\). But checking manually again, we find that the lemma holds for \(3\) also. ∎

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

(AKS theorem, strong version.) Let \(n \ge 2\), \(r\) be relatively prime to \(n\) with \(o_r(n) \gt \lg^2 n\), and \(M \gt \sqrt{\varphi(r)} \lg n\). Furthermore, let \(n\) have no prime factor less than \(M\) and let \[ (X + a)^n \equiv X^n + a \; (\mathrm{mod} \; X^r - 1, n)\text{.} \] for \(0 \le a \lt M\). Then \(n\) is prime.

Proof. From Lemma 1, we know that \(a^n = a \; (\mathrm{mod} \; n)\) for \(1 \le a \lt M\). Since \(M \gt \lfloor \sqrt{t} \rfloor \lg n \gt \lfloor \lg n \rfloor \cdot \lg n \gt \ln^2 n\) by Lemma 2, we can apply Lenstra's squarefree test to show that \(n\) is squarefree. From the weak version of the AKS theorem, we also know that \(n\) is a prime power. But since \(n\) is squarefree, it must have distinct prime factors, which immediately implies that \(n\) is prime. ∎

6. Finding a suitable \(r\)

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:

(Upper bound for \(r\).) Let \(n \ge 2\). Then there exists some \(r \le \max(3, \lceil \lg n \rceil^5)\) such that \(o_r(n) \gt \lceil \lg n \rceil^2\).[2]

(Proof.) Let's first prove the following lemma:

(Lemma 3.) Let \(n \ge 9\) and \(b = \lceil \lg n \rceil\). Then for \(m \ge 1\), there exists some \(r \le b^{2m + 1}\) such that \(o_r(n) \gt b^m\).

(Proof.) Let \[ N = n \cdot (n - 1) \cdot (n^2 - 1) \cdot \ldots \cdot (n^{b^m} - 1)\text{.} \] Note that \(r\) divides \(N\) if and only if \(o_r(n) \le b^m\). So it suffices to find some \(r\) that does not divide \(N\).

We can see that: \[ \begin{aligned} N &= n \cdot (n - 1) \cdot (n^2 - 1) \cdot \ldots \cdot (n^{b^m} - 1) \\ &\lt n \cdot n \cdot n^2 \cdot \ldots \cdot n^{b^m} \\ &= n^{1 + 1 + 2 + 3 + \ldots + b^m} \\ &= n^{1 + b^m (b^m + 1) / 2} \\ &= n^{\frac{1}{2} b^{2m} + \frac{1}{2} b^m + 1}\text{.} \end{aligned} \] Furthermore, we can upper-bound the exponent of \(n\): \[ \begin{aligned} b^{2m} &\gt \frac{1}{2} b^{2m} + \frac{1}{2} b^m + 1 \\ \frac{1}{2} b^{2m} - \frac{1}{2} b^m - 1 &\gt 0 \\ b^{2m} - b^m - 2 &\gt 0 \\ (b^m - 2) \cdot (b^m + 1) &\gt 0\text{.} \end{aligned} \] The last statement holds when \(b^m \gt 2\), which is always since \(b \ge 4\) and \(m \ge 1\).

Applying the upper bound, \[ \begin{aligned} N &\lt n^{\frac{1}{2} b^{2m} + \frac{1}{2} b^m + 1} \\ &\lt n^{b^{2m}} \\ &\le 2^{b^{2m + 1}}\text{.} \end{aligned} \]

We can then use the following theorem, which we'll state without proof:

(Primorial lower bound.) For \(x \ge 31\), the product of primes \(\le x\) exceeds \(2^x\).[3] That is, \[ x\# = \prod_{p \le x\text{, }p\text{ is prime}} p \gt 2^x\text{.} \]

Since \(b \ge 4\) and \(m \ge 1\), \(b^{2m + 1} \ge 31\), and so \(2^{b^{2m + 1}} \lt (b^{2m + 1})\#\). Therefore, \[ N \lt 2^{b^{2m + 1}} \lt (b^{2m + 1})\#\text{.} \] But that implies that there is some prime number \(p_0 \le b^{2m + 1}\) that does not divide \(N\); if they all did, then \(N\) would be at least their product \((b^{2m + 1})\#\), contradicting the inequality above. Therefore, \(o_{p_0}(n) \gt b^m\). ∎

We can then prove our theorem: for \(n \ge 9\), apply Lemma 3 with \(m = 2\). Here are explicit values for the rest: for \(n = 2\), \(r = 3\); \(n = 3\), \(r = 7\); \(n \in \{ 4, 6, 7, 8\}\), \(r = 11\); and for \(n = 5\), \(r = 17\). ∎

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

(Tight upper bound for some \(r\).) Let \(n \equiv \pm 3 \; (\mathrm{mod} \; 8)\). Then there exists some \(r \lt 8 \lceil \lg n \rceil^2\) such that \(o_r(n) \gt \lceil \lg n \rceil^2\).[4]

7. The AKS algorithm (simple version)

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”.

As we've showed in the previous section, there always exists an \(r\) such that \(o_r(n) \gt \lceil \lg n \rceil^2\), so step 1 will terminate. All other steps are bounded, so the entire algorithm will always terminate.

In step 2, since \(\varphi(r) \le r - 1\), the value of \(M\) that we compute is always greater than \(\sqrt{\varphi(r)} \lceil \lg n \rceil\). Step 4 checks if \((X + a)^n \equiv X^n + a \; (\mathrm{mod} \; (X^r - 1, n))\) holds. Therefore, By the strong AKS theorem, if the algorithm returns “prime”, then \(n\) is prime. Furthermore, by the weak version of Fermat's little theorem for polynomials, if the algorithm returns “composite”, then \(n\) is composite.

Since the algorithm always terminates and it returns the correct answer when it terminates, it is totally correct.

As shown in the previous section, we have to test \(O(\lg^5 n)\) values to find a suitable \(r\). Assuming a straightforward algorithm to compute the multiplicative order that bails out once \(\lfloor \lg n \rfloor^2\) is reached, and assuming we use the division-based Euclidean algorithm for computing the greatest common divisor, testing each value takes \(O(\lg^2 n)\) multiplies and \(O(\lg r) = O(\lg \lg n)\) divisions of \(O(\lg r)\)-bit numbers. Let \(M(b)\) be the cost to multiply two \(b\)-bit numbers. The complexity of division is asymptotically the same as multiplication, so the total cost of step 1 is \(O(\lg^5 n \cdot (\lg^2 n + \lg \lg n) \cdot M(\lg \lg n)) = O(\lg^7 n \cdot M(\lg \lg n))\), assuming \(M(O(b)) = O(M(b))\).

Step 2 involves one square root, one multiplication, and one increment, all involving \(O(\lg \lg n)\)-bit numbers. The complexity of taking the square root is asymptotically the same as multiplication, so the total cost of step 2 is \(O(M(\lg \lg n))\).

Step 3 takes a square root and tests \(M = O(\lg^{7/2} n)\) numbers, and each test involves dividing two \(O(\lg M)\)-bit numbers, so the total cost of step 3 is \(O(\lg^{7/2} n \cdot M(\lg \lg n))\).

Steps 4 and 5 test \(O(\lg^{7/2} n)\) polynomials. Testing each polynomial involves exponentiating it by \(n\), reducing power mod \(r\) and coefficients mod \(n\) at each step, which requires \(O(\lg n)\) multiplications of polynomials with \(O(r)\) coefficients each of size \(O(\lg n)\). The cost of multiplying two polynomials with \(s\) coefficients of size \(b\) is \(M(s) \cdot M(b)\), so the total cost of steps 4 and 5 is \(O(\lg^{9/2} n \cdot M(\lg^5 n \cdot \lg \lg n))\), assuming \(M(a) \cdot M(b) = M(a \cdot b)\).

If long multiplication is used, then it costs \(M(b) = b^2\), which gives a total cost of \(O(\lg^{29/2} n \cdot \lg^2 \lg n) = O(\lg^{15} n)\) for the whole algorithm. More complicated multiplication methods cost only \(M(b) = b \lg b\), which gives a total cost of \(O(\lg^{10} n)\) for the whole algorithm. Either way, the AKS primality test is shown to be implementable in polynomial time.

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 \(\varphi(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);
}

Let n = .

(To-do: Have an interactive box to demonstrate how the per-\(a\) AKS test works.)

8. The AKS algorithm (improved version)

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 \(\varphi(r)\) using \(r\)'s prime factors.
    6. Compute \(M = \lfloor \sqrt{\varphi(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”.

The logic of steps 1 to 3 of the simple version is essentially merged together to form steps 1 and 2 of this version; since each \(r\) has to be checked for co-primality with \(n\), that effectively also checks if \(r\) is a prime factor of \(n\), so we only have to check for prime factors of \(n\) up to the lower bound of \(r\). Furthermore, both the multiplicative order as well as the totient function can be computed more quickly given a complete prime factorization, so we can compute that for each \(r\). Third, we use \(\varphi(r)\) instead of \(r - 1\) to give a tighter bound for \(M\). Finally, the last two steps are the same as in the simple version.

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

(To-do: Wrap up and lead into what will be shown in part 3.)

[1] This is a version of Theorem 2 from Lenstra's paper Miller's Primality Test.

[2] We work with \(\lceil \lg n \rceil^2\) instead of \(\lceil \lg^2 n \rceil\) or \(\lg^2 n\) as it's easier to work with in an actual implementation.

[3] This is exercise 1.27 from Prime Numbers: A Computational Perspective.

[4] This is an adapted from section 8.4 of Granville's It is Easy to Determine Whether a Given Number is Prime.

[5] The SNat class used is the same as in my previous article, An Introduction to Primality Testing.