A Thread-Safe, Arbitrary Precision Computation Package
David H. Bailey
http://www.davidhbailey.com
Lawrence Berkeley National Lab (retired)
University of California, Davis

March 17, 2015

Some applications where high precision is essential
1. Planetary orbit calculations (32 digits).
2. Supernova simulations (32–64 digits).
3. Climate modeling (32 digits).
4. Optimization problems in biology and other fields (32 digits).
5. Coulomb n-body atomic system simulations (32–120 digits).
6. Schrodinger solutions for lithium and helium atoms (32 digits).
7. Electromagnetic scattering theory (32–100 digits).
8. Scattering amplitudes of fundamental particles (32 digits).
9. Discrete dynamical systems (32 digits).
10. The Taylor algorithm for ODEs (100–600 digits).
11. Ising integrals from mathematical physics (100–1000 digits).
12. Problems in experimental mathematics (100–50,000 digits).

Coulomb n-body atomic system simulations
I Alexei Frolov of Queen’s University in Canada has used high-precision arithmetic

to solve a generalized eigenvalue problem that arises in Coulomb n-bodyinteractions.
I Matrices are typically 5,000 x 5,000 and are very nearly singular.
I Computations typically involve massive cancellation, and high-precision arithmetic

must be employed to obtain numerically reproducible results.
I Frolov has also computed elements of the Hamiltonian matrix and the overlap

matrix in four- and five-body systems.
I These computations typically require 120-digit arithmetic.

Frolov: “We can consider and solve the bound state few-body problems ... beyond ourimagination even four years ago.”
I A. M. Frolov and D. H. Bailey, “Highly accurate evaluation of the few-body auxiliary functions andfour-body integrals,” Journal of Physics B, vol. 36, no. 9 (14 May 2003), pg. 1857–1867.

Taylor’s method for ODEs with high-precision arithmetic
1 period - TIDES (16 digits)
16 periods -TIDES (300 digits)

45
40
35

z

30
25
20
15
10
5
0
−10

x

0

10

−20

0

y

20

First point TIDES (16 digits)First-Last point TIDES (300 digits)
Last point TIDES (16 digits)

Numerical integration of the L25-R25 unstable periodic orbit for the Lorenz model
during 16 time periods using TIDES with 300 digits, versus 1 time period using DP.

Experimental mathematics: Discovering new mathematical results
Methodology:
1. Compute various mathematical entities (limits, infinite series sums, definite
integrals, etc.) to high precision, typically 100–10,000 digits.
2. Use algorithms such as PSLQ to recognize these numerical values in terms ofwell-known mathematical constants.
3. When results are found experimentally, seek formal mathematical proofs of the
discovered relations.
Many results have recently been found using this methodology, both in puremathematics and in mathematical physics.
“If mathematics describes an objective world just like physics, there is no reason why
inductive methods should not be applied in mathematics just the same as inphysics.” – Kurt Godel

The PSLQ integer relation algorithm
Let (xn ) be a given vector of real numbers. An integer relation algorithm either finds
integers (an ) such thata1 x1 + a2 x2 + · · · + an xn = 0(to within the “epsilon” of the arithmetic being used), or else finds bounds within
which no relation can exist.
I The “PSLQ” algorithm of mathematician-sculptor Helaman Ferguson is the mostwidely used integer relation algorithm.
I Integer relation detection requires very high precision (at least n × d digits, whered is the size in digits of the largest ak ), both in the input data and the algorithm.
I The multipair PSLQ algorithm is parallelizable, but only at a fairly low level. Only
low-level, thread-based parallelism is efficient for this algorithm.
1. H. R. P. Ferguson, D. H. Bailey and S. Arno, “Analysis of PSLQ, An Integer Relation Finding Algorithm,”
Mathematics of Computation, vol. 68, no. 225 (Jan 1999), pg. 351–369.
2. D. H. Bailey and D. J. Broadhurst, “Parallel Integer Relation Detection: Techniques and Applications,”
Mathematics of Computation, vol. 70, no. 236 (Oct 2000), pg. 1719–1736.

High-precision numerical integration
Given f (x) defined on (−1, 1), define g (t) = tanh(π/2 sinh t). Then setting x = g (t)yields
Z 1
Z ∞
N
X
0
f (x) dx =f (g (t))g (t) dt ≈ h
wj f (xj ),
−1

−∞

j=−N

where xj = g (hj ) and wj = g 0 (hj ). Since g 0 (t) goes to zero very rapidly for large t, theproduct f (g (t))g 0 (t) typically is a nice bell-shaped function, so that the simple
summation above converges very rapidly. Reducing h by half typically doubles thenumber of correct digits.
We have found that tanh-sinh is the best general-purpose integration scheme forfunctions with vertical derivatives or singularities at endpoints, or for any function atvery high precision (> 1000 digits). Otherwise we use Gaussian quadrature.
1. D. H. Bailey, X. S. Li and K. Jeyabalan, “A Comparison of Three High-Precision Quadrature Schemes,”
Experimental Mathematics, vol. 14 (2005), no. 3, pg. 317–329.
2. H. Takahasi and M. Mori, “Double Exponential Formulas for Numerical Integration,” Publications of
RIMS, Kyoto University, vol. 9 (1974), pg. 721-741.

Ising integrals from mathematical physics
We applied our methods to study three classes of integrals: Cn are connected toquantum field theory, Dn arise in the Ising theory of mathematical physics, while the
En integrands are derived from Dn:
Z
Z ∞
du1
1
dun
4 ∞
···
···
Cn :=
2
P
n! 0
u
un
n
1
0
j=1 (uj + 1/uj )Q  ui −uj 2
Z ∞Z ∞i<j ui +uj
4
dun
du1
Dn :=
···
···
P
2
n! 0
u1
un
n
0
j=1 (uj + 1/uj )

2
Z 1
Z 1
Y uk − uj dt2 dt3 · · · dtn
En = 2···
uk + uj
0
0
1≤j<k≤nwhere in the last line uk = t1 t2 · · · tk .
D. H. Bailey, J. M. Borwein and R. E. Crandall, “Integrals of the Ising class,” Journal of Physics A:
Mathematical and General, vol. 39 (2006), pg. 12271–12302.

Limiting value of Cn : What is this number?
Key observation: The Cn integrals can be converted to one- dimensional integrals
involving the modified Bessel function K0 (t):Z
2n ∞ nCn =tK0 (t) dt
n! 0
High-precision numerical values, computed using this formula and tanh-sinhquadrature, approach a limit. For example:
C1024 = 0.6304735033743867961220401927108789043545870787 . . .
What is this number? We copied the first 50 digits into the online Inverse Symbolic
Calculator (ISC) at http://carma-lx1.newcastle.edu.au:8087. The result was:lim Cn = 2e −2γ .

n→∞where γ denotes Euler’s constant. This is now proven.

Other Ising integral evaluations found using PSLQ
D3 = 8 + 4π 2 /3 − 27 L−3 (2)D4 = 4π 2 /9 − 1/6 − 7ζ(3)/2E2 = 6 − 8 log 2E3 = 10 − 2π 2 − 8 log 2 + 32 log2 2
E4 = 22 − 82ζ(3) − 24 log 2 + 176 log2 2 − 256(log3 2)/3+16π 2 log 2 − 22π 2 /3
E5 = 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 ζ is the Riemann zeta function and Lin (x) is the polylog function. E5 remained a“numerical conjecture” for several years, but was proven in March 2014 by Erik Panzer.
E5 was reduced to a 3-D integral of a 60-line integrand, which was evaluated usingtanh-sinh quadrature to 250-digit arithmetic using over 1000 CPU-hours on a highly
parallel system. The PSLQ calculation required only seconds.

Algebraic numbers in Poisson potential functions
Lattice sums arising from the Poisson equation have been studied widely inmathematical physics and image processing. We numerically discovered, and thenproved, that for rational (x, y ), the 2-D Poisson potential function satisfies
1 X cos(mπx) cos(nπy )1log αφ2 (x, y ) ==π2
m2 + n2πm,n odd

