AKS test for primes

From Rosetta Code
Jump to: navigation, search
Task
AKS test for primes
You are encouraged to solve this task according to the task description, using any language you may know.

The AKS algorithm for testing whether a number is prime is a polynomial-time algorithm based on an elementary theorem about Pascal triangles.

The theorem on which the test is based can be stated as follows:

  • a number p is prime if and only if all the coefficients of the polynomial expansion of
(x − 1)p − (xp − 1)

are divisible by p.

For example, trying p = 3:

(x − 1)3 − (x3 − 1) = (x3 − 3x2 + 3x − 1) − (x3 − 1) = − 3x2 + 3x
And all the coefficients are divisible by 3 so 3 is prime.
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 (x − 1)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 theorem.
  4. Use your 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).
Note

The task given here is related to the elementary theorem, not the actual AKS algorithm. Using the elementary theorem directly as a way of testing for primes is interesting as an exercise but impractical.

References

Contents

[edit] AutoHotkey

Works with: AutoHotkey L
; 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
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.

[edit] 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.

( (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
);

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:

( 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
)
)
);

Output:

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

[edit] 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;
}

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

[edit] Clojure

The *' function is an arbitrary precision multiplication.

(defn c 
"kth coefficient of (x - 1)^n"
[n k]
(/ (apply *' (range n (- n k) -1))
(apply *' (range k 0 -1))
(if (and (even? k) (< k n)) -1 1)))
 
(defn cs
"coefficient series for (x - 1)^n, k=[0..n]"
[n]
(map #(c n %) (range (inc n))))
 
(defn aks? [p] (->> (cs p) rest butlast (every? #(-> % (mod p) zero?))))
 
(println "coefficient series n (k[0] .. k[n])")
(doseq [n (range 10)] (println n (cs n)))
(println)
(println "primes < 50 per AKS:" (filter aks? (range 2 50)))
Output:
coefficient series n (k[0] .. k[n])
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)
8 (-1 8 -28 56 -70 56 -28 8 1)
9 (-1 9 -36 84 -126 126 -84 36 -9 1)

primes < 50 per AKS: (2 3 5 7 11 13 17 19 23 29 31 37 41 43 47)

[edit] CoffeeScript

pascal = () ->
a = []
return () ->
if a.length is 0 then a = [1]
else
b = (a[i] + a[i+1] for i in [0 ... a.length - 1])
a = [1].concat(b).concat [1]
 
show = (a) ->
show_x = (e) ->
switch e
when 0 then ""
when 1 then "x"
else "x^#{e}"
 
degree = a.length - 1
str = "(x - 1)^#{degree} ="
sgn = 1
 
for i in [0...a.length]
str += ' ' + (if sgn > 0 then "+" else "-") + ' ' + a[i] + show_x(degree - i)
sgn = -sgn
 
return str
 
primerow = (row) ->
degree = row.length - 1
row[1 ... degree].every (x) -> x % degree is 0
 
p = pascal()
console.log show p() for i in [0..7]
 
p = pascal()
p(); p() # skip 0 and 1
 
primes = (i+1 for i in [1..49] when primerow p())
 
console.log ""
console.log "The primes upto 50 are: #{primes}"
Output:
(x - 1)^0 = + 1
(x - 1)^1 = + 1x - 1
(x - 1)^2 = + 1x^2 - 2x + 1
(x - 1)^3 = + 1x^3 - 3x^2 + 3x - 1
(x - 1)^4 = + 1x^4 - 4x^3 + 6x^2 - 4x + 1
(x - 1)^5 = + 1x^5 - 5x^4 + 10x^3 - 10x^2 + 5x - 1
(x - 1)^6 = + 1x^6 - 6x^5 + 15x^4 - 20x^3 + 15x^2 - 6x + 1
(x - 1)^7 = + 1x^7 - 7x^6 + 21x^5 - 35x^4 + 35x^3 - 21x^2 + 7x - 1

The primes upto 50 are: 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47

[edit] Common 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 "~%"))
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 

[edit] D

Translation of: Python
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;
}
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]

[edit] Erlang

Translation of: CoffeeScript

The Erlang io module can print out lists of characters with any level of nesting as a flat string. (e.g. ["Er", ["la", ["n"]], "g"] prints as "Erlang") which is useful when constructing the strings to print out for the binomial expansions. The program also shows how lazy lists can be implemented in Erlang.

#! /usr/bin/escript
 
-import(lists, [all/2, seq/2, zip/2]).
 
iterate(F, X) -> fun() -> [X | iterate(F, F(X))] end.
 
take(0, _lazy) -> [];
take(N, Lazy) ->
[Value | Next] = Lazy(),
[Value | take(N-1, Next)].
 
 
pascal() -> iterate(fun (Row) -> [1 | sum_adj(Row)] end, [1]).
 
sum_adj([_] = L) -> L;
sum_adj([A, B | _] = Row) -> [A+B | sum_adj(tl(Row))].
 
 
show_binomial(Row) ->
Degree = length(Row) - 1,
["(x - 1)^", integer_to_list(Degree), " =", binomial_rhs(Row, 1, Degree)].
 
show_x(0) -> "";
show_x(1) -> "x";
show_x(N) -> [$x, $^ | integer_to_list(N)].
 
binomial_rhs([], _, _) -> [];
binomial_rhs([Coef | Coefs], Sgn, Exp) ->
SignChar = if Sgn > 0 -> $+; true -> $- end,
[$ , SignChar, $ , integer_to_list(Coef), show_x(Exp) | binomial_rhs(Coefs, -Sgn, Exp-1)].
 
 
primerow(Row, N) -> all(fun (Coef) -> (Coef =:= 1) or (Coef rem N =:= 0) end, Row).
 
main(_) ->
[io:format("~s~n", [show_binomial(Row)]) || Row <- take(8, pascal())],
io:format("~nThe primes upto 50: ~p~n",
[[N || {Row, N} <- zip(tl(tl(take(51, pascal()))), seq(2, 50)),
primerow(Row, N)]]).
 
Output:
(x - 1)^0 = + 1
(x - 1)^1 = + 1x - 1
(x - 1)^2 = + 1x^2 - 2x + 1
(x - 1)^3 = + 1x^3 - 3x^2 + 3x - 1
(x - 1)^4 = + 1x^4 - 4x^3 + 6x^2 - 4x + 1
(x - 1)^5 = + 1x^5 - 5x^4 + 10x^3 - 10x^2 + 5x - 1
(x - 1)^6 = + 1x^6 - 6x^5 + 15x^4 - 20x^3 + 15x^2 - 6x + 1
(x - 1)^7 = + 1x^7 - 7x^6 + 21x^5 - 35x^4 + 35x^3 - 21x^2 + 7x - 1

The primes upto 50: [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47]

[edit] Go

package main
 
import "fmt"
 
func bc(p int) []int64 {
c := make([]int64, p+1)
r := int64(1)
for i, half := 0, p/2; i <= half; i++ {
c[i] = r
c[p-i] = r
r = r * int64(p-i) / int64(i+1)
}
for i := p - 1; i >= 0; i -= 2 {
c[i] = -c[i]
}
return c
}
 
func main() {
for p := 0; p <= 7; p++ {
fmt.Printf("%d:  %s\n", p, pp(bc(p)))
}
for p := 2; p < 50; p++ {
if aks(p) {
fmt.Print(p, " ")
}
}
fmt.Println()
}
 
var e = []rune("²³⁴⁵⁶⁷")
 
func pp(c []int64) (s string) {
if len(c) == 1 {
return fmt.Sprint(c[0])
}
p := len(c) - 1
if c[p] != 1 {
s = fmt.Sprint(c[p])
}
for i := p; i > 0; i-- {
s += "x"
if i != 1 {
s += string(e[i-2])
}
if d := c[i-1]; d < 0 {
s += fmt.Sprintf(" - %d", -d)
} else {
s += fmt.Sprintf(" + %d", d)
}
}
return
}
 
func aks(p int) bool {
c := bc(p)
c[p]--
c[0]++
for _, d := range c {
if d%int64(p) != 0 {
return false
}
}
return true
}
Output:
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
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 

[edit] freebasic

 
 
'METHOD -- Use the Pascal triangle to retrieve the coefficients
'UPPER LIMIT OF FREEBASIC ULONGINT GETS PRIMES UP TO 70
Sub string_split(s_in As String,char As String,result() As String)
Dim As String s=s_in,var1,var2
Dim As Integer n,pst
#macro split(stri,char,var1,var2)
pst=Instr(stri,char)
var1="":var2=""
If pst<>0 Then
var1=Mid(stri,1,pst-1)
var2=Mid(stri,pst+1)
Else
var1=stri
End If
Redim Preserve result(1 To 1+n-((Len(var1)>0)+(Len(var2)>0)))
result(n+1)=var1
#endmacro
Do
split(s,char,var1,var2):n=n+1:s=var2
Loop Until var2=""
Redim Preserve result(1 To Ubound(result)-1)
End Sub
 
'Get Pascal triangle components
Function pasc(n As Integer,flag As Integer=0) As String
n+=1
Dim As Ulongint V(n):V(1)=1ul
Dim As String s,sign
For r As Integer= 2 To n
s=""
For i As Integer = r To 1 Step -1
V(i) += V(i-1)
If i Mod 2=1 Then sign="" Else sign="-"
s+=sign+Str(V(i))+","
Next i
Next r
If flag Then 'formatted output
Dim As String i,i2,i3,g
Redim As String a(0)
string_split(s,",",a())
For n1 As Integer=1 To Ubound(a)
If Left(a(n1),1)="-" Then sign="" Else sign="+"
If n1=Ubound(a) Then i2="" Else i2=a(n1)
If n1=2 Then i3="x" Else i3="x^"+Str(n1-1)
If n1=1 Then i="":sign=" " Else i=i3
g+=sign+i2+i+" "
Next n1
g="(x-1)^"+Str(n-1)+" = "+g
Return g
End If
Return s
End Function
 
Function isprime(num As Integer) As Integer
Redim As String a(0)
string_split(pasc(num),",",a())
For n As Integer=Lbound(a)+1 To Ubound(a)-1
If (Valulng(Ltrim(a(n),"-"))) Mod num<>0 Then Return 0
Next n
Return -1
End Function
'====================================
'Formatted output
For n As Integer=1 To 9
Print pasc(n,1)
Next n
 
Print
'Limit of Freebasic Ulongint sets about 70 max
Print "Primes up to 70:"
For n As Integer=2 To 70
If isprime(n) Then Print n;
Next n
 
Sleep
 
(x-1)^1 =  -1 +x
(x-1)^2 =  1 -2x +x^2
(x-1)^3 =  -1 +3x -3x^2 +x^3
(x-1)^4 =  1 -4x +6x^2 -4x^3 +x^4
(x-1)^5 =  -1 +5x -10x^2 +10x^3 -5x^4 +x^5
(x-1)^6 =  1 -6x +15x^2 -20x^3 +15x^4 -6x^5 +x^6
(x-1)^7 =  -1 +7x -21x^2 +35x^3 -35x^4 +21x^5 -7x^6 +x^7
(x-1)^8 =  1 -8x +28x^2 -56x^3 +70x^4 -56x^5 +28x^6 -8x^7 +x^8
(x-1)^9 =  -1 +9x -36x^2 +84x^3 -126x^4 +126x^5 -84x^6 +36x^7 -9x^8 +x^9

Primes up to 70:
 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67


[edit] Haskell

expand p = scanl (\z i -> z * (p-i+1) `div` i) 1 [1..p]
 
 
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 <- [l-1,l-2..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])
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]

