Wilson primes of order n

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

Task
Wilson primes of order n
You are encouraged to solve this task according to the task description, using any language you may know.
Definition
     (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.


Task

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.


Related task


ALGOL 68

Translation of: Visual Basic .NET

... but using a sieve for primeallity checking.
As with the various BASIC samples, all calculations are done MOD p2 so arbitrary precision integers are not needed.

BEGIN # find Wilson primes of order n, primes such that:                    #
      #                    ( ( n - 1 )! x ( p - n )! - (-1)^n ) mod p^2 = 0 #
    PR read "primes.incl.a68" PR                  # include prime utilities #
    []BOOL primes = PRIMESIEVE 11 000;         # sieve the primes to 11 500 #
    # returns TRUE if p is an nth order Wilson prime                        #
    PROC is wilson = ( INT n, p )BOOL:
         IF p < n THEN FALSE
         ELSE
            LONG INT p2    = p * p;
            LONG INT prod := 1;
            FOR i TO n - 1 DO # prod := ( n - 1 )! MOD p2                    #
                prod := ( prod * i ) MOD p2
            OD;
            FOR i TO p - n DO # prod := ( ( p - n )! * ( n - 1 )! ) MOD p2   #
                prod := ( prod * i ) MOD p2
            OD;
            0 = ( p2 + prod + IF ODD n THEN 1 ELSE -1 FI ) MOD p2
         FI # is wilson # ;
    # find the Wilson primes #
    print( ( " n: Wilson primes", newline ) );
    print( ( "-----------------", newline ) );
    FOR n TO 11 DO
        print( ( whole( n, -2 ), ":" ) ); 
        IF is wilson( n, 2 ) THEN print( ( " 2" ) ) FI;
        FOR p FROM 3 BY 2 TO UPB primes DO
            IF primes[ p ] THEN
                IF is wilson( n, p ) THEN print( ( " ", whole( p, 0 ) ) ) FI
            FI
        OD;
        print( ( newline ) )
    OD
END
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

BASIC

Applesoft BASIC

Translation of: Chipmunk Basic
100 home
110 print "n:     Wilson primes"
120 print "--------------------"
130 for n = 1 to 11
140  print n;chr$(9);
150  for p = 2 to 18
160   gosub 240
170   if pt = 0 then goto 200
180   gosub 340
190   if wnpt = 1 then print p,
200  next p
210  print
220 next n
230 end
240 rem tests if the number P is prime
250 rem result is stored in PT
260 pt = 1
270 if p = 2 then return
280 if p * 2 - int(p / 2) = 0 then pt = 0 : return
290 j = 3
300 if j*j > p then return
310 if p * j - int(p / j) = 0 then pt = 0 : return
320 j = j+2
330 goto 300
340 rem tests if the prime p is a Wilson prime of order n
350 rem make sure it actually is prime first
360 rem result is stored in wnpt
370 wnpt = 0
380 if p = 2 and n = 2 then wnpt = 1 : return
390 if n > p then wnpt = 0 : return
400 prod = 1 : p2 = p*p
410 for i = 1 to n-1
420  prod = (prod*i) : gosub 500
430 next i
440 for i = 1 to p-n
450  prod = (prod*i) : gosub 500
460 next i
470 prod = (p2+prod-(-1)^n) : gosub 500
480 if prod = 0 then wnpt = 1 : return
490 wnpt = 0 : return
500 rem prod mod p2 fails if prod > 32767 so brew our own modulus function
510 prod = prod-int(prod/p2)*p2
520 return

BASIC256

Translation of: FreeBASIC
function isPrime(v)
	if v <= 1 then return False
	for i = 2 To int(sqr(v))
		if v % i = 0 then return False
	next i
	return True
end function

function isWilson(n, p)
	if p < n then return false
	prod = 1
	p2 = p*p #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 false
end function

print " n:      Wilson primes"
print "----------------------"
for n = 1 to 11
	print n;" :  ";
	for p = 3 to 10499 step 2
		if isPrime(p) and isWilson(n, p) then print p; " ";
	next p
	print
next n
end

Chipmunk Basic

Works with: Chipmunk Basic version 3.6.4
Works with: GW-BASIC
Works with: MSX_BASIC
Works with: PC-BASIC version any
Works with: QBasic
Translation of: GW-BASIC
100 cls
110 print "n:     Wilson primes"
120 print "--------------------"
130 for n = 1 to 11
140  print n;chr$(9);
150  for p = 2 to 18
160   gosub 240
170   if pt = 0 then goto 200
180   gosub 340
190   if wnpt = 1 then print p,
200  next p
210  print
220 next n
230 end
240 rem tests if the number P is prime
250 rem result is stored in PT
260 pt = 1
270 if p = 2 then return
280 if p mod 2 = 0 then pt = 0 : return
290 j = 3
300 if j*j > p then return
310 if p mod j = 0 then pt = 0 : return
320 j = j+2
330 goto 300
340 rem tests if the prime p is a Wilson prime of order n
350 rem make sure it actually is prime first
360 rem result is stored in wnpt
370 wnpt = 0
380 if p = 2 and n = 2 then wnpt = 1 : return
390 if n > p then wnpt = 0 : return
400 prod = 1 : p2 = p*p
410 for i = 1 to n-1
420  prod = (prod*i) : gosub 500
430 next i
440 for i = 1 to p-n
450  prod = (prod*i) : gosub 500
460 next i
470 prod = (p2+prod-(-1)^n) : gosub 500
480 if prod = 0 then wnpt = 1 : return
490 wnpt = 0 : return
500 rem prod mod p2 fails if prod > 32767 so brew our own modulus function
510 prod = prod-int(prod/p2)*p2
520 return

QBasic

Works with: QBasic
Works with: QuickBasic
Translation of: FreeBASIC
FUNCTION isPrime (ValorEval)
    IF ValorEval < 2 THEN isPrime = False
    IF ValorEval MOD 2 = 0 THEN isPrime = 2
    IF ValorEval MOD 3 = 0 THEN isPrime = 3
    d = 5
    WHILE d * d <= ValorEval
        IF ValorEval MOD d = 0 THEN isPrime = False ELSE d = d + 2
    WEND
    isPrime = True
END FUNCTION

FUNCTION isWilson (n, p)
    IF p < n THEN isWilson = False
    prod = 1
    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 isWilson = True ELSE isWilson = False
END FUNCTION

PRINT " n:      Wilson primes"
PRINT "----------------------"
FOR n = 1 TO 11
    PRINT USING "##:      "; n;
    FOR p = 3 TO 10099 STEP 2
        If isPrime(p) AND isWilson(n, p) Then Print p; " ";
    NEXT p
    PRINT
NEXT n
END

MSX Basic

Both the GW-BASIC and Chipmunk Basic solutions work without change.

Visual Basic .NET

Translation of: QBasic

...but includes 2 and the 4th order Wilson Prime.

Option Strict On
Option Explicit On

Module WilsonPrimes

    Function isPrime(p As Integer) As Boolean
        If p < 2 Then Return False
        If p Mod 2 = 0 Then Return p = 2
        IF p Mod 3 = 0 Then Return p = 3
        Dim d As Integer = 5
        Do While d * d <= p
            If p Mod d = 0 Then
                Return False
            Else
                d = d + 2
            End If
        Loop
        Return True
    End Function

    Function isWilson(n As Integer, p As Integer) As Boolean
        If p < n Then Return False
        Dim prod As Long = 1
        Dim p2 As Long = p * p
        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 - If(n Mod 2 = 0, 1, -1)) Mod p2
        Return prod = 0
    End Function

    Sub Main()
        Console.Out.WriteLine(" n:      Wilson primes")
        Console.Out.WriteLine("----------------------")
        For n = 1 To 11
            Console.Out.Write(n.ToString.PadLeft(2) & ":      ")
            If isWilson(n, 2) Then Console.Out.Write("2 ")
            For p = 3 TO 10499 Step 2
                If isPrime(p) And isWilson(n, p) Then Console.Out.Write(p & " ")
            Next p
            Console.Out.WriteLine()
        Next n
    End Sub

