Prime test and triangular numbers​
​by Bill Gosper
I don’t think anyone I showed this figured out exactly why this works
In[]:=
oddPrimeQ[x_]:=x+1SeriesCoefficient[Sum[q^Binomial[n,2],{n,∞}]^4,{q,0,(x-1)/2}]
In[]:=
tim[xp_]:=(Print[#[[1]],",",Length[#[[2]]]];​​#[[2]])&@AbsoluteTiming[xp]
In[]:=
SetAttributes[tim,HoldAll];
In[]:=
Select[Range@239,oddPrimeQ]//tim
2.74415,51
Out[]=
{3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239}
It’s slow! But why does it work at all? This is so obscure that I’m spoiling it for you up front. The number of ways to partition n into four triangular numbers, counting 0s and permutations, equals the sum of the divisors of 2n+1. The generating function for the triangular numbers is
In[]:=
Sum[q^Binomial[n,2],{n,∞}]
Out[]=
EllipticTheta2,0,
q

2
1/8
q
Then
In[]:=
EllipticTheta[2,0,Sqrt[q]]^4/(16Sqrt[q])Sum[(1+2n)q^n/(1-q^(1+2n)),{n,0,∞}]Sum[DivisorSigma[1,2n+1]q^n,{n,0,∞}]
Out[]=
4
EllipticTheta2,0,
q

16
q

∞
∑
n=0
(1+2n)
n
q
1-
1+2n
q

∞
∑
n=0
n
q
DivisorSigma[1,1+2n]
In[]:=
Series[%/.∞11,{q,0,9}]
Out[]=
1+4q+6
2
q
+8
3
q
+13
4
q
+12
5
q
+14
6
q
+24
7
q
+18
8
q
+20
9
q
+
73/8
O[q]
1+4q+6
2
q
+8
3
q
+13
4
q
+12
5
q
+14
6
q
+24
7
q
+18
8
q
+20
9
q
+
10
O[q]
1+4q+6
2
q
+8
3
q
+13
4
q
+12
5
q
+14
6
q
+24
7
q
+18
8
q
+20
9
q
+
10
O[q]
Slash dot means substitute. It is amazing that Series still can’t figure out how many terms it needs of an infinite sum. So 19 is prime because the coefficient of q^9 is 20, the number of ways to express 9 as the sum of four triangular numbers:
In[]:=
Sum[(#q)^Binomial[n,2],{n,∞}]&/@(abcd)
Out[]=
EllipticTheta2,0,
aq
EllipticTheta2,0,
bq
EllipticTheta2,0,
cq
EllipticTheta2,0,
dq
16
1/8
(aq)
1/8
(bq)
1/8
(cq)
1/8
(dq)

In[]:=
Expand@SeriesCoefficient[%,{q,0,9}]
Out[]=
6
a
3
b
+
3
a
6
b
+
6
a
3
c
+
3
a
3
b
3
c
+
6
b
3
c
+
3
a
6
c
+
3
b
6
c
+
6
a
bcd+a
6
b
cd+ab
6
c
d+
6
a
3
d
+
3
a
3
b
3
d
+
6
b
3
d
+
3
a
3
c
3
d
+
3
b
3
c
3
d
+
6
c
3
d
+
3
a
6
d
+
3
b
6
d
+abc
6
d
+
3
c
6
d
E.g., a^6 d^3 says 9 = 6+0+0+3. 0, 1, 3, and 6 are triangular numbers.
In[]:=
Length@%
Out[]=
20