[edit] J

Solution:
   binomialExpansion =:  (!~ * _1 ^ 2 | ]) i.&.:<:         NB. 1) Create a function that gives the coefficients of (x-1)^p.
testAKS =: 0 *./ .= ] | binomialExpansion NB. 3) Use that function to create another which determines whether p is prime using AKS.
Examples:
   binomialExpansion&.> i. 8   NB.  2) show the polynomial expansions p in the range 0 to at 7 inclusive.
+-++--+----+-------+-----------+---------------+------------------+
|0||_2|_3 3|_4 6 _4|_5 10 _10 5|_6 15 _20 15 _6|_7 21 _35 35 _21 7|
+-++--+----+-------+-----------+---------------+------------------+
(#~ testAKS&> ) 2+i. 35 NB. 4) Generate a list of all primes under 35.
2 3 5 7 11 13 17 19 23 29 31
(#~ testAKS&> ) 2+i. 50 NB. 5) [stretch] Generate all primes under 50
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47
i.&.:(_1&p:) 50 NB. Double-check our results using built-in prime filter.
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47

[edit] Maple

Maple handles algebraic manipulation of polynomials natively.

> for xpr in seq( expand( (x-1)^p ), p = 0 .. 7 ) do print( xpr ) end:
1
 
x - 1
 
2
x - 2 x + 1
 
3 2
x - 3 x + 3 x - 1
 
4 3 2
x - 4 x + 6 x - 4 x + 1
 
5 4 3 2
x - 5 x + 10 x - 10 x + 5 x - 1
 
6 5 4 3 2
x - 6 x + 15 x - 20 x + 15 x - 6 x + 1
 
7 6 5 4 3 2
x - 7 x + 21 x - 35 x + 35 x - 21 x + 7 x - 1
 

To implement the primality test, we write the following procedure that uses the (built-in) polynomial expansion to generate a list of coefficients of the expanded polynomial.

 polc := p -> [coeffs]( expand( (x-1)^p - (x^p-1) ) ):

Use polc to implement prime? which does the primality test.

prime? := n -> n > 1 and {op}( map( modp, polc( n ), n ) ) = {0}

Of course, rather than calling polc, we can inline it, just for the sake of making the whole thing a one-liner (while adding argument type-checking for good measure):

prime? := (n::posint) -> n > 1 and {op}( map( modp, [coeffs]( expand( (x-1)^n - (x^n-1) ) ), n ) ) = {0}

This agrees with the built-in primality test isprime:

> evalb( seq( prime?(i), i = 1 .. 1000 ) = seq( isprime( i ), i = 1 .. 1000 ) );
true
 

Use prime? with the built-in Maple select procedure to pick off the primes up to 50:

> select( prime?, [seq](1..50) );
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
 


[edit] Mathematica

Algebraic manipulation is built into Mathematica, so there's no need to create a function to do (x-1)^p

Print["powers of (x-1)"]
(x - 1)^( Range[0, 7]) // Expand // TableForm
Print["primes under 50"]
poly[p_] := (x - 1)^p - (x^p - 1) // Expand;
coefflist[p_Integer] := Coefficient[poly[p], x, #] & /@ Range[0, p - 1];
AKSPrimeQ[p_Integer] := (Mod[coefflist[p] , p] // Union) == {0};
Select[Range[1, 50], AKSPrimeQ]
Output:
powers of (x-1)
1
-1+x
1-2 x+x^2
-1+3 x-3 x^2+x^3
1-4 x+6 x^2-4 x^3+x^4
-1+5 x-10 x^2+10 x^3-5 x^4+x^5
1-6 x+15 x^2-20 x^3+15 x^4-6 x^5+x^6
-1+7 x-21 x^2+35 x^3-35 x^4+21 x^5-7 x^6+x^7

primes under 50
{1, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47}

[edit] PARI/GP

getPoly(n)=('x-1)^n;
vector(8,n,getPoly(n-1))
AKS_slow(n)=my(P=getPoly(n));for(i=1,n-1,if(polcoeff(P,i)%n,return(0))); 1;
AKS(n)=my(X=('x-1)*Mod(1,n));X^n=='x^n-1;
select(AKS, [1..50])
Output:
 [1, x - 1, x^2 - 2*x + 1, x^3 - 3*x^2 + 3*x - 1, x^4 - 4*x^3 + 6*x^2 - 4*x + 1, x^5 - 5*x^4 + 10*x^3 - 10*x^2 + 5*x - 1, x^6 - 6*x^5 + 15*x^4 - 20*x^3 + 15*x^2 - 6*x + 1, x^7 - 7*x^6 + 21*x^5 - 35*x^4 + 35*x^3 - 21*x^2 + 7*x - 1]
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]

[edit] Perl

use strict;
use warnings;
# Select one of these lines. Math::BigInt is in core, but slower.
use Math::BigInt; sub binomial { Math::BigInt->new(shift)->bnok(shift) }
# use Math::Pari "binomial";
# use Math::Prime::Util "binomial";
 
sub binprime {
my $p = shift;
do { return 0 if binomial($p,$_) % $p } for 1 .. $p-1;
1;
}
sub coef { # For prettier printing
my($n,$e) = @_;
return $n unless $e;
$n = "" if $n==1;
$e==1 ? "${n}x" : "${n}x^$e";
}
sub binpoly {
my $p = shift;
join(" ", coef(1,$p),
map { join("",("+","-")[($p-$_)&1]," ",coef(binomial($p,$_),$_)) }
reverse 0..$p-1 );
}
print "expansions of (x-1)^p:\n";
print binpoly($_),"\n" for 0..9;
print "Primes to 80: [", join(",", grep { binprime($_) } 2..80), "]\n";
Output:
expansions of (x-1)^p:
1
x - 1
x^2 - 2x + 1
x^3 - 3x^2 + 3x - 1
x^4 - 4x^3 + 6x^2 - 4x + 1
x^5 - 5x^4 + 10x^3 - 10x^2 + 5x - 1
x^6 - 6x^5 + 15x^4 - 20x^3 + 15x^2 - 6x + 1
x^7 - 7x^6 + 21x^5 - 35x^4 + 35x^3 - 21x^2 + 7x - 1
x^8 - 8x^7 + 28x^6 - 56x^5 + 70x^4 - 56x^3 + 28x^2 - 8x + 1
x^9 - 9x^8 + 36x^7 - 84x^6 + 126x^5 - 126x^4 + 84x^3 - 36x^2 + 9x - 1
Primes to 80: [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79]

[edit] Perl 6

constant expansions = [1], [1,-1], -> @prior { [@prior,0 Z- 0,@prior] } ... *;
 
sub polyprime($p where 2..*) { so expansions[$p].[1 ..^ */2].all %% $p }

The expansions are generated similarly to how most FP languages generate sequences that resemble Pascal's triangle, using a zipwith meta-operator (Z) with subtraction, applied between two lists that add a 0 on either end to the prior list. Here we define a constant infinite sequence using the ... sequence operator with a "whatever" endpoint. In fact, the second term [1,-1] could have been generated from the first term, but we put it in there for documentation so the reader can see what direction things are going.

The polyprime function pretty much reads like the original description. Is it "so" that the p'th expansion's coefficients are all divisible by p? The .[1 ..^ */2] slice is done simply to weed out divisions by 1 or by factors we've already tested (since the coefficients are symmetrical in terms of divisibility). If we wanted to write polyprime even more idiomatically, we could have made it another infinite constant list that is just a mapping of the first list, but we decided that would just be showing off. :-)

Showing the expansions:

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;
[~] 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;
}
}
)
}
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

And testing the function:

print "\nPrimes up to 100:\n  { grep &polyprime, 2..100 }\n";
Output:
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

[edit] PicoLisp

(de pascal (N)
(let D 1
(make
(for X (inc N)
(link D)
(setq D
(*/ D (- (inc N) X) (- X)) ) ) ) ) )
 
(for (X 0 (> 10 X) (inc X))
(println X '-> (pascal X) ) )
 
(println
(filter
'((X)
(fully
'((Y) (=0 (% Y X)))
(cdr (head -1 (pascal X))) ) )
(range 2 50) ) )
 
(bye)
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)
8 -> (1 -8 28 -56 70 -56 28 -8 1)
9 -> (1 -9 36 -84 126 -126 84 -36 9 -1)
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47)

[edit] Prolog

[edit] Prolog(ue)

The theorem as stated ties together two elementary concepts in mathematics: prime numbers and the Pascal triangle. The simplicity of the connection can be expressed directly in Prolog by the following prime number generator:

 
prime(P) :-
pascal([1,P|Xs]),
append(Xs, [1], Rest),
forall( member(X,Xs), 0 is X mod P).
 

where pascal/1 is a generator of rows of the Pascal triangle, for example as defined below; the other predicates used above are standard.

This solution to the Rosetta Code problems will accordingly focus on the Pascal triangle, but to illustrate a number of points, we shall exploit its symmetry by representing each of its rows by the longest initial non-decreasing segment of that row, as illustrated in the third column of the following table:

Row   Pascal Row   optpascal
1       1           [1]
2      1 1          [1, 1]
3     1 2 1         [1, 2]
4    1 3 3 1        [1, 3, 3]

We shall refer to this condensed representation of a row as an "optpascal list". Using it, we can simplify and improve the above prime number generator by defining it as follows:

 prime(N) :- optpascal([1,N|Xs]), forall( member(X,Xs), 0 is X mod N).

Using SWI-Prolog without modifying any of the memory management parameters, this prime number generator was used to generate all primes up to and including 75,659.

Since Pascal triangles are the foundation of our approach to addressing the specific Rosetta Code problems, we being by defining the generator pascal/2 that is required by the first problem, but we do so by defining it in terms of an efficient generator, optpascal/1.

[edit] Pascal Triangle Generator

 
% To generate the n-th row of a Pascal triangle
% pascal(+N, Row)
pascal(0, [1]).
pascal(N, Row) :-
N > 0, optpascal( [1, N|Xs] ),
!,
pascalize( [1, N|Xs], Row ).
 
pascalize( Opt, Row ) :-
% if Opt ends in a pair, then peel off the pair:
( append(X, [R,R], Opt) -> true ; append(X, [R], Opt) ),
reverse(X, Rs),
append( Opt, Rs, Row ).
 
% optpascal(-X) generates optpascal lines:
optpascal(X) :-
optpascal_successor( [], X).
 
% optpascal_successor(+P, -Q) is true if Q is an optpascal list beneath the optpascal list P:
optpascal_successor(P, Q) :-
optpascal(P, NextP),
(Q = NextP ; optpascal_successor(NextP, Q)).
 
% optpascal(+Row, NextRow) is true if Row and NextRow are adjacent rows in the Pascal triangle.
% optpascal(+Row, NextRow) where the optpascal representation is used
optpascal(X, [1|Y]) :-
add_pairs(X, Y).
 
% add_pairs(+OptPascal, NextOptPascal) is a helper function for optpascal/2.
% Given one OptPascal list, it generates the next by adding adjacent
% items, but if the last two items are unequal, then their sum is
% repeated. This is intended to be a deterministic predicate, and to
% avoid a probable compiler limitation, we therefore use one cut.
add_pairs([], []).
add_pairs([X], [X]).
add_pairs([X,Y], Ans) :-
S is X + Y,
(X = Y -> Ans=[S] ; Ans=[S,S]),
!. % To overcome potential limitation of compiler
 
add_pairs( [X1, X2, X3|Xs], [S|Ys]) :-
S is X1 + X2,
add_pairs( [X2, X3|Xs], Ys).
 

[edit] Solutions

Solutions with output from SWI-Prolog:

 
%%% Task 1: "A method to generate the coefficients of (1-X)^p"
 
coefficients(N, Coefficients) :-
pascal(N, X),
alternate_signs(X, Coefficients).
 
alternate_signs( [], [] ).
alternate_signs( [A], [A] ).
alternate_signs( [A,B | X], [A, MB | Y] ) :-
MB is -B,
alternate_signs(X,Y).
 
%%% Task 2. "Show here the polynomial expansions of (x − 1)p for p in the range 0 to at least 7, inclusive."
 
coefficients(Coefficients) :-
optpascal( Opt),
pascalize( Opt, Row ),
alternate_signs(Row, Coefficients).
 
 
% As required by the problem statement, but necessarily very inefficient:
:- between(0, 7, N), coefficients(N, Coefficients), writeln(Coefficients), fail ; true.
 
[1]
[1,-1]
[1,-2,1]
[1,-3,3,-1]
[1,-4,6,-4,1]
[1,-5,10,-10,5,-1]
[1,-6,15,-20,15,-6,1]
[1,-7,21,-35,35,-21,7,-1]

The following would be more efficient because backtracking saves recomputation:

:- coefficients(Coefficients),
   writeln(Coefficients),
   Coefficients = [_,N|_], N = -7.
 
%%% Task 3. Use the previous function in creating [sic]
%%% another function that when given p returns whether p is prime
%%% using the AKS test.
 
% Even for testing whether a given number, N, is prime,
% this approach is inefficient, but here is a Prolog implementation:
 
prime_test_per_requirements(N) :-
coefficients(N, [1|Coefficients]),
append(Cs, [_], Coefficients),
forall( member(C, Cs), 0 is C mod N).
 

The following is more efficient (because it relies on optpascal lists rather than the full array of coefficients), and more flexible (because it can be used to generate primes without requiring recomputation):

 
prime(N) :- optpascal([1,N|Xs]), forall( member(X,Xs), 0 is X mod N).
 
 
%%% Task 4. Use your AKS test to generate a list of all primes under 35.
 
:- prime(N), (N < 35 -> write(N), write(' '), fail ; nl).
 
% Output: 1 2 3 5 7 11 13 17 19 23 29 31
 
%%% Task 5. As a stretch goal, generate all primes under 50.
 
:- prime(N), (N < 50 -> write(N), write(' '), fail ; nl).
 
% Output: 1 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47
 

[edit] Python

def expand_x_1(p):
ex = [1]
for i in range(p):
ex.append(ex[-1] * -(p-i) / (i+1))
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)])
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]

