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

From Rosetta Code
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

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.


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[edit]

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 ) )
OD
END
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++[edit]

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#[edit]

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[edit]

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[edit]

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[edit]

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[edit]

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

Java[edit]

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[edit]

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[edit]

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[edit]

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[edit]

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

Perl[edit]

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[edit]

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[edit]

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[edit]

# 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[edit]

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 [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[edit]

// [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[edit]

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[edit]

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