The PSLQ Algorithm:
Techniques for Efficient Computation
David H Bailey
Lawrence Berkeley National Laboratory
http://crd.lbl.gov/~dhbailey

1

LBNL’s high-precision software
QD: double-double (31 digits) and quad-double (62 digits).
ARPREC: arbitrary precision.
Low-level routines written in C++.
C++ and Fortran-90 translation modules permit use with existing C++ andFortran-90 programs -- only minor code changes are required.
 Includes many common functions: sqrt, cos, exp, gamma, etc.
 PSLQ, root finding, numerical integration.
Available at: http://www.experimentalmath.info
Authors: Xiaoye Li, Yozo Hida, Brandon Thompson and DHB.
Currently being maintained with the able assistance of Alex Kaiser of LBNL.

2

The PSLQ integer relation algorithm
Let (xn) be a given vector of real numbers. An integer relation algorithm
finds integers (an) such thata1 x1 + a2 x2 + · · · + an xn= 0(or within “epsilon” of zero, where epsilon = 10-p and p is the precision).
At the present time the “PSLQ” algorithm of mathematician-sculptor
Helaman Ferguson is the most widely used integer relation algorithm. It
was named one of ten “algorithms of the century” by Computing in Scienceand Engineering.

1. H. R. P. Ferguson, DHB and S. Arno, “Analysis of PSLQ, an integer relation finding algorithm,”
Mathematics of Computation, vol. 68, no. 225 (Jan 1999), pg. 351-369.
2. DHB and D. J. Broadhurst, “Parallel integer relation detection: Techniques and applications,” Mathematics
of Computation, vol. 70, no. 236 (Oct 2000), pg. 1719-1736.

3

PSLQ, continued
 PSLQ constructs a sequence of integer-valued matrices Bn that reduces
the vector y = x Bn, until either the relation is found (as one of the columnsof Bn), or else precision is exhausted.
 At the same time, PSLQ generates a steadily growing bound on the sizeof any possible relation.
 When a relation is found, the size of smallest entry of the vector y abruptly
drops to roughly “epsilon” (i.e. 10-p, where p is the number of digits ofprecision).
 The size of this drop can be viewed as a “confidence level” that therelation is real and not merely a numerical artifact -- a drop of 20+ ordersof magnitude almost always indicates a real relation.
 PSLQ (or any other integer relation scheme) requires very high precision
arithmetic (at least nd digits, where d is the size in digits of the largest ak),
both in the input data and in the operation of the algorithm.

4Decrease of log10(mink |yk|) as a functionof iteration number in a typical PSLQ run

5

Methodology for using PSLQ torecognize an unknown constant α	 Calculate α to high precision – typically 100 - 1000 digits. This is often themost computationally expensive part of the entire process.
 Based on experience with similar constants or relations, make a list of
possible terms on the right-hand side (RHS) of a linear formula for α, thencalculate each of the n RHS terms to the same precision as α.
 If you suspect α is algebraic of degree n (the root of a degree-n polynomialwith integer coefficients), compute the vector (1, α, α2, α3, …, αn).
 Apply PSLQ to the (n+1)-long vector, using the same numeric precision asα, but with a detection threshold a few orders of magnitude larger than“epsilon”– e.g., 10-480 instead of 10-500 for 500-digit arithmetic.
 When PSLQ runs, look for a detection following a drop in the size of the
reduced y vector by at least 20 orders of magnitude, to value near epsilon.
 If no credible relation is found, try expanding the list of RHS terms.
 Another possibility is to search for multiplicative relations (i.e., monomial
expressions), which can be done by taking logarithms of α and constants.
6

Bifurcation points in chaos theory:
The first “real” application of PSLQ
Let t = B3 = the smallest r suchthat the “logistic iteration”

xn+1= rxn (1 − xn )exhibits 8-way periodicity insteadof 4-way periodicity.
By means of a sequential
approximation scheme, one canobtain the numerical value of t toany desired precision:
3.544090359551922853615965986604804540583099845444573675457812
5303058429428588630122562585664248917999626…
Applying PSLQ to (1, t, t2, t3, …, t12), we obtained the result that t is a root of:

0 = 4913 + 2108t2 − 604t3 − 977t4 + 8t5 + 44t6 + 392t7−193t8 − 40t9 + 48t10 − 12t11 + t12
7

Ising integrals
We recently applied our methods to study three classes of integrals thatarise in the Ising theory of mathematical physics – Dn and two others:

Cn

:=

Dn

:=

En

=

� ∞

� ∞

