Search for the “best” arctan relation

           Jörg Arndt  <arndt@jjj.de>

What we are looking for

We search relations of the form

M_{1} arctan(1/X_{1})+M_{2} arctan(1/X_{2})+…+M_{n} arctan(1/X_{n}) = k π/4

where M_{j} ∈ ℤ, X_{j} ∈ ℕ, and k ∈ ℤ.

An example (a five-term relation):

π/4 = +88 arctan(1/192) + 39 arctan(1/239) + 100 arctan(1/515) −

−32 arctan(1/1068) − 56 arctan(1/173932)

Shorthand notation:

+88[192] +39[239] +100[515] -32[1068] -56[173932] == 1*Pi/4

Sort 6-term relations (all ==1*Pi/4) according to arguments:

+322[577]  +76[682] +139[1393] +156[12943] +132[32807] +44[1049433]
+122[319]  +61[378]  +115[557]   +29[1068]   +22[3458]   +44[27493]
+100[319] +127[378]   +71[557]   -15[1068]   +66[2943]  +44[478707]
+337[307] -193[463] +151[4193]  +305[4246] -122[39307]  -83[390112]

A relation is considered “better” than another if it precedes it. The “best” one is the first in our list of n-term relations.

Finding one relation

A 5-term relation:

+88[192] +39[239] +100[515] -32[1068] -56[173932] == 1 * Pi/4

For all (inverse) arguments X_{j}, factor X_{j}^{2} + 1:

192^2+1 == 36865 == 5 73 101
239^2+1 == 57122 == 2 13 13 13 13
515^2+1 == 265226 == 2 13 101 101
1068^2+1 == 1140625 == 5 5 5 5 5 5 73
173932^2+1 == 30252340625 == 5 5 5 5 5 13 73 101 101

All odd prime factors are {5, 13, 73, 101}.

Given a set of arguments, {X_{1}, …, X_{n}} so that all odd prime factors X_{j}^{2} + 1 are in a set {p_{1}, …, p_{n−1}}, how can we find the relation?

Finding one relation [cont.]

Set up a matrix whose entries are the signed exponents in the factorizations:

  X_j \ p_k 5  13 73 101
  173932  [-5, -1, 1, -2]
    1068  [ 6,  0, 1,  0]  (1068^2+1==5^6*73)
     515  [ 0,  1, 0, -2]
     239  [ 0, -4, 0,  0]  (239^2+1==2*13^4)
     192  [-1,  0, 1,  1]

Sign with argument X_{j} and prime p_{k} is minus if p_{k} < p_{j}/2 (mod p_{k}), else plus.

To find the multipliers M_{j}, compute the nullspace of the matrix:

? K=matkerint(M)   \\ pari/gp
[56] [32] [-100] [-39] [-88]

Thereby:

+88[192] +39[239] +100[515] -32[1068] -56[173932]   == k * Pi/4

A Floating point computation of the lhs. gives the value of k. Here it is +1.

Finding one n-term relation

Select a set 𝕊 of n − 1 odd primes 4k + 1.

Find numbers X that factor into (2 and) S by testing whether power products 2^{e_{0}} p_{1}^{e_{1}} p_{2}^{e_{2}} …, p_{n−1}^{e_{n−1}} are of the form X^{2} + 1.

Gauss used the set of the 8 smallest primes {5, 13, 17, 29, 37, 41, 53, 61} and found the 9-term relation

 +2805[5257] -398[9466] +1950[12943] +1850[34208] +2021[44179]
 +2097[85353] +1484[114669] +1389[330182] +808[485298]== 1 * Pi/4

The “best” 9-term relation known today is

 +3286[34208] +9852[39307] +5280[41688] +7794[44179] +7608[60443]
 +4357[275807] -1484[390112] -1882[619858] +776[976283]   == 1 * Pi/4

S = {5, 13, 17, 29, 41, 53, 97, 269}

Finding all n-term relations

Strategy: Find all X so that the odd prime factors of X^{2} +1 are in the set {5, 13, …, 761} (the first 64 odd primes of the form 4k+1). Let 𝕏 be the set of “candidate” values X. For n = 2, 3, 4, … find n-term relations as follows: For all (n − 1)-subsets S of the prime-set 𝕊: Select the subset of 𝕏 so that that all X^{2} + 1 factor into (2 and) S and compute the relation.

About “all”:

• Primes: 64 first primes for infinitely many.

