I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

# Wilson primes of order n

Wilson primes of order n
You are encouraged to solve this task according to the task description, using any language you may know.
Definition

A Wilson prime of order n is a prime number   p   such that   p2   exactly divides:

```     (n − 1)! × (p − n)! − (− 1)n
```

If   n   is   1,   the latter formula reduces to the more familiar:   (p - n)! + 1   where the only known examples for   p   are   5,   13,   and   563.

Calculate and show on this page the Wilson primes, if any, for orders n = 1 to 11 inclusive and for primes p < 18   or,
if your language supports big integers, for p < 11,000.

## ALGOL 68

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32
Translation of: Nim
which is
Translation of: Go
which is
Translation of: Wren

Algol 68G supports long integers, however the precision must be specified.
As the memory required for a limit of 11 000 would exceed he maximum supported by Algol 68G under WIndows, a limit of 5 500 is used here, which is sufficient to find all but the 4th order Wilson prime.

`BEGIN # find Wilson primes of order n, primes such that:                    #      #                    ( ( n - 1 )! x ( p - n )! - (-1)^n ) mod p^2 = 0 #    INT limit = 5 508; # max prime to consider #     # Build list of primes. #    []INT primes =        BEGIN            # sieve the primes to limit^2 which will hopefully be enough for primes #            [ 1 : limit * limit ]BOOL prime;            prime[ 1 ] := FALSE; prime[ 2 ] := TRUE;            FOR i FROM 3 BY 2 TO UPB prime DO prime[ i ] := TRUE  OD;            FOR i FROM 4 BY 2 TO UPB prime DO prime[ i ] := FALSE OD;            FOR i FROM 3 BY 2 TO ENTIER sqrt( UPB prime ) DO                IF prime[ i ] THEN FOR s FROM i * i BY i + i TO UPB prime DO prime[ s ] := FALSE OD FI            OD;            # count the primes up to the limit #            INT p count := 0; FOR i TO limit DO IF prime[ i ] THEN p count +:= 1 FI OD;            # construct a list of the primes #            [ 1 : p count ]INT primes;            INT p pos := 0;            FOR i WHILE p pos < UPB primes DO IF prime[ i ] THEN primes[ p pos +:= 1 ] := i FI OD;            primes        END;     # Build list of factorials. #    PR precision 20000 PR # set the number of digits for a LONG LONG INT #    [ 0 : primes[ UPB primes ] ]LONG LONG INT facts;    facts[ 0 ] := 1; FOR i TO UPB facts DO facts[ i ] := facts[ i - 1 ] * i OD;     # find the Wilson primes #    INT sign := 1;    print( ( " n: Wilson primes", newline ) );    print( ( "-----------------", newline ) );    FOR n TO 11 DO        print( ( whole( n, -2 ), ":" ) );         sign := - sign;        LONG LONG INT f n minus 1 = facts[ n - 1 ];        FOR p pos FROM LWB primes TO UPB primes DO            INT p = primes[ p pos ];            IF p >= n THEN                LONG LONG INT f = f n minus 1 * facts[ p - n ] - sign;                IF f MOD ( p * p ) = 0 THEN print( ( " ", whole( p, 0 ) ) ) FI            FI        OD;        print( ( newline ) )    ODEND`
Output:
``` n: Wilson primes
-----------------
1: 5 13 563
2: 2 3 11 107 4931
3: 7
4:
5: 5 7 47
6: 11
7: 17
8:
9: 541
10: 11 1109
11: 17 2713
```

## C++

Library: GMP
`#include <iomanip>#include <iostream>#include <vector>#include <gmpxx.h> std::vector<int> generate_primes(int limit) {    std::vector<bool> sieve(limit >> 1, true);    for (int p = 3, s = 9; s < limit; p += 2) {        if (sieve[p >> 1]) {            for (int q = s; q < limit; q += p << 1)                sieve[q >> 1] = false;        }        s += (p + 1) << 2;    }    std::vector<int> primes;    if (limit > 2)        primes.push_back(2);    for (int i = 1; i < sieve.size(); ++i) {        if (sieve[i])            primes.push_back((i << 1) + 1);    }    return primes;} int main() {    using big_int = mpz_class;    const int limit = 11000;    std::vector<big_int> f{1};    f.reserve(limit);    big_int factorial = 1;    for (int i = 1; i < limit; ++i) {        factorial *= i;        f.push_back(factorial);    }    std::vector<int> primes = generate_primes(limit);    std::cout << " n | Wilson primes\n--------------------\n";    for (int n = 1, s = -1; n <= 11; ++n, s = -s) {        std::cout << std::setw(2) << n << " |";        for (int p : primes) {            if (p >= n && (f[n - 1] * f[p - n] - s) % (p * p) == 0)                std::cout << ' ' << p;        }        std::cout << '\n';    }}`
Output:
``` n | Wilson primes
--------------------
1 | 5 13 563
2 | 2 3 11 107 4931
3 | 7
4 | 10429
5 | 5 7 47
6 | 11
7 | 17
8 |
9 | 541
10 | 11 1109
11 | 17 2713
```

