Wilson primes of order n
A Wilson prime of order n is a prime number p such that p2 exactly divides:
- Definition
(n − 1)! × (p − n)! − (− 1)n
If n is 1, the latter formula reduces to the more familiar: (p - n)! + 1 where the only known examples for p are 5, 13, and 563.
- Task
Calculate and show on this page the Wilson primes, if any, for orders n = 1 to 11 inclusive and for primes p < 18 or, if your language supports big integers, for p < 11,000.
- Related task
Go
<lang 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()
}
}</lang>
- 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
Julia
<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
if (facts[n < 2 ? 1 : n - 1] * facts[p - n] - sgn) % p^2 == 0
print("$p ")
end
end
end
println()
end
end
wilsonprimes() </lang> Output: Same as Wren example.
Nim
As in Nim there is not (not yet?) a standard module to deal with big numbers, we use the third party module “bignum”. <lang Nim>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(" ")</lang>
- 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
<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
}</lang>
- 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
<lang perl6># 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
}</lang>
- 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
<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.*/ parse arg oLO oHI hip . /*obtain optional arguments 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 /*search for Wilson primes of order N.*/
nmf= !(n-1); pom = -1**n /*precalculate a factorial and +|- */
$=
lim= #; if n==1 then lim= 103 /*limit to known primes for 1st order. */
do j=1 for lim; p= @.j /*search through allowable primes. */
x= nmf * !(p-n) - pom /*calculate (n-1)! * (p-n)! - (-1)**n */
if x//sq.j \==0 then iterate /*is X ÷ prime squared? No, then skip.*/ /* ◄■■■■■■■ 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 to hip-1 /*find odd primes from here on. */
parse var j -1 _; if _==5 then iterate /*J divisible by 5? (right dig)*/
if j// 3==0 then iterate /*" " " 3? */
if j// 7==0 then iterate /*" " " 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</lang>
- 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 ───────┴─────────────────────────────────────────────────────────────
Wren
<lang ecmascript>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()
}</lang>
- 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