AKS test for primes

AKS test for primes
You are encouraged to solve this task according to the task description, using any language you may know.

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 result)))))

(defun primep (p)

 (cond
   ((< p 2) nil)

   (t (let ((c (coefficients p)))
        (decf (elt c 0))
        (loop for i from 0 upto (/ (length c) 2)
              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 1000
       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]

Perl 6

<lang perl6>constant coeff = [1], [1,-1], -> @p { [@p,0 Z- 0,@p] } ... *;

sub aks($p) { $p > 1 and so coeff[$p][1 ..^ */2].all %% $p }

say ' p: (x-1)ᵖ'; say '-----------';

for ^13 -> $p {

   print $p.fmt('%2i: ');
   $_ = ~reverse gather
       for coeff[$p].reverse.kv -> $e, $n {
           my $c = $n < 0 ?? "- {abs $n}" !! "+ $n";
           given $e {
               when 0  { take $c }
               when 1  { take $c ~ 'x' }
               default { take $c ~ 'x' ~ $e.trans('0123456789' => '⁰¹²³⁴⁵⁶⁷⁸⁹') }
           }
       }
   s/^ '+ ' //;
   s/^ '1x' /x/;
   .say;

} say ;

say "Primes up to 100:"; say " ", grep { aks $_ }, 2..100;</lang>

Output:

 p: (x-1)ᵖ
-----------
 0: 1
 1: x - 1
 2: x² - 2x + 1
 3: x³ - 3x² + 3x - 1
 4: x⁴ - 4x³ + 6x² - 4x + 1
 5: x⁵ - 5x⁴ + 10x³ - 10x² + 5x - 1
 6: x⁶ - 6x⁵ + 15x⁴ - 20x³ + 15x² - 6x + 1
 7: x⁷ - 7x⁶ + 21x⁵ - 35x⁴ + 35x³ - 21x² + 7x - 1
 8: x⁸ - 8x⁷ + 28x⁶ - 56x⁵ + 70x⁴ - 56x³ + 28x² - 8x + 1
 9: x⁹ - 9x⁸ + 36x⁷ - 84x⁶ + 126x⁵ - 126x⁴ + 84x³ - 36x² + 9x - 1
10: x¹⁰ - 10x⁹ + 45x⁸ - 120x⁷ + 210x⁶ - 252x⁵ + 210x⁴ - 120x³ + 45x² - 10x + 1
11: x¹¹ - 11x¹⁰ + 55x⁹ - 165x⁸ + 330x⁷ - 462x⁶ + 462x⁵ - 330x⁴ + 165x³ - 55x² + 11x - 1
12: x¹² - 12x¹¹ + 66x¹⁰ - 220x⁹ + 495x⁸ - 792x⁷ + 924x⁶ - 792x⁵ + 495x⁴ - 220x³ + 66x² - 12x + 1

Primes up to 100:
  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]

Python: Output formatted for wiki

Using a wikitable and math features with the following additional code produces better formatted polynomial output:

<lang python>print(

for p in range(12): print('! $\%i$\n| $\%s$\n|-'  % (p, ' '.join('%s%s' % (('%+i' % e) if (e != 1 or not p or (p and not n) ) else '+', for n,e in enumerate(expand_x_1(p))))) print('|}')</lang>

Extra output rendered by wiki

Polynomial Expansions

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

Polynomial Expansions

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

Racket

This example is incomplete. Needs to show the polynomial expansion rather than simply the coeffiients. Please ensure that it meets all task requirements and remove this message.

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)

Rust

This example is incomplete. Needs to show the polynomial expansion rather than simply the coefficients. Please ensure that it meets all task requirements and remove this message.

<lang rust>//Rust 0.9

fn coefficients(p: uint) -> ~[i64] {

   if p==0 {
       ~[1i64]
   } else {
       let mut result: ~[i64] = ~[1, -1];
       let zero = Some(0i64);
       for _ in range(1, p) {
           result = {
               let a = result.iter().chain(zero.iter());
               let b = zero.iter().chain(result.iter());
               a.zip(b).map(|(x, &y)| x-y).to_owned_vec()
           };
       }
       result
   }

}

fn is_prime(p: uint) -> bool {

   if p<2 {
       false
   } else {
       let mut c = coefficients(p);
       c[0] -= 1;
       for i in range(0, c.iter().len()/2) {
           if (c[i] % (p as i64)) != 0 {
               return false
           }
       }
       true
   }

}

fn main() {

   for p in range(0, 8) {
       print(p.to_str() + ": ");
       println(coefficients(p as uint).to_str());
   }

   for p in range(1, 51) {
       if is_prime(p as uint) {
           print(p.to_str() + " ");
       }
   }
   println("");

}</lang>

Output:

0: [1]
1: [1, -1]
2: [1, -2, 1]
3: [1, -3, 3, -1]
4: [1, -4, 6, -4, 1]
5: [1, -5, 10, -10, 5, -1]
6: [1, -6, 15, -20, 15, -6, 1]
7: [1, -7, 21, -35, 35, -21, 7, -1]
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47