My son gave me the following recurrence formula for $x_n$ ($n\ge2$):

$$(n+1)(n-2)x_{n+1}=n(n^2-n-1)x_n-(n-1)^3x_{n-1}\tag{1}$$ $$x_2=x_3=1$$

The task I got from him:

  • The sequence has an interesting property, find it out.
  • Make a conjecture and prove it.

Obviously I had to start with a few values and calculating them by hand turned out to be difficult. So I used Mathematica and defined the sequence as follows:

b[n_] := b[n] = n (n^2 - n - 1) a[n] - (n - 1)^3 a[n - 1];
a[n_] := a[n] = b[n - 1]/(n (n - 3));
a[2] = 1;
a[3] = 1;

And I got the following results:

$$ a_4=\frac{7}{4}, \ a_5=5, \ a_6=\frac{121}{6}, \ a_7=103, \ a_8=\frac{5041}{8}, \ a_9=\frac{40321}{9}, \\ a_{10}=\frac{362881}{10}, \ a_{11}=329891, \ a_{12}=\frac{39916801}{12}, \ a_{13}=36846277, \ a_{14}=\frac{6227020801}{14}\dots$$

Numbers don't make any sense but it's strange that the sequence produces integer values from time to time. It's not something that I expected from a pretty complex definition like (1).

So I decided to find the values of $n$ producing integer values of $a_n$. I did an experiment for $2\le n \le 100$:

