# Compute the Wilson numbers

Given a positive integer n, compute the nth Wilson number W(n) where and e = 1 if n has a primitive root modulo n, otherwise e = -1. In other words, n has a primitive root if there does not exist an integer x where 1 < x < n-1 and x2 = 1 mod n.

• This is so create the shortest code for a function or program that computes the nth Wilson number for an input integer n > 0.
• You may use either 1-based or 0-based indexing. You may also choose to output the first n Wilson numbers.
• This is the OEIS sequence A157249.

## Test Cases

n  W(n)
1  2
2  1
3  1
4  1
5  5
6  1
7  103
8  13
9  249
10 19
11 329891
12 32
13 36846277
14 1379
15 59793
16 126689
17 1230752346353
18 4727
19 336967037143579
20 436486
21 2252263619
22 56815333
23 48869596859895986087
24 1549256
25 1654529071288638505

• Also, Oeis divides by n afterwards – H.PWiz Sep 4 '17 at 19:27
• @EriktheOutgolfer I added what is meant by having a primitive root. – miles Sep 4 '17 at 19:30
• Are we supposed to divide by n? – Leaky Nun Sep 4 '17 at 19:39
• As far as I'm aware, if k = 1 and e = -1, the result of the product would be 0. (sorry asking many questions but I need clarifications for my answer :p) – Erik the Outgolfer Sep 4 '17 at 19:41
• These numbers are called Wilson quotients. A Wilson number is an integer that divides its Wilson quotient evenly. For example, 13 is a Wilson number since 13 | 36846277. Also, W(n) usually excludes the denominator. – Dennis Sep 5 '17 at 1:00

# Jelly, 8 7 bytes

1 byte thanks to Dennis.

gRỊTP‘:


Try it online!

You don't really have to compute e since you need to divide anyway.

• gRỊT saves a byte. – Dennis Sep 5 '17 at 1:00
• Dennis getting down to the gRỊTty details of Jelly... – corsiKa Sep 5 '17 at 15:52

# Husk, 11 bytes

S÷ȯ→Π§foε⌋ḣ


Try it online!

### Explanation

          ḣ   Range from 1 to input
§foε⌋    Keep only those whose gcd with the input is 1
Π         Product
ȯ→          Plus 1
S÷            Integer division with input

• Please add explanation? I think you've got a nifty algo there... – Erik the Outgolfer Sep 4 '17 at 19:57

# Mathematica, 91 bytes