du1
dun
·
·
·
�2
��
u1
un
n
0
0
(u
+
1/u
)
j
j=1 j
�2
�
�
ui −uj� ∞� ∞i<j ui +uj
du1
4
dun
···
···
�2
��
n! 0
u1
un
n
0
(u
+
1/u
)
j
j=1 j
2

� 1
� 1
� uk − uj dt2 dt3 · · · dtn
2···
uk + uj
0
0
4
n!

···

1

1≤j<k≤nwhere in the last line uk = t1 t2 … tk.

DHB, J. M. Borwein and R. E. Crandall, “Integrals of the Ising class,” Journal of Physics A: Mathematicaland General, vol. 39 (2006), pg. 12271-12302.
8

Ising integral evaluations
D2

=

1/3

D3

=

D4

=

8 + 4π 2 /3 − 27 L−3 (2)

E2

=

E3

=

E4

=

E5

?

=

4π 2 /9 − 1/6 − 7ζ(3)/2
6 − 8 log 2

10 − 2π 2 − 8 log 2 + 32 log2 2

22 − 82ζ(3) − 24 log 2 + 176 log2 2 − 256(log3 2)/3+16π 2 log 2 − 22π 2 /3

42 − 1984 Li4 (1/2) + 189π 4 /10 − 74ζ(3) − 1272ζ(3) log 2+40π 2 log2 2 − 62π 2 /3 + 40(π 2 log 2)/3 + 88 log4 2+464 log2 2 − 40 log 2

where Lin(x) is the polylog function. D2, D3 and D4 were originally provided
to us by mathematical physicist Craig Tracy, who hoped that our tools could
help identify D5.
9Recursions in Ising integrals
Consider the 2-parameter class of Ising integrals (which arises in QFT for odd k):

Cn,k

=

4
n!

� ∞
0

···

� ∞
0

du1
dun
···
�k+1
��
u1
un
n
(u
+
1/u
)
j
j=1 j
1

After computing 1000-digit numerical values for all n up to 36 and all k up to 75(performed on a highly parallel computer system), we discovered (using PSLQ)
linear relations in the rows of this array. For example, when n = 3:

0 = C3,0 − 84C3,2 + 216C3,4

0 = 2C3,1 − 69C3,3 + 135C3,5
0 = C3,2 − 24C3,4 + 40C3,6

0 = 32C3,3 − 630C3,5 + 945C3,7

0 = 125C3,4 − 2172C3,6 + 3024C3,8

Similar, but more complicated, recursions have been found for all n.
DHB, D. Borwein, J.M. Borwein and R.E. Crandall, “Hypergeometric forms for Ising-class integrals,”
Experimental Mathematics, vol. 16 (2007), pg. 257-276.
J. M. Borwein and B. Salvy, “A proof of a recursion for Bessel moments,” Experimental Mathematics, vol. 17(2008), pg. 223-230.
10

Recent result (18 Jan 2009)∆3 (−1) =
=

√
� ∞
2
2
(−1 + e−u + π u erf(u))3
√
du
u6
π 0
√
√
√
√ �
1 �
6 + 6 2 − 12 3 − 10π + 30 log(1 + 2) + 30 log(2 + 3)
15

As in many of the previous results, this was found by first computing the
integral to high precision (250 to 1000 digits), conjecturing possible termson the right-hand side, then applying PSLQ to look for a relation. We now
have proven this result.
Dozens of similar results have since been found (see next few viewgraphs),raising hope that all box integrals eventually will be evaluated in closed
form.
DHB, J. M. Borwein and R. E. Crandall, “Advances in the theory of box integrals,” manuscript, Mar 2009,
available at http://crd.lbl.gov/~dhbailey/dhbpapers/BoxII.pdf.
11

Recent evaluations of box integrals
n
any
1
2
2
2
2
2
2
3
3
3
3
3
3

s
even s ≥ 0
s �= −1
-4
-3
-1
1
3
s �= −2
-5
-4
-2
-1
1
3

Bn (s)
rational, e.g., : B2 (2) = 2/3
1
s+1
1
−4√
− π8

− 2√
2 log(1 + 2)√
√
1
1
2
+
3√
3 log(1 + √2)
73
2+
log(1 + 2)�
5
20
�
2
1
s 3
F
,
−
2
1
2+s
2
2 ; 2 ; −1
√
1
1
−√
3
−
6
12 π
− 32 2 arctan √12√√3−3G + 2 π log(1 + 2)�+ 3 √Ti2�(3 − 2 2)
1
3
−
π
+
3√ �
4
2 log 2 +
√
�
1
1
1
3
−
π
+
4√
24
2 log �2 + √3 �
1
7
2
3
−
π
−
3
5
60
20 log 2 +Here F is hypergeometric function; G is Catalan; Ti is Lewin’s inverse-tan function.
12

