Primalty by Wilson's theorem

From Rosetta Code
Primalty by Wilson's theorem is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
Task

Write a boolean function that tells whether a given integer is prime.

Remember that   1   and all non-positive numbers are not prime.

Use Wilson's theorem.

A number is prime if p divides (p - 1)! + 1.


Ref

Factor[edit]

Works with: Factor version 0.99 2019-10-06
USING: formatting grouping io kernel lists lists.lazy math
math.functions memoize prettyprint sequences ;
 
MEMO: factorial ( m -- n )  ! memoize factorial function
[ 1 ] [ [ 1 - factorial ] [ * ] bi ] if-zero ;
 
: wilson ( n -- ? ) [ 1 - factorial 1 + ] [ divisor? ] bi ;
: prime? ( n -- ? ) dup 2 < [ drop f ] [ wilson ] if ;
: primes ( -- list ) 1 lfrom [ prime? ] lfilter ;
 
"n prime?\n--- -----" print
{ 2 3 9 15 29 37 47 57 67 77 87 97 237 409 659 }
[ dup prime? "%-3d  %u\n" printf ] each nl
 
"First 120 primes via Wilson's theorem:" print
120 primes ltake list>array 20 group simple-table. nl
 
"1000th through 1015th primes:" print
16 primes 999 [ cdr ] times ltake list>array
[ pprint bl ] each nl
Output:
n    prime?
---  -----
2    t
3    t
9    f
15   f
29   t
37   t
47   t
57   f
67   t
77   f
87   f
97   t
237  f
409  t
659  t

First 120 primes via Wilson's theorem:
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 503 509 521 523 541
547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659

1000th through 1015th primes:
7919 7927 7933 7937 7949 7951 7963 7993 8009 8011 8017 8039 8053 8059 8069 8081

Forth[edit]

 
: fac-mod ( n m -- r )
>r 1 swap
begin dup 0 > while
dup rot * [email protected] mod swap 1-
repeat drop rdrop ;
 
: ?prime ( n -- f )
dup 1- tuck swap fac-mod = ;
 
: .primes ( n -- )
cr 2 ?do i ?prime if i . then loop ;
 
Output:
128 .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  ok

Go[edit]

Needless to say, Wilson's theorem is an extremely inefficient way of testing for primalty with 'big integer' arithmetic being needed to compute factorials greater than 20.

Presumably we're not allowed to make any trial divisions here except by the number two where all even positive integers, except two itself, are obviously composite.

package main
 
import (
"fmt"
"math/big"
)
 
var (
zero = big.NewInt(0)
one = big.NewInt(1)
prev = big.NewInt(factorial(20))
)
 
// Only usable for n <= 20.
func factorial(n int64) int64 {
res := int64(1)
for k := n; k > 1; k-- {
res *= k
}
return res
}
 
// If memo == true, stores previous sequential
// factorial calculation for odd n > 21.
func wilson(n int64, memo bool) bool {
if n <= 1 || (n%2 == 0 && n != 2) {
return false
}
if n <= 21 {
return (factorial(n-1)+1)%n == 0
}
b := big.NewInt(n)
r := big.NewInt(0)
z := big.NewInt(0)
if !memo {
z.MulRange(2, n-1) // computes factorial from scratch
} else {
prev.Mul(prev, r.MulRange(n-2, n-1)) // uses previous calculation
z.Set(prev)
}
z.Add(z, one)
return r.Rem(z, b).Cmp(zero) == 0
}
 
func main() {
numbers := []int64{2, 3, 9, 15, 29, 37, 47, 57, 67, 77, 87, 97, 237, 409, 659}
fmt.Println(" n prime")
fmt.Println("--- -----")
for _, n := range numbers {
fmt.Printf("%3d  %t\n", n, wilson(n, false))
}
 
// sequential memoized calculation
fmt.Println("\nThe first 120 prime numbers are:")
for i, count := int64(2), 0; count < 1015; i += 2 {
if wilson(i, true) {
count++
if count <= 120 {
fmt.Printf("%3d ", i)
if count%20 == 0 {
fmt.Println()
}
} else if count >= 1000 {
if count == 1000 {
fmt.Println("\nThe 1,000th to 1,015th prime numbers are:")
}
fmt.Printf("%4d ", i)
}
}
if i == 2 {
i--
}
}
fmt.Println()
}
Output:
  n  prime
---  -----
  2  true
  3  true
  9  false
 15  false
 29  true
 37  true
 47  true
 57  false
 67  true
 77  false
 87  false
 97  true
237  false
409  true
659  true

The first 120 prime numbers are:
  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 503 509 521 523 541 
547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659

The 1,000th to 1,015th prime numbers are:
7919 7927 7933 7937 7949 7951 7963 7993 8009 8011 8017 8039 8053 8059 8069 8081 

Haskell[edit]

