AKS test for primes

AKS test for primes is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

The AKS test for primes states that a number $p$ is prime if all the coefficients of the polynomial expansion of

$(x-1)^{p}-(x^{p}-1)$

are divisible by $p$.

For example, trying $p=3$:

$(x-1)^{3}-(x^{3}-1)=(X^{3}-3X^{2}+3X-1)-(x^{3}-1)=-3X^{2}+3X$

And all the coefficients are divisible by 3 so 3 is prime by the AKS test.

The task

1. Create a function/subroutine/method that given p generates the coefficients of the expanded polynomial representation of $(x-1)^{p}$.
2. Use the function to show here the polynomial expansions of p for p in the range 0 to at least 7, inclusive.
3. Use the previous function in creating another function that when given p returns whether p is prime using the AKS test.
4. Use your AKS test to generate a list of all primes under 35.
5. As a stretch goal, generate all primes under 50 (Needs greater than 31 bit integers).

Reference

• Fool-Proof Test for Primes - Numberphile (Video).

C

<lang c>#include <stdio.h>

1. include <stdlib.h>

long long c[100];

void coef(int n) { int i, j;

if (n < 0 || n > 63) abort(); // gracefully deal with range issue

for (c[i=0] = 1; i < n; c[0] = -c[0], i++) for (c[1 + (j=i)] = 1; j > 0; j--) c[j] = c[j-1] - c[j]; }

int is_prime(int n) { int i;

coef(n); c[0] += 1, c[i=n] -= 1; while (i-- && !(c[i] % n));

return i < 0; }

void show(int n) { do printf("%+lldx^%d", c[n], n); while (n--); }

int main(void) { int n;

for (n = 0; n < 10; n++) { coef(n); printf("(x-1)^%d = ", n); show(n); putchar('\n'); }

printf("\nprimes (never mind the 1):"); for (n = 1; n <= 63; n++) if (is_prime(n)) printf(" %d", n);

putchar('\n'); return 0; }</lang>

The ugly output:

(x-1)^0 = +1x^0
(x-1)^1 = +1x^1-1x^0
(x-1)^2 = +1x^2-2x^1+1x^0
(x-1)^3 = +1x^3-3x^2+3x^1-1x^0
(x-1)^4 = +1x^4-4x^3+6x^2-4x^1+1x^0
(x-1)^5 = +1x^5-5x^4+10x^3-10x^2+5x^1-1x^0
(x-1)^6 = +1x^6-6x^5+15x^4-20x^3+15x^2-6x^1+1x^0
(x-1)^7 = +1x^7-7x^6+21x^5-35x^4+35x^3-21x^2+7x^1-1x^0
(x-1)^8 = +1x^8-8x^7+28x^6-56x^5+70x^4-56x^3+28x^2-8x^1+1x^0
(x-1)^9 = +1x^9-9x^8+36x^7-84x^6+126x^5-126x^4+84x^3-36x^2+9x^1-1x^0

primes (never mind the 1): 1 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61

Common Lisp