Elliptic Integrals
Recent research in integrals of elliptic functions have revealed hundreds ofheretofore unknown identities, for instance:

2

� 1
0

K 2 (x) dx − 4

� 1
0

−2

−2
K(x)E(x) dx + 3
� 1

� 1

� 1
0

xK(x) dx + 30

E 2 (x) dx −K 3 (x)K � (x)E � (x) dx +

0

� 1

� 1

� 1

xE(x) dx = 00K � (x)E � (x) dx = 0

0

E(x)K �3 (x)E � (x) dx = 00These studies involved computing thousands of individual definite
integrals, each to at least 1600-digit precision, then searching for
relations among them using PSLQ.

13

Statement of the PSLQ algorithm
Initialization:
1. For j := 1 to n: for i := 1 to n: if i = j then set Aij := 1 and Bij := 1
else set Aij := 0 and Bij := 0; endfor; endfor.
��n2
2. For k := 1 to n: set sk :=j=k xj ; endfor. Set t = 1/s1 . For k := 1 ton: set yk := txk ; sk := tsk ; endfor.
3. Initial H: For j := 1 to n − 1: for i := 1 to j − 1: set Hij := 0; endfor;
set Hjj := sj+1 /sj ; for i := j + 1 to n: set Hij := −yi yj /(sj sj+1 ); endfor;
endfor.
4. Reduce H: For i := 2 to n: for j := i − 1 to 1 step −1: set t :=nint(Hij /Hjj ); and yj := yj + tyi ; for k := 1 to j: set Hik := Hik − tHjk ;endfor; for k := 1 to n: set Aik := Aik − tAjk and Bkj := Bkj + tBki ;endfor; endfor; endfor.

14

Statement of the PSLQ algorithm,
continued
Iteration:
1. Select m such that γ i |Hii | is maximal when i = m.
2. Exchange the entries of y indexed m and m + 1, the corresponding rowsof A and H, and the corresponding columns of B.
�22
3. Remove corner on H diagonal: If m ≤ n−2 then set t0 := Hmm+ Hm,m+1,
t1 := Hmm /t0 and t2 := Hm,m+1 /t0 ; for i := m to n: set t3 := Him ,t4 := Hi,m+1 , Him := t1 t3 + t2 t4 and Hi,m+1 := −t2 t3 + t1 t4 ; endfor;
endif.

4. Reduce H: For i := m + 1 to n: for j := min(i − 1, m + 1) to 1 step
−1: set t := nint(Hij /Hjj ) and yj := yj + tyi ; for k := 1 to j: set
Hik := Hik − tHjk ; endfor; for k := 1 to n: set Aik := Aik − tAjk andBkj := Bkj + tBki ; endfor; endfor; endfor.
5. Norm bound: Compute M := 1/ maxj |Hjj |. Then there can exist no
relation vector whose Euclidean norm is less than M .
6. Termination test: If the largest entry of A exceeds the level of numericprecision used, then precision is exhausted. If the smallest entry of the yvector is less than the detection threshold (see below), a relation has been
detected and is given in the corresponding column of B.
15

Multi-level implementations of PSLQ In spite of the effectiveness of PSLQ, computation time grows cubicallywith the vector size n and also increases sharply with the level ofnumerical precision employed.
 Huge savings in run time can be achieved by employing a “two-level”
implementation of PSLQ: Perform most iterations with simple doubleprecision arithmetic, periodically updating the multiprecision arrays.
 Because double-precision arithmetic is so much faster, savings of up to
100X can be achieved.
 For very large n and very high precision, additional savings can beachieved by employing a “three-level” scheme: double precision,“intermediate precision” (typically 120 digits), and full multiprecision.

16

Details of two-level PSLQ
1. Perform standard PSLQ initialization.
2. Perform DP “re-initialization”: set A and B = I, y = y (scaled so max =1.0), and H = H (here underscore denotes DP approximations).
3. Perform LQ matrix factorization of H, and replace H with lower-diagonal.
4. Perform PSLQ iterations using DP arrays. When the largest entry of A
and B exceeds 1013, or if a very small entry is detected in the y vector,
update full-precision y, A, B, H arrays by matrix multiplication.
5. Check full-precision y vector for an entry that is zero (or smaller thansome acceptable epsilon). If so, a relation has been detected.
6. Compute norm bound.
7. Repeat beginning with step 2.
Notes:
 Computation of full-precision A matrix may be omitted.
 Full-precision H matrix will not be in lower-triangular form, and thus