If[(k=#)==1,2,(Times@@Select[Range@k,CoprimeQ[k,#]&]+If[IntegerQ@PrimitiveRoot@#,1,-1])/#]&

• @BillSteihn please do not directly edit others' answers (relevant meta discussion). If you have a golfing suggestion, please leave a comment instead! – JungHwan Min Sep 4 '17 at 21:23
• @JungHwanMin Yes, I noticed that edit! thanks for helping new users with the rules – ZaMoC Sep 4 '17 at 21:28

# Pyth, 11 bytes

/h*Ff>2iTQS


Try it here!

# How?

• /h*Ff>2iTQS - Full program.

• S - Generate the inclusive range [1, input]

• f - Filter-keep those:

• iTQ - whose GCD with the input.

• >2 - Is less than two (can be replaced by either of the following: q1, !t)

• *F - Apply multiplication repeatedly. In other words, the product of the list.

• h - Increment the product by 1.

• / - Floor division with the input.

TL;DR: Get all the coprimes to the input in the range [1, input], get their product, increment it and divide it by the input.

# Python 2, 62 bytes

n=input();k=r=1
exec'r*=k/n*k%n!=1or k/n;k+=1;'*n*n
print-~r/n


Try it online!

# J, 33 bytes

3 :'<.%&y>:*/(#~1&=@(+.&y))1+i.y'


This one is more of a request to see an improvement than anything else. I tried a tacit solution first, but it was longer than this.

# explanation

This is fairly straightforward translation of Mr. Xcoder's solution into J.

Try it online!

# 05AB1E, 8 bytes

Lʒ¿}P>I÷


Try it online!

# R, 82 bytes

function(n)(prod((1:n)[g(n,1:n)<2])+1)%/%n
g=function(a,b)ifelse(o<-a%%b,g(b,o),b)


Uses integer division rather than figuring out e like many answers here, although I did work out that e=2*any((1:n)^2%%n==1%%n)-1 including the edge case of n=1 which I thought was pretty neat.

Try it online!

# Pari/GP, 36 bytes

n->(prod(i=1,n,i^(gcd(i,n)==1))+1)\n


Try it online!

## JavaScript (ES6), 7270 68 bytes

f=(n,p=1,i=n,a=n,b=i)=>i?f(n,b|a-1?p:p*i,i-=!b,b||n,b?a%b:i):-~p/n|0
<input type=number min=1 oninput=o.textContent=f(+this.value)><pre id=o>

Integer division strikes again. Edit: Saved 2 bytes thanks to @Shaggy. Saved a further 2 bytes by making it much more recursive, so it may fail for smaller values than it used to.

• 70 bytes (although I haven't had a chance to run a full set of tests on it yet): f=(n,i=n,p=1,g=(a,b)=>b?g(b,a%b):a)=>--i?f(n,i,g(n,i)-1?p:p*i):-~p/n|0 – Shaggy Sep 6 '17 at 14:28
• I went back to the recursive solution I was working on before I decided to try mapping an array instead and got it down to 70 bytes too. It's a bit of a mess but you might be able to salvage something from it to help get your solution down below 70: (n,x=n)=>(g=s=>--x?g(s*(h=(y,z)=>z?h(z,y%z):--y?1:x)(n,x)):++s)(1)/n|0 – Shaggy Sep 6 '17 at 15:45
• @Shaggy Well, I was inspired to have another look at it, but I'm not sure that's what you were expecting... – Neil Sep 7 '17 at 13:54

# Haskell, 42 bytes

f n=div(product[x|x<-[1..n],gcd x n<2]+1)n


Try it online!

Uses the integer division trick as all the other answers.
Uses 1-based indices.

Explanation

f n=                                       -- function
div                                  n -- integer division of next arg by n
(product                            -- multiply all entries in the following list
[x|                         -- return all x with ...
x<-[1..n],               -- ... 1 <= x <= n and ...
gcd x n<2]     -- ... gcd(x,n)==1
+1)  -- fix e=1


# Japt, 11 bytes

õ fjU ×Ä zU


Try it

## Explanation

Implicit input of integer U.

õ


Generate an array of integers from 1 to U.

fjU


Filter (f) co-primes of U.

×


Reduce by multiplication.

Ä


Add 1.

zU


Divide by U, floor the result and implicitly output.

• for n=25 it return 1654529071288638400 and it would be wrong because it would be 1654529071288638505 – user58988 Sep 6 '17 at 5:41
• @RosLuP: As confirmed by the challenge author, we don't need to handle numbers of more than 32-bits. – Shaggy Sep 6 '17 at 8:06

# Axiom, 121 bytes

f(n)==(e:=p:=1;for i in 1..n repeat(if gcd(i,n)=1 then p:=p*i;e=1 and i>1 and i<n-1 and(i*i)rem n=1=>(e:=-1));(p+e)quo n)


add some type, ungolf that and result

w(n:PI):PI==
e:INT:=p:=1
for i in 1..n repeat
if gcd(i,n)=1 then p:=p*i
e=1 and i>1 and i<n-1 and (i*i)rem n=1=>(e:=-1)
(p+e)quo n

(5) -> [[i,f(i)] for i in 1..25]
(5)
[[1,2], [2,1], [3,1], [4,1], [5,5], [6,1], [7,103], [8,13], [9,249],
[10,19], [11,329891], [12,32], [13,36846277], [14,1379], [15,59793],
[16,126689], [17,1230752346353], [18,4727], [19,336967037143579],
[20,436486], [21,2252263619], [22,56815333], [23,48869596859895986087],
[24,1549256], [25,1654529071288638505]]
Type: List List Integer

(8) -> f 101
(8)
9240219350885559671455370183788782226803561214295210046395342959922534652795_
041149400144948134308741213237417903685520618929228803649900990099009900990_
09901
Type: PositiveInteger


# JavaScript (ES6), 8381807876 68 bytes

My first pass at this was a few bytes longer than Neil's solution, which is why I originally ditched it in favour of the array reduction solution below. I've since golfed it down to tie with Neil.

n=>(g=s=>--x?g(s*(h=(y,z)=>z?h(z,y%z):--y?1:x)(n,x)):++s)(1,x=n)/n|0


### Try it

o.innerText=(f=
n=>(g=s=>--x?g(s*(h=(y,z)=>z?h(z,y%z):--y?1:x)(n,x)):++s)(1,x=n)/n|0
)(i.value=8);oninput=_=>o.innerText=f(+i.value)
<input id=i type=number><pre id=o>

# Non-recursive, 76 bytes

I wanted to give a non-recursive solution a try to see how it would turn out - not as bad as I expected.

n=>-~[...Array(x=n)].reduce(s=>s*(g=(y,z)=>z?g(z,y%z):y<2?x:1)(--x,n),1)/n|0


## Try it

o.innerText=(f=
n=>-~[...Array(x=n)].reduce(s=>s*(g=(y,z)=>z?g(z,y%z):y<2?x:1)(--x,n),1)/n|0
)(i.value=8);oninput=_=>o.innerText=f(+i.value)
<input id=i type=number><pre id=o>