[edit] Python: Output formatted for wiki

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

print('''
{| class="wikitable" style="text-align:left;"
|+ Polynomial Expansions and AKS prime test
|-
! <math>p</math>
! <math>(x-1)^p</math>
|-'''
)
for p in range(12):
print('! <math>%i</math>\n| <math>%s</math>\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 > 1 else 'x') if n else '')
for n,e in enumerate(expand_x_1(p))),
aks_test(p)))
print('|}')
Output:
Polynomial Expansions and AKS prime test
p (x − 1)p Prime(p)?
0 + 1 False
1 − 1 + x False
2 + 1 − 2x + x2 True
3 − 1 + 3x − 3x2 + x3 True
4 + 1 − 4x + 6x2 − 4x3 + x4 False
5 − 1 + 5x − 10x2 + 10x3 − 5x4 + x5 True
6 + 1 − 6x + 15x2 − 20x3 + 15x4 − 6x5 + x6 False
7 − 1 + 7x − 21x2 + 35x3 − 35x4 + 21x5 − 7x6 + x7 True
8 + 1 − 8x + 28x2 − 56x3 + 70x4 − 56x5 + 28x6 − 8x7 + x8 False
9 − 1 + 9x − 36x2 + 84x3 − 126x4 + 126x5 − 84x6 + 36x7 − 9x8 + x9 False
10 + 1 − 10x + 45x2 − 120x3 + 210x4 − 252x5 + 210x6 − 120x7 + 45x8 − 10x9 + x10 False
11 − 1 + 11x − 55x2 + 165x3 − 330x4 + 462x5 − 462x6 + 330x7 − 165x8 + 55x9 − 11x10 + x11 True