cannot be used to compute norm bounds, but H can be used for this.
17

Some potential difficulties
 Occasionally a very large entry (> 253 = 9 x 1015 approx.) is produced in
the double precision A or B matrix, in spite of the 1013 cutoff.
 One must also check if any intermediate value produced in thecomputation of any entry of A or B exceeds 253.
 When this happens, it is necessary to:
1.
2.
3.
4.
5.

Abandon the current set of DP iterations.
Retreat to stored values of y, A, B and H at some previous iteration.
Update the multiprecision arrays using these stored values.
Perform a full LQ matrix factorization on H.
Perform multiprecision iterations, checking periodically (typically every ten
iterations) to see if the dynamic range of the y vector has reduced to thepoint that DP iterations are safe again.

 A three-level implementation saves time on very large problems above a
two-level scheme, mostly because intermediate precision can be used tohandle these special cases instead of full precision.

18Performance results for 1-2-3 level PSLQ:
Recover polynomial for alpha = 3^1/r – 2^1/s

r, s
5,5
5,6
6,6
6,7
7,7
7,8
8,8
8,9
9,9
9,10
10,10

n
26
31
37
43
50
57
65
73
82
91
101

Iterations
5143
9357
15217
25361
36947
60817
86684
124521
174140245443342931One-level
Digits
Time
180
32.37
240
105.48310
298.85420942.66500 2363.71

19

Two-level
Digits
Time
190
1.29
250
3.16
320
7.19
420
17.22
510
36.29
680
90.08
850
195.191050
425.671310
934.961620 2032.69
2000 4968.64

Three-level
Digits
Time

910
1120
1370
1680
2060

233.48460.34922.901780.65
3366.92

The “multi-pair” PSLQ algorithm
 The standard PSLQ algorithm (even the basic one-level scheme) is poorly
suited to parallel processing, because of recursions in the inner loop of
the reduction step.
 In attempt to devise a scheme better suited for parallel processing, DHB
developed the “multi-pair” PSLQ algorithm.
 The multi-pair PSLQ algorithm selects certain adjacent pairs of indicesthat can be reduced in one step, independently of the other pairs, thus
providing an opportunity for parallel processing.
 As it turns out, the multi-pair PSLQ algorithm runs faster even on a single
processor.
 Two-level and three-level multi-pair PSLQ schemes have been devised.
Parallel versions of these have also been devised, albeit it with great
difficulty.
 DHB has used “multi-pair” PSLQ (mostly 1- and 2-level) exclusively inrecent research work.
20

Statement of the multi-pair PSLQ
algorithm
Initialize:
1. For j := 1 to n: for i := 1 to n: if i = j then set Aij := 1 and Bij := 1
else set Aij := 0 and Bij := 0; endfor; endfor.
��n
2 ; endfor; set t = 1/s ; for k := 1 to
2. For k := 1 to n: set sk :=
x
1
j=k j
n: set yk := txk ; sk := tsk ; endfor.
3. Initial H: For j := 1 to n − 1: for i := 1 to j − 1: set Hij := 0; endfor;
set Hjj := sj+1 /sj ; for i := j + 1 to n: set Hij := −yi yj /(sj sj+1 ); endfor;
endfor.

21

Statement of the multi-pair PSLQ
algorithm, continued
Iterate.
1. Sort the entries of the (n − 1)-long vector {γ i |Hii |} in decreasing order,
producing the sort indices.
2. Beginning at the sort index m1 corresponding to the largest γ i |Hii |, selectpairs of indices (mi , mi + 1), where mi is the sort index. If at any stepeither mi or mi + 1 has already been selected, pass to the next index in
the list. Continue until either βn pairs have been selected, or the list
is exhausted. Let p denote the number of pairs actually selected in thismanner.
3. For i := 1 to p, exchange the entries of y indexed mi and mi + 1, and thecorresponding rows of A, B and H; endfor.
4. Remove� corners on H diagonal: For i := 1 to p: if mi ≤ n − 2 then set
22t0 := Hm
+ Hm
, t1 := Hmi ,mi /t0 and t2 := Hmi ,mi +1 /t0 ; fori ,mii ,mi +1
i := mi to n: set t3 := Hi,mi ; t4 := Hi,mi +1 ; Hi,mi := t1 t3 + t2 t4 ; andHi,mi +1 := −t2 t3 + t1 t4 ; endfor; endif; endfor.
22