<lang lisp>(defun coefficients (p)

 (cond
   ((= p 0) #(1))

   (t (loop for i from 1 upto p
            for result = #(1 -1) then (map 'vector
                                           #'-
                                           (concatenate 'vector result #(0))
                                           (concatenate 'vector #(0) result))
            finally (return (reverse result))))))

(defun primep (p)

 (cond
   ((< p 2) nil)

   (t (let ((c (coefficients p)))
        (incf (elt c 0))
        (loop for i from 2 upto (length c)
              for x across c
              never (/= (mod x p) 0))))))

(defun main ()

 (format t "# p: (x-1)^p for small p:~%")
 (loop for p from 0 upto 7
       do (format t "~D: " p)
          (loop for i from 0
                for x across (coefficients p)
                do (when (>= x 0) (format t "+"))
                   (format t "~D" x)
                   (if (> i 0)
                       (format t "X^~D " i)
                       (format t " ")))
          (format t "~%"))
 (loop for i from 0 to 50
       do (when (primep i) (format t "~D " i)))
 (format t "~%"))</lang>

Output:

# p: (x-1)^p for small p:
0: +1 
1: -1 +1X^1 
2: +1 -2X^1 +1X^2 
3: -1 +3X^1 -3X^2 +1X^3 
4: +1 -4X^1 +6X^2 -4X^3 +1X^4 
5: -1 +5X^1 -10X^2 +10X^3 -5X^4 +1X^5 
6: +1 -6X^1 +15X^2 -20X^3 +15X^4 -6X^5 +1X^6 
7: -1 +7X^1 -21X^2 +35X^3 -35X^4 +21X^5 -7X^6 +1X^7 
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 

D

Translation of: Python

<lang d>import std.stdio, std.range, std.algorithm, std.string, std.bigint;

BigInt[] expandX1(in uint p) pure /*nothrow*/ {

   if (p == 0) return [1.BigInt];
   typeof(return) r = [BigInt(1), BigInt(-1)];
   foreach (immutable _; 1 .. p)
       r = zip(r~0.BigInt, 0.BigInt~r).map!(xy => xy[0]-xy[1]).array;
   r.reverse();
   return r;

}

bool aksTest(in uint p) pure /*nothrow*/ {

   if (p < 2) return false;
   auto ex = p.expandX1;
   ex[0]++;
   return !ex[0 .. $ - 1].any!(mult => mult % p);

}

void main() {

   "# p: (x-1)^p for small p:".writeln;
   foreach (immutable p; 0 .. 12)
       writefln("%3d: %s", p, p.expandX1.zip(iota(p + 1))
           .map!q{"%+dx^%d ".format(a[])}.join.replace("x^0", ""));

   "\nSmall primes using the AKS test:".writeln;
   101.iota.filter!aksTest.writeln;

}</lang>

Output:

# p: (x-1)^p for small p:
  0: +1 
  1: -1 +1x^1 
  2: +1 -2x^1 +1x^2 
  3: -1 +3x^1 -3x^2 +1x^3 
  4: +1 -4x^1 +6x^2 -4x^3 +1x^4 
  5: -1 +5x^1 -10x^2 +10x^3 -5x^4 +1x^5 
  6: +1 -6x^1 +15x^2 -20x^3 +15x^4 -6x^5 +1x^6 
  7: -1 +7x^1 -21x^2 +35x^3 -35x^4 +21x^5 -7x^6 +1x^7 
  8: +1 -8x^1 +28x^2 -56x^3 +70x^4 -56x^5 +28x^6 -8x^7 +1x^8 
  9: -1 +9x^1 -36x^2 +84x^3 -126x^4 +126x^5 -84x^6 +36x^7 -9x^8 +1x^9 
 10: +1 -10x^1 +45x^2 -120x^3 +210x^4 -252x^5 +210x^6 -120x^7 +45x^8 -10x^9 +1x^10 
 11: -1 +11x^1 -55x^2 +165x^3 -330x^4 +462x^5 -462x^6 +330x^7 -165x^8 +55x^9 -11x^10 +1x^11 

Small primes using the AKS test:
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

Python

<lang python>def expand_x_1(p):

   if p == 0: return [1]
   ex = [1, -1]
   for i in range(1, p):
       ex = [x - y for x,y in zip(ex+[0], [0]+ex)]
   return ex[::-1]

def aks_test(p):

   if p < 2: return False
   ex = expand_x_1(p)
   ex[0] += 1
   return not any(mult % p for mult in ex[0:-1])
   
   

print('# p: (x-1)^p for small p') for p in range(12):

   print('%3i: %s' % (p, ' '.join('%+i%s' % (e, ('X^%i' % n) if n else )
                                  for n,e in enumerate(expand_x_1(p)))))

print('\n# small primes using the aks test') print([p for p in range(101) if aks_test(p)])</lang>

Output:

# p: (x-1)^p for small p
  0: +1
  1: -1 +1X^1
  2: +1 -2X^1 +1X^2
  3: -1 +3X^1 -3X^2 +1X^3
  4: +1 -4X^1 +6X^2 -4X^3 +1X^4
  5: -1 +5X^1 -10X^2 +10X^3 -5X^4 +1X^5
  6: +1 -6X^1 +15X^2 -20X^3 +15X^4 -6X^5 +1X^6
  7: -1 +7X^1 -21X^2 +35X^3 -35X^4 +21X^5 -7X^6 +1X^7
  8: +1 -8X^1 +28X^2 -56X^3 +70X^4 -56X^5 +28X^6 -8X^7 +1X^8
  9: -1 +9X^1 -36X^2 +84X^3 -126X^4 +126X^5 -84X^6 +36X^7 -9X^8 +1X^9
 10: +1 -10X^1 +45X^2 -120X^3 +210X^4 -252X^5 +210X^6 -120X^7 +45X^8 -10X^9 +1X^10
 11: -1 +11X^1 -55X^2 +165X^3 -330X^4 +462X^5 -462X^6 +330X^7 -165X^8 +55X^9 -11X^10 +1X^11

# small primes using the aks test
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

Racket

With copious use of the math/number-theory library...

<lang racket>#lang racket (require math/number-theory)

1. coefficients of expanded polynomial (x-1)^p

produces a vector because in-vector can provide a start

and stop (of 1 and p) which allow us to drop the (-1)^p

and the x^p terms, respectively.

(vector-ref (coefficients p) e) is the coefficient for p^e

(define (coefficients p)

 (for/vector ((exponent (in-range 0 (add1 p))))
   (define sign (expt -1 (- p exponent)))
   (* sign (binomial p exponent))))

2. Show the polynomial expansions from p=0 .. 7 (inclusive)

(for ((p (in-range 0 (add1 7)))) (printf "p=~a: ~a~%" p (coefficients p)))

3. AKS primeality test

(define (prime?/AKS p)

 (for/and
     ((coefficient (in-vector (coefficients p) 1 p)))
   (divides? p coefficient)))

4. list of numbers < 35 that are prime (note that 1 is prime

by the definition of the AKS test for primes)

(displayln (for/list ((i (in-range 1 35)) #:when (prime?/AKS i)) i))

5. stretch goal

all prime numbers under 50

(displayln (for/list ((i (in-range 1 50)) #:when (prime?/AKS i)) i))</lang>

Output:

p=0: #(1)
p=1: #(-1 1)
p=2: #(1 -2 1)
p=3: #(-1 3 -3 1)
p=4: #(1 -4 6 -4 1)
p=5: #(-1 5 -10 10 -5 1)
p=6: #(1 -6 15 -20 15 -6 1)
p=7: #(-1 7 -21 35 -35 21 -7 1)
(1 2 3 5 7 11 13 17 19 23 29 31)
(1 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47)