import qualified Data.Text as T
import Data.List
 
main = do
putStrLn $ showTable True ' ' '-' ' ' $ ["p","isPrime"]:map (\p -> [show p, show $ isPrime p]) numbers
putStrLn $ "The first 120 prime numbers are:"
putStrLn $ see 20 $ take 120 primes
putStrLn "The 1,000th to 1,015th prime numbers are:"
putStrLn $ see 16.take 16 $ drop 999 primes
 
 
numbers = [2,3,9,15,29,37,47,57,67,77,87,97,237,409,659]
 
primes = [p | p <- 2:[3,5..], isPrime p]
 
isPrime :: Integer -> Bool
isPrime p = if p < 2 then False else 0 == mod (succ $ product [1..pred p]) p
 
bagOf :: Int -> [a] -> [[a]]
bagOf _ [] = []
bagOf n xs = let (us,vs) = splitAt n xs in us : bagOf n vs
 
see :: Show a => Int -> [a] -> String
see n = unlines.map unwords.bagOf n.map (T.unpack.T.justifyRight 3 ' '.T.pack.show)
 
showTable::Bool -> Char -> Char -> Char -> [[String]] -> String
showTable _ _ _ _ [] = []
showTable header ver hor sep contents = unlines $ hr:(if header then z:hr:zs else intersperse hr zss) ++ [hr]
where
vss = map (map length) $ contents
ms = map maximum $ transpose vss ::[Int]
hr = concatMap (\ n -> sep : replicate n hor) ms ++ [sep]
top = replicate (length hr) hor
bss = map (\ps -> map (flip replicate ' ') $ zipWith (-) ms ps) $ vss
zss@(z:zs) = zipWith (\us bs -> (concat $ zipWith (\x y -> (ver:x) ++ y) us bs) ++ [ver]) contents bss
Output:
 --- ------- 
 p   isPrime 
 --- ------- 
 2   True    
 3   True    
 9   False   
 15  False   
 29  True    
 37  True    
 47  True    
 57  False   
 67  True    
 77  False   
 87  False   
 97  True    
 237 False   
 409 True    
 659 True    
 --- ------- 

The first 120 prime numbers are:
  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 503 509 521 523 541
547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659

The 1,000th to 1,015th prime numbers are:
7919 7927 7933 7937 7949 7951 7963 7993 8009 8011 8017 8039 8053 8059 8069 8081

Java[edit]

Wilson's theorem is an extremely inefficient way of testing for primality. As a result, optimizations such as caching factorials not performed.

 
import java.math.BigInteger;
 
public class PrimaltyByWilsonsTheorem {
 
public static void main(String[] args) {
System.out.printf("Primes less than 100 testing by Wilson's Theorem%n");
for ( int i = 0 ; i <= 100 ; i++ ) {
if ( isPrime(i) ) {
System.out.printf("%d ", i);
}
}
}
 
 
private static boolean isPrime(long p) {
if ( p <= 1) {
return false;
}
return fact(p-1).add(BigInteger.ONE).mod(BigInteger.valueOf(p)).compareTo(BigInteger.ZERO) == 0;
}
 
private static BigInteger fact(long n) {
BigInteger fact = BigInteger.ONE;
for ( int i = 2 ; i <= n ; i++ ) {
fact = fact.multiply(BigInteger.valueOf(i));
}
return fact;
}
 
}
 
Output:
Primes less than 100 testing by Wilson's Theorem
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 

Julia[edit]

iswilsonprime(p) = (p < 2 || (p > 2 && iseven(p))) ? false : foldr((x, y) -> (x * y) % p, 1:p - 1) == p - 1
 
wilsonprimesbetween(n, m) = [i for i in n:m if iswilsonprime(i)]
 
println("First 120 Wilson primes: ", wilsonprimesbetween(1, 1000)[1:120])
println("\nThe first 40 Wilson primes above 7900 are: ", wilsonprimesbetween(7900, 9000)[1:40])
 
Output:
First 120 Wilson 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, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659]

The first 40 Wilson primes above 7900 are: [7901, 7907, 7919, 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017, 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219, 8221, 8231, 8233, 8237, 8243, 8263, 8269]

Perl[edit]

Library: ntheory
use strict;
use warnings;
use feature 'say';
use ntheory qw(factorial);
 
my($ends_in_7, $ends_in_3);
 
sub is_wilson_prime {
my($n) = @_;
$n > 1 or return 0;
(factorial($n-1) % $n) == ($n-1) ? 1 : 0;
}
 
for (0..50) {
my $m = 3 + 10 * $_;
$ends_in_3 .= "$m " if is_wilson_prime($m);
my $n = 7 + 10 * $_;
$ends_in_7 .= "$n " if is_wilson_prime($n);
}
 