End Module
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

Yabasic

Translation of: FreeBASIC
print "n:      Wilson primes"
print "---------------------"
for n = 1 to 11
    print n, ":",
    for p = 3 to 10099 step 2
        if isPrime(p) and isWilson(n, p) then print p, " ", : fi
    next p
    print
next n
end

sub isPrime(v)
    if v < 2 then return False : fi
    if mod(v, 2) = 0 then return v = 2 : fi
    if mod(v, 3) = 0 then return v = 3 : fi
    d = 5
    while d * d <= v
        if mod(v, d) = 0 then return False else d = d + 2 : fi
    end while
    return True
end sub

sub isWilson(n, p)
    if p < n then return False : fi
    prod = 1
	p2 = p**2
    for i = 1 to n-1
        prod = mod((prod*i), p2)
    next i
    for i = 1 to p-n
        prod = mod((prod*i), p2)
    next i
    prod = mod((p2 + prod - (-1)**n), p2)
    if prod = 0 then return True else return False : fi
end sub

C

Translation of: Wren
Library: GMP
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <gmp.h>

bool *sieve(int limit) {
    int i, p;
    limit++;
    // True denotes composite, false denotes prime.
    bool *c = calloc(limit, sizeof(bool)); // all false by default
    c[0] = true;
    c[1] = true;
    for (i = 4; i < limit; i += 2) c[i] = true;
    p = 3; // Start from 3.
    while (true) {
        int p2 = p * p;
        if (p2 >= limit) break;
        for (i = p2; i < limit; i += 2 * p) c[i] = true;
        while (true) {
            p += 2;
            if (!c[p]) break;
        }
    }
    return c;
}

