AKS test for primes: Difference between revisions
m →{{header|J}}: minor formatting & alignment |
m →{{header|J}}: whitespace in code |
||
Line 777: | Line 777: | ||
'''Solution''':<lang j> binomialExpansion =: (!~ * _1 ^ 2 | ]) i.&.:<: NB. 1) Create a function that gives the coefficients of (x-1)^p. |
'''Solution''':<lang j> binomialExpansion =: (!~ * _1 ^ 2 | ]) i.&.:<: NB. 1) Create a function that gives the coefficients of (x-1)^p. |
||
testAKS =: 0 *. |
testAKS =: 0 *./ .= ] | binomialExpansion NB. 3) Use that function to create another which determines whether p is prime using AKS.</lang> |
||
'''Examples''':<lang j> binomialExpansion&.> i. 8 NB. 2) show the polynomial expansions p in the range 0 to at 7 inclusive. |
'''Examples''':<lang j> binomialExpansion&.> i. 8 NB. 2) show the polynomial expansions p in the range 0 to at 7 inclusive. |
Revision as of 16:41, 15 May 2014
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
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 can be stated as follows:
- a number is prime if and only 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.
- 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 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 theorem.
- Use your 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.
BASIC
<lang ubasic>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</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 (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
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
Clojure
The *' function is an arbitrary precision multiplication. <lang clojure>(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)))</lang>
- 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)
CoffeeScript
<lang 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}"</lang>
- 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
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 = [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]
Erlang
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.
<lang 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", | {Row, N} <- zip(tl(tl(take(51, pascal()))), seq(2, 50)), primerow(Row, N)).
</lang>
- 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]
Go
<lang 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
}</lang>
- 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
freebasic
<lang 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 </lang>
(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
Haskell
<lang 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])</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]
J
Solution:<lang j> 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.</lang>
Examples:<lang j> 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</lang>
Mathematica
Algebraic manipulation is built into Mathematica, so there's no need to create a function to do (x-1)^p <lang Mathematica>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]</lang>
- 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}
PARI/GP
<lang parigp>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])</lang>
- 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]
Perl
<lang 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";</lang>
- 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]
Perl 6
<lang perl6>constant expansions = [1], [1,-1], -> @prior { [@prior,0 Z- 0,@prior] } ... *;
sub polyprime($p where 2..*) { so expansions[$p].[1 ..^ */2].all %% $p }</lang>
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: <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; [~] 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; } } )
}</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
And testing the function: <lang perl6>print "\nPrimes up to 100:\n { grep &polyprime, 2..100 }\n";</lang>
- 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
Prolog
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: <lang prolog>
prime(P) :- pascal([1,P|Xs]), append(Xs, [1], Rest), forall( member(X,Xs), 0 is X mod P).
</lang> 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.
Pascal Triangle Generator
<lang prolog> % 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).
</lang>
Solutions
Solutions with output from SWI-Prolog:
<lang 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.
</lang>
[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.
<lang prolog> %%% 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).
</lang>
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):
<lang prolog>
prime(N) :- optpascal([1,N|Xs]), forall( member(X,Xs), 0 is X mod N).
</lang>
<lang prolog> %%% 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 </lang>
Python
<lang 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)])</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| %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('|}')</lang>- Output:
Prime(p)? | ||
---|---|---|
False | ||
False | ||
True | ||
True | ||
False | ||
True | ||
False | ||
True | ||
False | ||
False | ||
False | ||
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
- 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' </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 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
Ruby
Using the `polynomial` Rubygem, this can be written directly from the definition in the description:
<lang ruby>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(',')</lang>
Or without the dependency:
<lang ruby>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(',')</lang>
- 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
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
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
zkl
<lang zkl>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(*));</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 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)