Statement of the multi-pair PSLQ,
continued
6. Reduce H: For i := 2 to n: for j := 1 to n − i + 1: set l := i + j − 1; fork := j+1 to l−1: set Hlj := Hlj −Tlk Hkj ; endfor; set Tlj := nint(Hlj /Hjj )and Hlj := Hlj − Tlj Hjj ; endfor; endfor.
7. Update y: For j := 1 to n − 1: for i := j + 1 to n: set yj := yj + Tij yi ;endfor; endfor.
8. Update A and B: For k := 1 to n: for j := 1 to n − 1: for i := j + 1 to n:
set Aik := Aik − Tij Ajk and Bjk := Bjk + Tij Bik ; endfor; endfor; endfor.
9. Norm bound: Compute M := 1/ maxj |Hjj |. Then there can exist no
relation vector whose Euclidean norm is less than M .
10. Termination test: If the largest entry of A exceeds the level of numericprecision used, then precision is exhausted. If the smallest entry of the yvector is less than the detection threshold (see section 2), a relation has
been detected and is given in the corresponding row of B.

23

Performance results for multi-pair PSLQ:
Recover polynomial for alpha = 3^1/r – 2^1/s

r, s
5,5
5,6
6,6
6,7
7,7
7,8
8,8
8,9
9,9
9,10
10,10

n
26
31
37
43
50
57
65
73
82
91
101

Iterations
558
840
1136
1625
2071
2410
3723
4943
6169
7850
10017

One-level
Digits
Time
180
26.08
230
70.71
310
189.27400
479.07500 1130.85

24

Two-level
Digits
Time
180
1.48
240
3.43
310
7.84
410
17.22
500
35.64
660
69.39
800
169.621010
358.071260
744.201560 1556.37
1890 3283.08

Three-level
Digits
Time

880
1100
1320
1600
1950

214.66427.29804.511450.292747.12

Large test problems
1. “Fibonacci”:
Studies a conjecture of Broadhurst on the dimension of certain Euler sums.
n = 145.
5,000 digits.

2. “Bifurcation”:


Finds the polynomial satisfies by the 8-16 bifurcation in the logistic iteration.
n = 120.
10,000 digits.

3. “S(20)”:

Finds an analytic evaluation of certain constants involving Clausen values.
n = 118.
5,000 digits.

4. “Ladder”:

Finds a relation among constants associated with a root of Lehmer’s poly.
n = 125.
50,000 digits.
25

Parallel results on three test problems

Processors
1
2
4
8
16
32
48

Fibonacci
Time Speedup
47788
1.00
24665
1.94
12945
3.69
7076
6.75
4180
11.43
2994
15.96
2463
19.40

B4
Time Speedup
90855
1.00
46134
1.97
23966
3.79
12924
7.03
7424
12.24
4865
18.68
4049
22.44

26

S(20)
Time Speedup
23208
1.00
11973
1.94
6305
3.68
3470
6.69
2126
10.92
1548
14.99
1303
17.81

Work to be done These results are nearly 10 years old, employed an old all-Fortran
multiprecision library and the OpenMP parallel model.
 The ARPREC package does not lend itself to OpenMP parallelization –each individual “thread” requires too much context.
 The MPI model can be used, but only for very large problems, because of
greater data transfer overhead.
 The most likely parallel target for a parallel implementation is a sharedmemory, multicore system (such as recent Apple workstations).
 The UPC and Co-Array Fortran are more appropriate for such anenvironment.

27

Work to be done, continued
In addition, perhaps the whole approach to high-precision computation forexperimental math applications needs to be re-thought: One of the key developers of the QD and ARPREC libraries, is no longer
interesting in doing development or maintenance.
 Recent progress by the GMP and MPFR packages show good promise –they out-performs ARPREC in certain precision levels, particularly above
1000 digits.
 The ARPREC high-level C++ and Fortran-90 translation modules
continue to be very effective (high-level modules for GMP Can we “marry” the high-level ARPREC modules with the low-level GMP
modules?
Outside the box:
 Are there other, completely novel programming environments that weshould consider for this type of computation?
 Can we devise PSLQ-like algorithms for more general spaces?
28Summary The emerging “experimental” methodology in mathematics andmathematical physics often requires hundreds or even thousands ofdigits of precision.
 High-precision evaluation of integrals, followed by constant-recognition
techniques, has been a particularly fruitful area of recent research, withmany new results in pure math and mathematical physics.
 The PSLQ algorithm continues to be extremely effective in identifying
constants and in finding relationships between constants.
 Multi-level and multi-pair variants of PSLQ are much faster and more
effective than standard PSLQ.
 Much work needs to be done in developing new, parallel PSLQ software,and in improving multiprecision software in general.

29