AKS test for primes: Difference between revisions

Content added Content deleted

Inline

Revision as of 19:16, 10 February 2014

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

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

Reference

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

AutoHotkey

Works with: AutoHotkey L

<lang autohotkey>; 1. Create a function/subroutine/method that given p generates the coefficients of the expanded polynomial representation of (x-1)^p.

Function modified from http://rosettacode.org/wiki/Pascal%27s_triangle#AutoHotkey

pascalstriangle(n=8) ; n rows of Pascal's triangle { p := Object(), z:=Object() Loop, % n Loop, % row := A_Index col := A_Index , p[row, col] := row = 1 and col = 1 ? 1 : (p[row-1, col-1] = "" ; math operations on blanks return blanks; I want to assume zero ? 0 : p[row-1, col-1]) - (p[row-1, col] = "" ? 0 : p[row-1, col]) Return p }

2. Use the function to show here the polynomial expansions of p for p in the range 0 to at least 7, inclusive.

For k, v in pascalstriangle() { s .= "`n(x-1)^" k-1 . "=" For k, w in v s .= "+" w "x^" k-1 } s := RegExReplace(s, "\+-", "-") s := RegExReplace(s, "x\^0", "") s := RegExReplace(s, "x\^1", "x") Msgbox % clipboard := s

3. Use the previous function in creating another function that when given p returns whether p is prime using the AKS test.