table = Table[{i, a[i]}, {i, 2, 100}];
integers = Select[table, (IntegerQ[#[[2]]]) &];
itegerIndexes = Map[#[[1]] &, integers]

...and the output was:

{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 
59, 61, 67, 71, 73, 79, 83, 89, 97}

Conjecture (pretty amazing, at least to me):

$a_n$ is an integer if and only if $n$ is prime.

Interesting primality test, isn't it? The trick is to prove that it's correct. I have started with the substitution:

$$y_n=n x_n$$

...which simplifies (1) a bit:


...but I did not get much further (the next step, I guess, should be rearrangement).

Myunghyun Song
  • 21,403
  • 2
  • 22
  • 60
  • 15,427
  • 1
  • 18
  • 54
  • 2
    The fact that this sequence starts with the first prime index looks like a huge hint, seeing as mathematicians would start with n= 1 and software devs with n = 0 :-) – Carl Witthoft Mar 15 '19 at 15:51
  • 2
    If you're interested in returning the favor to your son, there are [Diophantine Equations](https://en.wikipedia.org/wiki/Formula_for_primes#Formula_based_on_a_system_of_Diophantine_equations) which enumerate all primes. The result is an inequality such that, if this polynomial function (whose inputs are 26 whole numbers labeled `a` through `z`) is greater than zero, then `k`, the 11th argument to the function, is prime, and *every* prime on the entire numberline is enumerated this way. – Cort Ammon Mar 15 '19 at 15:56
  • 2
    Looking at OEIS reveals origins of this problem, numerators are https://oeis.org/A005450, denominators are https://oeis.org/A005451, and in the references in the second one you find this is from Problem 10578 in The American Mathematical Monthly Vol. 104, No. 3 (Mar., 1997), page 270, see https://www.jstor.org/stable/2974795?read-now=1&seq=1#page_scan_tab_contents. – Sil Mar 16 '19 at 18:05

4 Answers4


The given difference equation can be solved in the following way. We have for $n\ge 3$, $$ (n-2)(y_{n+1}-y_n) = (n-1)^2(y_n-y_{n-1}),\\ \frac{y_{n+1}-y_n}{n-1}=(n-1)\frac{y_{n}-y_{n-1}}{n-2}. $$ If we let $\displaystyle z_n=\frac{y_{n}-y_{n-1}}{n-2}$, it follows that $$ z_{n+1}=(n-1)z_n=(n-1)(n-2)z_{n-1}=\cdots =(n-1)!z_3=(n-1)!\frac{3x_3-2x_2}{1}=(n-1)! $$ This gives $$ y_{n+1}-y_n=(n-1)(n-1)!=n!-(n-1)!, $$ hence $y_n = nx_n = (n-1)!+c$ for some $c$. Plugging $n=2$ yields $2=1!+c$, thus we have that $$\displaystyle x_n =\frac{(n-1)!+1}{n}.$$

Myunghyun Song
  • 21,403
  • 2
  • 22
  • 60

The $n^{\rm{th}}$ term of the sequence is $\dfrac{(n-1)! + 1}{n}$, which is an integer if and only if $n$ is prime (according to Wilson's Theorem).

  • 4,207
  • 1
  • 14
  • 25

Too long for a comment:

Numbers don't make any sense

Actually, they do ! :-) Just take a closer look at the sequence's composite-index denominators, and notice the following:

$$ a_{\color{blue}4}=\frac{7}{\color{blue}4},\quad a_{5}=5,\quad a_{\color{blue}6}=\frac{121}{\color{blue}6},\quad a_{7}=103,\quad a_{\color{blue}8}=\frac{5041}{\color{blue}8},\quad a_{\color{blue}9}=\frac{40321}{\color{blue}9},\quad \\~\\~\\ a_{\color{blue}{10}}=\frac{362881}{\color{blue}{10}},\quad a_{11}=329891,\quad a_{\color{blue}{12}}=\frac{39916801}{\color{blue}{12}},\quad a_{13}=36846277,\quad\dots $$

Let us now rewrite the sequence's prime-indexed elements in a similar manner, for a more uniform approach:

$$ a_{\color{blue}4}=\frac{7}{\color{blue}4},\quad a_{\color{blue}5}=\frac{25}{\color{blue}5},\quad a_{\color{blue}6}=\frac{121}{\color{blue}6},\quad a_{\color{blue}7}=\frac{721}{\color{blue}7},\quad a_{\color{blue}8}=\frac{5041}{\color{blue}8},\quad a_{\color{blue}9}=\frac{40321}{\color{blue}9},\quad \\~\\~\\ a_{\color{blue}{10}}=\frac{362881}{\color{blue}{10}}, a_{\color{blue}{11}}=\frac{3628801}{\color{blue}{11}},\quad a_{\color{blue}{12}}=\frac{39916801}{\color{blue}{12}},\quad a_{\color{blue}{13}}=\frac{479001601}{\color{blue}{13}},\quad\dots $$

At this point, we might be able to cheat, and use OEIS to identify the afferent sequence $\color{blue}{b_n=n\cdot a_n}$ by its first few elements, yielding three possible suspects: but let's say that our mathematical virtue and intellectual integrity will prevail over our base urges of rampant curiosity, and we might resist the temptation to do so. Could we then, without any outside aide, still manage to deduce an expression for $b_n$ ? Indeed, even a superficial glance will undoubtedly reveal the growth to resemble what one might otherwise expect to see in a geometric progression, rather than an arithmetic one. We then have:

$$ r_{\color{blue}4}=\frac{25}{7}\simeq\color{blue}4,\quad r_{\color{blue}5}=\frac{121}{25}\simeq\color{blue}5,\quad r_{\color{blue}6}=\frac{721}{121}\simeq\color{blue}6,\quad r_{\color{blue}7}=\frac{5041}{721}\simeq\color{blue}7,\quad\dots $$

In short, $\color{blue}{r_n=\dfrac{b_{n+1}}{b_n}\simeq n}.~$ A type of “modified geometric progression”, as it were, with a variable ratio equal to the index itself: Does this sound in any way familiar ? If no, there is still no need to worry: We'll get to the bottom of it soon enough, anyway. More precisely, we have $\color{blue}{b_{n+1}=n\cdot b_n-(n-1)}.~$ Now we are left with showing that, for prime values of the index p, $~\color{blue}{a_p=\dfrac{b_p}p\in\mathbb N},~$ starting from $\color{blue}{b_2=2}.~$ Could mathematical induction work in this case, or does the seemingly random distribution of primes throw a wrench into such an approach ? And, if so, could we then maybe slightly modify it, to fit the new conditions ? What would you say ? :-)

  • 46,778
  • 1
  • 77
  • 148
  • I'm confused. What does this add to the existing answers which already give the explicit formula $x_n=((n-1)!+1)/n$? – YiFan Mar 16 '19 at 03:22
  • 5
    @YiFan: Clarity, perhaps ? Maybe a bit of intuition ? – Lucian Mar 16 '19 at 03:23
  • 5
    +1. These types of answers are by far the most helpful IMO. They encourage discovery and intuitive thinking rather than just giving the answer. – Nico A Mar 17 '19 at 01:25

My process to find the solution probably involves too far a leap of faith, but I will post it anyway.

First, after writing out at least 8 terms, I noticed that, apart from integers, all terms are of the form $a_n=?/n$. So naturally I turned towards the sequence $b_n\equiv a_n\cdot n$, which goes like $3, 7, 25, 121, 721, 5041, 40321, \dots$.

Now my intuition senses the number sequence $504*, 4032*$, which resembles the factorial of 7&8, moreover, a lot of them seems to end in $1$. So I made a further conjecture that


A quick check of Wilson's theorem shows that it is compatible with the original conjecture. Thus, we seem to be on the right track. This can be easily proved by induction.

  • 3,044
  • 1
  • 7
  • 27