AKS test for primes: Difference between revisions

* [http://www.youtube.com/watch?v=HvMSRWTE2mI Fool-Proof Test for Primes] - Numberphile (Video). The accuracy of this video is disputed -- at best it is an oversimplification.

<br><br>

 

=={{header|8th}}==

with: a

 

"The primes upto 50 are (via AKS): " . .primes-via-aks cr

 

{{out}}

<pre>

 

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

</pre>

 

=={{header|Ada}}==

 

procedure Test_For_Primes is

Show_Pascal_Triangles;

Show_Primes;

 

{{out}}

The code below uses Algol 68 Genie which provides arbitrary precision arithmetic for LONG LONG modes.

 

BEGIN

COMMENT

print (newline)

END

{{out}}

<pre>

=={{header|ARM Assembly}}==

{{works with|as|Raspberry Pi <br> or android 32 bits with application Termux}}

/* ARM assembly Raspberry PI or android 32 bits */

/* program AKS.s */

/***************************************************/

.include "../affichage.inc"

{{Output}}

<pre>

=={{header|AutoHotkey}}==

{{works with|AutoHotkey L}}

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

pascalstriangle(n=8) ; n rows of Pascal's triangle

}

Msgbox % t

{{out}}

<pre>(x-1)^0=+1

 

The primality test uses a pattern that looks for a fractional factor. If such a factor is found, the test fails. Otherwise it succeeds.

& (expandx-1P=.forceExpansion$((x+-1)^!arg))

& ( isPrime

& out$"2 <= Primes <= 50:"

& out$!primes

Output:

<pre>Polynomial representations of (x-1)^p for p <= 7 :

2 3 5 7 11 13 17 19 23 29 31 37 41 43 47</pre>

The AKS test kan be written more concisely than the task describes. This prints the primes between 980 and 1000:

& 980:?n

& whl

)

)

Output:

<pre>Primes between 980 and 1000, short version:

 

=={{header|C}}==

#include <stdlib.h>

 

putchar('\n');

return 0;

 

The ugly output:

=={{header|C sharp|C#}}==

{{trans|C}}

using System;

public class AksTest

}

}

 

=={{header|C++}}==

{{trans|Pascal}}

#include <iomanip>

#include <iostream>

cout << endl;

}

{{out}}

<pre>

 

=={{header|Clojure}}==

The *' function is an arbitrary precision multiplication.

"kth coefficient of (x - 1)^n"

[n k]

(doseq [n (range 10)] (println n (cs n)))

(println)

{{output}}

<pre>coefficient series n (k[0] .. k[n])

 

=={{header|CoffeeScript}}==

a = []

return () ->

 

console.log ""

{{out}}

<pre>(x - 1)^0 = + 1

 

=={{header|Common Lisp}}==