where α is algebraic, i.e., the root of an integer polynomial
0 = a0 + a1 α + a2 α 2 + · · · + an α n
The minimal polynomials for these α were found by PSLQ calculations, with the(n + 1)-long vector (1, α, α2 , · · · , αn ) as input, where α = exp(πφ2 (x, y )). PSLQ
returns the vector of integer coefficients (a0 , a1 , a2 , · · · , an ) as output.
1. D. H. Bailey, J. M. Borwein, R. E. Crandall and J. Zucker, “Lattice sums arising from the Poisson
equation,” Journal of Physics A: Mathematical and Theoretical, vol. 46 (2013), pg. 115201,
http://www.davidhbailey.com/dhbpapers/PoissonLattice.pdf.
2. D. H. Bailey and J. M. Borwein, “Compressed lattice sums arising from the Poisson equation: Dedicatedto Professor Hari Sirvastava,” Boundary Value Problems, vol. 75 (2013), DOI:
10.1186/1687-2770-2013-75, http://www.boundaryvalueproblems.com/content/2013/1/75.

Samples of minimal polynomials found by PSLQ
k
5
6
7
8
9

10

Minimal polynomial for exp(8πφ2 (1/k, 1/k))
1 + 52α − 26α2 − 12α3 + α4
1 − 28α + 6α2 − 28α3 + α4
−1 − 196α + 1302α2 − 14756α3 + 15673α4 + 42168α5 − 111916α6 + 82264α7−35231α8 + 19852α9 − 2954α10 − 308α11 + 7α12
1 − 88α + 92α2 − 872α3 + 1990α4 − 872α5 + 92α6 − 88α7 + α8
−1 − 534α + 10923α2 − 342864α3 + 2304684α4 − 7820712α5 + 13729068α6−22321584α7 + 39775986α8 − 44431044α9 + 19899882α10 + 3546576α11−8458020α12 + 4009176α13 − 273348α14 + 121392α15−11385α16 − 342α17 + 3α18
1 − 216α + 860α2 − 744α3 + 454α4 − 744α5 + 860α6 − 216α7 + α8

The minimal polynomial for exp(8πφ2 (1/32, 1/32)) has degree 128, with individual
coefficients ranging from 1 to over 1056 (see next viewgraph). This computation
required 10,000-digit precision. Most of the time was spent on PSLQ, not quadrature.
Other polynomials required up to 50,000-digit precision, with runs extending to several
days. The ARPREC software failed on at least one of these problems.

Degree-128 minimal polynomial found for k = 32

What is needed for high-precision floating-point software?

1. A rock-solid-reliable arithmetic engine, with precision scalable to 1,000,000 digits
or more — 1,000 or even 10,000 digits is no longer enough.
2. Highly efficient algorithms, so single-processor performance is as fast as possible,
given constraints of portability and ease of installation.
3. A thread-safe design to facilitate parallel processing, both thread-based (low-level)
and message-passing (high-level) parallelism.
4. A comprehensive library of tuned transcendentals and special functions:sin, cos, exp, gamma, incomplete gamma, polylogarithms, Bessel functions, etc.
5. High-level language interfaces for C++, Fortran-90 and other languages.
6. Interoperability with Maple and Mathematica.

Free software for high-precision computation
1. ARPREC: Arbitrary precision, with numerous algebraic and transcendental
functions. High-level interfaces for C++ and Fortran-90.
http://crd.lbl.gov/~dhbailey/mpdist.
2. GFORTRAN: Now provides full REAL*16 (IEEE 128-bit, or 34-digit) support.
3. GMP: Produced by a volunteer effort and distributed under the GNU license.
http://gmplib.org.
4. MPFR: C library for multiple-precision floating-point computations with exact
rounding, based on GMP. http://www.mpfr.org.
5. MPFR++: High-level C++ interface to MPFR.
http://perso.ens-lyon.fr/nathalie.revol/software.html.
6. MPFUN90: Similar to ARPREC, but is written entirely in Fortran-90 and provides
only a Fortran-90 interface. http://crd.lbl.gov/~dhbailey/mpdist.
7. QD: Performs “double-double” (31 digits) and “quad-double” (62 digits)
arithmetic. High-level interfaces for C++ and Fortran-90.
http://crd.lbl.gov/~dhbailey/mpdist.