[edit] Racket

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

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

[edit] REXX

[edit] version 1

/* REXX ---------------------------------------------------------------
* 09.02.2014 Walter Pachl
* 22.02.2014 WP fix 'accounting' problem (courtesy GS)
*--------------------------------------------------------------------*/

c.=1
Numeric Digits 100
limit=200
pl=''
mmm=0
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
mmm=max(mmm,c.p.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 ' '
Say 'Primes:' subword(pl,2,27)
Say ' ' subword(pl,29)
Say 'Largest coefficient:' mmm
Say 'This has' length(mmm) 'digits'
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
        107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199
Largest coefficient: 45274257328051640582702088538742081937252294837706668420660
This has 59 digits    

[edit] version 2

This REXX version is an optimized version (of version 1) and modified to address each of the requirements.
The program determines programmatically the required number of digits (precision) for the large coefficients.

/*REXX pgm calculates primes via the Agrawal-Kayal-Saxena primality test*/
parse arg top .; if top=='' then top=200 /*Not specified? Use default.*/
oTop=top; tell=top<0; top=abs(top) /*TOP negative? Show expression.*/
numeric digits max(9,top%3) /*define a dynamic number of digs*/
@.=1; big=1 /*set all coefficients to unity. */
#= /*define a list of prime numbers.*/
do p=3 for top; pm=p-1; pp=p+1 /*PM & PP: used as a convenience.*/
do m=2 to pp%2; mm=m-1 /*calc. coefficients for a power.*/
@.p.m=@.pm.mm + @.pm.m; mh=pp-m /*calculate left side of coeff.*/
@.p.mh=@.p.m /* " right " " " */
if @.p.m>big then big=@.p.m /*This coefficient the biggest? */
end /*m*/ /* [↑] The M DO loop does both*/
end /*p*/ /* sides in the same loop, */
/* saving a bunch of time. */
if tell then say '(x-1)^0: 1' /*maybe show the first expression*/
$.0='-'; $.1="+" /*$.x is the sign to be used.*/
/* [↓] test for primality by ÷ */
do n=2 for top; nh=n%2; nm=n-1 /*create expressions/find primes.*/
do k=3 to nh until @.n.k//nm\==0 /*coefficients divisible by N-1 ?*/
end /*k*/ /* [↑] skip the 1st & 2nd coeff.*/
/* [↓] search for a good coeff. */
if k>nh & nm\==1 & n\==5 then #=# nm /*add a number to the prime list.*/
if \tell then iterate /*¬tell? Don't show expressions.*/
s=1 /*S: is the sign indicator. */
y='(x-1)^'nm": " /*define first part of expression*/
/* [↓] create higher powers 1st.*/
do j=n to 2 by -1; jm=j-1 /*JM is used as a convenience.*/
if j==2 then exp='x' /*if power=1, don't show the pow.*/
else exp='x^'jm /* ··· else show the power with ^*/
if j==n then y=y exp /*no sign for the 1st expression.*/
else y=y $.s @.n.j'∙'exp /*build the expression with sign.*/
s=\s /*flip the sign in the expression*/
end /*j*/ /* [↑] the sign (now) is 0 | 1,*/
/* and is shown as - | + */
say y $.s 1 /*just show first N expressions, */
end /*n*/ /* [↑] ··· but only for neg TOP.*/
say /* [↓] Has TOP a leading +? Show*/
if left(oTop,1)=='+' then say top is 'prime.' /*tell is/isn't. */
else say 'primes:' # /*display prime # list. */
say /* [↓] size of big 'un.*/
say 'Found ' words(#) ' primes and the largest coefficient has' ,
length(big) "decimal digits." /*stick a fork in it, we're done.*/

output for requirement #2, showing twenty expressions using as input:   -20

(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
(x-1)^9:  x^9 - 9∙x^8 + 36∙x^7 - 84∙x^6 + 126∙x^5 - 126∙x^4 + 84∙x^3 - 36∙x^2 + 9∙x - 1
(x-1)^10:  x^10 - 10∙x^9 + 45∙x^8 - 120∙x^7 + 210∙x^6 - 252∙x^5 + 210∙x^4 - 120∙x^3 + 45∙x^2 - 10∙x + 1
(x-1)^11:  x^11 - 11∙x^10 + 55∙x^9 - 165∙x^8 + 330∙x^7 - 462∙x^6 + 462∙x^5 - 330∙x^4 + 165∙x^3 - 55∙x^2 + 11∙x - 1
(x-1)^12:  x^12 - 12∙x^11 + 66∙x^10 - 220∙x^9 + 495∙x^8 - 792∙x^7 + 924∙x^6 - 792∙x^5 + 495∙x^4 - 220∙x^3 + 66∙x^2 - 12∙x + 1
(x-1)^13:  x^13 - 13∙x^12 + 78∙x^11 - 286∙x^10 + 715∙x^9 - 1287∙x^8 + 1716∙x^7 - 1716∙x^6 + 1287∙x^5 - 715∙x^4 + 286∙x^3 - 78∙x^2 + 13∙x - 1
(x-1)^14:  x^14 - 14∙x^13 + 91∙x^12 - 364∙x^11 + 1001∙x^10 - 2002∙x^9 + 3003∙x^8 - 3432∙x^7 + 3003∙x^6 - 2002∙x^5 + 1001∙x^4 - 364∙x^3 + 91∙x^2 - 14∙x + 1
(x-1)^15:  x^15 - 15∙x^14 + 105∙x^13 - 455∙x^12 + 1365∙x^11 - 3003∙x^10 + 5005∙x^9 - 6435∙x^8 + 6435∙x^7 - 5005∙x^6 + 3003∙x^5 - 1365∙x^4 + 455∙x^3 - 105∙x^2 + 15∙x - 1
(x-1)^16:  x^16 - 16∙x^15 + 120∙x^14 - 560∙x^13 + 1820∙x^12 - 4368∙x^11 + 8008∙x^10 - 11440∙x^9 + 12870∙x^8 - 11440∙x^7 + 8008∙x^6 - 4368∙x^5 + 1820∙x^4 - 560∙x^3 + 120∙x^2 - 16∙x + 1
(x-1)^17:  x^17 - 17∙x^16 + 136∙x^15 - 680∙x^14 + 2380∙x^13 - 6188∙x^12 + 12376∙x^11 - 19448∙x^10 + 24310∙x^9 - 24310∙x^8 + 19448∙x^7 - 12376∙x^6 + 6188∙x^5 - 2380∙x^4 + 680∙x^3 - 136∙x^2 + 17∙x - 1
(x-1)^18:  x^18 - 18∙x^17 + 153∙x^16 - 816∙x^15 + 3060∙x^14 - 8568∙x^13 + 18564∙x^12 - 31824∙x^11 + 43758∙x^10 - 48620∙x^9 + 43758∙x^8 - 31824∙x^7 + 18564∙x^6 - 8568∙x^5 + 3060∙x^4 - 816∙x^3 + 153∙x^2 - 18∙x + 1
(x-1)^19:  x^19 - 19∙x^18 + 171∙x^17 - 969∙x^16 + 3876∙x^15 - 11628∙x^14 + 27132∙x^13 - 50388∙x^12 + 75582∙x^11 - 92378∙x^10 + 92378∙x^9 - 75582∙x^8 + 50388∙x^7 - 27132∙x^6 + 11628∙x^5 - 3876∙x^4 + 969∙x^3 - 171∙x^2 + 19∙x - 1
(x-1)^20:  x^20 - 20∙x^19 + 190∙x^18 - 1140∙x^17 + 4845∙x^16 - 15504∙x^15 + 38760∙x^14 - 77520∙x^13 + 125970∙x^12 - 167960∙x^11 + 184756∙x^10 - 167960∙x^9 + 125970∙x^8 - 77520∙x^7 + 38760∙x^6 - 15504∙x^5 + 4845∙x^4 - 1140∙x^3 + 190∙x^2 - 20∙x + 1

primes:  2 3 5 7 11 13 17 19

Found  8  primes and the largest coefficient has 6 decimal digits.

output for requirement #3, showing if 2221 is prime (or not) using for input:   +2221

(Output note:   this number is really pushing at the limits of REXX's use of virtual memory;   the version of
Regina REXX used herein has a limit of around 2 Gbytes.)

2221 is prime.

Found  331  primes and the largest coefficient has 668 decimal digits.

output for requirement #4, showing all primes under 35 using for input:   35

primes:  2 3 5 7 11 13 17 19 23 29 31

Found  11  primes and the largest coefficient has 10 decimal digits.

output for requirement #5 (stretch goal), showing all primes under 50 using for input:   50

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

Found  15  primes and the largest coefficient has 15 decimal digits.

output when using for input:   500

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 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499

Found  95  primes and the largest coefficient has 150 decimal digits.

[edit] Ruby

Using the `polynomial` Rubygem, this can be written directly from the definition in the description:

require 'polynomial'
 
def x_minus_1_to_the(p)
return Polynomial.new(-1,1)**p
end
 
def prime?(p)
return false if p < 2
(x_minus_1_to_the(p) - Polynomial.from_string("x**#{p}-1")).coefs.all?{|n| n%p==0}
end
 
8.times do |n|
# the default Polynomial#to_s would be OK here; the substitutions just make the
# output match the other version below.
puts "(x-1)^#{n} = #{x_minus_1_to_the(n).to_s.gsub(/\*\*/,'^').gsub(/\*/,'')}"
end
 
puts "\nPrimes below 50:", 50.times.select {|n| prime? n}.join(',')

Or without the dependency:

def x_minus_1_to_the(p)
p.times.inject([1]) do |ex, _|
([0] + ex).zip(ex + [0]).map { |x,y| x - y }
end
end
 
def prime?(p)
return false if p < 2
coeff = x_minus_1_to_the(p)
coeff[0] += coeff.pop
coeff.all?{|n| n%p==0}
end
 
8.times do |n|
puts "(x-1)^#{n} = " +
x_minus_1_to_the(n).
each_with_index.
map { |c, p|
if p.zero? then c.to_s
else
(c<0 ? " - " : " + ") + (c.abs==1 ? "x" : "#{c.abs}x") + (p==1 ? "" : "^#{p}")
end
}.join
end
 
puts "\nPrimes below 50:", 50.times.select {|n| prime? n}.join(',')
Output:
(x-1)^0 = 1
(x-1)^1 = -1 + x
(x-1)^2 = 1 - 2x + x^2
(x-1)^3 = -1 + 3x - 3x^2 + x^3
(x-1)^4 = 1 - 4x + 6x^2 - 4x^3 + x^4
(x-1)^5 = -1 + 5x - 10x^2 + 10x^3 - 5x^4 + x^5
(x-1)^6 = 1 - 6x + 15x^2 - 20x^3 + 15x^4 - 6x^5 + x^6
(x-1)^7 = -1 + 7x - 21x^2 + 35x^3 - 35x^4 + 21x^5 - 7x^6 + x^7

Primes below 50:
2,3,5,7,11,13,17,19,23,29,31,37,41,43,47

[edit] Rust

//Rust 0.12
fn aks_coefficients(k: uint) -> Vec<i64> {
let mut coefficients = Vec::from_elem(k + 1, 0i64);
coefficients[0] = 1;
for i in range(1, k + 1) {
coefficients[i] = -range(1, i).fold(coefficients[0], |prev, j|{
let old = coefficients[j];
coefficients[j] = old - prev;
old
});
}
coefficients
}
 
fn is_prime(p: uint) -> bool {
if p < 2 {
false
} else {
let c = aks_coefficients(p);
std::iter::range_inclusive(1, (c.len() - 1) / 2)
.all(|i| (c[i] % (p as i64)) == 0)
}
}
 
fn main() {
for i in range(0u, 8) {
println!("{}: {}", i, aks_coefficients(i));
}
for i in range(1u, 51).filter(|i| is_prime(*i)) {
print!("{} ", i);
}
}
 
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 

An alternative version which computes the coefficients in a more functional but less efficient way.

 
fn aks_coefficients(k: uint) -> Vec<i64> {
if k == 0 {
vec![1i64]
} else {
let zero = Some(0i64);
range(1, k).fold(vec![1i64, -1], |r, _| {
let a = r.iter().chain(zero.iter());
let b = zero.iter().chain(r.iter());
a.zip(b).map(|(x, &y)| x-y).collect()
})
}
}
 

[edit] Scilab

 
clear
xdel(winsid())
 
stacksize('max')
sz=stacksize();
 
n=7; //For the expansion up to power of n
g=50; //For test of primes up to g
 
function X = pascal(g) //Pascal´s triangle
X(1,1)=1; //Zeroth power
X(2,1)=1; //First power
X(2,2)=1;
for q=3:1:g+1 //From second power use this loop
X(q,1)=1;
X(q,q)=1;
for p=2:1:q-1
X(q,p)=X(q-1,p-1)+X(q-1,p);
end
end
endfunction
 
Z=pascal(g); //Generate Pascal's triangle up to g
 
Q(0+1)="(x-1)^0 = 1"; //For nicer display
Q(1+1)="(x-1)^1 = x^1-1"; //For nicer display
 
disp(Q(1))
disp(Q(2))
 
function cf=coef(Z,q,p) //Return coeffiecents for nicer display of expansion without "ones"
if Z(q,p)==1 then
cf="";
else
cf=string(Z(q,p));
end
endfunction
 
for q=3:n+1 //Generate and display the expansions
Q(q)=strcat(["(x-1)^",string(q-1)," = "]);
sing=""; //Sign of coeff.
for p=1:q-1 //Number of coefficients equals power minus 1
Q(q)=strcat([Q(q),sing,coef(Z,q,p),"x^",string(q-p)]);
if sing=="-" then sing="+"; else sing="-"; end
end
Q(q)=strcat([Q(q),sing,string(1)]);
disp(Q(q))
clear Q
end
 
function prime=prime(Z,g)
prime="true";
for p=2:g
if abs(floor(Z(g+1,p)/g)-Z(g+1,p)/g)>0 then
prime="false";
break;
end
end
endfunction
 
R="2"; //For nicer display
for r=3:g
if prime(Z,r)=="true" then
R=strcat([R, ", ",string(r)]);
end
end
disp(R)
 
Output:
(x-1)^0 = 1   
 
 (x-1)^1 = x^1-1   
 
 (x-1)^2 = x^2-2x^1+1   
 
 (x-1)^3 = x^3-3x^2+3x^1-1   
 
 (x-1)^4 = x^4-4x^3+6x^2-4x^1+1   
 
 (x-1)^5 = x^5-5x^4+10x^3-10x^2+5x^1-1   
 
 (x-1)^6 = x^6-6x^5+15x^4-20x^3+15x^2-6x^1+1   
 
 (x-1)^7 = x^7-7x^6+21x^5-35x^4+35x^3-21x^2+7x^1-1   
 
 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47 

[edit] Seed7

$ include "seed7_05.s7i";
 
const func array integer: expand_x_1 (in integer: p) is func
result
var array integer: ex is [] (1);
local
var integer: i is 0;
begin
for i range 0 to p - 1 do
ex := [] (ex[1] * -(p - i) div (i + 1)) & ex;
end for;
end func;
 
const func boolean: aks_test (in integer: p) is func
result
var boolean: aks_test is FALSE;
local
var array integer: ex is 0 times 0;
var integer: idx is 0;
begin
if p >= 2 then
ex := expand_x_1(p);
ex[1] +:= 1;
for idx range 1 to pred(length(ex)) until ex[idx] rem p <> 0 do
noop;
end for;
aks_test := idx = length(ex);
end if;
end func;
 
const proc: main is func
local
var integer: p is 0;
var integer: n is 0;
var integer: e is 0;
begin
writeln("# p: (x-1)^p for small p");
for p range 0 to 11 do
write(p lpad 3 <& ": ");
for n key e range expand_x_1(p) do
write(" ");
if n >= 0 then
write("+");
end if;
write(n);
if e > 1 then
write("x^" <& pred(e));
end if;
end for;
writeln;
end for;
writeln;
writeln("# small primes using the aks test");
for p range 0 to 61 do
if aks_test(p) then
write(p <& " ");
end if;
end for;
writeln;
end func;
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 

[edit] Tcl

A recursive method with memoization would be more efficient, but this is sufficient for small-scale work.

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 ,]"
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

[edit] uBasic/4tH

For n = 0 To 9
Push n : Gosub _coef : Gosub _drop
Print "(x-1)^";n;" = ";
Push n : Gosub _show
Print
Next
 
Print
Print "primes (never mind the 1):";
 
For n = 1 To 34
Push n : Gosub _isprime
If Pop() Then Print " ";n;
Next
 
Print
End
 
' show polynomial expansions
_show ' ( n --)
Do
If @(Tos()) > -1 Then Print "+";
Print @(Tos());"x^";Tos();
While (Tos())
Push Pop() - 1
Loop
 
Gosub _drop
Return
 
' test whether number is a prime
_isprime ' ( n --)
Gosub _coef
 
i = Tos()
@(0) = @(0) + 1
@(i) = @(i) - 1
 
 
Do While (i) * ((@(i) % Tos()) = 0)
i = i - 1
Loop
 
Gosub _drop
Push (i = 0)
Return
 
' generate coefficients
_coef ' ( n -- n)
If (Tos() < 0) + (Tos() > 34) Then End
' gracefully deal with range issue
i = 0
@(i) = 1
 
Do While i < Tos()
j = i
@(j+1) = 1
 
Do While j > 0
@(j) = @(j-1) - @(j)
j = j - 1
Loop
 
@(0) = -@(0)
i = i + 1
Loop
Return
 
' drop a value from the stack
_drop ' ( n --)
If Pop() Endif
Return
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

[edit] zkl

Translation of: Python
var BN=Import("zklBigNum");
fcn expand_x_1(p){
ex := L(BN(1));
foreach i in (p){ ex.append(ex[-1] * -(p-i) / (i+1)) }
return(ex.reverse())
}
fcn aks_test(p){
if (p < 2) return(False);
ex := expand_x_1(p);
ex[0] = ex[0] + 1;
return(not ex[0,-1].filter('%.fp1(p)));
}
println("# p: (x-1)^p for small p");
foreach p in (12){
println("%3d: ".fmt(p),expand_x_1(p).enumerate()
.pump(String,fcn([(n,e)]){"%+d%s ".fmt(e,n and "x^%d".fmt(n) or "")}));
}
 
println("\n# small primes using the aks test");
println([0..110].filter(aks_test).toString(*));
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
L(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,107,109)
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox