Wilson primes of order n: Difference between revisions
m (Repaired Wren entry.) |
m (→{{header|RPL}}: removed a debug instruction) |
||
(47 intermediate revisions by 22 users not shown) | |||
Line 1: | Line 1: | ||
{{ |
{{task|Prime Numbers}} |
||
;Definition |
;Definition |
||
Line 12: | Line 12: | ||
;Task |
;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, |
Calculate and show on this page the Wilson primes, if any, for orders '''n = 1 to 11 inclusive''' and for primes '''p < 18''' or, |
||
<br>if your language supports ''big integers'', for '''p < 11,000'''. |
|||
Line 18: | Line 19: | ||
:* [[Primality by Wilson's theorem]] |
:* [[Primality by Wilson's theorem]] |
||
=={{header|ALGOL 68}}== |
|||
{{Trans|Visual Basic .NET}} |
|||
... but using a sieve for primeallity checking. |
|||
<br>As with the various BASIC samples, all calculations are done MOD p2 so arbitrary precision integers are not needed. |
|||
<syntaxhighlight lang="algol68"> |
|||
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 |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|BASIC}}== |
|||
==={{header|Applesoft BASIC}}=== |
|||
{{trans|Chipmunk Basic}} |
|||
<syntaxhighlight lang="qbasic">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</syntaxhighlight> |
|||
==={{header|BASIC256}}=== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="basic256">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</syntaxhighlight> |
|||
==={{header|Chipmunk Basic}}=== |
|||
{{works with|Chipmunk Basic|3.6.4}} |
|||
{{works with|GW-BASIC}} |
|||
{{works with|MSX_BASIC}} |
|||
{{works with|PC-BASIC|any}} |
|||
{{works with|QBasic}} |
|||
{{trans|GW-BASIC}} |
|||
<syntaxhighlight lang="qbasic">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</syntaxhighlight> |
|||
==={{header|QBasic}}=== |
|||
{{works with|QBasic}} |
|||
{{works with|QuickBasic}} |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="qbasic">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</syntaxhighlight> |
|||
==={{header|MSX Basic}}=== |
|||
Both the [[#GW-BASIC|GW-BASIC]] and [[#Chipmunk_Basic|Chipmunk Basic]] solutions work without change. |
|||
==={{header|Visual Basic .NET}}=== |
|||
{{Trans|QBasic}} |
|||
...but includes 2 and the 4th order Wilson Prime. |
|||
<syntaxhighlight lang="vbnet"> |
|||
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 |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
==={{header|Yabasic}}=== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="yabasic">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</syntaxhighlight> |
|||
=={{header|C}}== |
|||
{{trans|Wren}} |
|||
{{libheader|GMP}} |
|||
<syntaxhighlight lang="c">#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; |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|C++}}== |
|||
{{libheader|GMP}} |
|||
<syntaxhighlight lang="cpp">#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'; |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|EasyLang}}== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight> |
|||
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 "" |
|||
. |
|||
</syntaxhighlight> |
|||
=={{header|F_Sharp|F#}}== |
|||
This task uses [http://www.rosettacode.org/wiki/Extensible_prime_generator#The_functions Extensible Prime Generator (F#)] |
|||
<syntaxhighlight lang="fsharp"> |
|||
// 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 "") |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|Factor}}== |
|||
{{works with|Factor|0.99 2021-06-02}} |
|||
<syntaxhighlight lang="factor">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</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|FreeBASIC}}== |
|||
This excludes the trivial case p=n=2. |
|||
<syntaxhighlight lang="freebasic">#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 |
|||
</syntaxhighlight> |
|||
{{out}}<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|Go}}== |
=={{header|Go}}== |
||
{{trans|Wren}} |
{{trans|Wren}} |
||
{{libheader|Go-rcu}} |
{{libheader|Go-rcu}} |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 61: | Line 711: | ||
fmt.Println() |
fmt.Println() |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 79: | Line 729: | ||
11: 17 2713 |
11: 17 2713 |
||
</pre> |
</pre> |
||
=={{header|GW-BASIC}}== |
|||
<syntaxhighlight lang="gwbasic">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</syntaxhighlight> |
|||
=={{header|J}}== |
|||
<syntaxhighlight lang="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 │ |
|||
├──┼───────────────┤ |
|||
│10│11 1109 │ |
|||
├──┼───────────────┤ |
|||
│11│17 2713 │ |
|||
└──┴───────────────┘</syntaxhighlight> |
|||
=={{header|Java}}== |
|||
<syntaxhighlight lang="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; |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|jq}}== |
|||
'''Works with [[jq]]''' (*) |
|||
'''Works with gojq, the Go implementation of jq''' |
|||
See e.g. [[Erd%C5%91s-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''' |
|||
<syntaxhighlight lang="jq">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));</syntaxhighlight> |
|||
'''Wilson primes''' |
|||
<syntaxhighlight lang="jq"># 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</syntaxhighlight> |
|||
{{out}} |
|||
gojq -ncr -f rc-wilson-primes.jq |
|||
<pre> |
|||
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] |
|||
</pre> |
|||
=={{header|Julia}}== |
|||
{{trans|Wren}} |
|||
<syntaxhighlight lang="julia">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() |
|||
</syntaxhighlight> |
|||
Output: Same as Wren example. |
|||
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
|||
<syntaxhighlight lang="mathematica">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} |
|||
]</syntaxhighlight> |
|||
{{out}} |
|||
<pre>{5,13,563} |
|||
{2,3,11,107,4931} |
|||
{7} |
|||
{10429} |
|||
{5,7,47} |
|||
{11} |
|||
{17} |
|||
{} |
|||
{541} |
|||
{11,1109} |
|||
{17,2713}</pre> |
|||
=={{header|Nim}}== |
=={{header|Nim}}== |
||
Line 84: | Line 982: | ||
{{libheader|bignum}} |
{{libheader|bignum}} |
||
As in Nim there is not (not yet?) a standard module to deal with big numbers, we use the third party module “bignum”. |
As in Nim there is not (not yet?) a standard module to deal with big numbers, we use the third party module “bignum”. |
||
< |
<syntaxhighlight lang="nim">import strformat, strutils |
||
import bignum |
import bignum |
||
Line 113: | Line 1,011: | ||
if f mod (p * p) == 0: |
if f mod (p * p) == 0: |
||
wilson.add p |
wilson.add p |
||
echo &"{n:2}: ", wilson.join(" ")</ |
echo &"{n:2}: ", wilson.join(" ")</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 129: | Line 1,027: | ||
10: 11 1109 |
10: 11 1109 |
||
11: 17 2713</pre> |
11: 17 2713</pre> |
||
=={{header|PARI/GP}}== |
|||
{{trans|Julia}} |
|||
<syntaxhighlight lang="PARI/GP"> |
|||
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(); |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|Perl}}== |
|||
{{libheader|ntheory}} |
|||
<syntaxhighlight lang="perl">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 |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> 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</pre> |
|||
=={{header|Phix}}== |
|||
{{trans|Wren}} |
|||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">limit</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">11000</span> |
|||
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span> |
|||
<span style="color: #004080;">mpz</span> <span style="color: #000000;">f</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</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: #7060A8;">get_primes_le</span><span style="color: #0000FF;">(</span><span style="color: #000000;">limit</span><span style="color: #0000FF;">),</span> |
|||
<span style="color: #000000;">facts</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_inits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">limit</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- (nb 0!==1!, same slot)</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;">limit</span> <span style="color: #008080;">do</span> <span style="color: #7060A8;">mpz_mul_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">facts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #000000;">facts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span><span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">sgn</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</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;">" n: Wilson primes\n"</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;">"--------------------\n"</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">11</span> <span style="color: #008080;">do</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;">"%2d: "</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">sgn</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">sgn</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</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: #008080;">do</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">primes</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">>=</span><span style="color: #000000;">n</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">facts</span><span style="color: #0000FF;">[</span><span style="color: #7060A8;">max</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)],</span><span style="color: #000000;">facts</span><span style="color: #0000FF;">[</span><span style="color: #7060A8;">max</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">-</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)])</span> |
|||
<span style="color: #7060A8;">mpz_sub_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">sgn</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpz_divisible_ui_p</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</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: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%d "</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</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;">"\n"</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<!--</syntaxhighlight>--> |
|||
Output: Same as Wren example. |
|||
=={{header|Prolog}}== |
|||
{{works with|SWI Prolog}} |
|||
<syntaxhighlight lang="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).</syntaxhighlight> |
|||
Module for finding prime numbers up to some limit: |
|||
<syntaxhighlight lang="prolog">:- 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).</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|Python}}== |
|||
<syntaxhighlight lang="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)}") |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|Racket}}== |
|||
<syntaxhighlight lang="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)))</syntaxhighlight> |
|||
{{out}} |
|||
<pre>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)</pre> |
|||
=={{header|Raku}}== |
=={{header|Raku}}== |
||
<syntaxhighlight lang="raku" line># Factorial |
|||
<lang perl6>sub postfix:<!> (Int $n) { (constant f = 1, |[\×] 1..*)[$n] } |
|||
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; |
my @primes = ^1.1e4 .grep: *.is-prime; |
||
Line 139: | Line 1,337: | ||
────────────────────'; |
────────────────────'; |
||
.say for (1..40).hyper(:1batch).map: -> \𝒏 { |
|||
-> \n { printf "%3d: %s\n", n, @primes.grep( ->\p { (p ≥ n) && ((n - 1)! × (p - n)! - (-1) ** n) %% p² } ).Str } for 1..11;</lang> |
|||
sprintf "%3d: %s", 𝒏, @primes.grep( -> \𝒑 { (𝒑 ≥ 𝒏) && ((𝒏 - 1)!(𝒑 - 𝒏)! - (-1) ** 𝒏) %% 𝒑² } ).Str |
|||
}</syntaxhighlight> |
|||
{{out}} |
{{out}} |
||
<pre> n: Wilson primes |
<pre> n: Wilson primes |
||
Line 153: | Line 1,353: | ||
9: 541 |
9: 541 |
||
10: 11 1109 |
10: 11 1109 |
||
11: 17 2713 |
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</pre> |
|||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
For more (extended) results, see this task's discussion page. |
|||
<lang rexx>/*REXX program finds and displays Wilson primes: a prime P such that P**2 divides:*/ |
|||
<syntaxhighlight lang="rexx">/*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.*/ |
/*────────────────── (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.*/ |
parse arg oLO oHI hip . /*obtain optional argument from the CL.*/ |
||
Line 172: | Line 1,402: | ||
say ' order │'center(title, w ) |
say ' order │'center(title, w ) |
||
say '───────┼'center("" , w, '─') |
say '───────┼'center("" , w, '─') |
||
do n=oLO to oHI |
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. */ |
|||
$= |
|||
lim= # |
else lim= # /* " " all " " orders ≥ 2.*/ |
||
$= /*$: a line (output) of Wilson primes.*/ |
|||
do j=1 for lim; p= @.j /*search through some generated primes.*/ |
|||
if |
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.*/ |
$= $ ' ' commas(p) /*add a commatized prime ──► $ list.*/ |
||
end /*p*/ |
end /*p*/ |
||
Line 195: | Line 1,425: | ||
!.=0; !.2=1; !.3=1; !.5=1; !.7=1; !.11=1 /* " " " " semaphores. */ |
!.=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.*/ |
sq.1=4; sq.2=9; sq.3= 25; sq.4= 49; #= 5; sq.#= @.#**2 /*squares of low primes.*/ |
||
do j=@.#+2 by 2 |
do j=@.#+2 by 2 for max(0, hip%2-@.#%2-1) /*find odd primes from here on. */ |
||
parse var j '' -1 _; if |
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? */ |
|||
if j// 7==0 then iterate /*" " " 7? */ |
|||
do k=5 while sq.k<=j /* [↓] divide by the known odd primes.*/ |
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. ___ */ |
if j // @.k == 0 then iterate j /*Is J ÷ X? Then not prime. ___ */ |
||
end /*k*/ /* [↑] only process numbers ≤ √ J */ |
end /*k*/ /* [↑] only process numbers ≤ √ J */ |
||
#= #+1; @.#= j; sq.#= j*j; !.j= 1 /*bump # of Ps; assign next P; P²; P# */ |
#= #+1; @.#= j; sq.#= j*j; !.j= 1 /*bump # of Ps; assign next P; P²; P# */ |
||
end /*j*/; return</ |
end /*j*/; return</syntaxhighlight> |
||
{{out|output|text= when using the default inputs:}} |
{{out|output|text= when using the default inputs:}} |
||
<pre> |
<pre> |
||
Line 220: | Line 1,449: | ||
11 │ 17 2,713 |
11 │ 17 2,713 |
||
───────┴───────────────────────────────────────────────────────────── |
───────┴───────────────────────────────────────────────────────────── |
||
</pre> |
|||
=={{header|RPL}}== |
|||
{{works with|RPL|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''' |
|||
» » '<span style="color:blue">TASK</span>' STO |
|||
{{out}} |
|||
<pre> |
|||
1: { { 5 13 } { 2 3 11 } { 7 } { } { 5 7 } { 11 } { 17 } { } { } { 11 } { 17 } } |
|||
</pre> |
|||
=={{header|Ruby}}== |
|||
<syntaxhighlight lang="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 |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> 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] |
|||
</pre> |
|||
=={{header|Rust}}== |
|||
<syntaxhighlight lang="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; |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|Sidef}}== |
|||
<syntaxhighlight lang="ruby">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) }) |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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] |
|||
</pre> |
</pre> |
||
Line 226: | Line 1,623: | ||
{{libheader|Wren-big}} |
{{libheader|Wren-big}} |
||
{{libheader|Wren-fmt}} |
{{libheader|Wren-fmt}} |
||
< |
<syntaxhighlight lang="wren">import "./math" for Int |
||
import "/big" for BigInt |
import "./big" for BigInt |
||
import "/fmt" for Fmt |
import "./fmt" for Fmt |
||
var limit = 11000 |
var limit = 11000 |
||
Line 247: | Line 1,644: | ||
} |
} |
||
System.print() |
System.print() |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
Latest revision as of 09:39, 10 March 2024
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
... 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
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
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
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
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
...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
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
#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++
#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
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
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
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 │
├──┼───────────────┤
│10│11 1109 │
├──┼───────────────┤
│11│17 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
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
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
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
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
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
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
« → 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
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
- Programming Tasks
- Prime Numbers
- ALGOL 68
- BASIC
- Applesoft BASIC
- BASIC256
- Chipmunk Basic
- QBasic
- MSX Basic
- Visual Basic .NET
- Yabasic
- C
- GMP
- C++
- EasyLang
- F Sharp
- Factor
- FreeBASIC
- Go
- Go-rcu
- GW-BASIC
- J
- Java
- Jq
- Julia
- Mathematica
- Wolfram Language
- Nim
- Bignum
- PARI/GP
- Perl
- Ntheory
- Phix
- Prolog
- Python
- Racket
- Raku
- REXX
- RPL
- Ruby
- Rust
- Sidef
- Wren
- Wren-math
- Wren-big
- Wren-fmt