Lessons learned from my own packages: ARPREC and QD
I Complicated configure files and install scripts are a problem — many users are

reluctant to use these packages for this reason. Many support issues are for these
items, not the underlying multiprecision software.
I My older package MPFUN90 (all Fortran-90) continues to be very popular

because of its simplicity — downloads and citations still exceed those of ARPREC.
I Support is a nagging issue, given that Yozo Hida and Brandon Thompson are nolonger available to help, and Xiaoye Li is very busy with other tasks.
I Special system settings and system-dependent features are problematic.
I In some cases, language features are not supported by certain systems.
I Double precision constants and expressions in code, e.g., 0.1d0, are notautomatically converted to high precision, often producing unexpectedly poor
accuracy in multiprecision results. Can something be done here?

Thread safety: A major challenge for arbitrary precision software
Since run times are greatly magnified by high-precision arithmetic, parallel processing isincreasingly essential. However, available arbitrary-precision packages are notthread-safe at the application user level, and this kills thread-based parallelization, themost efficient form of parallelization:
I All high-level packages I am aware of use global variables, and these variables aredeadly for thread safety.
I Most of these global variables can be dispensed with, but not the working

precision level, which is changed frequently within the package itself and often byusers also.
I The only reasonable design for user programming is via operator overloading,available in languages such as Fortran-90 and C++, which do not permit theworking precision level to be passed.
What can be done?

MPFUN2015: A NEW thread-safe arbitrary precision package
I 100% thread safe. NO global read/write global variables.
I NO initialization required (unless over 20,000 digits required).
I NO system-dependent features or settings — the code only assumes 64-bit IEEE

format (it does not depend on rounding mode).
I NO reliance on sophisticated language features that might not be supported.
I Installation is a snap — one compile line.
I Precision is scalable to over billion digits; largest size is over 1064,000,000 .
I A full package of basic arithmetic functions. Higher-level routines forbinary-decimal conversion, I/O and transcendentals.
I A full-featured high-level language interface for Fortran-90 is provided, so thatmost users need only make minor changes to existing double precision code.
I Extensive error checking and a solution to the double precision accuracy problem.
I Currently only Fortran-90 version; C++ version is being developed in Australia.
I Available at: http://www.davidhbailey.com.

The MPFUN2015 multiprecision data structure
Each multiprecision (MP) datum consists of a vector of double-precision (DP) floats:
1. Word 0: Total space allocated for array, in words.
2. Word 1: Working precision level (in words) associated with this data.
3. Word 2: Number of mantissa words N; sign of word 2 is sign of value.
4. Word 3: Exponent, base 248 .
5. Word 4 to N + 3: Mantissa (whole numbers between 0 and 248 − 1).
6. Word N + 4 and N + 5: scratch; set to zero in output arrays.
Example: (12., 6., −2., 0., 3., 70368744177664., 0., 0., 0., 0., 0., 0.) represents −3.25(since 70368744177664 = 0.25 × 248 ). The working precision associated with this data
is 6 (as specified in word 1). Note that N mantissa words require at least N + 6 totalarray space.

Solution to the thread safety problem
Note that word 1 of each multiprecision datum is the working precision level associatedwith that datum (variable or array element). This solves the thread safety problem,
although it requires a somewhat different programming style at the user level:
I The user sets a default (maximum) working precision level in a compile-time datastatement.
I The result of any operation involving MP variables or array elements “inherits”
the working precision level of the input operands.
I In the case of an operation with two or more operands with different working
precision levels, the largest is chosen for the result.
I The built-in function mpreal allows one to convert a MP variable or arrayelement to one with a different working precision level.
I The built-in function mpwprec returns the working precision level currently
assigned to any MP variable or array element.
I Complications: Assignments to DP variables or constants, reading values fromfile, etc. (See later viewgraphs.)Low-level routine module
The following low-level routines are provided:
I Addition, subtraction, multiplication and division of MP arguments.
I DP-to-MP and MP-to-DP conversion; MP times DP and MP divided by DP.
I Square root routine employs a Newton iteration with a dynamically varying