• Candidates: X ≤ 10^{14} for X ≤ some-large-unknown-value.

• Subsets: \binom{64}{20} = 19, 619, 725, 782, 651, 120, i.e. too many.

A sieve for X so that X^{2} + 1 is smooth

For each prime p in S = {5, 13, …, 761} compute the two roots of the equation x^{2} + 1 ≡ 0 (mod p), (mod p^{2}), (mod p^{3}), … for p^{e} small enough.

Fill priority queue with events occurring at times x where x^{2} +1 ≡ 0 (mod p_{k}^{e}) (for all p_{k} and e). The corresponding event is to add log(p_{k}) at bucket x. Re-schedule event at x + p^{e}.

When all events up to time x are processed then test whether the sum in bucket x equals log(x^{2} + 1). If so, a candidate was found.

A sieve for X so that X^{2} + 1 is smooth [cont.]

Optimization: Let x_{0} be the candidate found most recently, and L_{0} = log(x_{0}^{2} + 1) then for the next candidate x_{1} we must have L_{1} > L_{0} (where L_{1} = log(x_{1}^{2} + 1)). That is, compute a logarithm only if the current sum is strictly greater than L_{0}. In practice a logarithm is computed exactly with each candidate found (one in 2 · 10^{9} tests).

We need about 240 CPU cycles per bucket. On average 1.5 additions happen with each bucket. Hmmmmmm … Overhead of priority queue may be a waste with frequent events (small p_{k}^{e}).

Use priority queues only for rare events … Use priority queues only for very rare events … Don’t use priority queues at all … Just over 11 CPU cycles per bucket. Happiness!

A sieve for X so that X^{2} + 1 is smooth [cont.]

Within 8 days all candidates X < 10^{14} were found (almost 44 thousand). The last few are

99501239756693
99627378461772
99759820688082
99849755159917
99950583525307
99955223464153

One can extend the list, e.g. using the relation

[x] = [x+d] + [x + (x^2+1)/d]  for all d dividing x^2+1

Almost 49 thousand values < 2^{64} are obtained

6541485364253438682
6775933831149960693
// 1e19
12637340271431925782
15674579038778108145
15689622322199725018
//276914859479857813947 > 2^64

Searching for n-term relations

Discard all X where X^{2}+1 has more than n−1 odd prime factors.

Discard small X.

Setup array of words whose bit-pattern correspond to odd primes occurring in factorizations of x^{2} + 1.

Search through all (n − 1) subsets of primes using only bit-and, bit-or and tests whether zero. Generate bit-combinations (used as masks) via

static inline ulong next_colex_comb(ulong x)
// Return smallest integer greater than x
// with the same number of bits set.
{
ulong r = x & -x; // lowest set bit
x += r;           // replace lowest block by a one left to it
if ( 0==x )  return 0;  // input was last combination
ulong l = x & -x; // first zero beyond lowest block
l -= r;           // lowest block  (cf. lowest_block())
while ( 0==(l&1) )  { l >>= 1; }  // move block to low end of word
return  x | (l>>1);  // need one bit less of low block
}

(cf. HAKMEM Memo 239)

Searching for n-term relations [cont.]

Whole array is accessed with each scan for with a new combination of primes. We want better memory locality (cache!).

Partition set of N = 64 primes into the 20 smallest and b = 64 − 20 = 44 “big” ones.

\binom{64}{q} = \binom{b}{0} \binom{N − b}{q} + \binom{b}{1} \binom{N − b}{q − 1} + \binom{b}{2} \binom{N − b}{q − 2} + …

= ∑_{j=0}^{q} \binom{b}{j} \binom{N − b}{q − j}

Now select arguments in two stages. Speedup by a factor of 25.

Small primes seem to be preferred in best relations:

     n prime set of best relation
     2 {13}
     3 {5, 13}
     4 {5, 13, 61}
     5 {5, 13, 73, 101}
     6 {5, 13, 61, 89, 197}
     7 {5, 13, 17, 29, 97, 433}
     8 {5, 13, 29, 37, 61, 97, 337}
     9 {5, 13, 17, 29, 41, 53, 97, 269}
    10 {5, 13, 17, 41, 53, 73, 97, 101, 157}
    11 {5, 13, 17, 29, 37, 53, 61, 89, 97, 101}

Omit large values of j in the sum with n > 12.

Best relations for 2 ≤ n ≤ 12