say $ends_in_3;
say $ends_in_7;
Output:
3 13 23 43 53 73 83 103 113 163 173 193 223 233 263 283 293 313 353 373 383 433 443 463 503
7 17 37 47 67 97 107 127 137 157 167 197 227 257 277 307 317 337 347 367 397 457 467 487

Phix[edit]

include mpfr.e
 
function wilson(integer p)
mpz f = mpz_init()
mpz_fac_ui(f,p-1)
mpz_add_ui(f,f,1)
return mpz_fdiv_ui(f,p)=0
end function
 
atom t0 = time()
sequence primes = {}
integer p = 2
while length(primes)<1015 do
if wilson(p) then
primes &= p
end if
p += 1
end while
printf(1,"The first 25 primes: %v\n",{primes[1..25]})
printf(1," '' builtin: %v\n",{get_primes(-25)})
printf(1,"primes[1000..1015]: %v\n",{primes[1000..1015]})
printf(1," '' builtin: %v\n",{get_primes(-1015)[1000..1015]})
?elapsed(time()-t0) -- this method really is hideously slow!
Output:
The first 25 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}
         '' builtin: {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}
primes[1000..1015]: {7919,7927,7933,7937,7949,7951,7963,7993,8009,8011,8017,8039,8053,8059,8069,8081}
        '' builtin: {7919,7927,7933,7937,7949,7951,7963,7993,8009,8011,8017,8039,8053,8059,8069,8081}
"3.4s"

Python[edit]

No attempt is made to optimise this as this method is a very poor primality test.

from math import factorial
 
def is_wprime(n):
return n > 1 and bool(n == 2 or
(n % 2 and (factorial(n - 1) + 1) % n == 0))
 
if __name__ == '__main__':
c = 100
print(f"Primes under {c}:", end='\n ')
print([n for n in range(c) if is_wprime(n)])
Output:
Primes under 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]

Raku[edit]

(formerly Perl 6)

Works with: Rakudo version 2019.11

Not a particularly recommended way to test for primality, especially for larger numbers. It works, but is slow and memory intensive.

sub postfix:<!> (Int $n) { (constant f = 1, |[\*] 1..*)[$n] }
 
sub is-wilson-prime (Int $p where * > 1) { (($p - 1)! + 1) %% $p }
 
# Pre initialize factorial routine (not thread safe)
9000!;
 
# Testing
put ' p prime?';
printf("%4d  %s\n", $_, .&is-wilson-prime) for 2, 3, 9, 15, 29, 37, 47, 57, 67, 77, 87, 97, 237, 409, 659;
 
my $wilsons = (2,3,*+2*).hyper.grep: *.&is-wilson-prime;
 
put "\nFirst 120 primes:";
put $wilsons[^120].rotor(20)».fmt('%3d').join: "\n";
 
put "\n1000th through 1015th primes:";
put $wilsons[999..1014];
Output:
   p  prime?
   2  True
   3  True
   9  False
  15  False
  29  True
  37  True
  47  True
  57  False
  67  True
  77  False
  87  False
  97  True
 237  False
 409  True
 659  True

First 120 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 503 509 521 523 541
547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659

1000th through 1015th primes:
7919 7927 7933 7937 7949 7951 7963 7993 8009 8011 8017 8039 8053 8059 8069 8081

REXX[edit]

Some effort was made to optimize the factorial computation by using memoization and also minimize the size of the
decimal digit precision     (NUMERIC DIGITS expression).

Also, a "pretty print" was used to align the displaying of a list.

/*REXX pgm tests for primality via Wilson's theorem: a # is prime if p divides (p-1)! +1*/
parse arg LO zz /*obtain optional arguments from the CL*/
if LO=='' | LO=="," then LO= 120 /*Not specified? Then use the default.*/
if zz ='' | zz ="," then zz=2 3 9 15 29 37 47 57 67 77 87 97 237 409 659 /*use default?*/
sw= linesize() - 1; if sw<1 then sw= 79 /*obtain the terminal's screen width. */
digs = digits() /*the current number of decimal digits.*/
#= 0 /*number of (LO) primes found so far.*/
!.= 1 /*placeholder for factorial memoization*/
$= /* " to hold a list of primes.*/
do p=1 until #=LO; oDigs= digs /*remember the number of decimal digits*/
 ?= isPrimeW(p) /*test primality using Wilson's theorem*/
if digs>Odigs then numeric digits digs /*use larger number for decimal digits?*/
if \? then iterate /*if not prime, then ignore this number*/
#= # + 1; $= $ p /*bump prime counter; add prime to list*/
end /*p*/
 
call show 'The first ' LO " prime numbers are:"
w= max( length(LO), length(word(reverse(zz),1))) /*used to align the number being tested*/
@is.0= " isn't"; @is.1= 'is' /*2 literals used for display: is/ain't*/
say
do z=1 for words(zz); oDigs= digs /*remember the number of decimal digits*/
p= word(zz, z) /*get a number for user-supplied list. */
 ?= isPrimeW(p) /*test primality using Wilson's theorem*/