precision level (roughly doubling with each iteration).
I n-th powers and n-roots.
I Nearest integer.
I Rounding and normalization.

Binary-to-decimal and I/O routine module

I Convert MP to decimal character string in En1.n2 format.
I Convert MP to decimal character string in Fn1.n2 format.
I Convert decimal character string to MP (may be in E or F format).
I Input decimal value from input file (may be in E or F format).
I Output MP variable or array element to file in E format.

Transcendental constant and function module
I log(2) and π to arbitrary precision. For requested precision levels less than 20,000

digits, the binary value is merely copied from a compile-time data statement.
I exp(x) and log(x).
I cos(x) and sin(x).
I cosh(x) and sinh(x).
I Conversion of (x, y ) to angle in (−π, π).
I Routines use somewhat tricky variations of relatively conventional algorithms.

None of these routines requires any precalculated tables. This preserves thread safety.

Thread-safe algorithms for exp(x) and sin(x)Arbitrary precision packages require an initialization step (not thread safe!) to
precalculate data for exp and cos/sin. MPFUN2015 does not:I The exp routine first reduces the input argument to within the interval(− log(2)/2, log(2)/2). Then it divides this result by 210 = 1024, producing a verysmall value, which is then input to the Taylor’s series for exp(x). When complete,
the result is squared ten times, and then corrected for the initial reduction.
I The sine routine first reduces the result to within (−π, π), and then to the nearest

quadrant. The result is then divided by 256, producing a very small value, whichis thenq input to the Taylor’s series for sin(x). When complete, cos(x) is computed
as 1 − sin2 (x), and then the double-angle formulas for sin(x) and cos(x) areapplied eight times. The result is adjusted for quadrant.

Note that neither scheme requires any precalculated data. Run time performanceappears to be within a factor of two of ARPREC schemes, which precalculate extensive
tables in an (thread-unsafe) initialization step.

High-level language interface moduleThis module contains the Fortran language interfaces that connect to all lower-level
functions from user source code:
I Binary operations (+, −, ×, ÷).
I Intrinsics (sin, cos, exp, etc.).
I I/O facilities.
I Some special functions and subroutines provided as a part this package.

Easy compilation and linking

The package is simple to compile. Using the gfortran compiler, for instance, type:gfortran -O2 -c mpmodule.f90 mpfuna.f90 mpfunb.f90 mpfunc.f90 mpfund.f90 mpfunf.f90

You may have to repeat this command to fully resolve the modules. Thereafter, tocompile and link an application program prog, type:gfortran -O2 prog.f90 mpmodule.o mpfuna.o mpfunb.o mpfunc.o mpfund.o mpfunf.oBasic Fortran application programming instructions
Almost any working Fortran-90 code using double precision (DP) arithmetic can beeasily converted to use the package:
I To invoke the package, place this line in every subprogram that contains a

multiprecision variable or array:use mpmodule
I To designate a variable or array as MP, use the Fortran-90 type statement:

type (mp real) a, b(m), c(m,n)
I All basic operations are supported just as in conventional code, e.g.:d = a + b(i) * sqrt(3.d0 - c(i,j))
I The following standard Fortran-90 intrinsics are supported with MP arguments:

abs, acos, aint, anint, asin, atan, atan2, cos, cosh, dble, exp, log,
max, min, sign, sin, sqrt, tan, tanhHandling double precision constants and expressions
I Problem: The code a= 3.14159d0, where the LHS is MP, does NOT produce
the full-precision equivalent of 3.14159, since by standard rules of precedence, theRHS is evaluated in DP, then converted (by zero extension) to LHS.

I Solution: Write this as a

= ’3.14159’. In MPFUN2015, this forces the constant
to be evaluated to full precision.

I Related problem: The code b= a + 3.d0 * sqrt (2.d0) does NOT produce afully accurate result, since 3.d0 * sqrt (2.d0) is performed in DP.

