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 is prime if all the coefficients of the polynomial expansion of
are divisible by .
For example, trying :
- 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 .
- Use the function to show here the polynomial expansions of 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
<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.
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
<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('! \n| \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
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)
Rust
<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