(cond

((= p 0) #(1))

(loop for i from 0 to 50

do (when (primep i) (format t "~D " i)))

{{out}}

<pre># p: (x-1)^p for small p:

=={{header|Crystal}}==

{{trans|Ruby}}

p.times.reduce([1]) do |ex, _|

end

end

 

puts "\nPrimes below 50:", 50.times.select { |n| prime? n }.join(',')

 

{{out}}

=={{header|D}}==

{{trans|Python}}

 

BigInt[] expandX1(in uint p) pure /*nothrow*/ {

"\nSmall primes using the AKS test:".writeln;

101.iota.filter!aksTest.writeln;

{{out}}

<pre># p: (x-1)^p for small p:

=={{header|EchoLisp}}==

We use the math.lib library and the poly functions to compute and display the required polynomials. A polynomial P(x) = a0 +a1*x + .. an*x^n is a list of coefficients (a0 a1 .... an).

(lib 'math.lib)

;; 1 - x^p : P = (1 0 0 0 ... 0 -1)

(test (lambda(a) (zero? (modulo a p))))) ;; p divides a[i] ?

(apply and (map test P)))) ;; returns #t if true for all a[i]

{{Output}}

(show-them 13) →

p 1 x -1

 

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

 

=={{header|Elena}}==

{{trans|C#}}

singleton AksTest

c[i] := 1l;

c[1 + i] := 1l;

c[j] := c[j - 1] - c[j]

};

public program()

{

AksTest.coef(n);

console.print("Primes:");

if (AksTest.is_prime(n))

{

console.printLine().readChar()

{{out}}

<pre>

=={{header|Elixir}}==

{{trans|Erlang}}

def iterate(f, x), do: fn -> [x | iterate(f, f.(x))] end

end

 

 

{{out}}

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.

 

 

-import(lists, [all/2, seq/2, zip/2]).

[[N || {Row, N} <- zip(tl(tl(take(51, pascal()))), seq(2, 50)),

primerow(Row, N)]]).

{{out}}

<pre>(x - 1)^0 = + 1

 

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

</pre>

 

=={{header|Factor}}==

math.polynomials prettyprint sequences ;

IN: rosetta-code.aks-test

: aks-test ( -- ) part1 nl part2 ;

 

{{out}}

<pre>

 

=={{header|Forth}}==

1 swap 1+ dup 1 ?do over over i - i */ negate swap loop drop ;

 

11 0 ?do i . ." : " i .poly cr loop cr

50 1 ?do i prime? if i . then loop

{{out}}

<pre>

 

=={{header|Fortran}}==

program aks

implicit none

end function

end program aks

{{out}}

<pre>

 

=={{header|FreeBASIC}}==

'UPPER LIMIT OF FREEBASIC ULONGINT GETS PRIMES UP TO 70

Sub string_split(s_in As String,char As String,result() As String)

Next n

 

{{out}}

<pre>(x-1)^1 = -1 +x

 

=={{header|Go}}==

 

import "fmt"

}

return true

{{out}}

<pre>

 

=={{header|Haskell}}==

 

 

putStrLn $ unlines [show i ++ ": " ++ printPoly (expand i) | i <- [0..10]]

putStrLn "-- Primes up to 100:"

{{out}}

<pre>-- p: (x-1)^p for small p

 

=={{header|Idris}}==

 

-- Computes Binomial Coefficients

printExpansions 10

putStrLn "\n-- Primes Up To 100:"

{{out}}

<pre>-- p: (x-1)^p for small p

=={{header|J}}==

 

 

+-++--+----+-------+-----------+---------------+------------------+

|0||_2|_3 3|_4 6 _4|_5 10 _10 5|_6 15 _20 15 _6|_7 21 _35 35 _21 7|

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.

 

=={{header|Java}}==

 

{{trans|C}}

private static final long[] c = new long[64];

 

System.out.println();

}

Output:

<pre>

=={{header|JavaScript}}==

{{trans|CoffeeScript}}

 

pascal = function() {

console.log("");

 

{{out}}

<pre>(x - 1)^0 = + 1

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

Reviewed (ES6):

var cs = []; if (n) while (n--) coef(); return coef

function coef() {

 

document.write('<br>Primes: ');

{{output}}

(x-1)⁰ = +1

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

{{trans|C}}

for (var c=[1], i=0; i<n; c[0]=-c[0], i+=1) {

c[i+1]=1; for (var j=i; j; j-=1) c[j] = c[j-1]-c[j]

document.write('<br>Primes: ');

{{output}}

(x-1)⁰ = +1

easily be adapted to use a BigInt library such as the one at

https://github.com/joelpurra/jq-bigint

# Input: an OptPascal array

# Output: the next OptPascal array (obtained by adding adjacent items,

def pascal: nth(.; pascals);

 

 

'''Task 1:''' "A method to generate the coefficients of (x-1)^p"

def alternate_signs: . as $in

| reduce range(0; length) as $i ([]; . + [$in[$i] * (if $i % 2 == 0 then 1 else -1 end )]);

'''Task 2:''' "Show here the polynomial expansions of (x − 1)^p for p in the range 0 to at least 7, inclusive."

{{out}}

Coefficients for (x - 1)^1: [1,-1]

Coefficients for (x - 1)^2: [1,-2,1]

Coefficients for (x - 1)^5: [1,-5,10,-10,5,-1]

Coefficients for (x - 1)^6: [1,-6,15,-20,15,-6,1]

 

'''Task 3:''' Prime Number Test

 

For brevity, we show here only the relatively efficient solution based on optpascal/0:

. as $N

| if . < 2 then false

else (1+.) | optpascal

| all( .[2:][]; . % $N == 0 )

 

'''Task 4:''' "Use your AKS test to generate a list of all primes under 35."

{{out}}

3

5

23

29

 

'''Task 5:''' "As a stretch goal, generate all primes under 50."

{{out}}

 

=={{header|Julia}}==

'''Task 1'''

 

function polycoefs(n::Int64)

pc = typeof(n)[]

return pc

end

 

Perhaps this should be done with a comprehension, but properly accounting for the sign is tricky in that case.

'''Task 2'''

 

 

function stringpoly(n::Int64)

return st

end

 

Of course this could be simpler, but this produces a nice payoff in typeset equations that do on include extraneous characters (leading pluses and coefficients of 1).

'''Task 3'''

 

function isaksprime(n::Int64)

if n < 2

return true

end

 

'''Task 4'''

 

println("<math>")

println("\\begin{array}{lcl}")

end

println()

 

{{out}}

 

=={{header|Kotlin}}==

 

fun binomial(n: Int, k: Int): Long = when {

println("\nThe prime numbers under 62 are:")

println(primes)

 

{{out}}

The prime numbers under 62 are:

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61]

</pre>

=={{header|Lua}}==

local function coefs(n)

local list = {[0]=1}

if (isprimeaks(i)) then primes[#primes+1] = i end

end

{{out}}

<pre>(x-1)^0 : 1

 

=={{header|Lambdatalk}}==

{require lib_BN} // for big numbers

 

-> 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 . . .

 

=={{header|Liberty BASIC}}==

{{trans|Pascal}}

{{works with|Just BASIC|any}}

global pasTriMax

pasTriMax = 61

wend

end sub

{{out}}

<pre>

=={{header|Maple}}==

Maple handles algebraic manipulation of polynomials natively.

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.

Use <code>polc</code> to implement <code>prime?</code> which does the primality test.

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

This agrees with the built-in primality test <code>isprime</code>:

true

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

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]

 

=={{header|Mathematica}} / {{header|Wolfram Language}}==

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

(x - 1)^( Range[0, 7]) // Expand // TableForm

Print["primes under 50"]

coefflist[p_Integer] := Coefficient[poly[p], x, #] & /@ Range[0, p - 1];

AKSPrimeQ[p_Integer] := (Mod[coefflist[p] , p] // Union) == {0};

{{out}}

<pre>powers of (x-1)

primes under 50

{1, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47}</pre>

 

=={{header|Nim}}==

from math import binom

import strutils

primes.addInt(p)

echo "\nPrimes under 50: ", primes

 

{{out}}

=={{header|Objeck}}==

{{trans|Java}}

@c : static : Int[];

} while (n-- <> 0);

}

 

Output:

{{trans|Clojure}}

Uses [http://github.com/c-cube/gen gen] library for lazy streams and [https://forge.ocamlcore.org/projects/zarith zarith] for arbitrarily sized integers. Runs as is through the [https://github.com/diml/utop utop] REPL.

#require "zarith"

open Z

print_endline ("primes < 50 per AKS: " ^

(Gen.filter aks (range (of_int 2) (of_int 50)) |>

 

{{output}}

 

=={{header|Oforth}}==

 

: nextCoef( prev -- [] )

0 10 for: i [ System.Out "(x-1)^" << i << " = " << coefs( i ) << cr ]

50 seq filter( #prime? ) apply(#.) printcr

 

{{out}}

 

=={{header|PARI/GP}}==

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;

{{out}}

<pre> [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]

=={{header|Pascal}}==

{{works with|Free Pascal}}

const

pasTriMax = 61;

Write(n:3);

WriteLn;

;output:

<pre>

 

=={{header|Perl}}==

use warnings;

# Select one of these lines. Math::BigInt is in core, but quite slow.

print "expansions of (x-1)^p:\n";

print binpoly($_),"\n" for 0..9;

{{out}}

<pre>expansions of (x-1)^p:

 

{{libheader|ntheory}}

# Uncomment next line to see the r and s values used. Set to 2 for more detail.

# prime_set_config(verbose => 1);

{{out}}

<pre>1000000007 1000000009 1000000021 1000000033 1000000087 1000000093 1000000097</pre>

 

=={{header|Phix}}==

{{out}}

<pre>

 

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

</pre>

 

=={{header|PicoLisp}}==

(let D 1

(make

(range 2 50) ) )

 

{{out}}

 

=={{header|PL/I}}==

AKS: procedure options (main, reorder); /* 16 September 2015, derived from Fortran */

 

 

end AKS;

Results obtained:

<pre>

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

 

===Pascal Triangle Generator===

% To generate the n-th row of a Pascal triangle

% pascal(+N, Row)

S is X1 + X2,

add_pairs( [X2, X3|Xs], Ys).

===Solutions===

 

Solutions with output from SWI-Prolog:

 

%%% Task 1: "A method to generate the coefficients of (1-X)^p"

 

% As required by the problem statement, but necessarily very inefficient:

:- between(0, 7, N), coefficients(N, Coefficients), writeln(Coefficients), fail ; true.

<pre>

[1]

</pre>

 

%%% Task 3. Use the previous function in creating [sic]

%%% another function that when given p returns whether p is prime

append(Cs, [_], Coefficients),

forall( member(C, Cs), 0 is C mod N).

 

The following is more efficient (because it relies on optpascal

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.

 

 

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

 

=={{header|PureBasic}}==

Define vzr.b = -1, vzc.b = ~vzr, nMAX.i = 10, n.i , k.i

 

If isPrime(n,pd()) : Print(Str(n)+Space(2)) : EndIf

Next

{{out}}

<pre> 0: + 1

=={{header|Python}}==

 

# This version uses a generator and thus less computations

c =1

# we stop without computing all possible solutions

return False

or equivalently:

if p==2:return True

c=1

c=c*(p-i)//(i+1)

if c%p:return False

alternatively:

ex = [1]

for i in range(p):

 

print('\n# small primes using the aks test')

 

{{out}}

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

 

{| class="wikitable" style="text-align:left;"

|+ Polynomial Expansions and AKS prime test

for n,e in enumerate(expand_x_1(p))),

aks_test(p)))

 

{{out}}

|-

|}

 

=={{header|R}}==

 

i<-2:p-1

l<-unique(factorial(p) / (factorial(p-i) * factorial(i)))

print(noquote("It isn't prime."))

}

 

=={{header|Racket}}==

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

 

(require math/number-theory)

 

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

 

{{out}}

The <tt>polyprime</tt> function pretty much reads like the original description. Is it "so" that the p'th expansion's coefficients are all divisible by p? The <tt>.[1 ..^ */2]</tt> 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 <tt>polyprime</tt> 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. <tt>:-)</tt>

 

 

sub polyprime($p where 2..*) { so expansions[$p].[1 ..^ */2].all %% $p }

# And testing the function:

 

 

{{out}}

=={{header|REXX}}==

===version 1===

* 09.02.2014 Walter Pachl

* 22.02.2014 WP fix 'accounting' problem (courtesy GS)

Say ' ' subword(pl,29)

Say 'Largest coefficient:' mmm

{{out}}

<pre>(x-1)**0 = 1

 

The program determines programmatically the required number of decimal digits (precision) for the large coefficients.

parse arg Z . /*obtain optional argument from the CL.*/

if Z=='' | Z=="," then Z= 200 /*Not specified? Then use the default.*/

if #=='' then #= "none"; return # /*if null, return "none"; else return #*/

/*──────────────────────────────────────────────────────────────────────────────────────*/

{{out|output|text= &nbsp; for task requirement #2, showing thirty-one expressions using as input: &nbsp; <tt> -31 </tt>}}

<br>(Shown at five-sixth size.)

 

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

</pre>

 

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

 

 

def x_minus_1_to_the(p)

end

 

 

Or without the dependency:

 

p.times.inject([1]) do |ex, _|

([0] + ex).zip(ex + [0]).map { |x,y| x - y }

end

 

 

{{out}}

 

=={{header|Rust}}==

let mut coefficients = vec![0i64; k + 1];

coefficients[0] = 1;

print!("{} ", i);

}

{{out}}

<pre>0: [1]

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

 

fn aks_coefficients(k: usize) -> Vec<i64> {

if k == 0 {

}

}

 

=={{header|Scala}}==

 

val pascal = (( Vector(Vector(BigInt(1))) /: (1 to 50)) { (rows, i) =>

val primes = (2 to 50).filter(isPrime)

println

 

{{out}}

=={{header|Scheme}}==

 

;; implement mod m arithmetic with polnomials in x

;; as lists of coefficients, x^0 first.

;; > (filter rosetta-aks-test (iota 50))

;; '(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47)

 

=={{header|Scilab}}==

 

clear

xdel(winsid())

end

disp(R)

 

{{out}}

 

=={{header|Seed7}}==

 

const func array integer: expand_x_1 (in integer: p) is func

ex := expand_x_1(p);

ex[1] +:= 1;

noop;

end for;

end for;

writeln;

 

{{out}}

=={{header|Sidef}}==

{{trans|Perl}}

p >= 2 || return false

for i in (1 .. p>>1) {

say "expansions of (x-1)^p:"

for i in ^10 { say binpoly(i) }

{{out}}

<pre>

Using the [https://ideas.repec.org/c/boc/bocode/s455001.html moremata] library to print the polynomial coefficients. They are in decreasing degree order. To install moremata, type '''ssc install moremata''' in Stata. Since Stata is using double precision floating-point instead of 32 bit integers, the polynomials are exact up to p=54.

 

function pol(n) {

a=J(1,n+1,1)

 

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

 

=={{header|Swift}}==

var result = [Int](count : n+1, repeatedValue : 0)

}

}

 

{{out}}

=={{header|Tcl}}==

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

set clist 1

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

}

}

{{out}}

<pre>

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

</pre>

 

=={{header|uBasic/4tH}}==

Push n : Gosub _coef : Gosub _drop

Print "(x-1)^";n;" = ";

_drop ' ( n --)

If Pop() Endif

{{out}}

<pre>

=={{header|VBA}}==

{{trans|Phix}}

'-- Does not work for primes above 97, which is actually beyond the original task anyway.

'-- Translated from the C version, just about everything is (working) out-by-1, what fun.

Next n

End Sub

(x-1)^0 = 1

(x-1)^1 = x-1

primes (<= 99 ):

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

</pre>

 

=={{header|Wren}}==

{{trans|Go}}

var c = List.filled(p+1, 0)

var r = 1

System.print("\nAll primes under 50:")

for (p in 2..49) if (aks.call(p)) System.write("%(p) ")

 

{{out}}

=={{header|Yabasic}}==

{{trans|Phix}}

// Translated from the C version, just about everything is (working) out-by-1, what fun.

 

end sub

 

 

=={{header|Zig}}==

Uses Zig's compile-time interpreter to create Pascal's triangle at compile-time.

const std = @import("std");

const assert = std.debug.assert;

 

pub fn main() !void {

var i: u6 = 0;

while (i < 8) : (i += 1)

 

try stdout.print("\nThe primes upto 50 (via AKS) are: ", .{});

}

 

const row = binomial(n).?;

var sign: u8 = '+';

var exp = row.len;

for (row) |coef| {

if (exp != row.len)

exp -= 1;

if (coef != 1 or exp == 0)

if (exp >= 1) {

if (exp > 1)

}

sign = if (sign == '+') '-' else '+';

}

}

 

const rmax = 68;

 

const pascal = build: {

@setEvalBranchQuota(100_000);

break :build coefficients;

};

{{Out}}

<pre>

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

</pre>

=={{header|zkl}}==

{{trans|Python}}

fcn expand_x_1(p){

ex := L(BN(1));

println("\n# small primes using the aks test");

{{out}}

<pre>