I Solution: Write this as b

= a + 3.d0 * sqrt (mpreal (2.d0)). This forces
all operations to be done using MP routines.

Run-time solution to DP accuracy problems: The MPFUN2015 software checks all
DP entities (constants, variables and expression values) at execution time to see if they
have more than 40 significant bits. If so, it is flagged as an error. This feature catches
99.99% of the precision loss problems due to DP usage.
If you are sure the DP value is OK as is, use special routines (see next viewgraph).

Functions and subroutine provided by the packageNamempcsshf
mpcssnf
mplog2
mpeform
mpfform
mpnrtfmppi
mpprodd
mpquotd
mpreadmprealmpreald
mpwprec
mpwriteType
Subroutine
Subroutine
MP function
Subroutine
Subroutine
MP function
MP function
MP function
MP function
Subroutine
MP function
MP function
Int function
Subroutine

Description
Returns both cosh and sinh of MP argument.
Returns both cos and sin of MP argument.
Returns log(2).
Returns a character string with MP argument in E format.
Returns a character string with MP argument in F format.
Returns the n-th root of MP argument.
Returns π.
Returns the product of MP and DP args (without checking).
Returns the quotient of MP and DP args (without checking).
Reads MP argument from file.
Converts the argument (DP, integer or string) to MP.
Converts the argument (DP) to MP (without checking).
Returns working precision (in words) assigned to MP argument.
Writes MP argument to output unit in E format.

Dynamically changing the working precision
There are two versions of the language interface module:
1. Version 1 (in file mpfunf1.f90): Precision specifications are optional. Intended for
applications with only occasional changes in the working precision level, or none.
2. Version 2: (in file mpfunf2.f90): Precision specifications are required. Intended for
applications that frequently change the working precision level.
In particular, with Version 1:
I Assignments of the form a = x, where a is MP and x is DP, integer or literal, arepermitted (a is converted to MP with the default precision level).
I The functions mpreal (which converts a DP, integer or literal argument to MP),mpread (which inputs MP data from file), mplog2 (which returns log(2)), andmppi (which returns π) have an (optional) integer working precision argument.
With Version 2:
I Assignments of the form a = x, where a is MP and x is DP, integer or literal, are
not permitted. You must use mpreal with precision level parameter instead.
I The optional precision level arguments in the above routines are required.

Working applications
One-, two- and three-level multipair PSLQ. These are workhorses of my own research.
They are very unforgiving of even minor arithmetic errors:
I One-level multipair PSLQ (907 lines): Constant precision level.
I Two-level multipair PSLQ (1682 lines): Changes between MP and DP.
I Three-level multipair PSLQ (2103 lines): Changes between full precision (1000
digits), medium precision (120 digits) and DP; uses Version 2 of MPFUNF.
The first two were easy to convert — the third was more work.
The tanh-sinh quadrature test suite (911 lines):
I Tests the tanh-sinh quadrature routine on a suite of 15 integrals, producing (inmost cases) 400-digit results.
I Many of the integrands are complicated expressions involving varioustranscendental functions.
I Dynamically changes the working precision level between 400 and 1000 digits:uses version 2 of MPFUNF.
I All tests appear to work satisfactorily.

Remaining work

I Additional testing and documentation.
I Support for the complex MP datatype.
I Implementation of routines for extra-high precision computation (e.g., FFT-based

multiplication, quadratically convergent transcendentals, etc.).
I High-level language interfaces for other languages: C++, etc.
I Explore parallel implementations of high-level applications using the package.

Try it (and please send me bug reports)
The MPFUN2015 software is available at:
http://crd-legacy.lbl.gov/~dhbailey/mpdist.
The download zip file currently includes:
I All MPFUN2015 Fortran modules.
I Three PSLQ codes and one tanh-sinh code, together with sample output files forcomparison with your own runs.
I A README.txt file with detailed instructions for installation and coding.
I A full technical paper will be available shortly.

This talk is available at
http://www.davidhbailey.com/dhbtalks/dhb-mpfun2015.pdf.