if digs>Odigs then numeric digits digs /*use larger number for decimal digits?*/
say right(p, max(w,length(p) ) ) @is.? "prime."
end /*z*/
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
isPrimeW: procedure expose !. digs; parse arg x '' -1 last;  != 1; xm= x - 1
if x<2 then return 0 /*is the number too small to be prime? */
if x==2 | x==5 then return 1 /*is the number a two or a five? */
if last//2==0 | last==5 then return 0 /*is the last decimal digit even or 5? */
if !.xm\==1 then != !.xm /*has the factorial been pre-computed? */
else do; if xm>!.0 then do; base= !.0+1; _= !.0;  != !._; end
else base= 2 /* [↑] use shortcut.*/
do j=!.0+1 to xm;  != ! * j /*compute factorial.*/
if pos(., !)\==0 then do; parse var ! 'E' expon
numeric digits expon +99
digs = digits()
end /* [↑] has exponent,*/
end /*j*/ /*bump numeric digs.*/
if xm<999 then do; !.xm=!; !.0=xm; end /*assign factorial. */
end /*only save small #s*/
if (!+1)//x==0 then return 1 /*X is a prime.*/
return 0 /*" isn't " " */
/*──────────────────────────────────────────────────────────────────────────────────────*/
show: parse arg header,oo; say header /*display header for the first N primes*/
w= length( word($, LO) ) /*used to align prime numbers in $ list*/
do k=1 for LO; _= right( word($, k), w) /*build list for displaying the primes.*/
if length(oo _)>sw then do; say substr(oo,2); oo=; end /*a line overflowed?*/
oo= oo _ /*display a line. */
end /*k*/ /*does pretty print.*/
if oo\='' then say substr(oo, 2); return /*display residual (if any overflowed).*/

Programming note:   This REXX program makes use of   LINESIZE   REXX program   (or BIF)   which is used to determine the screen width
(or linesize)   of the terminal (console).   Some REXXes don't have this BIF.

The   LINESIZE.REX   REXX program is included here   ───►   LINESIZE.REX.


output   when using the default inputs:
The first  120  prime numbers are:
  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 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659

  2 is prime.
  3 is prime.
  9             isn't prime.
 15             isn't prime.
 29 is prime.
 37 is prime.
 47 is prime.
 57             isn't prime.
 67 is prime.
 77             isn't prime.
 87             isn't prime.
 97 is prime.
237             isn't prime.
409 is prime.
659 is prime.

Sidef[edit]

func is_wilson_prime_slow(n) {
n > 1 || return false
(n-1)! % n == n-1
}
 
func is_wilson_prime_fast(n) {
n > 1 || return false
factorialmod(n-1, n) == n-1
}
 
say 25.by(is_wilson_prime_slow) #=> [2, 3, 5, ..., 83, 89, 97]
say 25.by(is_wilson_prime_fast) #=> [2, 3, 5, ..., 83, 89, 97]
 
say is_wilson_prime_fast(2**43 - 1) #=> false
say is_wilson_prime_fast(2**61 - 1) #=> true

zkl[edit]

Library: GMP
GNU Multiple Precision Arithmetic Library and primes
var [const] BI=Import("zklBigNum");  // libGMP
fcn isWilsonPrime(p){
if(p<=1 or (p%2==0 and p!=2)) return(False);
BI(p-1).factorial().add(1).mod(p) == 0
}
fcn wPrimesW{ [2..].tweak(fcn(n){ isWilsonPrime(n) and n or Void.Skip }) }
numbers:=T(2, 3, 9, 15, 29, 37, 47, 57, 67, 77, 87, 97, 237, 409, 659);
println(" n prime");
println("--- -----");
foreach n in (numbers){ println("%3d  %s".fmt(n, isWilsonPrime(n))) }
 
println("\nFirst 120 primes via Wilson's theorem:");
wPrimesW().walk(120).pump(Void, T(Void.Read,15,False),
fcn(ns){ vm.arglist.apply("%4d".fmt).concat(" ").println() });
 
println("\nThe 1,000th to 1,015th prime numbers are:");
wPrimesW().drop(999).walk(15).concat(" ").println();
Output:
  n  prime
---  -----
  2  True
  3  True
  9  False
 15  False
 29  True
 37  True
 47  True
 57  False
 67  True
 77  False
 87  False
 97  True
237  False
409  True
659  True

First 120 primes via Wilson's theorem:
   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  503
 509  521  523  541  547  557  563  569  571  577  587  593  599  601  607  613
 617  619  631  641  643  647  653  659

The 1,000th to 1,015th prime numbers are:
7919 7927 7933 7937 7949 7951 7963 7993 8009 8011 8017 8039 8053 8059 8069