int main() {
    const int limit = 11000;
    int i, j, n, pc = 0;
    unsigned long p;
    bool *c = sieve(limit);
    for (i = 0; i < limit; ++i) {
        if (!c[i]) ++pc;
    }
    unsigned long *primes = (unsigned long *)malloc(pc * sizeof(unsigned long));
    for (i = 0, j = 0; i < limit; ++i) {
        if (!c[i]) primes[j++] = i;
    }
    mpz_t *facts = (mpz_t *)malloc(limit *sizeof(mpz_t));
    for (i = 0; i < limit; ++i) mpz_init(facts[i]);
    mpz_set_ui(facts[0], 1);
    for (i = 1; i < limit; ++i) mpz_mul_ui(facts[i], facts[i-1], i);
    mpz_t f, sign;
    mpz_init(f);
    mpz_init_set_ui(sign, 1);
    printf(" n:  Wilson primes\n");
    printf("--------------------\n");
    for (n = 1; n < 12; ++n) {
        printf("%2d:  ", n);
        mpz_neg(sign, sign);
        for (i = 0; i < pc; ++i) {
            p = primes[i];
            if (p < n) continue;
            mpz_mul(f, facts[n-1], facts[p-n]);
            mpz_sub(f, f, sign);
            if (mpz_divisible_ui_p(f, p*p)) printf("%ld ", p);
        }
        printf("\n");
    }
    free(c);
    free(primes);
    for (i = 0; i < limit; ++i) mpz_clear(facts[i]);
    free(facts);
    return 0;
}
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 

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

EasyLang

Translation of: FreeBASIC
func isprim num .
   i = 2
   while i <= sqrt num
      if num mod i = 0
         return 0
      .
      i += 1
   .
   return 1
.
func is_wilson n p .
   if p < n
      return 0
   .
   prod = 1
   p2 = p * p
   for i = 1 to n - 1
      prod = prod * i mod p2
   .
   for i = 1 to p - n
      prod = prod * i mod p2
   .
   prod = (p2 + prod - pow -1 n) mod p2
   if prod = 0
      return 1
   .
   return 0
.
print "n:  Wilson primes"
print "-----------------"
for n = 1 to 11
   write n & "   "
   for p = 3 step 2 to 10099
      if isprim p = 1 and is_wilson n p = 1
         write p & " "
      .
   .
   print ""
.

F#