aks(n) { isnotprime := False For k, v in pascalstriangle(n+1)[n+1] (k != 1 and k != n+1) ? isnotprime |= !(v // n = v / n) ; if any is not divisible, returns true Return !isnotprime }

4. Use your AKS test to generate a list of all primes under 35.

i := 49 p := pascalstriangle(i+1) Loop, % i { n := A_Index isnotprime := False For k, v in p[n+1] (k != 1 and k != n+1) ? isnotprime |= !(v // n = v / n) ; if any is not divisible, returns true t .= isnotprime ? "" : A_Index " " } Msgbox % t Return</lang>

Output:

(x-1)^0=+1
(x-1)^1=-1+1x
(x-1)^2=+1-2x+1x^2
(x-1)^3=-1+3x-3x^2+1x^3
(x-1)^4=+1-4x+6x^2-4x^3+1x^4
(x-1)^5=-1+5x-10x^2+10x^3-5x^4+1x^5
(x-1)^6=+1-6x+15x^2-20x^3+15x^4-6x^5+1x^6
(x-1)^7=-1+7x-21x^2+35x^3-35x^4+21x^5-7x^6+1x^7

1 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47

Function maxes out at i = 61 as AutoHotkey supports up to 64-bit signed integers.

Bracmat

Bracmat automatically normalizes symbolic expressions with the algebraic binary operators +, *, ^ and \L (logartithm). It can differentiate such expressions using the \D binary operator. (These operators were implemented in Bracmat before all other operators!). Some algebraic values can exist in two evaluated forms. The equivalent x*(a+b) and x*a+x*b are both considered "normal", but x*(a+b)+-1 is not, and therefore expanded to -1+a*x+b*x. This is used in the forceExpansion function to convert e.g. x*(a+b) to x*a+x*b.

The primality test uses a pattern that looks for a fractional factor. If such a factor is found, the test fails. Otherwise it succeeds. <lang bracmat>( (forceExpansion=.1+!arg+-1) & (expandx-1P=.forceExpansion$((x+-1)^!arg)) & ( isPrime

 =
   .         forceExpansion
           $ (!arg^-1*(expandx-1P$!arg+-1*(x^!arg+-1)))
         : ?+/*?+?
       & ~`
     |
 )

& out$"Polynomial representations of (x-1)^p for p <= 7 :" & -1:?n & whl

 ' ( 1+!n:~>7:?n
   & out$(str$("n=" !n ":") expandx-1P$!n)
   )

& 1:?n & :?primes & whl

 ' ( 1+!n:~>50:?n
   & ( isPrime$!n&!primes !n:?primes
     |
     )
   )

& out$"2 <= Primes <= 50:" & out$!primes );</lang> Output:

Polynomial representations of (x-1)^p for p <= 7 :
n=0: 1
n=1: -1+x
n=2: 1+-2*x+x^2
n=3: -1+3*x+-3*x^2+x^3
n=4: 1+-4*x+6*x^2+-4*x^3+x^4
n=5: -1+5*x+-10*x^2+10*x^3+-5*x^4+x^5
n=6: 1+-6*x+15*x^2+-20*x^3+15*x^4+-6*x^5+x^6
  n=7:
    -1
  + 7*x
  + -21*x^2
  + 35*x^3
  + -35*x^4
  + 21*x^5
  + -7*x^6
  + x^7
2 <= Primes <= 50:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47

The AKS test kan be written more concisely than the task describes. This prints the primes between 980 and 1000: <lang bracmat>( out$"Primes between 980 and 1000, short version:" & 980:?n & whl

 ' ( !n+1:<1000:?n
   & ( 1+!n^-1*((x+-1)^!n+-1*(x^!n+-1))+-1:?+/*?+?
     | out$!n
     )
   )

);</lang> Output:

Primes between 980 and 1000, short version:
983
991
997

C

<lang c>#include <stdio.h>

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 (reverse (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 = [1.BigInt, 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)).retro
                .map!q{"%+dx^%d ".format(a[])}.join.replace("x^0", "")
                .replace("^1 ", " ").replace("+", "+ ")
                .replace("-", "- ").replace(" 1x", " x")[2 .. $]);

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

}</lang>

Output:

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

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]

Haskell

<lang haskell>expand 0 = [1] expand p = zipWith (+) (r++[0]) (0:r)

   where r = expand (p-1)

test p | p < 2 = False

      | otherwise = and [mod n p == 0 | n <- init . tail $ expand p]

printPoly [1] = "1" printPoly p = concat [ unwords [pow i, sgn (l-i), show (p!!(i-1))]

                      | i <- reverse [1..l-1] ] where
   l = length p
   sgn i = if even i then "+" else "-"
   pow i = take i "x^" ++ if i > 1 then show i else ""

main = do

   putStrLn "-- p: (x-1)^p for small p"
   putStrLn $ unlines [show i ++ ": " ++ printPoly (expand i) | i <- [0..10]]
   putStrLn "-- Primes up to 100:"
   print (filter test [1..100])</lang>

Output:

-- p: (x-1)^p for small p
0: 1
1: x - 1
2: x^2 - 2x + 1
3: x^3 - 3x^2 + 3x - 1
4: x^4 - 4x^3 + 6x^2 - 4x + 1
5: x^5 - 5x^4 + 10x^3 - 10x^2 + 5x - 1
6: x^6 - 6x^5 + 15x^4 - 20x^3 + 15x^2 - 6x + 1
7: x^7 - 7x^6 + 21x^5 - 35x^4 + 35x^3 - 21x^2 + 7x - 1
8: x^8 - 8x^7 + 28x^6 - 56x^5 + 70x^4 - 56x^3 + 28x^2 - 8x + 1
9: x^9 - 9x^8 + 36x^7 - 84x^6 + 126x^5 - 126x^4 + 84x^3 - 36x^2 + 9x - 1
10: x^10 - 10x^9 + 45x^8 - 120x^7 + 210x^6 - 252x^5 + 210x^4 - 120x^3 + 45x^2 - 10x + 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]

Perl 6

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

sub aks($p where 2..*) { so expansions[$p].[1 ..^ */2].all %% $p }</lang> Testing: <lang perl6>say ' p: (x-1)ᵖ'; say '-----------';

sub super ($n) {

   $n.trans: '0123456789'
          => '⁰¹²³⁴⁵⁶⁷⁸⁹';

}

for ^13 -> $d {

   say $d.fmt('%2i: '), (
       expansions[$d].kv.map: -> $i, $n {
           my $p = $d - $i;
           join , gather {
               take < + - >[$n < 0] ~ ' ' unless $p == $d;
               take $n.abs                unless $p == $d > 0;
               take 'x'                   if $p > 0;
               take super $p - $i         if $p > 1;
           }
       }
   )

}

print "\nPrimes up to 100:\n { grep &aks, 2..100 }\n";</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| %r\n|-' % (p, ' '.join('%s%s' % (('%+i' % e) if (e != 1 or not p or (p and not n) ) else '+', ('x^{%i}' % n) if n else ) for n,e in enumerate(expand_x_1(p))), aks_test(p))) print('|}')</lang>

Extra output rendered by wiki

Polynomial Expansions and AKS prime test
$p$	$(x-1)^{p}$

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

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 ((e (in-range 0 (add1 p))))
   (define sign (expt -1 (- p e)))
   (* sign (binomial p e))))

2. Show the polynomial expansions from p=0 .. 7 (inclusive) (it's possible some of these can be merged...)

(define (format-coefficient c e leftmost?)

 (define (format-c.x^e c e)
   (define +c (abs c))
   (match* (+c e)
     [(_ 0) (format "~a" +c)]
     [(1 _) (format "x^~a" e)]
     [(_ _) (format "~ax^~a" +c e)]))
 (define +/- (if (negative? c) "-" "+"))
 (define +c.x^e (format-c.x^e c e))
 (match* (c e leftmost?)
   [(0 _ _) ""]
   [((? negative?) _ #t) (format "-~a" +c.x^e)]
   [(_ _ #t) +c.x^e]
   [(_ _ _) (format " ~a ~a" +/- +c.x^e)]))

(define (format-polynomial cs)

 (define cs-length (sequence-length cs))
 (apply
  string-append
  (reverse ; convention is to display highest exponent first
   (for/list ((c cs) (e (in-naturals)))
     (format-coefficient c e (= e (sub1 cs-length)))))))

(for ((p (in-range 0 (add1 11))))

 (printf "p=~a: ~a~%" p (format-polynomial (coefficients p))))

3. AKS primeality test

(define (prime?/AKS p)

 (define cs (coefficients p))
 (and
  (or (= (vector-ref cs 0) -1) ; c_0 = -1 -> c_0 - (-1) = 0
      (divides? p 2))        ; c_0 = 1 -> c_0 - (-1) = 2 -> divides?
  (for/and ((c (in-vector cs 1 p))) (divides? p c))))

there is some discussion (see Discussion) about what to do with the perennial "1" case. This is my way of saying that I'm ignoring it

(define lowest-tested-number 2)

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 lowest-tested-number 35)) #:when (prime?/AKS i)) i))

5. stretch goal

all prime numbers under 50

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

Output:

p=0: 1
p=1: x^1 - 1
p=2: x^2 - 2x^1 + 1
p=3: x^3 - 3x^2 + 3x^1 - 1
p=4: x^4 - 4x^3 + 6x^2 - 4x^1 + 1
p=5: x^5 - 5x^4 + 10x^3 - 10x^2 + 5x^1 - 1
p=6: x^6 - 6x^5 + 15x^4 - 20x^3 + 15x^2 - 6x^1 + 1
p=7: x^7 - 7x^6 + 21x^5 - 35x^4 + 35x^3 - 21x^2 + 7x^1 - 1
p=8: x^8 - 8x^7 + 28x^6 - 56x^5 + 70x^4 - 56x^3 + 28x^2 - 8x^1 + 1
p=9: x^9 - 9x^8 + 36x^7 - 84x^6 + 126x^5 - 126x^4 + 84x^3 - 36x^2 + 9x^1 - 1
p=10: x^10 - 10x^9 + 45x^8 - 120x^7 + 210x^6 - 252x^5 + 210x^4 - 120x^3 + 45x^2 - 10x^1 + 1
p=11: x^11 - 11x^10 + 55x^9 - 165x^8 + 330x^7 - 462x^6 + 462x^5 - 330x^4 + 165x^3 - 55x^2 + 11x^1 - 1
(2 3 5 7 11 13 17 19 23 29 31)
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47)
(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)

REXX

<lang rexx>/* REXX ---------------------------------------------------------------

09.02.2014 Walter Pachl
--------------------------------------------------------------------*/

c.=1 Numeric Digits 100 limit=200 pl= Do p=3 To limit

 pm1=p-1
 c.p.1=1
 c.p.p=1
 Do j=2 To p-1
   jm1=j-1
   c.p.j=c.pm1.jm1+c.pm1.j
   End
 End

Say '(x-1)**0 = 1' do i=2 To limit

 im1=i-1
 sign='+'
 ol='(x-1)^'im1 '='
 Do j=i to 2 by -1
   If j=2 Then
     term='x  '
   Else
     term='x^'||(j-1)
   If j=i Then
     ol=ol term
   Else
     ol=ol sign c.i.j'*'term
   sign=translate(sign,'+-','-+')
   End
 If i<10 then
   Say ol sign 1
 Do j=2 To i-1
   If c.i.j//(i-1)>0 Then
     Leave
   End
 If j>i-1 Then
   pl=pl (i-1)
 End

Say 'Primes:' subword(pl,2,27) Say ' ' subword(pl,28) </lang>

Output:

(x-1)**0 = 1
(x-1)^1 = x   - 1
(x-1)^2 = x^2 - 2*x   + 1
(x-1)^3 = x^3 - 3*x^2 + 3*x   - 1
(x-1)^4 = x^4 - 4*x^3 + 6*x^2 - 4*x   + 1
(x-1)^5 = x^5 - 5*x^4 + 10*x^3 - 10*x^2 + 5*x   - 1
(x-1)^6 = x^6 - 6*x^5 + 15*x^4 - 20*x^3 + 15*x^2 - 6*x   + 1
(x-1)^7 = x^7 - 7*x^6 + 21*x^5 - 35*x^4 + 35*x^3 - 21*x^2 + 7*x   - 1
(x-1)^8 = x^8 - 8*x^7 + 28*x^6 - 56*x^5 + 70*x^4 - 56*x^3 + 28*x^2 - 8*x   + 1
Primes: 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 101 103
        103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199

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

Tcl

A recursive method with memoization would be more efficient, but this is sufficient for small-scale work. <lang tcl>proc coeffs {p {signs 1}} {

   set clist 1
   for {set i 0} {$i < $p} {incr i} {

set clist [lmap x [list 0 {*}$clist] y [list {*}$clist 0] { expr {$x + $y} }]

   }
   if {$signs} {

set s -1 set clist [lmap c $clist {expr {[set s [expr {-$s}]] * $c}}]

   }
   return $clist

} proc aksprime {p} {

   if {$p < 2} {

return false

   }
   foreach c [coeffs $p 0] {

if {$c == 1} continue if {$c % $p} { return false }

   }
   return true

}

for {set i 0} {$i <= 7} {incr i} {

   puts -nonewline "(x-1)^$i ="
   set j $i
   foreach c [coeffs $i] {

puts -nonewline [format " %+dx^%d" $c $j] incr j -1

   }
   puts ""

}

set sub35primes {} for {set i 1} {$i < 35} {incr i} {

   if {[aksprime $i]} {

lappend sub35primes $i

} puts "primes under 35: [join $sub35primes ,]"

set sub50primes {} for {set i 1} {$i < 50} {incr i} {

   if {[aksprime $i]} {

lappend sub50primes $i

} puts "primes under 50: [join $sub50primes ,]"</lang>

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
primes under 35: 2,3,5,7,11,13,17,19,23,29,31
primes under 50: 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47