## F#

This task uses Extensible Prime Generator (F#)

` // Wilson primes. Nigel Galloway: July 31st., 2021let rec fN g=function n when n<2I->g |n->fN(n*g)(n-1I) let fG (n:int)(p:int)=let g,p=bigint n,bigint p in (((fN 1I (g-1I))*(fN 1I (p-g))-(-1I)**n)%(p*p))=0I[1..11]|>List.iter(fun n->printf "%2d -> " n; let fG=fG n in pCache|>Seq.skipWhile((>)n)|>Seq.takeWhile((>)11000)|>Seq.filter fG|>Seq.iter(printf "%d "); printfn "") `
Output:
``` 1 -> 5 13 563
2 -> 2 3 11 107 4931
3 -> 7
4 -> 10429
5 -> 5 7 47
6 -> 11
7 -> 17
8 ->
9 -> 541
10 -> 11 1109
11 -> 17 2713
```

## Factor

Works with: Factor version 0.99 2021-06-02
`USING: formatting infix io kernel literals math math.functionsmath.primes math.ranges prettyprint sequences sequences.extras ; << CONSTANT: limit 11,000 >> CONSTANT: primes \$[ limit primes-upto ] CONSTANT: factorials    \$[ limit [1,b] 1 [ * ] accumulate* 1 prefix ] : factorial ( n -- n! ) factorials nth ; inline INFIX:: fn ( p n -- m )    factorial(n-1) * factorial(p-n) - -1**n ; : wilson? ( p n -- ? ) [ fn ] keepd sq divisor? ; inline : order ( n -- seq )    primes swap [ [ < ] curry drop-while ] keep    [ wilson? ] curry filter ; : order. ( n -- )    dup "%2d:  " printf order [ pprint bl ] each nl ; " n:  Wilson primes\n--------------------" print11 [1,b] [ order. ] each`
Output:
``` n:  Wilson primes
--------------------
1:  5 13 563
2:  2 3 11 107 4931
3:  7
4:  10429
5:  5 7 47
6:  11
7:  17
8:
9:  541
10:  11 1109
11:  17 2713
```

## FreeBASIC

This excludes the trivial case p=n=2.

`#include "isprime.bas" function is_wilson( n as uinteger, p as uinteger ) as boolean    'tests if p^2 divides (n-1)!(p-n)! - (-1)^n    'does NOT test the primality of p; do that first before you call this!    'using mods no big nums are required    if p<n then return false    dim as uinteger prod = 1, i, p2 = p^2    for i = 1 to n-1        prod = (prod*i) mod p2    next i    for i = 1 to p-n        prod = (prod*i) mod p2    next i    prod = (p2 + prod - (-1)^n) mod p2    if prod = 0 then return true else return falseend function print "n:      Wilson primes"print "--------------------"for n as uinteger = 1 to 11    print using "##      ";n;        for p as uinteger = 3 to 10099 step 2        if isprime(p) andalso is_wilson(n, p) then print p;" ";    next p    printnext n `
Output:
```
n:      Wilson primes
--------------------
1      5 13 563
2      3 11 107 4931
3      7
4
5      5 7 47
6      11
7      17
8
9      541
10      11 1109
11      17 2713

```

## Go

Translation of: Wren
Library: Go-rcu
`package main import (    "fmt"    "math/big"    "rcu") func main() {    const LIMIT = 11000    primes := rcu.Primes(LIMIT)    facts := make([]*big.Int, LIMIT)    facts[0] = big.NewInt(1)    for i := int64(1); i < LIMIT; i++ {        facts[i] = new(big.Int)        facts[i].Mul(facts[i-1], big.NewInt(i))    }    sign := int64(1)    f := new(big.Int)    zero := new(big.Int)    fmt.Println(" n:  Wilson primes")    fmt.Println("--------------------")    for n := 1; n < 12; n++ {        fmt.Printf("%2d:  ", n)        sign = -sign        for _, p := range primes {            if p < n {                continue            }            f.Mul(facts[n-1], facts[p-n])            f.Sub(f, big.NewInt(sign))            p2 := int64(p * p)            bp2 := big.NewInt(p2)            if f.Rem(f, bp2).Cmp(zero) == 0 {                fmt.Printf("%d ", p)            }        }        fmt.Println()    }}`
Output:
``` n:  Wilson primes
--------------------
1:  5 13 563
2:  2 3 11 107 4931
3:  7
4:  10429
5:  5 7 47
6:  11
7:  17
8:
9:  541
10:  11 1109
11:  17 2713
```

## GW-BASIC

`10 PRINT "n:     Wilson primes"20 PRINT "--------------------"30 FOR N = 1 TO 1140 PRINT USING "##";N;50 FOR P=2 TO 1860 GOSUB 14070 IF PT=0 THEN GOTO 10080 GOSUB 23090 IF WNPT=1 THEN PRINT P;100 NEXT P110 PRINT120 NEXT N130 END140 REM tests if the number P is prime150 REM result is stored in PT160 PT = 1170 IF P=2 THEN RETURN175 IF P MOD 2 = 0 THEN PT=0:RETURN180 J=3190 IF J*J>P THEN RETURN200 IF P MOD J = 0 THEN PT = 0: RETURN210 J = J + 2220 GOTO 190230 REM tests if the prime P is a Wilson prime of order N240 REM make sure it actually is prime first250 REM RESULT is stored in WNPT260 WNPT=0270 IF P=2 AND N=2 THEN WNPT = 1: RETURN280 IF N>P THEN WNPT=0: RETURN290 PROD# = 1 : P2 = P*P300 FOR I = 1 TO N-1310 PROD# = (PROD#*I) : GOSUB 3000320 NEXT I330 FOR I = 1 TO P-N340 PROD# = (PROD#*I) : GOSUB 3000350 NEXT I360 PROD# = (P2+PROD#-(-1)^N) : GOSUB 3000370 IF PROD# = 0 THEN WNPT = 1: RETURN380 WNPT = 0: RETURN3000 REM    PROD# MOD P2 fails if PROD#>32767 so brew our own modulus function3010 PROD# = PROD# - INT(PROD#/P2)*P23020 RETURN`

## Java

`import java.math.BigInteger;import java.util.*; public class WilsonPrimes {    public static void main(String[] args) {        final int limit = 11000;        BigInteger[] f = new BigInteger[limit];        f[0] = BigInteger.ONE;        BigInteger factorial = BigInteger.ONE;        for (int i = 1; i < limit; ++i) {            factorial = factorial.multiply(BigInteger.valueOf(i));            f[i] = factorial;        }        List<Integer> primes = generatePrimes(limit);        System.out.printf(" n | Wilson primes\n--------------------\n");        BigInteger s = BigInteger.valueOf(-1);        for (int n = 1; n <= 11; ++n) {            System.out.printf("%2d |", n);            for (int p : primes) {                if (p >= n && f[n - 1].multiply(f[p - n]).subtract(s)                        .mod(BigInteger.valueOf(p * p))                        .equals(BigInteger.ZERO))                    System.out.printf(" %d", p);            }            s = s.negate();            System.out.println();        }    }     private static List<Integer> generatePrimes(int limit) {        boolean[] sieve = new boolean[limit >> 1];        Arrays.fill(sieve, true);        for (int p = 3, s = 9; s < limit; p += 2) {            if (sieve[p >> 1]) {                for (int q = s; q < limit; q += p << 1)                    sieve[q >> 1] = false;            }            s += (p + 1) << 2;        }        List<Integer> primes = new ArrayList<>();        if (limit > 2)            primes.add(2);        for (int i = 1; i < sieve.length; ++i) {            if (sieve[i])                primes.add((i << 1) + 1);        }         return primes;    }}`
Output:
``` n | Wilson primes
--------------------
1 | 5 13 563
2 | 2 3 11 107 4931
3 | 7
4 | 10429
5 | 5 7 47
6 | 11
7 | 17
8 |
9 | 541
10 | 11 1109
11 | 17 2713
```

## jq

Works with jq (*)

Works with gojq, the Go implementation of jq

See e.g. Erdős-primes#jq for a suitable implementation of `is_prime` as used here.

(*) The C implementation of jq lacks the precision for handling the case p<11,000 so the output below is based on a run of gojq.

Preliminaries

`def emit_until(cond; stream): label \$out | stream | if cond then break \$out else . end; # For 0 <= \$n <= ., factorials[\$n] is \$n !def factorials:   reduce range(1; .+1) as \$n ([1];    .[\$n] = \$n * .[\$n-1]); def lpad(\$len): tostring | (\$len - length) as \$l | (" " * \$l)[:\$l] + .; def primes: 2, (range(3; infinite; 2) | select(is_prime));`

Wilson primes

`# Input: the limit of \$pdef wilson_primes:  def sgn: if . % 2 == 0 then 1 else -1 end;   . as \$limit  | factorials as \$facts  | " n:  Wilson primes\n--------------------",    (range(1;12) as \$n     | "\(\$n|lpad(2)) :  \(       [emit_until( . >= \$limit; primes)        | select(. as \$p            | \$p >= \$n and              ((\$facts[\$n - 1] * \$facts[\$p - \$n] - (\$n|sgn))                % (\$p*\$p) == 0 )) ])" ); 11000 | wilson_primes`
Output:

gojq -ncr -f rc-wilson-primes.jq

``` n:  Wilson primes
--------------------
1 :  [5,13,563]
2 :  [2,3,11,107,4931]
3 :  [7]
4 :  [10429]
5 :  [5,7,47]
6 :  [11]
7 :  [17]
8 :  []
9 :  [541]
10 :  [11,1109]
11 :  [17,2713]
```

## Julia

Translation of: Wren
`using Primes function wilsonprimes(limit = 11000)    sgn, facts = 1, accumulate(*, 1:limit, init = big"1")    println(" n:  Wilson primes\n--------------------")    for n in 1:11        print(lpad(n, 2), ":  ")        sgn = -sgn        for p in primes(limit)            if p > n && (facts[n < 2 ? 1 : n - 1] * facts[p - n] - sgn) % p^2 == 0                print("\$p ")            end        end        println()    endend wilsonprimes() `

Output: Same as Wren example.

## Mathematica/Wolfram Language

`ClearAll[WilsonPrime]WilsonPrime[n_Integer] := Module[{primes, out},  primes = Prime[Range[PrimePi[11000]]];  out = [email protected][     If[Divisible[((n - 1)!) ((p - n)!) - (-1)^n, p^2], Sow[p]]     ,     {p, primes}     ];  First[out[[2]], {}] ]Do[ Print[WilsonPrime[n]] , {n, 1, 11}]`
Output:
```{5,13,563}
{2,3,11,107,4931}
{7}
{10429}
{5,7,47}
{11}
{17}
{}
{541}
{11,1109}
{17,2713}```

## Nim

Translation of: Go
Library: bignum

As in Nim there is not (not yet?) a standard module to deal with big numbers, we use the third party module “bignum”.

`import strformat, strutilsimport bignum const Limit = 11_000 # Build list of primes using "nextPrime" function from "bignum".var primes: seq[int]var p = newInt(2)while p < Limit:  primes.add p.toInt  p = p.nextPrime() # Build list of factorials.var facts: array[Limit, Int]facts[0] = newInt(1)for i in 1..<Limit:  facts[i] = facts[i - 1] * i var sign = 1echo " n: Wilson primes"echo "—————————————————"for n in 1..11:  sign = -sign  var wilson: seq[int]  for p in primes:    if p < n: continue    let f = facts[n - 1] * facts[p - n] - sign    if f mod (p * p) == 0:      wilson.add p  echo &"{n:2}:  ", wilson.join(" ")`
Output:
```n: Wilson primes
—————————————————
1:  5 13 563
2:  2 3 11 107 4931
3:  7
4:  10429
5:  5 7 47
6:  11
7:  17
8:
9:  541
10:  11 1109
11:  17 2713```

## Perl

Library: ntheory
`use strict;use warnings;use ntheory <primes factorial>; my @primes = @{primes( 10500 )}; for my \$n (1..11) {    printf "%3d: %s\n", \$n, join ' ', grep { \$_ >= \$n && 0 == (factorial(\$n-1) * factorial(\$_-\$n) - (-1)**\$n) % \$_**2 } @primes}`
Output:
```  1: 5 13 563
2: 2 3 11 107 4931
3: 7
4: 10429
5: 5 7 47
6: 11
7: 17
8:
9: 541
10: 11 1109
11: 17 2713```

## Phix

Translation of: Wren
```with javascript_semantics
constant limit = 11000
include mpfr.e
mpz f = mpz_init()
sequence primes = get_primes_le(limit),
facts = mpz_inits(limit,1) -- (nb 0!==1!, same slot)
for i=2 to limit do mpz_mul_si(facts[i],facts[i-1],i) end for
integer sgn = 1
printf(1," n:  Wilson primes\n")
printf(1,"--------------------\n")
for n=1 to 11 do
printf(1,"%2d:  ", n)
sgn = -sgn
for i=1 to length(primes) do
integer p = primes[i]
if p>=n then
mpz_mul(f,facts[max(n-1,1)],facts[max(p-n,1)])
mpz_sub_si(f,f,sgn)
if mpz_divisible_ui_p(f,p*p) then
printf(1,"%d ", p)
end if
end if
end for
printf(1,"\n")
end for
```

Output: Same as Wren example.

## Prolog

Works with: SWI Prolog
`main:-    wilson_primes(11000). wilson_primes(Limit):-    writeln('  n | Wilson primes\n---------------------'),    make_factorials(Limit),    find_prime_numbers(Limit),    wilson_primes(1, 12, -1). wilson_primes(N, N, _):-!.wilson_primes(N, M, S):-    wilson_primes(N, S),    S1 is -S,    N1 is N + 1,    wilson_primes(N1, M, S1). wilson_primes(N, S):-    writef('%3r |', [N]),    N1 is N - 1,    factorial(N1, F1),    is_prime(P),    P >= N,    PN is P - N,    factorial(PN, F2),    0 is (F1 * F2 - S) mod (P * P),    writef(' %w', [P]),    fail.wilson_primes(_, _):-    nl. make_factorials(N):-    retractall(factorial(_, _)),    make_factorials(N, 0, 1). make_factorials(N, N, F):-    assert(factorial(N, F)),    !.make_factorials(N, M, F):-    assert(factorial(M, F)),    M1 is M + 1,    F1 is F * M1,    make_factorials(N, M1, F1).`

Module for finding prime numbers up to some limit:

`:- module(prime_numbers, [find_prime_numbers/1, is_prime/1]).:- dynamic is_prime/1. find_prime_numbers(N):-    retractall(is_prime(_)),    assertz(is_prime(2)),    init_sieve(N, 3),    sieve(N, 3). init_sieve(N, P):-    P > N,    !.init_sieve(N, P):-    assertz(is_prime(P)),    Q is P + 2,    init_sieve(N, Q). sieve(N, P):-    P * P > N,    !.sieve(N, P):-    is_prime(P),    !,    S is P * P,    cross_out(S, N, P),    Q is P + 2,    sieve(N, Q).sieve(N, P):-    Q is P + 2,    sieve(N, Q). cross_out(S, N, _):-    S > N,    !.cross_out(S, N, P):-    retract(is_prime(S)),    !,    Q is S + 2 * P,    cross_out(Q, N, P).cross_out(S, N, P):-    Q is S + 2 * P,    cross_out(Q, N, P).`
Output:
```  n | Wilson primes
---------------------
1 | 5 13 563
2 | 2 3 11 107 4931
3 | 7
4 | 10429
5 | 5 7 47
6 | 11
7 | 17
8 |
9 | 541
10 | 11 1109
11 | 17 2713
```

## Raku

`# Factorialsub postfix:<!> (Int \$n) { (constant f = 1, |[\×] 1..*)[\$n] } # Invisible timessub infix:<⁢> is tighter(&infix:<**>) { \$^a * \$^b }; # Prime the iterator for thread safetysink 11000!; my @primes = ^1.1e4 .grep: *.is-prime; say'  n: Wilson primes────────────────────'; .say for (1..40).hyper(:1batch).map: -> \𝒏 {     sprintf "%3d: %s", 𝒏, @primes.grep( -> \𝒑 { (𝒑 ≥ 𝒏) && ((𝒏 - 1)!⁢(𝒑 - 𝒏)! - (-1) ** 𝒏) %% 𝒑² } ).Str}`
Output:
```  n: Wilson primes
────────────────────
1: 5 13 563
2: 2 3 11 107 4931
3: 7
4: 10429
5: 5 7 47
6: 11
7: 17
8:
9: 541
10: 11 1109
11: 17 2713
12:
13: 13
14:
15: 349
16: 31
17: 61 251 479
18:
19: 71
20: 59 499
21:
22:
23:
24: 47 3163
25:
26:
27: 53
28: 347
29:
30: 137 1109 5179
31:
32: 71
33: 823 1181 2927
34: 149
35: 71
36:
37: 71 1889
38:
39: 491
40: 59 71 1171```

## REXX

For more (extended) results,   see this task's discussion page.

`/*REXX program finds and displays Wilson primes:  a prime   P   such that  P**2 divides:*//*────────────────── (n-1)! * (P-n)!  -  (-1)**n   where  n  is 1 ──◄ 11,   and  P < 18.*/parse arg oLO oHI hip .                          /*obtain optional argument from the CL.*/if  oLO=='' |  oLO==","  then  oLO=      1       /*Not specified?  Then use the default.*/if  oHI=='' |  oHI==","  then  oHI=     11       /* "      "         "   "   "     "    */if  hip=='' |  hip==","  then  hip=  11000       /* "      "         "   "   "     "    */call genP                                        /*build array of semaphores for primes.*/!!.= .                                           /*define the  default  for factorials. */bignum= !(hip)                                   /*calculate a ginormous factorial prod.*/parse value bignum 'E0' with ex 'E' ex .         /*obtain possible exponent of factorial*/numeric digits (max(9, ex+2) )                   /*calculate max # of dec. digits needed*/call facts hip                                   /*go & calculate a number of factorials*/title= ' Wilson primes  P  of order '  oLO " ──► " oHI',  where  P < '  commas(hip)w= length(title) + 1                             /*width of columns of possible numbers.*/say ' order │'center(title, w     )say '───────┼'center(""   , w, '─')      do n=oLO  to oHI;  nf= !(n-1)              /*precalculate a factorial product.    */      z= -1**n                                   /*     "       " plus or minus (+1│-1).*/      if n==1   then lim= 103                    /*limit to known primes for 1st order. */                else lim=   #                    /*  "    "  all     "    "  orders ≥ 2.*/      \$=                                         /*\$:  a line (output) of Wilson primes.*/         do j=1  for lim;    p= @.j              /*search through some generated primes.*/         if (nf*!(p-n)-z)//sq.j\==0 then iterate /*expression ~ q.j ?  No, then skip it.*/       /* ◄■■■■■■■ the filter.*/         \$= \$  ' '  commas(p)                    /*add a commatized prime  ──►  \$  list.*/         end   /*p*/       if \$==''  then \$= '        (none found within the range specified)'      say center(n, 7)'│'  substr(\$, 2)          /*display what Wilson primes we found. */      end   /*n*/say '───────┴'center(""   , w, '─')exit 0                                           /*stick a fork in it,  we're all done. *//*──────────────────────────────────────────────────────────────────────────────────────*/!:      arg x; if !!.x\==.  then return !!.x;  a=1;  do f=1  for x;   a=a*f; end; return acommas: parse arg ?;  do jc=length(?)-3  to 1  by -3; ?=insert(',', ?, jc); end;  return ?facts:  !!.= 1;   x= 1;  do  f=1  for hip;   x= x * f;   !!.f= x;           end;  return/*──────────────────────────────────────────────────────────────────────────────────────*/genP:        @.1=2; @.2=3; @.3=5; @.4=7;  @.5=11 /*define some low primes.              */      !.=0;  !.2=1; !.3=1; !.5=1; !.7=1;  !.11=1 /*   "     "   "    "     semaphores.  */      sq.1=4; sq.2=9; sq.3= 25; sq.4= 49; #= 5;  sq.#= @.#**2   /*squares of low primes.*/        do [email protected].#+2  by 2  for max(0, hip%[email protected].#%2-1)     /*find odd primes from here on. */        parse var  j   ''  -1  _;  if    _==5  then iterate    /*J ÷ 5?   (right digit).*/        if j//3==0  then iterate;  if j//7==0  then iterate    /*" " 3?    Is J ÷ by 7? */               do k=5  while sq.k<=j             /* [↓]  divide by the known odd primes.*/               if j // @.k == 0  then iterate j  /*Is  J ÷ X?  Then not prime.     ___  */               end   /*k*/                       /* [↑]  only process numbers  ≤  √ J   */        #= #+1;    @.#= j;    sq.#= j*j;  !.j= 1 /*bump # of Ps; assign next P;  P²; P# */        end          /*j*/;               return`
output   when using the default inputs:
``` order │ Wilson primes  P  of order  1  ──►  11,  where  P <  11,000
───────┼─────────────────────────────────────────────────────────────
1   │   5   13   563
2   │   2   3   11   107   4,931
3   │   7
4   │   10,429
5   │   5   7   47
6   │   11
7   │   17
8   │        (none found within the range specified)
9   │   541
10   │   11   1,109
11   │   17   2,713
───────┴─────────────────────────────────────────────────────────────
```

## Rust

`// [dependencies]// rug = "1.13.0" use rug::Integer; fn generate_primes(limit: usize) -> Vec<usize> {    let mut sieve = vec![true; limit >> 1];    let mut p = 3;    let mut sq = p * p;    while sq < limit {        if sieve[p >> 1] {            let mut q = sq;            while q < limit {                sieve[q >> 1] = false;                q += p << 1;            }        }        sq += (p + 1) << 2;        p += 2;    }    let mut primes = Vec::new();    if limit > 2 {        primes.push(2);    }    for i in 1..sieve.len() {        if sieve[i] {            primes.push((i << 1) + 1);        }    }    primes} fn factorials(limit: usize) -> Vec<Integer> {    let mut f = vec![Integer::from(1)];    let mut factorial = Integer::from(1);    f.reserve(limit);    for i in 1..limit {        factorial *= i as u64;        f.push(factorial.clone());    }    f} fn main() {    let limit = 11000;    let f = factorials(limit);    let primes = generate_primes(limit);    println!(" n | Wilson primes\n--------------------");    let mut s = -1;    for n in 1..=11 {        print!("{:2} |", n);        for p in &primes {            if *p >= n {                let mut num = Integer::from(&f[n - 1] * &f[*p - n]);                num -= s;                if num % ((p * p) as u64) == 0 {                    print!(" {}", p);                }            }        }        println!();        s = -s;    }}`
Output:
``` n | Wilson primes
--------------------
1 | 5 13 563
2 | 2 3 11 107 4931
3 | 7
4 | 10429
5 | 5 7 47
6 | 11
7 | 17
8 |
9 | 541
10 | 11 1109
11 | 17 2713
```

## Sidef

`func is_wilson_prime(p, n = 1) {    var m = p*p    (factorialmod(n-1, m) * factorialmod(p-n, m) - (-1)**n) % m == 0} var primes = 1.1e4.primes say "  n: Wilson primes\n────────────────────" for n in (1..11) {    printf("%3d: %s\n", n, primes.grep {|p| is_wilson_prime(p, n) })}`
Output:
```  n: Wilson primes
────────────────────
1: [5, 13, 563]
2: [2, 3, 11, 107, 4931]
3: [7]
4: [10429]
5: [5, 7, 47]
6: [11]
7: [17]
8: []
9: [541]
10: [11, 1109]
11: [17, 2713]
```

## Wren

Library: Wren-math
Library: Wren-big
Library: Wren-fmt
`import "/math" for Intimport "/big" for BigIntimport "/fmt" for Fmt var limit = 11000var primes = Int.primeSieve(limit)var facts = List.filled(limit, null)facts[0] = BigInt.onefor (i in 1...11000) facts[i] = facts[i-1] * ivar sign = 1System.print(" n:  Wilson primes")System.print("--------------------")for (n in 1..11) {    Fmt.write("\$2d:  ", n)    sign = -sign    for (p in primes) {        if (p < n) continue        var f = facts[n-1] * facts[p-n] - sign        if (f.isDivisibleBy(p*p)) Fmt.write("%(p) ", p)    }    System.print()}`
Output:
``` n:  Wilson primes
--------------------
1:  5 13 563
2:  2 3 11 107 4931
3:  7
4:  10429
5:  5 7 47
6:  11
7:  17
8:
9:  541
10:  11 1109
11:  17 2713
```