+4[5] -1[239]   == 1 * Pi/4
+12[18] +8[57] -5[239]   == 1 * Pi/4
+44[57] +7[239] -12[682] +24[12943]   == 1 * Pi/4
+88[192] +39[239] +100[515] -32[1068] -56[173932]   == 1 * Pi/4
+322[577] +76[682] +139[1393] +156[12943] +132[32807] +44[1049433]   == 1 * Pi/4
+1587[2852] +295[4193] +593[4246] +359[39307]
 +481[55603] +625[211050] -708[390112]   == 1 * Pi/4
+2192[5357] +2097[5507] -227[9466] +832[12943]
 +537[34522] -2287[39307] -171[106007] -708[1115618]   == 1 * Pi/4
+3286[34208] +9852[39307] +5280[41688] +7794[44179]
 +7608[60443] +4357[275807] -1484[390112] -1882[619858] +776[976283]   == 1 * Pi/4
+1106[54193] -30569[78629] -28687[88733] -13882[173932]
 +9127[390112] -9852[478707] -24840[1131527] +4357[3014557]
 +21852[5982670] +23407[201229582]   == -1 * Pi/4  <--=  ***NOTE***
+36462[390112] +135908[485298] +274509[683982] -39581[1984933]
 +178477[2478328] -114569[3449051] -146571[18975991] +61914[22709274]
 -69044[24208144] -89431[201229582] -43938[2189376182]   == 1 * Pi/4
+893758[1049433] +655711[1264557] +310971[1706203] +503625[1984933]
 -192064[2478328] -229138[3449051] -875929[18975991] -616556[21638297]
 -187143[22709274] -171857[24208144] -251786[201229582] -432616[2189376182]
 == 2 * Pi/4  <--=  ***NOTE***

Best relations, first arguments

n-terms       min-arg
 2                5  Machin (1706)   <--=  +4[5] -1[239]
 3               18  Gauss (YY?)     <--=  +12[18] +8[57] -5[239]
 4               57  Stormer (1896)  <--=  +44[57] +7[239] -12[682] +24[12943]
 5              192  JJ (1993), prev: Stormer (1896) 172
 6              577  JJ (1993)
 7            2,852  JJ (1993)
 8            5,357  JJ (2006), prev: JJ (1993) 4,246
 9           34,208  JJ (2006), prev: JJ (1993) 12,943, prev: Gauss (Y?) 5,257
10           54,193  JJ (2006), prev: JJ (1993) 51,387
11          390,112  JJ (1993)
12        1,049,433  JJ (2006), prev: JJ (1993) 683,982
13        3,449,051  JJ (2006), prev: JJ (1993) 1,984,933
14        6,826,318  JJ (2006)
15       20,942,043  HCL (1997), prev: MRW (1997) 18.975,991
16       53,141,564  JJ (2006)
17      201,229,582  JJ (2006)
18      299,252,491  JJ (2006)
19      778,401,733  JJ (2006)
20    2,674,664,693  JJ (2006)
21    5,513,160,193  JJ (2006)
22   17,249,711,432  JJ (2006), prev: 16,077,395,443 MRW (27-Jan-2003)
23   58,482,499,557  JJ (2006)
24  102,416,588,812  JJ (2006)
25  160,422,360,532  JJ (2006)
26  392,943,720,343  JJ (2006)
27  970,522,492,753  JJ (2006)
MRW := Michael Roby Wetherfield
HCL := Hwang Chien-lih
JJ := Joerg Arndt

Open questions Way. Too. Many.

More...

More details are given in the online draft of my book,

online at http://www.jjj.de/fxt/#fxtbook
All relations are given at http://www.jjj.de/arctan/

The following 11-term relation is still my favorite (found 1993):

\frac{π}{4} = 36462 arctan (\frac{1}{390112}) + 135908 arctan (\frac{1}{485298})

+274509 arctan (\frac{1}{683982}) − 39581 arctan (\frac{1}{1984933}) +178477 arctan (\frac{1}{2478328}) − 114569 arctan (\frac{1}{3449051}) −146571 arctan (\frac{1}{18975991}) + 61914 arctan (\frac{1}{22709274})

−69044 arctan (\frac{1}{24208144}) − 89431 arctan (\frac{1}{201229582})

−43938 arctan (\frac{1}{2189376182})

Primes: {5, 13, 17, 29, 37, 53, 61, 89, 97, 101}

            Jörg Arndt  <arndt@jjj.de>