Primality by Wilson's theorem
From Rosetta Code
You are encouraged to solve this task according to the task description, using any language you may know.
- Task
Write a boolean function that tells whether a given integer is prime using Wilson's theorem.
By Wilson's theorem, a number p is prime if and only if p divides (p - 1)! + 1.
Remember that 1 and all non-positive integers are not prime.
- See also
Contents
- 1 11l
- 2 8086 Assembly
- 3 Ada
- 4 ALGOL 68
- 5 ALGOL W
- 6 APL
- 7 AppleScript
- 8 Arturo
- 9 AWK
- 10 C
- 11 C#
- 12 C++
- 13 CLU
- 14 Common Lisp
- 15 Cowgol
- 16 D
- 17 Erlang
- 18 F#
- 19 Factor
- 20 Fermat
- 21 Forth
- 22 FreeBASIC
- 23 Fōrmulæ
- 24 GAP
- 25 Go
- 26 Haskell
- 27 J
- 28 Java
- 29 jq
- 30 Julia
- 31 Lua
- 32 Mathematica/Wolfram Language
- 33 Nim
- 34 PARI/GP
- 35 Perl
- 36 Phix
- 37 Plain English
- 38 PL/I
- 39 PL/M
- 40 Polyglot:PL/I and PL/M
- 41 Python
- 42 Quackery
- 43 Raku
- 44 REXX
- 45 Ring
- 46 Ruby
- 47 Rust
- 48 Sidef
- 49 Swift
- 50 Tiny BASIC
- 51 Wren
- 52 zkl
11l[edit]
F is_wprime(Int64 n)
R n > 1 & (n == 2 | (n % 2 & (factorial(n - 1) + 1) % n == 0))
V c = 20
print(‘Primes under #.:’.format(c), end' "\n ")
print((0 .< c).filter(n -> is_wprime(n)))
- Output:
Primes under 20: [2, 3, 5, 7, 11, 13, 17, 19]
8086 Assembly[edit]
cpu 8086
org 100h
section .text
jmp demo
;;; Wilson primality test of CX.
;;; Zero flag set if CX prime. Destroys AX, BX, DX.
wilson: xor ax,ax ; AX will hold intermediate fac-mod value
inc ax
mov bx,cx ; BX = factorial loop counter
dec bx
.loop: mul bx ; DX:AX = AX*BX
div cx ; modulus goes in DX
mov ax,dx
dec bx ; Next value
jnz .loop ; If not zero yet, go again
inc ax ; fac-mod + 1 equal to input?
cmp ax,cx ; Set flags according to result
ret
;;; Demo: print primes under 256
demo: mov cx,2
.loop: call wilson ; Is it prime?
jnz .next ; If not, try next number
mov ax,cx
call print ; Otherwise, print the number
.next: inc cl ; Next number.
jnz .loop ; If <256, try next number
ret
;;; Print value in AX using DOS syscall
print: mov bp,10 ; Divisor
mov bx,numbuf ; Pointer to buffer
.digit: xor dx,dx
div bp ; Divide AX and get digit in DX
add dl,'0' ; Make ASCII
dec bx ; Store in buffer
mov [bx],dl
test ax,ax ; Done yet?
jnz .digit ; If not, get next digit
mov dx,bx ; Print buffer
mov ah,9 ; 9 = MS-DOS syscall to print a string
int 21h
ret
section .data
db '*****' ; Space to hold ASCII number for output
numbuf: db 13,10,'$'
- Output:
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
Ada[edit]
--
-- Determine primality using Wilon's theorem.
-- Uses the approach from Algol W
-- allowing large primes without the use of big numbers.
--
with Ada.Text_IO; use Ada.Text_IO;
procedure Main is
type u_64 is mod 2**64;
package u_64_io is new modular_io (u_64);
use u_64_io;
function Is_Prime (n : u_64) return Boolean is
fact_Mod_n : u_64 := 1;
begin
if n < 2 then
return False;
end if;
for i in 2 .. n - 1 loop
fact_Mod_n := (fact_Mod_n * i) rem n;
end loop;
return fact_Mod_n = n - 1;
end Is_Prime;
num : u_64 := 1;
type cols is mod 12;
count : cols := 0;
begin
while num < 500 loop
if Is_Prime (num) then
if count = 0 then
New_Line;
end if;
Put (Item => num, Width => 6);
count := count + 1;
end if;
num := num + 1;
end loop;
end Main;
- Output:
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
ALGOL 68[edit]
As with many samples on this page, applies the modulo operation at each step in calculating the factorial, to avoid needing large integeres.
BEGIN
# find primes using Wilson's theorem: #
# p is prime if ( ( p - 1 )! + 1 ) mod p = 0 #
# returns true if p is a prime by Wilson's theorem, false otherwise #
# computes the factorial mod p at each stage, so as to #
# allow numbers whose factorial won't fit in 32 bits #
PROC is wilson prime = ( INT p )BOOL:
BEGIN
INT factorial mod p := 1;
FOR i FROM 2 TO p - 1 DO factorial mod p *:= i MODAB p OD;
factorial mod p = p - 1
END # is wilson prime # ;
FOR i TO 100 DO IF is wilson prime( i ) THEN print( ( " ", whole( i, 0 ) ) ) FI OD
END
- Output:
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
ALGOL W[edit]
As with the APL, Tiny BASIC and other samples, this computes the factorials mod p at each multiplication to avoid needing numbers larger than the 32 bit limit.
begin
% find primes using Wilson's theorem: %
% p is prime if ( ( p - 1 )! + 1 ) mod p = 0 %
% returns true if n is a prime by Wilson's theorem, false otherwise %
% computes the factorial mod p at each stage, so as to %
% allow numbers whose factorial won't fit in 32 bits %
logical procedure isWilsonPrime ( integer value n ) ;
begin
integer factorialModN;
factorialModN := 1;
for i := 2 until n - 1 do factorialModN := ( factorialModN * i ) rem n;
factorialModN = n - 1
end isWilsonPrime ;
for i := 1 until 100 do if isWilsonPrime( i ) then writeon( i_w := 1, s_w := 0, " ", i );
end.
- Output:
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
APL[edit]
This version avoids huge intermediate values by calculating the modulus after each step of the factorial multiplication. This is necessary for the function to work correctly with more than the first few numbers.
wilson ← {⍵<2:0 ⋄ (⍵-1)=(⍵|×)/⍳⍵-1}
- Output:
wilson {(⍺⍺¨⍵)/⍵} ⍳200
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
The naive version (using APL's built-in factorial) looks like this:
naiveWilson ← {⍵<2:0 ⋄ 0=⍵|1+!⍵-1}
But due to loss of precision with large floating-point values, it only works correctly up to number 19 even with ⎕CT set to zero:
- Output:
⎕CT←0 ⋄ naiveWilson {(⍺⍺¨⍵)/⍵} ⍳20
2 3 5 7 11 13 17 19 20
AppleScript[edit]
Nominally, the AppleScript solution would be as follows, the 'mod n' at every stage of the factorial being to keep the numbers within the range the language can handle:
on isPrime(n)
if (n < 2) then return false
set f to n - 1
repeat with i from (n - 2) to 2 by -1
set f to f * i mod n
end repeat
return ((f + 1) mod n = 0)
end isPrime
local output, n
set output to {}
repeat with n from 0 to 500
if (isPrime(n)) then set end of output to n
end repeat
output
- Output:
{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}
In fact, though, the modding by n after every multiplication means there are only three possibilities for the final value of f: n - 1 (if n's a prime), 2 (if n's 4), or 0 (if n's any other non-prime). So the test at the end of the handler could be simplified. Another thing is that if f becomes 0 at some point in the repeat, it obviously stays that way for the remaining iterations, so quite a bit of time can be saved by testing for it and returning false immediately if it occurs. And if 2 and its multiples are caught before the repeat, any other non-prime will guarantee a jump out of the handler. Simply reaching the end will mean n's a prime.
It turns out too that false results only occur when multiplying numbers between √n and n - √n and that only multiplying numbers in this range still leads to the correct outcomes. And if this isn't abusing Wilson's theorem enough, multiples of 2 and 3 can be prechecked and omitted from the "factorial" process altogether, much as they can be skipped in tests for primality by trial division:
on isPrime(n)
-- Check for numbers < 2 and 2 & 3 and their multiples.
if (n < 4) then return (n > 1)
if ((n mod 2 = 0) or (n mod 3 = 0)) then return false
-- Only multiply numbers in the range √n -> n - √n that are 1 less and 1 more than multiples of 6,
-- starting with a number that's 1 less than a multiple of 6 and as close as practical to √n.
tell (n ^ 0.5 div 1) to set f to it - (it - 2) mod 6 + 3
repeat with i from f to (n - f - 6) by 6
set f to f * i mod n * (i + 2) mod n
if (f = 0) then return false
end repeat
return true
end isPrime
Arturo[edit]
factorial: function [x]-> product 1..x
wprime?: function [n][
if n < 2 -> return false
zero? mod add factorial sub n 1 1 n
]
print "Primes below 20 via Wilson's theorem:"
print select 1..20 => wprime?
- Output:
Primes below 20 via Wilson's theorem: 2 3 5 7 11 13 17 19
AWK[edit]
# syntax: GAWK -f PRIMALITY_BY_WILSONS_THEOREM.AWK
# converted from FreeBASIC
BEGIN {
start = 2
stop = 200
for (i=start; i<=stop; i++) {
if (is_wilson_prime(i)) {
printf("%5d%1s",i,++count%10?"":"\n")
}
}
printf("\nWilson primality test range %d-%d: %d\n",start,stop,count)
exit(0)
}
function is_wilson_prime(n, fct,i) {
fct = 1
for (i=2; i<=n-1; i++) {
# because (a mod n)*b = (ab mod n)
# it is not necessary to calculate the entire factorial
fct = (fct * i) % n
}
return(fct == n-1)
}
- Output:
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 Wilson primality test range 2-200: 46
C[edit]
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
uint64_t factorial(uint64_t n) {
uint64_t product = 1;
if (n < 2) {
return 1;
}
for (; n > 0; n--) {
uint64_t prev = product;
product *= n;
if (product < prev) {
fprintf(stderr, "Overflowed\n");
return product;
}
}
return product;
}
// uses wilson's theorem
bool isPrime(uint64_t n) {
uint64_t large = factorial(n - 1) + 1;
return (large % n) == 0;
}
int main() {
uint64_t n;
// Can check up to 21, more will require a big integer library
for (n = 2; n < 22; n++) {
printf("Is %llu prime: %d\n", n, isPrime(n));
}
return 0;
}
- Output:
Is 2 prime: 1 Is 3 prime: 1 Is 4 prime: 0 Is 5 prime: 1 Is 6 prime: 0 Is 7 prime: 1 Is 8 prime: 0 Is 9 prime: 0 Is 10 prime: 0 Is 11 prime: 1 Is 12 prime: 0 Is 13 prime: 1 Is 14 prime: 0 Is 15 prime: 0 Is 16 prime: 0 Is 17 prime: 1 Is 18 prime: 0 Is 19 prime: 1 Is 20 prime: 0 Is 21 prime: 0
C#[edit]
Performance comparison to Sieve of Eratosthenes.
using System;
using System.Linq;
using System.Collections;
using static System.Console;
using System.Collections.Generic;
using BI = System.Numerics.BigInteger;
class Program {
// initialization
const int fst = 120, skp = 1000, max = 1015; static double et1, et2; static DateTime st;
static string ms1 = "Wilson's theorem method", ms2 = "Sieve of Eratosthenes method",
fmt = "--- {0} ---\n\nThe first {1} primes are:", fm2 = "{0} prime thru the {1} prime:";
static List<int> lst = new List<int>();
// dumps a chunk of the prime list (lst)
static void Dump(int s, int t, string f) {
foreach (var item in lst.Skip(s).Take(t)) Write(f, item); WriteLine("\n"); }
// returns the ordinal string representation of a number
static string Ord(int x, string fmt = "{0:n0}") {
var y = x % 10; if ((x % 100) / 10 == 10 || y > 3) y = 0;
return string.Format(fmt, x) + "thstndrd".Substring(y << 1, 2); }
// shows the results of one type of prime tabulation
static void ShowOne(string title, ref double et) {
WriteLine(fmt, title, fst); Dump(0, fst, "{0,-3} ");
WriteLine(fm2, Ord(skp), Ord(max)); Dump(skp - 1, max - skp + 1, "{0,4} ");
WriteLine("Time taken: {0}ms\n", et = (DateTime.Now - st).TotalMilliseconds); }
// for stand-alone computation
static BI factorial(int n) { BI res = 1; if (n < 2) return res;
while (n > 0) res *= n--; return res; }
static bool WTisPrimeSA(int n) { return ((factorial(n - 1) + 1) % n) == 0; }
static BI[] facts;
static void initFacts(int n) {
facts = new BI[n]; facts[0] = facts[1] = 1;
for (int i = 1, j = 2; j < n; i = j++)
facts[j] = facts[i] * j; }
static bool WTisPrime(int n) { return ((facts[n - 1] + 1) % n) == 0; }
// end stand-alone
static void Main(string[] args) { st = DateTime.Now;
BI f = 1; for (int n = 2; lst.Count < max; f *= n++) if ((f + 1) % n == 0) lst.Add(n);
ShowOne(ms1, ref et1);
st = DateTime.Now; int lmt = lst.Last(); lst.Clear(); BitArray flags = new BitArray(lmt + 1);
for (int n = 2; n <= lmt; n+=n==2?1:2) if (!flags[n]) {
lst.Add(n); for (int k = n * n, n2=n<<1; k <= lmt; k += n2) flags[k] = true; }
ShowOne(ms2, ref et2);
WriteLine("{0} was {1:0.0} times slower than the {2}.", ms1, et1 / et2, ms2);
// stand-alone computation
WriteLine("\n" + ms1 + " stand-alone computation:");
WriteLine("factorial computed for each item");
st = DateTime.Now;
for (int x = lst[skp - 1]; x <= lst[max - 1]; x++) if (WTisPrimeSA(x)) Write("{0,4} ", x);
WriteLine(); WriteLine("\nTime taken: {0}ms\n", (DateTime.Now - st).TotalMilliseconds);
WriteLine("factorials precomputed up to highest item");
st = DateTime.Now; initFacts(lst[max - 1]);
for (int x = lst[skp - 1]; x <= lst[max - 1]; x++) if (WTisPrime(x)) Write("{0,4} ", x);
WriteLine(); WriteLine("\nTime taken: {0}ms\n", (DateTime.Now - st).TotalMilliseconds);
}
}
- Output @ Tio.run:
--- Wilson's theorem method --- The first 120 primes 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 1,000th prime thru the 1,015th prime: 7919 7927 7933 7937 7949 7951 7963 7993 8009 8011 8017 8039 8053 8059 8069 8081 Time taken: 340.901ms --- Sieve of Eratosthenes method --- The first 120 primes 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 1,000th prime thru the 1,015th prime: 7919 7927 7933 7937 7949 7951 7963 7993 8009 8011 8017 8039 8053 8059 8069 8081 Time taken: 2.118ms Wilson's theorem method was 161.0 times slower than the Sieve of Eratosthenes method. Wilson's theorem method stand-alone computation: factorial computed for each item 7919 7927 7933 7937 7949 7951 7963 7993 8009 8011 8017 8039 8053 8059 8069 8081 Time taken: 11265.2768ms factorials precomputed up to highest item 7919 7927 7933 7937 7949 7951 7963 7993 8009 8011 8017 8039 8053 8059 8069 8081 Time taken: 177.7401ms
The "slow" factor may be different on different processors and programming environments. For example, on Tio.run, the "slow" factor is anywhere between 120 and 180 times slower. Slowness most likely caused by the sluggish BigInteger library usage. The SoE method, although quicker, does consume some memory (due to the flags BitArray). The Wilson's theorem method may consume considerable memory due to the large factorials (the f variable) when computing larger primes.
The Wilson's theorem method is better suited for computing single primes, as the SoE method causes one to compute all the primes up to the desired item. In this C# implementation, a running factorial is maintained to help the Wilson's theorem method be a little more efficient. The stand-alone results show that when having to compute a BigInteger factorial for every primality test, the performance drops off considerably more. The last performance figure illustrates that memoizing the factorials can help when calculating nearby prime numbers.
C++[edit]
#include <iomanip>
#include <iostream>
int factorial_mod(int n, int p) {
int f = 1;
for (; n > 0 && f != 0; --n)
f = (f * n) % p;
return f;
}
bool is_prime(int p) {
return p > 1 && factorial_mod(p - 1, p) == p - 1;
}
int main() {
std::cout << " n | prime?\n------------\n";
std::cout << std::boolalpha;
for (int p : {2, 3, 9, 15, 29, 37, 47, 57, 67, 77, 87, 97, 237, 409, 659})
std::cout << std::setw(3) << p << " | " << is_prime(p) << '\n';
std::cout << "\nFirst 120 primes by Wilson's theorem:\n";
int n = 0, p = 1;
for (; n < 120; ++p) {
if (is_prime(p))
std::cout << std::setw(3) << p << (++n % 20 == 0 ? '\n' : ' ');
}
std::cout << "\n1000th through 1015th primes:\n";
for (int i = 0; n < 1015; ++p) {
if (is_prime(p)) {
if (++n >= 1000)
std::cout << std::setw(4) << p << (++i % 16 == 0 ? '\n' : ' ');
}
}
}
- 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 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 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
CLU[edit]
% Wilson primality test
wilson = proc (n: int) returns (bool)
if n<2 then return (false) end
fac_mod: int := 1
for i: int in int$from_to(2, n-1) do
fac_mod := fac_mod * i // n
end
return (fac_mod + 1 = n)
end wilson
% Print primes up to 100 using Wilson's theorem
start_up = proc ()
po: stream := stream$primary_output()
for i: int in int$from_to(1, 100) do
if wilson(i) then
stream$puts(po, int$unparse(i) || " ")
end
end
stream$putl(po, "")
end start_up
- Output:
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
Common Lisp[edit]
(defun factorial (n)
(if (< n 2) 1 (* n (factorial (1- n)))) )
(defun primep (n)
"Primality test using Wilson's Theorem"
(unless (zerop n)
(zerop (mod (1+ (factorial (1- n))) n)) ))
- Output:
;; Primes under 20: (dotimes (i 20) (when (primep i) (print i))) 1 2 3 5 7 11 13 17 19
Cowgol[edit]
include "cowgol.coh";
# Wilson primality test
sub wilson(n: uint32): (out: uint8) is
out := 0;
if n >= 2 then
var facmod: uint32 := 1;
var ct := n - 1;
while ct > 0 loop
facmod := (facmod * ct) % n;
ct := ct - 1;
end loop;
if facmod + 1 == n then
out := 1;
end if;
end if;
end sub;
# Print primes up to 100 according to Wilson
var i: uint32 := 1;
while i < 100 loop
if wilson(i) == 1 then
print_i32(i);
print_char(' ');
end if;
i := i + 1;
end loop;
print_nl();
- Output:
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
D[edit]
import std.bigint;
import std.stdio;
BigInt fact(long n) {
BigInt f = 1;
for (int i = 2; i <= n; i++) {
f *= i;
}
return f;
}
bool isPrime(long p) {
if (p <= 1) {
return false;
}
return (fact(p - 1) + 1) % p == 0;
}
void main() {
writeln("Primes less than 100 testing by Wilson's Theorem");
foreach (i; 0 .. 101) {
if (isPrime(i)) {
write(i, ' ');
}
}
writeln;
}
- 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
Erlang[edit]
#! /usr/bin/escript
isprime(N) when N < 2 -> false;
isprime(N) when N band 1 =:= 0 -> N =:= 2;
isprime(N) -> fac_mod(N - 1, N) =:= N - 1.
fac_mod(N, M) -> fac_mod(N, M, 1).
fac_mod(1, _, A) -> A;
fac_mod(N, M, A) -> fac_mod(N - 1, M, A*N rem M).
main(_) ->
io:format("The first few primes (via Wilson's theorem) are: ~n~p~n",
[[K || K <- lists:seq(1, 128), isprime(K)]]).
- Output:
The first few primes (via Wilson's theorem) 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]
F#[edit]
// Wilsons theorem. Nigel Galloway: August 11th., 2020
let wP(n,g)=(n+1I)%g=0I
let fN=Seq.unfold(fun(n,g)->Some((n,g),((n*g),(g+1I))))(1I,2I)|>Seq.filter wP
fN|>Seq.take 120|>Seq.iter(fun(_,n)->printf "%A " n);printfn "\n"
fN|>Seq.skip 999|>Seq.take 15|>Seq.iter(fun(_,n)->printf "%A " n);printfn ""
- Output:
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 7919 7927 7933 7937 7949 7951 7963 7993 8009 8011 8017 8039 8053 8059 8069
Factor[edit]
USING: formatting grouping io kernel lists lists.lazy math
math.factorials math.functions prettyprint sequences ;
: 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
Fermat[edit]
Func Wilson(n) = if ((n-1)!+1)|n = 0 then 1 else 0 fi.;
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
FreeBASIC[edit]
function wilson_prime( n as uinteger ) as boolean
dim as uinteger fct=1, i
for i = 2 to n-1
'because (a mod n)*b = (ab mod n)
'it is not necessary to calculate the entire factorial
fct = (fct * i) mod n
next i
if fct = n-1 then return true else return false
end function
for i as uinteger = 2 to 100
if wilson_prime(i) then print i,
next i
- Output:
Primes below 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
Fōrmulæ[edit]
Fōrmulæ programs are not textual, visualization/edition of programs is done showing/manipulating structures but not text. Moreover, there can be multiple visual representations of the same program. Even though it is possible to have textual representation —i.e. XML, JSON— they are intended for storage and transfer purposes more than visualization and edition.
Programs in Fōrmulæ are created/edited online in its website, However they run on execution servers. By default remote servers are used, but they are limited in memory and processing power, since they are intended for demonstration and casual use. A local server can be downloaded and installed, it has no limitations (it runs in your own computer). Because of that, example programs can be fully visualized and edited, but some of them will not run if they require a moderate or heavy computation/memory resources, and no local server is being used.
In this page you can see the program(s) related to this task and their results.
GAP[edit]
# find primes using Wilson's theorem:
# p is prime if ( ( p - 1 )! + 1 ) mod p = 0
isWilsonPrime := function( p )
local fModP, i;
fModP := 1;
for i in [ 2 .. p - 1 ] do fModP := fModP * i; fModP := fModP mod p; od;
return fModP = p - 1;
end; # isWilsonPrime
prime := [];
for i in [ -4 .. 100 ] do if isWilsonPrime( i ) then Add( prime, i ); fi; od;
Display( prime );
- Output:
[ 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 ]
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
J[edit]
wilson=: 0 = (| !&.:<:)
(#~ wilson) x: 2 + i. 30
2 3 5 7 11 13 17 19 23 29 31
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
jq[edit]
Works with jq, subject to the limitations of IEEE 754 64-bit arithmetic.
Works with gojq, which supports unlimited-precision integer arithmetic.
'Adapted from Julia and Nim'
## Compute (n - 1)! mod m.
def facmod($n; $m):
reduce range(2; $n+1) as $k (1; (. * $k) % $m);
def isPrime: .>1 and (facmod(. - 1; .) + 1) % . == 0;
"Prime numbers between 2 and 100:",
[range(2;101) | select (isPrime)],
# Notice that `infinite` can be used as the second argument of `range`:
"First 10 primes after 7900:",
[limit(10; range(7900; infinite) | select(isPrime))]
- Output:
Prime numbers between 2 and 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]
First 10 primes after 7900:
[7901,7907,7919,7927,7933,7937,7949,7951,7963,7993]
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]
Lua[edit]
-- primality by Wilson's theorem
function isWilsonPrime( n )
local fmodp = 1
for i = 2, n - 1 do
fmodp = fmodp * i
fmodp = fmodp % n
end
return fmodp == n - 1
end
for n = -1, 100 do
if isWilsonPrime( n ) then
io.write( " " .. n )
end
end
- Output:
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
Mathematica/Wolfram Language[edit]
ClearAll[WilsonPrimeQ]
WilsonPrimeQ[1] = False;
WilsonPrimeQ[p_Integer] := Divisible[(p - 1)! + 1, p]
Select[Range[100], WilsonPrimeQ]
- Output:
Prime factors up to a 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}
Nim[edit]
import strutils, sugar
proc facmod(n, m: int): int =
## Compute (n - 1)! mod m.
result = 1
for k in 2..n:
result = (result * k) mod m
func isPrime(n: int): bool = (facmod(n - 1, n) + 1) mod n == 0
let primes = collect(newSeq):
for n in 2..100:
if n.isPrime: n
echo "Prime numbers between 2 and 100:"
echo primes.join(" ")
- Output:
Prime numbers between 2 and 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
PARI/GP[edit]
Wilson(n) = prod(i=1,n-1,Mod(i,n))==-1
Perl[edit]
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]
Uses the modulus method to avoid needing gmp, which was in fact about 7 times slower (when calculating the full factorials).
function wilson(integer n) integer facmod = 1 for i=2 to n-1 do facmod = remainder(facmod*i,n) end for return facmod+1=n 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)
- 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}
"0.5s"
Plain English[edit]
To run:
Start up.
Show some primes (via Wilson's theorem).
Wait for the escape key.
Shut down.
The maximum representable factorial is a number equal to 12. \32-bit signed
To show some primes (via Wilson's theorem):
If a counter is past the maximum representable factorial, exit.
If the counter is prime (via Wilson's theorem), write "" then the counter then " " on the console without advancing.
Repeat.
A prime is a number.
A factorial is a number.
To find a factorial of a number:
Put 1 into the factorial.
Loop.
If a counter is past the number, exit.
Multiply the factorial by the counter.
Repeat.
To decide if a number is prime (via Wilson's theorem):
If the number is less than 1, say no.
Find a factorial of the number minus 1. Bump the factorial.
If the factorial is evenly divisible by the number, say yes.
Say no.
- Output:
1 2 3 5 7 11
PL/I[edit]
/* primality by Wilson's theorem */
wilson: procedure options( main );
declare n binary(15)fixed;
isWilsonPrime: procedure( n )returns( bit(1) );
declare n binary(15)fixed;
declare ( fmodp, i ) binary(15)fixed;
fmodp = 1;
do i = 2 to n - 1;
fmodp = mod( fmodp * i, n );
end;
return ( fmodp = n - 1 );
end isWilsonPrime ;
do n = 1 to 100;
if isWilsonPrime( n ) then do;
put edit( n ) ( f(3) );
end;
end;
end wilson ;
- Output:
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
See also #Polyglot:PL/I and PL/M
PL/M[edit]
... under CP/M (or an emulator)100H: /* FIND PRIMES USING WILSON'S THEOREM: */
/* P IS PRIME IF ( ( P - 1 )! + 1 ) MOD P = 0 */
DECLARE FALSE LITERALLY '0';
BDOS: PROCEDURE( FN, ARG ); /* CP/M BDOS SYSTEM CALL */
DECLARE FN BYTE, ARG ADDRESS;
GOTO 5;
END BDOS;
PRINT$CHAR: PROCEDURE( C ); DECLARE C BYTE; CALL BDOS( 2, C ); END;
PRINT$STRING: PROCEDURE( S ); DECLARE S ADDRESS; CALL BDOS( 9, S ); END;
PRINT$NUMBER: PROCEDURE( N );
DECLARE N ADDRESS;
DECLARE V ADDRESS, N$STR( 6 ) BYTE, W BYTE;
V = N;
W = LAST( N$STR );
N$STR( W ) = '$';
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
DO WHILE( ( V := V / 10 ) > 0 );
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
END;
CALL PRINT$STRING( .N$STR( W ) );
END PRINT$NUMBER;
/* RETURNS TRUE IF P IS PRIME BY WILSON'S THEOREM, FALSE OTHERWISE */
/* COMPUTES THE FACTORIAL MOD P AT EACH STAGE, SO AS TO ALLOW */
/* FOR NUMBERS WHOSE FACTORIAL WON'T FIT IN 16 BITS */
IS$WILSON$PRIME: PROCEDURE( P )BYTE;
DECLARE P ADDRESS;
IF P < 2 THEN RETURN FALSE;
ELSE DO;
DECLARE ( I, FACTORIAL$MOD$P ) ADDRESS;
FACTORIAL$MOD$P = 1;
DO I = 2 TO P - 1;
FACTORIAL$MOD$P = ( FACTORIAL$MOD$P * I ) MOD P;
END;
RETURN FACTORIAL$MOD$P = P - 1;
END;
END IS$WILSON$PRIME;
DECLARE I ADDRESS;
DO I = 1 TO 100;
IF IS$WILSON$PRIME( I ) THEN DO;
CALL PRINT$CHAR( ' ' );
CALL PRINT$NUMBER( I );
END;
END;
EOF
- Output:
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
See also #Polyglot:PL/I and PL/M
Polyglot:PL/I and PL/M[edit]
The following Primality by Wilson's theorem solution will run under both PL/M and PL/I.
... under CP/M (or an emulator)Should work with many PL/I implementations.
The PL/I include file "pg.inc" can be found on the Polyglot:PL/I and PL/M page.
Note the use of text in column 81 onwards to hide the PL/I specifics from the PL/M compiler.
/* PRIMALITY BY WILSON'S THEOREM */
wilson_100H: procedure options (main);
/* PL/I DEFINITIONS */
%include 'pg.inc';
/* PL/M DEFINITIONS: CP/M BDOS SYSTEM CALL AND CONSOLE I/O ROUTINES, ETC. */ /*
DECLARE BINARY LITERALLY 'ADDRESS', CHARACTER LITERALLY 'BYTE';
DECLARE SADDR LITERALLY '.', BIT LITERALLY 'BYTE';
BDOSF: PROCEDURE( FN, ARG )BYTE;
DECLARE FN BYTE, ARG ADDRESS; GOTO 5; END;
BDOS: PROCEDURE( FN, ARG ); DECLARE FN BYTE, ARG ADDRESS; GOTO 5; END;
PRCHAR: PROCEDURE( C ); DECLARE C CHARACTER; CALL BDOS( 2, C ); END;
PRNL: PROCEDURE; CALL PRCHAR( 0DH ); CALL PRCHAR( 0AH ); END;
PRNUMBER: PROCEDURE( N );
DECLARE N ADDRESS;
DECLARE V ADDRESS, N$STR( 6 ) BYTE, W BYTE;
N$STR( W := LAST( N$STR ) ) = '$';
N$STR( W := W - 1 ) = '0' + ( ( V := N ) MOD 10 );
DO WHILE( ( V := V / 10 ) > 0 );
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
END;
CALL BDOS( 9, .N$STR( W ) );
END PRNUMBER;
MODF: PROCEDURE( A, B )ADDRESS;
DECLARE ( A, B )ADDRESS;
RETURN( A MOD B );
END MODF;
/* END LANGUAGE DEFINITIONS */
/* TASK */
DECLARE N BINARY;
ISWILSONPRIME: PROCEDURE( N )returns (
BIT )
;
DECLARE N BINARY;
DECLARE ( FMODP, I ) BINARY;
FMODP = 1;
DO I = 2 TO N - 1;
FMODP = MODF( FMODP * I, N );
END;
RETURN ( FMODP = N - 1 );
END ISWILSONPRIME ;
DO N = 1 TO 100;
IF ISWILSONPRIME( N ) THEN DO;
CALL PRCHAR( ' ' );
CALL PRNUMBER( N );
END;
END;
EOF: end wilson_100H ;
- Output:
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[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]
Quackery[edit]
[ 1 swap times [ i 1+ * ] ] is ! ( n --> n )
[ dup 2 < iff
[ drop false ] done
dup 1 - ! 1+
swap mod 0 = ] is prime ( n --> b )
say "Primes less than 500: "
500 times
[ i^ prime if
[ i^ echo sp ] ]
- Output:
Primes less than 500: 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
Raku[edit]
(formerly Perl 6)
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 from 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.
Ring[edit]
load "stdlib.ring"
decimals(0)
limit = 19
for n = 2 to limit
fact = factorial(n-1) + 1
see "Is " + n + " prime: "
if fact % n = 0
see "1" + nl
else
see "0" + nl
ok
next
Output:
Is 2 prime: 1 Is 3 prime: 1 Is 4 prime: 0 Is 5 prime: 1 Is 6 prime: 0 Is 7 prime: 1 Is 8 prime: 0 Is 9 prime: 0 Is 10 prime: 0 Is 11 prime: 1 Is 12 prime: 0 Is 13 prime: 1 Is 14 prime: 0 Is 15 prime: 0 Is 16 prime: 0 Is 17 prime: 1 Is 18 prime: 0 Is 19 prime: 1
Alternative version computing the factorials modulo n so as to avoid overflow.
# primality by Wilson's theorem
limit = 100
for n = 1 to limit
if isWilsonPrime( n )
see " " + n
ok
next n
func isWilsonPrime n
fmodp = 1
for i = 2 to n - 1
fmodp *= i
fmodp %= n
next i
return fmodp = n - 1
- Output:
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
Ruby[edit]
def w_prime?(i)
return false if i < 2
((1..i-1).inject(&:*) + 1) % i == 0
end
p (1..100).select{|n| w_prime?(n) }
- Output:
[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]
Rust[edit]
fn factorial_mod(mut n: u32, p: u32) -> u32 {
let mut f = 1;
while n != 0 && f != 0 {
f = (f * n) % p;
n -= 1;
}
f
}
fn is_prime(p: u32) -> bool {
p > 1 && factorial_mod(p - 1, p) == p - 1
}
fn main() {
println!(" n | prime?\n------------");
for p in vec![2, 3, 9, 15, 29, 37, 47, 57, 67, 77, 87, 97, 237, 409, 659] {
println!("{:>3} | {}", p, is_prime(p));
}
println!("\nFirst 120 primes by Wilson's theorem:");
let mut n = 0;
let mut p = 1;
while n < 120 {
if is_prime(p) {
n += 1;
print!("{:>3}{}", p, if n % 20 == 0 { '\n' } else { ' ' });
}
p += 1;
}
println!("\n1000th through 1015th primes:");
let mut i = 0;
while n < 1015 {
if is_prime(p) {
n += 1;
if n >= 1000 {
i += 1;
print!("{:>3}{}", p, if i % 16 == 0 { '\n' } else { ' ' });
}
}
p += 1;
}
}
- 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 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 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
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
Swift[edit]
Using a BigInt library.
import BigInt
func factorial<T: BinaryInteger>(_ n: T) -> T {
guard n != 0 else {
return 1
}
return stride(from: n, to: 0, by: -1).reduce(1, *)
}
func isWilsonPrime<T: BinaryInteger>(_ n: T) -> Bool {
guard n >= 2 else {
return false
}
return (factorial(n - 1) + 1) % n == 0
}
print((1...100).map({ BigInt($0) }).filter(isWilsonPrime))
- Output:
[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]
Tiny BASIC[edit]
PRINT "Number to test"
INPUT N
IF N < 0 THEN LET N = -N
IF N = 2 THEN GOTO 30
IF N < 2 THEN GOTO 40
LET F = 1
LET J = 1
10 LET J = J + 1
REM exploits the fact that (F mod N)*J = (F*J mod N)
REM to do the factorial without overflowing
LET F = F * J
GOSUB 20
IF J < N - 1 THEN GOTO 10
IF F = N - 1 THEN PRINT "It is prime"
IF F <> N - 1 THEN PRINT "It is not prime"
END
20 REM modulo by repeated subtraction
IF F < N THEN RETURN
LET F = F - N
GOTO 20
30 REM special case N=2
PRINT "It is prime"
END
40 REM zero and one are nonprimes by definition
PRINT "It is not prime"
END
Wren[edit]
Due to a limitation in the size of integers which Wren can handle (2^53-1) and lack of big integer support, we can only reliably demonstrate primality using Wilson's theorem for numbers up to 19.
import "/math" for Int
import "/fmt" for Fmt
var wilson = Fn.new { |p|
if (p < 2) return false
return (Int.factorial(p-1) + 1) % p == 0
}
for (p in 1..19) {
Fmt.print("$2d -> $s", p, wilson.call(p) ? "prime" : "not prime")
}
- Output:
1 -> not prime 2 -> prime 3 -> prime 4 -> not prime 5 -> prime 6 -> not prime 7 -> prime 8 -> not prime 9 -> not prime 10 -> not prime 11 -> prime 12 -> not prime 13 -> prime 14 -> not prime 15 -> not prime 16 -> not prime 17 -> prime 18 -> not prime 19 -> prime
zkl[edit]
GNU Multiple Precision Arithmetic Library and primesvar [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
Categories:
- Programming Tasks
- Prime Numbers
- 11l
- 8086 Assembly
- Ada
- ALGOL 68
- ALGOL W
- APL
- AppleScript
- Arturo
- AWK
- C
- C sharp
- System.Numerics
- C++
- CLU
- Common Lisp
- Cowgol
- D
- Erlang
- F Sharp
- Factor
- Fermat
- Forth
- FreeBASIC
- Fōrmulæ
- GAP
- Go
- Haskell
- J
- Java
- Jq
- Julia
- Lua
- Mathematica
- Wolfram Language
- Nim
- PARI/GP
- Perl
- Ntheory
- Phix
- Plain English
- PL/I
- PL/M
- Polyglot:PL/I and PL/M
- Python
- Quackery
- Raku
- REXX
- Ring
- Ruby
- Rust
- Sidef
- Swift
- Tiny BASIC
- Wren
- Wren-math
- Wren-fmt
- Zkl
- GMP