This task uses Extensible Prime Generator (F#)

// Wilson primes. Nigel Galloway: July 31st., 2021
let 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.functions
math.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--------------------" print
11 [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 false
end 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
    print
next 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 11
40 PRINT USING "##";N;
50 FOR P=2 TO 18
60 GOSUB 140
70 IF PT=0 THEN GOTO 100
80 GOSUB 230
90 IF WNPT=1 THEN PRINT P;
100 NEXT P
110 PRINT
120 NEXT N
130 END
140 REM tests if the number P is prime
150 REM result is stored in PT
160 PT = 1
170 IF P=2 THEN RETURN
175 IF P MOD 2 = 0 THEN PT=0:RETURN
180 J=3
190 IF J*J>P THEN RETURN
200 IF P MOD J = 0 THEN PT = 0: RETURN
210 J = J + 2
220 GOTO 190
230 REM tests if the prime P is a Wilson prime of order N
240 REM make sure it actually is prime first
250 REM RESULT is stored in WNPT
260 WNPT=0
270 IF P=2 AND N=2 THEN WNPT = 1: RETURN
280 IF N>P THEN WNPT=0: RETURN
290 PROD# = 1 : P2 = P*P
300 FOR I = 1 TO N-1
310 PROD# = (PROD#*I) : GOSUB 3000
320 NEXT I
330 FOR I = 1 TO P-N
340 PROD# = (PROD#*I) : GOSUB 3000
350 NEXT I
360 PROD# = (P2+PROD#-(-1)^N) : GOSUB 3000
370 IF PROD# = 0 THEN WNPT = 1: RETURN
380 WNPT = 0: RETURN
3000 REM    PROD# MOD P2 fails if PROD#>32767 so brew our own modulus function
3010 PROD# = PROD# - INT(PROD#/P2)*P2
3020 RETURN

J

   wilson=. 0 = (*:@] | _1&^@[ -~ -~ *&! <:@[)^:<:
   
   (>: i. 11x) ([ ;"0 wilson"0/ <@# ]) i.&.(p:inv) 11000
┌──┬───────────────┐
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            
├──┼───────────────┤
1011 1109        
├──┼───────────────┤
1117 2713        
└──┴───────────────┘

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 $p
def 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()
    end
end

wilsonprimes()

Output: Same as Wren example.

Mathematica /Wolfram Language

ClearAll[WilsonPrime]
WilsonPrime[n_Integer] := Module[{primes, out},
  primes = Prime[Range[PrimePi[11000]]];
  out = Reap@Do[
     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, strutils
import 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 = 1
echo " 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

PARI/GP

Translation of: Julia
default("parisizemax", "1024M");


\\ Define the function wilsonprimes with a default limit of 11000
wilsonprimes(limit) = {
    \\ Set the default limit if not specified
    my(limit = if(limit, limit, 11000));
    \\ Precompute factorial values up to the limit to save time
    my(facts = vector(limit, i, i!));
    \\ Sign variable for adjustment in the formula
    my(sgn = 1);
    print(" n:  Wilson primes\n--------------------");
    \\ Loop over the specified range (1 to 11 in the original code)
    for(n = 1, 11,
        print1(Str(" ", n, ":  "));
        sgn = -sgn; \\ Toggle the sign
        \\ Loop over all primes up to the limit
        forprime(p = 2, limit,
            \\ Check the Wilson prime condition modified for PARI/GP
            index=1;
            if(n<2,index=1,index=n-1);
            if(p > n && Mod(facts[index] * facts[p - n] - sgn, p^2) == 0,
                print1(Str(p, " "));
            )
        );
        print1("\n");
    );
}

\\ Execute the function with the default limit
wilsonprimes();
Output:
 n:  Wilson primes
--------------------
 1:  5 13 563 
 2:  3 11 107 4931 
 3:  7 
 4:  10429 
 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

Python

# wilson_prime.py by xing216
def sieve(n):
    multiples = []
    for i in range(2, n+1):
        if i not in multiples:
            yield i
            for j in range(i*i, n+1, i):
                multiples.append(j)
def intListToString(list):
    return " ".join([str(i) for i in list])
limit = 11000
primes = list(sieve(limit))
facs = [1]
for i in range(1,limit):
    facs.append(facs[-1]*i)
sign = 1
print(" n: Wilson primes")
print("—————————————————")

for n in range(1,12):
    sign = -sign
    wilson = []
    for p in primes:
        if p < n: continue
        f = facs[n-1] * facs[p-n] - sign
        if f % p**2 == 0: wilson.append(p)
    print(f"{n:2d}: {intListToString(wilson)}")
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

Racket

#lang racket

(require math/number-theory)

(define ((wilson-prime? n) p)
  (and (>= p n)
       (prime? p)
       (divides? (sqr p)
                 (- (* (factorial (- n 1))
                       (factorial (- p n)))
                    (expt -1 n)))))

(define primes<11000 (filter prime? (range 1 11000)))

(for ((n (in-range 1 (add1 11))))
  (printf "~a: ~a~%" n (filter (wilson-prime? n) primes<11000)))
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)

Raku

# Factorial
sub postfix:<!> (Int $n) { (constant f = 1, |[\×] 1..*)[$n] }

# Invisible times
sub infix:<⁢> is tighter(&infix:<**>) { $^a * $^b };

# Prime the iterator for thread safety
sink 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 a
commas: 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 j=@.#+2  by 2  for max(0, hip%2-@.#%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
───────┴─────────────────────────────────────────────────────────────

RPL

Works with: RPL version HP-49C
« → maxp
  « { } 
    1 11 FOR n 
      { } n
      IF DUP ISPRIME? NOT THEN NEXTPRIME END
      WHILE DUP maxp < REPEAT 
         n 1 - FACT OVER n - FACT * -1 n ^ -
         IF OVER SQ MOD NOT THEN SWAP OVER + SWAP END 
         NEXTPRIME
      END 
      DROP 1 →LIST +
    NEXT
» » 'TASK' STO
Output:
1: { { 5 13 } { 2 3 11 } { 7 } { } { 5 7 } { 11 } { 17 } { } { } { 11 } { 17 } }

Ruby

require "prime"
  
module Modulo
  refine Integer do
    def factorial_mod(m) = (1..self).inject(1){|prod, n| (prod *= n) % m }
  end
end

using Modulo
primes = Prime.each(11000).to_a

(1..11).each do |n| 
  res = primes.select do |pr|
    prpr = pr*pr
    ((n-1).factorial_mod(prpr) * (pr-n).factorial_mod(prpr) - (-1)**n) % (prpr) == 0
  end
  puts "#{n.to_s.rjust(2)}: #{res.inspect}"
end
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]

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 Int
import "./big" for BigInt
import "./fmt" for Fmt

var limit = 11000
var primes = Int.primeSieve(limit)
var facts = List.filled(limit, null)
facts[0] = BigInt.one
for (i in 1...11000) facts[i] = facts[i-1] * i
var sign = 1
System.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