View source for Primality by Wilson's theorem

{{task|Prime Numbers}}

;Task:
Write a boolean function that tells whether a given integer is prime using [http://www.cut-the-knot.org/blue/Wilson.shtml Wilson's theorem]. 

By Wilson's theorem, a number '''p''' is prime if and only if '''p''' divides <code>(p - 1)! + 1</code>.

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



;See also:
* [http://www.cut-the-knot.org/blue/Wilson.shtml Cut-the-knot: Wilson's theorem.]
* [[wp:Wilson's_theorem|Wikipedia: Wilson's theorem]]
<br><br>

=={{header|11l}}==
{{trans|Python}}

<syntaxhighlight lang="11l">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)))</syntaxhighlight>

{{out}}
<pre>
Primes under 20:
  [2, 3, 5, 7, 11, 13, 17, 19]
</pre>

=={{header|8086 Assembly}}==
<syntaxhighlight lang="asm">	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,'$'</syntaxhighlight>

{{out}}

<pre style="height: 50ex">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
</pre>

=={{header|Action!}}==
{{Trans|PL/M}}
<syntaxhighlight lang="action!">
;;; returns TRUE(1) if p is prime by Wilson's theorem, FALSE(0) otherwise
;;;         computes the factorial mod p at each stage, so as to allow
;;;         for numbers whose factorial won't fit in 16 bits
BYTE FUNC isWilsonPrime( CARD p )
  CARD i, factorial_mod_p
  BYTE result

  factorial_mod_p = 1
  FOR i = 2 TO p - 1 DO
    factorial_mod_p = ( factorial_mod_p * i ) MOD p
  OD

  IF factorial_mod_p = p - 1 THEN result = 1 ELSE result = 0 FI

RETURN( result )

PROC Main()
  CARD i

  FOR i = 1 TO 100 DO
    IF isWilsonPrime( i ) THEN
       Put(' ) PrintC( i )
    FI
  OD
RETURN
</syntaxhighlight>
{{out}}
<pre>
 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|Ada}}==
<syntaxhighlight lang="ada">--
-- 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;

</syntaxhighlight>
{{output}}
<pre>
     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
</pre>

=={{header|ALGOL 68}}==
{{Trans|ALGOL W}}
As with many samples on this page, applies the modulo operation at each step in calculating the factorial, to avoid needing large integeres.
<syntaxhighlight lang="algol68">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</syntaxhighlight>
{{out}}
<pre>
 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|ALGOL W}}==
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.
<syntaxhighlight lang="algolw">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.</syntaxhighlight>
{{out}}
<pre>
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|APL}}==

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.  

<syntaxhighlight lang="apl">wilson ← {⍵<2:0 ⋄ (⍵-1)=(⍵|×)/⍳⍵-1}</syntaxhighlight>

{{out}}

<pre>      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</pre>

The naive version (using APL's built-in factorial) looks like this:

<syntaxhighlight lang="apl">naiveWilson ← {⍵<2:0 ⋄ 0=⍵|1+!⍵-1}</syntaxhighlight>

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:

{{out}}

<pre>      ⎕CT←0 ⋄ naiveWilson {(⍺⍺¨⍵)/⍵} ⍳20
2 3 5 7 11 13 17 19 20</pre>

=={{header|AppleScript}}==
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:

<syntaxhighlight lang="applescript">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</syntaxhighlight>

{{output}}
<syntaxhighlight lang="applescript">{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}</syntaxhighlight>

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 <tt>false</tt> 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 <tt>false</tt> 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:

<syntaxhighlight lang="applescript">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</syntaxhighlight>

=={{header|Arturo}}==

<syntaxhighlight lang="rebol">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?</syntaxhighlight>

{{out}}

<pre>Primes below 20 via Wilson's theorem:
2 3 5 7 11 13 17 19</pre>
=={{header|AWK}}==
<syntaxhighlight lang="awk">
# 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)
}
</syntaxhighlight>
{{out}}
<pre>
    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
</pre>


=={{header|BASIC}}==
==={{header|Applesoft BASIC}}===
{{works with|Chipmunk Basic}}
{{works with|GW-BASIC}}
{{trans|FreeBASIC}}
<syntaxhighlight lang="qbasic">100 HOME : REM  100 CLS for Chipmunk Basic
110 PRINT "Primes below 100"+CHR$(10)
120 FOR n = 2 TO 100
130  GOSUB 160
140 NEXT n
150 END
160 rem FUNCTION WilsonPrime(n)
170  fct = 1
180  FOR i = 2 TO n-1
181   a = fct * i
190   fct = a - INT(a / n) * n
200  NEXT i
210  IF fct = n-1 THEN PRINT i;"  ";
220 RETURN</syntaxhighlight>

==={{header|BASIC256}}===
{{trans|FreeBASIC}}
<syntaxhighlight lang="freebasic">function wilson_prime(n)
    fct = 1
    for i = 2 to n-1
        fct = (fct * i) mod n
    next i
    if fct = n-1 then return True else return False
end function

print "Primes below 100" & Chr(10)
for i = 2 to 100
    if wilson_prime(i) then print i; "   ";
next i
end</syntaxhighlight>

==={{header|Chipmunk Basic}}===
{{works with|Chipmunk Basic|3.6.4}}
{{trans|FreeBASIC}}
<syntaxhighlight lang="qbasic">100 cls
110 print "Primes below 100"+chr$(10)
120 for i = 2 to 100
130  wilsonprime(i)
140 next i
150 end
160 function wilsonprime(n)
170  fct = 1
180  for i = 2 to n-1
190   fct = (fct*i) mod n
200  next i
210  if fct = n-1 then print i;
220 end function</syntaxhighlight>

==={{header|Craft Basic}}===
<syntaxhighlight lang="basic">for i = 2 to 100

	let f = 1

	for j = 2 to i - 1

		let f = (f * j) % i
		wait

	next j

	if f = i - 1 then

		print i

	endif

next i

end</syntaxhighlight>
{{out| Output}}<pre>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|GW-BASIC}}===
{{works with|Chipmunk Basic}}
{{works with|PC-BASIC|any}}
{{works with|MSX_Basic}}
{{works with|QBasic}}
{{trans|FreeBASIC}}
<syntaxhighlight lang="qbasic">100 CLS : REM  100 CLS for Chipmunk Basic
110 PRINT "Primes below 100"+CHR$(10)
120 FOR N = 2 TO 100
130  GOSUB 160
140 NEXT N
150 END
160 REM FUNCTION WilsonPrime(n)
170  FCT = 1
180  FOR I = 2 TO N-1
190   FCT = (FCT*I) MOD N
200  NEXT I
210  IF FCT = N-1 THEN PRINT I;"   ";
220 RETURN</syntaxhighlight>

==={{header|Minimal BASIC}}===
{{works with|QBasic}}
{{trans|FreeBASIC}}
<syntaxhighlight lang="qbasic">110 PRINT "Primes below 100"
120 FOR n = 2 TO 100
130  GOSUB 160
140 NEXT n
150 GOTO 250
160 rem FUNCTION WilsonPrime(n)
170  LET f = 1
180  FOR i = 2 TO n-1
181   LET a = f * i
190   LET f = a - INT(a / n) * n
200  NEXT i
210  IF f = n-1 THEN 230
220 RETURN
230 PRINT i
240 RETURN
250 END</syntaxhighlight>

==={{header|MSX Basic}}===
The [[#GW-BASIC|GW-BASIC]] solution works without any changes.

==={{header|QBasic}}===
{{works with|QBasic|1.1}}
{{works with|QuickBasic|4.5}}
{{works with|Run BASIC}}
{{trans|FreeBASIC}}
<syntaxhighlight lang="qbasic">FUNCTION wilsonprime(n)
    fct = 1
    FOR i = 2 TO n - 1
        fct = (fct * i) MOD n
    NEXT i
    IF fct = n - 1 THEN wilsonprime = 1 ELSE wilsonprime = 0
END FUNCTION

PRINT "Primes below 100"; CHR$(10)
FOR i = 2 TO 100
    IF wilsonprime(i) THEN PRINT i; "   ";
NEXT i</syntaxhighlight>

==={{header|Quite BASIC}}===
{{works with|QBasic}}
{{trans|FreeBASIC}}
<syntaxhighlight lang="qbasic">100 CLS
110 PRINT "Primes below 100": PRINT
120 FOR n = 2 TO 100
130  GOSUB 160
140 NEXT n
150 GOTO 250
160 rem FUNCTION WilsonPrime(n)
170  LET f = 1
180  FOR i = 2 TO n-1
181   LET a = f * i
190   LET f = a - INT(a / n) * n
200  NEXT i
210  IF f = n-1 THEN 230
220 RETURN
230 PRINT i;"  ";
240 RETURN
250 END</syntaxhighlight>

==={{header|PureBasic}}===
{{trans|FreeBASIC}}
<syntaxhighlight lang="purebasic">Procedure wilson_prime(n.i)
  fct.i = 1
  For i.i = 2 To n-1
    fct = (fct * i) % n
  Next i
  If fct = n-1 
    ProcedureReturn #True 
  Else 
    ProcedureReturn #False
  EndIf
EndProcedure

OpenConsole()
PrintN("Primes below 100")
For i = 2 To 100
  If wilson_prime(i) 
    Print(Str(i) + #TAB$)
  EndIf
Next i
PrintN("")
Input()
CloseConsole()</syntaxhighlight>

==={{header|Run BASIC}}===
{{works with|QBasic}}
{{trans|FreeBASIC}}
<syntaxhighlight lang="runbasic">print "Primes below 100"
for i = 2 to 100
    if wilsonprime(i) = 1 then print i; "   ";
next i
end

function wilsonprime(n)
    fct = 1
    for i = 2 to n-1
        fct = (fct * i) mod n
    next i
    if fct = n-1 then wilsonprime = 1 else wilsonprime = 0
end function</syntaxhighlight>

==={{header|True BASIC}}===
<syntaxhighlight lang="qbasic">FUNCTION wilsonprime(n)
    LET fct = 1
    FOR i = 2 TO n - 1
        LET fct = MOD((fct * i), n)
    NEXT i
    IF fct = n - 1 THEN LET wilsonprime = 1 ELSE LET wilsonprime = 0
END FUNCTION

PRINT "Primes below 100"; CHR$(10)
FOR i = 2 TO 100
    IF wilsonprime(i) = 1 THEN PRINT i; "   ";
NEXT i
END</syntaxhighlight>

==={{header|Yabasic}}===
{{trans|FreeBASIC}}
<syntaxhighlight lang="freebasic">print "Primes below 100\n"
for i = 2 to 100
    if wilson_prime(i)  print i, "   ";
next i

sub wilson_prime(n)
    local i, fct
    
    fct = 1
    for i = 2 to n-1
        fct = mod((fct * i), n)
    next i
    if fct = n-1 then return True else return False : fi
end sub</syntaxhighlight>


=={{header|BCPL}}==
<syntaxhighlight lang="bcpl">get "libhdr"

let wilson(n) = valof
$(  let f = n - 1
    if n < 2 then resultis false
    for i = n-2 to 2 by -1 do
        f := f*i rem n
    resultis (f+1) rem n = 0
$)

let start() be
    for i = 1 to 100 if wilson(i) do
        writef("%N*N", i)</syntaxhighlight>
{{out}}
<pre>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|C}}==
<syntaxhighlight lang="c">#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;
}</syntaxhighlight>
{{out}}
<pre>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</pre>

=={{header|C#|Csharp}}==
{{libheader|System.Numerics}}
Performance comparison to Sieve of Eratosthenes.
<syntaxhighlight lang="csharp">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);
    }
}</syntaxhighlight>
{{out|Output @ Tio.run}}
<pre style="white-space: pre-wrap;">--- 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

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

=={{header|C++}}==
<syntaxhighlight lang="cpp">#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' : ' ');
        }
    }
}</syntaxhighlight>

{{out}}
<pre>
  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
</pre>

=={{header|CLU}}==
<syntaxhighlight lang="clu">% 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 </syntaxhighlight>
{{out}}
<pre>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|Common Lisp}}==

<syntaxhighlight lang="lisp">
(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)) ))
</syntaxhighlight>

{{out}}
<pre>
;; Primes under 20:
(dotimes (i 20) (when (primep i) (print i)))

1 
2 
3 
5 
7 
11 
13 
17 
19 </pre>

=={{header|Cowgol}}==
<syntaxhighlight lang="cowgol">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();</syntaxhighlight>

{{out}}

<pre>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|D}}==
{{trans|Java}}
<syntaxhighlight lang="d">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;
}</syntaxhighlight>
{{out}}
<pre>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</pre>


=={{header|Dart}}==
{{trans|Swift}}
<syntaxhighlight lang="Dart">
BigInt factorial(BigInt n) {
  if (n == BigInt.zero) {
    return BigInt.one;
  }

  BigInt result = BigInt.one;
  for (BigInt i = n; i > BigInt.zero; i = i - BigInt.one) {
    result *= i;
  }

  return result;
}

bool isWilsonPrime(BigInt n) {
  if (n < BigInt.from(2)) {
    return false;
  }

  return (factorial(n - BigInt.one) + BigInt.one) % n == BigInt.zero;
}

void main() {
  var wilsonPrimes = [];
  for (var i = BigInt.one; i <= BigInt.from(100); i += BigInt.one) {
    if (isWilsonPrime(i)) {
      wilsonPrimes.add(i);
    }
  }

  print(wilsonPrimes);
}
</syntaxhighlight>
{{out}}
<pre>
[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|Draco}}==
<syntaxhighlight lang="draco">proc wilson(word n) bool:
    word f, i;
    if n<2 then
        false
    else
        f := n - 1;
        for i from n-2 downto 2 do
            f := (f*i) % n
        od;
        (f+1) % n = 0
    fi
corp

proc main() void:
    word i;
    for i from 1 upto 100 do
        if wilson(i) then
            write(i, ' ')
        fi
    od
corp</syntaxhighlight>
{{out}}
<pre>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|EasyLang}}==
{{trans|BASIC256}}
<syntaxhighlight>
func wilson_prime n .
   fct = 1
   for i = 2 to n - 1
      fct = fct * i mod n
   .
   return if fct = n - 1
.
for i = 2 to 100
   if wilson_prime i = 1
      write i & " "
   .
.
</syntaxhighlight>

{{out}}
<pre>
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|EDSAC order code}}==
{{trans|Pascal}}
A translation of the Pascal short-cut algorithm, for 17-bit) EDSAC signed integers. Finding primes in the range 65436..65536 took 80 EDSAC minutes, so there is not much point in implementing the unshortened algorithm or extending to 35-bit integers.
<syntaxhighlight lang="edsac">
[Primes by Wilson's Theoem, for Rosetta Code.]
[EDSAC program, Initial Orders 2.]
            T51K P64F  [address for G parameter: low-level subroutines]
            T47K P130F [M parameter: main routine + high-level subroutine]
   
[======== M parameter: Main routine + high-level subroutine ============]
            E25K TM GK 
[Editable range of integers to be tested for primality.]
[Integers are stored right-justified, so e.g. 1000 is P500F.]
      [0]   P500F   [lowest]
      [1]   P550F   [highest]
[Constants used with the M parameter]
      [2]   PD      [17-bit 1; also serves as letter P]
      [3]   K2048F  [set letters mode]
      [4]   #F      [set figures mode]
      [5]   RF      [letter R]
      [6]   IF      [letter I]
      [7]   MF      [letter M in letters mode, dot in figures mode]
      [8]   @F      [carriage return]
      [9]   &F      [line feed]
     [10]   !F      [space character]
     [11]   K4096F  [null character]

[Subroutine for testing whether 17-bit integer n is a prime,
  using Wilson's Theorem with short cut.]
[Input:  n in 6F.]
[Output: 0F holds 0 if n is prime, negative if n is not prime.]
     [12]   A3F T69@        [plant return link as usual         ]
            A6F S2F G68@    [acc := n - 2, exit if n < 2]
            A2@ T72@        [r := n - 1, clear acc]
            T7F             [extend n to 35 bits in 6D]
            A2@ U71@ U70@   [f := 1;  m := 1]
            A2F T73@        [m2inc := 3]
     [25]   A72@ S73@ G44@  [if r < m2inc jump to part 2]
            T72@            [dec( r, m2inc)]
            A70@ A2@ T70@   [inc( m)]
            H71@ V70@       [acc := f*m]
[Note that f and m are held as f/2^16 and m/2^16, so their product is (f*m)/2^32.
 We want to store the product as (f*m)/2^34, hence need to shift 2 right]
            R1F T4D         [shift product and pass to modulo subroutine]
     [36]   A36@ G38G       [call modulo subroutine]
            A4F T71@        [f := product modulo n]
            A73@ A2F T73@   [inc( m2inc, 2)]
            E25@            [always loop back]
[Part 2: Euclid's algorithm]
     [44]   TF              [clear acc]
            A6FT74@         [h := n]
     [47]   S71@ E63@       [if f = 0 then jump to test HCF]
            TF              [clear acc]
            A71@ T6F T7F    [f to 6F and extend to 35 bits in 6D]
            A74@ T4F T5F    [h to 4F and extend to 35 bits in 4D]
     [56]   A56@ G38G       [call subroutine, 4F := h modulo f]
            A71@ T74@       [h := f]
            A4F T71@        [f := (old h) modulo f]
            E47@            [always loop back]
[Here with acc = 0. Test for h = 1]
     [63]   A74@ S2@        [acc := h - 1]
            G68@            [return false if h = 0]
            TF SF           [acc := 1 - h]
     [68]   TF              [return result in 0F]
     [69]   ZF              [(planted) jump back to caller]
[Variables with names as in Pascal program]
     [70]   PF [m]
     [71]   PF [f]
     [72]   PF [r]
     [73]   PF [m2inc]
     [74]   PF [h]

[Subroutine for finding and printing primes between the passed-in limits]
[Input:  4F = minimum value, 5F = maximum value]
[Output: None. 4F and 5F are not preserved.]
     [75]   A3F T128@      [plant return link as usual]
           [Set letters mode, write 'PRIMES ', set figures mode]
            O3@ O2@ O5@ O6@ O7@ O124@ O104@ O10@ O4@
            A5F T130@       [store maximum value locally]
            A4F U129@       [store minimum value locally]
            TF              [pass minimum value to print subroutine]
            A11@ T1F        [pass null for leading zeros]
     [93]   A93@ GG         [call print subroutine]
            O7@ O7@         [print 2 dots for range]
            A130 @TF        [pass maximum value to print routine]
     [99]   A99@ GG         [call print subroutine]
            O8@ O9@         [print CRLF]
    [103]   A130@           [load n_max]
    [104]   S129@           [subtract n; also serves as letter S]
            G125@           [exit if n > n_max]
            TF              [clear acc]
            A129 @T6F       [pass current n to prime-testing subroutine]
    [109]   A109@ G12M      [call prime-testing subroutine]
            AF G120@        [load result, skip printing if n isn't prime]
            O10@            [print space]
            A129 @TF        [pass n to print subroutine]
            A11@ T1F        [pass null for leading zeros]
    [118]   A118@ GG        [call print subroutine]
    [120]   TF              [clear acc]
            A129@ A2@ T129@ [inc(n)]
    [124]   E103@           [always loop back; also serves as letter E]
    [125]   O8@ O9@         [print CRLF]
            TF              [clear acc before return (EDSAC convention)]
    [128]   ZF              [(planted) jump back to caller]
[Variables]
    [129]   PF [n]
    [130]   PF [n_max]

[Enter with acc = 0]
    [131]   A@ T4F          [pass lower limit to prime-finding subroutine]
            A1@ T5F         [pass upper limit to prime-finding subroutine]
    [135]   A135@ G75M      [call prime-finding subroutine]
            O11@            [print null to flush printer buffer]
            ZF              [stop]

[==================== G parameter: Low-level subroutines ====================]
            E25K TG 
[Subroutine to print non-negative 17-bit integer. Always prints 5 chars.]
[Caller specifies character for leading 0 (typically 0, space or null).]
[Parameters: 0F = integer to be printed (not preserved)]
            [1F = character for leading zero (preserved)]
[Workspace: 4F..7F, 38 locations]
      [0]   GKA3FT34@A1FT7FS35@T6FT4#FAFT4FH36@V4FRDA4#FR1024FH37@E23@O7FA2F
            T6FT5FV4#FYFL8FT4#FA5FL1024FUFA6FG16@OFTFT7FA6FG17@ZFP4FZ219DTF

[Subroutine to find X modulo M, where X and M are 35-bit integers.]
[Input:  X >= 0 in 4D, M > 0 in 6D.]
[Output: X modulo M in 4D, M preserved in 6D. Does not return the quotient.]
[Workspace: 0F.  27 locations.]
     [38]   GKA3FT26@A6DT8DA4DRDS8DG12@TFA8DLDE3@TF
            A4DS8DG17@T4DTFA6DS8DE26@TFA8DRDT8DE13@EF

[======== M parameter again ============]
            E25K TM GK
            E131Z   [define entry point]
            PF      [acc = 0 on entry]
[end]
</syntaxhighlight>
{{out}}
<pre>
PRIMES 1000..1100
 1009 1013 1019 1021 1031 1033 1039 1049 1051 1061 1063 1069 1087 1091 1093 1097
</pre>

=={{header|Erlang}}==
<syntaxhighlight lang="erlang">
#! /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)]]).
</syntaxhighlight>
{{Out}}
<pre>
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]
</pre>

=={{header|F_Sharp|F#}}==
<syntaxhighlight lang="fsharp">
// 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 ""</syntaxhighlight>
{{out}}
<pre>
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
</pre>

=={{header|Factor}}==
{{works with|Factor|0.99 2020-08-14}}
<syntaxhighlight lang="factor">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</syntaxhighlight>
{{out}}
<pre>
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
</pre>

=={{header|Fermat}}==
<syntaxhighlight lang="fermat">Func Wilson(n) = if ((n-1)!+1)|n = 0 then 1 else 0 fi.;</syntaxhighlight>

=={{header|Forth}}==
<syntaxhighlight lang="forth">
: fac-mod ( n m -- r )
    >r 1 swap
    begin dup 0> while
        dup rot * r@ 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 ;
</syntaxhighlight>
{{Out}}
<pre>
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
</pre>

=={{header|FreeBASIC}}==
<syntaxhighlight lang="freebasic">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</syntaxhighlight>
{{out}}
Primes below 100
<pre>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|FutureBasic}}==
<syntaxhighlight lang="futurebasic">
local fn WilsonPrime( n as long ) as BOOL
  long i, fct = 1
  BOOL result
  
  for i = 2 to n -1
    fct = (fct * i) mod n
  next i
  if fct == n - 1 then exit fn = YES else exit fn = NO
end fn = result

long i

print "Primes below 100:"

for i = 2 to 100
  if fn WilsonPrime(i) then print i
next

HandleEvents
</syntaxhighlight>
{{output}}
<pre>
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
</pre>

=={{header|Fōrmulæ}}==

{{FormulaeEntry|page=https://formulae.org/?script=examples/Primality_by_Wilson%27s_theorem}}

'''Solution'''

[[File:Fōrmulæ - Primality by Wilson's theorem 01.png]]

'''Test cases'''

[[File:Fōrmulæ - Primality by Wilson's theorem 02.png]]

[[File:Fōrmulæ - Primality by Wilson's theorem 03.png]]

=={{Header|GAP}}==
<syntaxhighlight lang="gap"># 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 );</syntaxhighlight>

{{out}}
<pre>
[ 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|Go}}==
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.
<syntaxhighlight lang="go">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()    
}</syntaxhighlight>

{{out}}
<pre>
  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 
</pre>

=={{header|Haskell}}==
<syntaxhighlight lang="haskell">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</syntaxhighlight>
{{out}}
<pre>
 --- ------- 
 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
</pre>

=={{header|J}}==
<syntaxhighlight lang="j">
   wilson=: 0 = (| !&.:<:)
   (#~ wilson) x: 2 + i. 30
2 3 5 7 11 13 17 19 23 29 31
</syntaxhighlight>

=={{header|Jakt}}==
<syntaxhighlight lang="jakt">
fn factorial_modulo<T>(anon n: T, modulus: T, accumulator: T = 1) throws -> T => match n {
    (..0) => { throw Error::from_string_literal("Negative factorial") }
    0 => accumulator
    else => factorial_modulo(n - 1, modulus, accumulator: (accumulator * n) % modulus)
}

fn is_prime(anon p: i64) throws -> bool => match p {
    (..1) => false
    else => factorial_modulo(p - 1, modulus: p) + 1 == p
}

fn main() {
    println("Primes under 100: ")
    for i in (-100)..100 {
        if is_prime(i) {
            print("{} ", i)
        }
    }
    println()
}
</syntaxhighlight>
{{out}}
<pre>
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 
</pre>

=={{header|Java}}==
Wilson's theorem is an ''extremely'' inefficient way of testing for primality.  As a result, optimizations such as caching factorials not performed.
<syntaxhighlight lang="java">
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;
    }

}
</syntaxhighlight>

{{out}}
<pre>
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 
</pre>

=={{header|jq}}==

'''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''''
<syntaxhighlight lang="jq">## 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))]</syntaxhighlight>
{{out}}
<syntaxhighlight lang="sh">
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]</syntaxhighlight>

=={{header|Julia}}==
<syntaxhighlight lang="julia">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])
</syntaxhighlight>{{out}}
<pre>
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]
</pre>

=={{header|Lua}}==
<syntaxhighlight lang="lua">-- 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</syntaxhighlight>
{{out}}
<pre>
 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|MAD}}==
<syntaxhighlight lang="mad">            NORMAL MODE IS INTEGER

            INTERNAL FUNCTION(N)
            ENTRY TO WILSON.
            WHENEVER N.L.2, FUNCTION RETURN 0B
            F = 1
            THROUGH FM, FOR I = N-1, -1, I.L.2
            F = F*I
FM          F = F-F/N*N
            FUNCTION RETURN N.E.F+1
            END OF FUNCTION

            PRINT COMMENT $ PRIMES UP TO 100$
            THROUGH TEST, FOR V=1, 1, V.G.100
TEST        WHENEVER WILSON.(V), PRINT FORMAT NUMF, V

            VECTOR VALUES NUMF = $I3*$
            END OF PROGRAM</syntaxhighlight>
{{out}}
<pre>PRIMES UP TO 100
  2
  3
  5
  7
 11
 13
 17
 19
 23
 29
 31
 37
 41
 43
 47
 53
 59
 61
 67
 71
 73
 79
 83
 89
 97</pre>

=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">ClearAll[WilsonPrimeQ]
WilsonPrimeQ[1] = False;
WilsonPrimeQ[p_Integer] := Divisible[(p - 1)! + 1, p]
Select[Range[100], WilsonPrimeQ]</syntaxhighlight>
{{out}}
Prime factors up to a 100:
<pre>{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|Miranda}}==
<syntaxhighlight lang="miranda">main :: [sys_message]
main = [Stdout (show (filter wilson [1..100]) ++ "\n")]

wilson :: num->bool
wilson n = False, if n<2
         = test (n-1) (n-2), otherwise
           where test f i = f+1 = n, if i<2
                          = test (f*i mod n) (i-1), otherwise</syntaxhighlight>
{{out}}
<pre>[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|Modula-2}}==
<syntaxhighlight lang="modula2">MODULE WilsonPrimes;
FROM InOut IMPORT WriteCard, WriteLn;

VAR i: CARDINAL;

PROCEDURE Wilson(n: CARDINAL): BOOLEAN;
VAR
    f, i: CARDINAL;
BEGIN
    IF n<2 THEN RETURN FALSE END;
    f := 1;
    FOR i := n-1 TO 2 BY -1 DO
        f := f*i MOD n
    END;
    RETURN f + 1 = n
END Wilson;

BEGIN
    FOR i := 1 TO 100 DO
        IF Wilson(i) THEN
            WriteCard(i, 3)
        END
    END;
    WriteLn
END WilsonPrimes.</syntaxhighlight>
{{out}}
<pre>  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|Nim}}==
<syntaxhighlight lang="nim">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(" ")</syntaxhighlight>

{{out}}
<pre>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</pre>

=={{header|PARI/GP}}==
<syntaxhighlight lang="parigp">Wilson(n) = prod(i=1,n-1,Mod(i,n))==-1
</syntaxhighlight>

=={{header|Pascal}}==
A console application in Free Pascal, created with the Lazarus IDE. Shows:

(1) Straightforward method that calculates (n - 1)! modulo n, where n is the number under test. Like most solutions on here, reduces the product modulo n at each step.

(2) Short cut, based on an observation in the AppleScript solution: if during the calculation of (n - 1)! a partial product is divisible by n, then n is not prime. In fact it suffices for a partial product and n to have a common factor greater than 1. Further, such a common factor must be present in s!, where s = floor(sqrt(n)). Having got s! modulo n we find its HCF with n by Euclid's algorithm; then n is prime if and only if the HCF is 1.
<syntaxhighlight lang="pascal">
program PrimesByWilson;
uses SysUtils;

(* Function to return whether 32-bit unsigned n is prime.
   Applies Wilson's theorem with full calculation of (n - 1)! modulo n. *)
function WilsonFullCalc( n : longword) : boolean;
var
  f, m : longword;
begin
  if n < 2 then begin
    result := false;  exit;
  end;
  f := 1;
  for m := 2 to n - 1 do begin
    f := (uint64(f) * uint64(m)) mod n; // typecast is needed
  end;
  result := (f = n - 1);
end;

(* Function to return whether 32-bit unsigned n is prime.
   Applies Wilson's theorem with a short cut. *)
function WilsonShortCut( n : longword) : boolean;
var
  f, g, h, m, m2inc, r : longword;
begin
  if n < 2 then begin
    result := false;  exit;
  end;
  (* Part 1: Factorial (modulo n) of floor(sqrt(n)) *)
  f := 1;
  m := 1;
  m2inc := 3; // (m + 1)^2 - m^2
  // Want to loop while m^2 <= n, but if n is close to 2^32 - 1 then least
  //   m^2 > n overflows 32 bits. Work round this by looking at r = n - m^2.
  r := n - 1;
  while r >= m2inc do begin
    inc(m);
    f := (uint64(f) * uint64(m)) mod n;
    dec( r, m2inc);
    inc( m2inc, 2);
  end;
 (* Part 2: Euclid's algorithm: at the end, h = HCF( f, n) *)
  h := n;
  while f <> 0 do begin
    g := h mod f;
    h := f;
    f := g;
  end;
  result := (h = 1);
end;

type TPrimalityTest = function( n : longword) : boolean;
procedure ShowPrimes( isPrime : TPrimalityTest;
                      minValue, maxValue : longword);
var
  n : longword;
begin
  WriteLn( 'Primes in ', minValue, '..', maxValue);
  for n := minValue to maxValue do
    if isPrime(n) then Write(' ', n);
  WriteLn;
end;

(* Main routine *)
begin
  WriteLn( 'By full calculation:');
  ShowPrimes( @WilsonFullCalc, 1, 100);
  ShowPrimes( @WilsonFullCalc, 1000, 1100);
  WriteLn; WriteLn( 'Using the short cut:');
  ShowPrimes( @WilsonShortCut, 1, 100);
  ShowPrimes( @WilsonShortCut, 1000, 1100);
  ShowPrimes( @WilsonShortCut, 4294967195, 4294967295 {= 2^32 - 1});
end.
</syntaxhighlight>
{{out}}
<pre>
By full calculation:
Primes in 1..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
Primes in 1000..1100
 1009 1013 1019 1021 1031 1033 1039 1049 1051 1061 1063 1069 1087 1091 1093 1097

Using the short cut:
Primes in 1..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
Primes in 1000..1100
 1009 1013 1019 1021 1031 1033 1039 1049 1051 1061 1063 1069 1087 1091 1093 1097
Primes in 4294967195..4294967295
 4294967197 4294967231 4294967279 4294967291
</pre>

=={{header|Perl}}==
{{libheader|ntheory}}
<syntaxhighlight lang="perl">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;</syntaxhighlight>
{{out}}
<pre>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</pre>

=={{header|Phix}}==
Uses the modulus method to avoid needing gmp, which was in fact about 7 times slower (when calculating the full factorials).
<!--<syntaxhighlight lang="phix">(phixonline)-->
 <span style="color: #008080;">function</span> <span style="color: #000000;">wilson</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
     <span style="color: #004080;">integer</span> <span style="color: #000000;">facmod</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span>
     <span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
         <span style="color: #000000;">facmod</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">remainder</span><span style="color: #0000FF;">(</span><span style="color: #000000;">facmod</span><span style="color: #0000FF;">*</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
     <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
     <span style="color: #008080;">return</span> <span style="color: #000000;">facmod</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">=</span><span style="color: #000000;">n</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
  
 <span style="color: #004080;">atom</span> <span style="color: #000000;">t0</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()</span>
 <span style="color: #004080;">sequence</span> <span style="color: #000000;">primes</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
 <span style="color: #004080;">integer</span> <span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">2</span> 
 <span style="color: #008080;">while</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">primes</span><span style="color: #0000FF;">)<</span><span style="color: #000000;">1015</span> <span style="color: #008080;">do</span>
     <span style="color: #008080;">if</span> <span style="color: #000000;">wilson</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
         <span style="color: #000000;">primes</span> <span style="color: #0000FF;">&=</span> <span style="color: #000000;">p</span>
     <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
     <span style="color: #000000;">p</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
 <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"The first 25 primes: %V\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">primes</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">25</span><span style="color: #0000FF;">]})</span>
 <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"         '' builtin: %V\n"</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">get_primes</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">25</span><span style="color: #0000FF;">)})</span>
 <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"primes[1000..1015]: %V\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">primes</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1000</span><span style="color: #0000FF;">..</span><span style="color: #000000;">1015</span><span style="color: #0000FF;">]})</span>
 <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"        '' builtin: %V\n"</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">get_primes</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">1015</span><span style="color: #0000FF;">)[</span><span style="color: #000000;">1000</span><span style="color: #0000FF;">..</span><span style="color: #000000;">1015</span><span style="color: #0000FF;">]})</span>
 <span style="color: #0000FF;">?</span><span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">)</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
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"
</pre>

=={{header|Plain English}}==
<syntaxhighlight lang="plainenglish">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.</syntaxhighlight>
{{out}}
<pre>
1 2 3 5 7 11
</pre>

=={{header|PL/I}}==
<syntaxhighlight lang="pli">/* 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 ;</syntaxhighlight>
{{out}}
<pre>
  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>

See also [[#Polyglot:PL/I and PL/M]]

=={{header|PL/M}}==
{{works with|8080 PL/M Compiler}} ... under CP/M (or an emulator)
<syntaxhighlight lang="pli">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</syntaxhighlight>
{{out}}
<pre>
 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>

See also [[#Polyglot:PL/I and PL/M]]

=={{header|Polyglot:PL/I and PL/M}}==
The following Primality by Wilson's theorem solution will run under both PL/M and PL/I.
{{works with|8080 PL/M Compiler}} ... under CP/M (or an emulator)
Should work with many PL/I implementations.
<br>
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.
<syntaxhighlight lang="pli">/* 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 ;</syntaxhighlight>
{{out}}
<pre>
 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|PROMAL}}==
<syntaxhighlight lang="promal">
;;; Find primes using Wilson's theorem:
;;;             p is prime if ( ( p - 1 )! + 1 ) mod p = 0

;;; returns TRUE(1) if p is prime by Wilson's theorem, FALSE(0) otherwise
;;;         computes the factorial mod p at each stage, so as to allow
;;;         for numbers whose factorial won't fit in 16 bits
PROGRAM wilson
INCLUDE library

FUNC BYTE isWilsonPrime
ARG WORD p
WORD i
WORD fModP
BYTE result
BEGIN
fModP = 1
IF p > 2
  FOR i = 2 TO p - 1
    fModP = ( fModP * i ) % p
IF fModP = p - 1
  result = 1
ELSE
  result = 0
RETURN result
END

WORD i
BEGIN
FOR i = 1 TO 100
  IF isWilsonPrime( i )
    OUTPUT " #W", i
END
</syntaxhighlight>
{{out}}
<pre>
 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|Python}}==
No attempt is made to optimise this as this method is a [https://en.wikipedia.org/wiki/Wilson%27s_theorem#Primality_tests very poor primality test].
<syntaxhighlight lang="python">from math import factorial

def is_wprime(n):
    return n == 2 or (
        n > 1
        and n % 2 != 0
        and (factorial(n - 1) + 1) % n == 0
    )

if __name__ == '__main__':
    c = int(input('Enter upper limit: '))
    print(f'Primes under {c}:')
    print([n for n in range(c) if is_wprime(n)])</syntaxhighlight>

{{out}}
<pre>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]</pre>

=={{header|Quackery}}==

<syntaxhighlight lang="quackery"> [ 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 ] ]</syntaxhighlight>

{{out}}

<pre>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 </pre>

=={{header|Raku}}==
(formerly Perl 6)
{{works with|Rakudo|2019.11}}
Not a particularly recommended way to test for primality, especially for larger numbers. It ''works'', but is slow and memory intensive.

<syntaxhighlight lang="raku" line>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];</syntaxhighlight>
{{out}}
<pre>   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</pre>

=={{header|REXX}}==
Some effort was made to optimize the factorial computation by using memoization and also minimize the size of the 
<br>decimal digit precision &nbsp; &nbsp; ('''NUMERIC DIGITS expression''').

Also, a "pretty print" was used to align the displaying of a list.
<syntaxhighlight lang="rexx">/*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).*/</syntaxhighlight>
Programming note:  &nbsp; This REXX program makes use of &nbsp; '''LINESIZE''' &nbsp; REXX program &nbsp; (or
BIF) &nbsp; which is used to determine the screen width
<br>(or linesize) &nbsp; of the terminal (console). &nbsp; Some
REXXes don't have this BIF.

The &nbsp; '''LINESIZE.REX''' &nbsp; REXX program is included
here &nbsp; ───► &nbsp; [[LINESIZE.REX]].


{{out|output|text=&nbsp; when using the default inputs:}}
<pre>
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.
</pre>

=={{header|Refal}}==
<syntaxhighlight lang="refal">$ENTRY Go {
    = <Prout <Filter Wilson <Iota 100>>>;
};

Wilson {
    s.N, <Compare s.N 2>: '-' = F;
    s.N = <Wilson s.N 1 <- s.N 1>>;
    s.N s.A 1, <- s.N 1>: { s.A = T; s.X = F; };
    s.N s.A s.C = <Wilson s.N <Mod <* s.A s.C> s.N> <- s.C 1>>;
};

Iota {
    s.N = <Iota 1 s.N>;
    s.N s.N = s.N;
    s.N s.M = s.N <Iota <+ 1 s.N> s.M>;
};

Filter {
    s.F = ;
    s.F t.I e.X, <Mu s.F t.I>: {
        T = t.I <Filter s.F e.X>;
        F = <Filter s.F e.X>;
    };
};</syntaxhighlight>
{{out}}
<pre>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|Ring}}==
<syntaxhighlight lang="ring">
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
</syntaxhighlight>
Output:
<pre>
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
</pre>

Alternative version computing the factorials modulo n so as to avoid overflow.
<syntaxhighlight lang="ring"># 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</syntaxhighlight>

{{out}}
<pre>
 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|RPL}}==
{{works with|HP|48S}}
 ≪ '''IF''' DUP 1 > '''THEN'''
       DUP → p j 
       ≪ 1 
          '''WHILE''' 'j' DECR 1 > '''REPEAT''' j * p MOD '''END'''
          1 + p MOD NOT
       ≫   
    '''ELSE''' DROP 0 '''END'''
 ≫ '<span style="color:blue">WILSON?<span>' STO

=={{header|Ruby}}==
<syntaxhighlight lang="ruby">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) }
</syntaxhighlight>
{{out}}
<pre>[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|Rust}}==
<syntaxhighlight lang="rust">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;
    }
}</syntaxhighlight>

{{out}}
<pre>
  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
</pre>

=={{header|Scala}}==
{{trans|Java}}
<syntaxhighlight lang="Scala">
import scala.math.BigInt

object PrimalityByWilsonsTheorem extends App {
  println("Primes less than 100 testing by Wilson's Theorem")
  (0 to 100).foreach(i => if (isPrime(i)) print(s"$i "))

  private def isPrime(p: Long): Boolean = {
    if (p <= 1) return false
    (fact(p - 1).+(BigInt(1))).mod(BigInt(p)) == BigInt(0)
  }

  private def fact(n: Long): BigInt = {
    (2 to n.toInt).foldLeft(BigInt(1))((fact, i) => fact * BigInt(i))
  }
}
</syntaxhighlight>
{{out}}
<pre>
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 
</pre>


=={{header|SETL}}==
<syntaxhighlight lang="setl">program wilsons_theorem;
    print({n : n in {1..100} | wilson n});

    op wilson(p);
        return p>1 and */{1..p-1} mod p = p-1;
    end op;
end program;</syntaxhighlight>
{{out}}
<pre>{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|Sidef}}==
<syntaxhighlight lang="ruby">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</syntaxhighlight>

=={{header|Swift}}==

Using a BigInt library.

<syntaxhighlight lang="swift">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))</syntaxhighlight>

{{out}}

<pre>[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|Tiny BASIC}}==
<syntaxhighlight lang="tinybasic">    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</syntaxhighlight>

=={{header|Wren}}==
{{libheader|Wren-gmp}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "./gmp" for Mpz
import "./fmt" for Fmt

var t = Mpz.new()

var wilson = Fn.new { |p|
    if (p < 2) return false
    return (t.factorial(p-1) + 1) % p == 0
}

var primes = [2]
var i = 3
while (primes.count < 1015) {
    if (wilson.call(i)) primes.add(i)
    i = i + 2
}

var candidates = [2, 3, 9, 15, 29, 37, 47, 57, 67, 77, 87, 97, 237, 409, 659]
System.print("  n | prime?\n------------")
for (cand in candidates) Fmt.print("$3d | $s", cand, wilson.call(cand))

System.print("\nThe first 120 prime numbers by Wilson's theorem are:")
Fmt.tprint("$3d", primes[0..119], 20)

System.print("\nThe 1,000th to 1,015th prime numbers are:")
System.print(primes[-16..-1].join(" "))</syntaxhighlight>

{{out}}
<pre>
  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 by 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 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
</pre>

=={{header|XPL0}}==
{{trans|ALGOL W}}
<syntaxhighlight lang "XPL0">    \ 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
    function IsWilsonPrime; integer N ;
    integer  FactorialModN, I;
    begin
        FactorialModN := 1;
        for I := 2 to N - 1 do FactorialModN := rem( FactorialModN * I / N );
        return FactorialModN = N - 1
    end \isWilsonPrime\ ;

    integer I;
    for I := 1 to 100 do if IsWilsonPrime( I ) then [IntOut(0, I); ChOut(0, ^ )]</syntaxhighlight>
{{out}}
<pre>
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|zkl}}==
{{libheader|GMP}} GNU Multiple Precision Arithmetic Library and primes
<syntaxhighlight lang="zkl">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 }) }</syntaxhighlight>
<syntaxhighlight lang="zkl">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();</syntaxhighlight>
{{out}}
<pre>
  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
</pre>