MPFUN2020: A thread-safe arbitrary
precision package with special functions(Full Documentation)David H. Bailey∗
May 30, 2024

Abstract
Numerous research studies have arisen, particularly in mathematical physics and experimental mathematics, that require extremely
high numeric precision — hundreds or thousands of digits in somecases. This paper describes a new arbitrary precision software package (“MPFUN2020”) to support such applications. It comes in two
versions: a self-contained all-Fortran version that is easy to install,and a version based on the MPFR package that runs faster although
it is significantly more complicated to install.
Both versions feature: (a) a completely thread-safe design, so user
codes can be converted for parallel execution at the application level;(b) a full-featured high-level Fortran interface, so that most applications can be converted to multiprecision with only minor changes tosource code; (c) full support for both real and complex datatypes; (d)an extensive set of transcendental and special functions; (e) run-time
checking to guard against problems with converting double precision
data; (f) a medium precision datatype that improves performance onlarge variable precision applications; and (g) interoperability, so thatwith a simple restriction, application codes written for one version canbe run using the other without any code changes.
∗

Lawrence Berkeley National Laboratory (retired), 1 Cyclotron Road, Berkeley, CA
94720, USA. E-mail: dhbailey@lbl.gov.

1Contents
1 Applications of high-precision computation
1.1 The PSLQ integer relation algorithm . . . . . . . . . . . . . .
1.2 High-precision numerical integration . . . . . . . . . . . . . .
1.3 Ising integrals . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Algebraic numbers in Poisson potential functions . . . . . . .
2 Current high-precision software
2.1 Thread safety and parallel implementations

3
3
4
4
6

8
. . . . . . . . . . 10

3 MPFUN2020: A new thread-safe package
3.1 Data structure . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Modules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 The MPFUN2020 solution to thread safety . . . . . . . . . . .

11
11
12
13

4 Installation, compilation and linking
14
4.1 The two variants . . . . . . . . . . . . . . . . . . . . . . . . . 14
5 Coding Fortran applications
5.1 Basic instructions for using the package . . . . . . . . . . . . .
5.2 Extra-high precision . . . . . . . . . . . . . . . . . . . . . . .
5.3 Functions and subroutines . . . . . . . . . . . . . . . . . . . .
5.4 Special functions . . . . . . . . . . . . . . . . . . . . . . . . .
5.5 Polynomial evaluation and solution . . . . . . . . . . . . . . .
5.6 Input and output of multiprecision data . . . . . . . . . . . .
5.7 Handling double precision values . . . . . . . . . . . . . . . .
5.8 Support for quad precision . . . . . . . . . . . . . . . . . . . .
5.9 Dynamically changing the working precision . . . . . . . . . .
5.10 Medium precision datatype . . . . . . . . . . . . . . . . . . . .

15
15
16
18
18
19
19
23
26
27
30

6 Sample applications and performance31
6.1 Timings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
7 Appendix: Numerical algorithms
7.1 Algorithms for basic arithmetic . . . . . . . . . . . . . . . . .
7.2 FFT-based multiplication . . . . . . . . . . . . . . . . . . . .
7.3 Transcendental constants and functions . . . . . . . . . . . . .
7.4 Special functions . . . . . . . . . . . . . . . . . . . . . . . . .
2

36
36
37
38
41

1

Applications of high-precision computation

For many scientific calculations, particularly those that employ empirical
data, IEEE 32-bit floating-point arithmetic is sufficiently accurate, and ispreferred since it saves memory, run time and energy usage. For other calculations, 64-bit floating-point arithmetic is required to produce results ofsufficient accuracy; still others switch between 32-bit and 64-bit.
Another group of applications, particularly in the fields of mathematicalphysics and experimental mathematics, require even higher precision — tens,
hundreds or even thousands of digits:
1. Supernova simulations (32–64 digits).
2. Optimization problems in biology and other fields (32–64 digits).
3. Coulomb n-body atomic system simulations (32–120 digits).
4. Electromagnetic scattering theory (32–100 digits).
5. The Taylor algorithm for ODEs (100–600 digits).
6. Ising integrals from mathematical physics (100–1000 digits).
7. Problems in experimental mathematics (100–50,000 digits and higher).
These applications are described in greater detail in [7, 5], which provides
detailed references. Here is a brief overview of several key applications:

1.1

The PSLQ integer relation algorithm

Very high-precision floating-point arithmetic is now considered an indispensable tool in experimental mathematics and mathematical physics [7]. Manyof these computations involve variants of Ferguson’s PSLQ integer relation
detection algorithm [25, 15]. Suppose one is given an n-long vector (xi ) ofreal or complex numbers (presented as a vector of high-precision values). The
PSLQ algorithm attempts to find nontrivial integer coefficients (ai ), if they
exist, such thata1 x 1 + a2 x 2 + · · · + an x n = 0
to available precision. Some researchers adapt variations of the LenstraLenstra-Lovasz (LLL) lattice basis reduction algorithm [32], or the “HJLS”
3

algorithm, which is based on LLL, for integer relation detection [21]. Whicheveralgorithm is used, integer relation detection requires very high precision—atleast n · maxi log10 |ai | digits, or else the true relation will be lost in a sea ofnumerical artifacts.

1.2

High-precision numerical integration

One of the most fruitful applications of the experimental methodology and
the PSLQ integer relation algorithm has been to identify classes of definite
integrals, based on very high-precision numerical values, in terms of simple
analytic expressions.
These studies typically employ either Gaussian quadrature or the tanhsinh quadrature scheme of Takahasi and Mori [38, 4]. The tanh-sinh quadrature algorithm approximates the integral of a function on (−1, 1) as
Z 1
f (x) dx ≈ h−1NX

wj f (xj ),

(1)

j=−N

where the abscissas xj and weights wj are given by
xj = tanh (π/2 · sinh(hj))
wj = π/2 · cosh(hj)/ cosh (π/2 · sinh(hj))2 ,(2)and where N is chosen large enough that summation terms in (1) beyond N(positive and negative) are smaller than the “epsilon” of the numeric precision
being used. Full details are given in [4]. An overview of applications of highprecision integration in experimental mathematics is given in [9].

1.3Ising integrals

In one study, tanh-sinh quadrature and PSLQ were employed to study thefollowing classes of integrals [11]. The Cn are connected to quantum field
theory, the Dn integrals arise in the Ising theory of mathematical physics,

4

while the En integrands are derived from Dn :
Z
Z ∞
du1
dun
4 ∞
1
···
Cn =
···
P
2
n! 0
u1
un
n
0
j=1 (uj + 1/uj )Q  ui −uj 2
Z ∞Z ∞i<j ui +uj
4
du1
dun
···
Dn =
···
P
2
n! 0
u1
un
n
0
j=1 (uj + 1/uj )
!2
Z 1
Z 1
Y uk − uj
En = 2···dt2 dt3 · · · dtn .
uk + uj
0
0
1≤j<k≤nQ
In the last line uk = ki=1 ti .
In general, it is very difficult to compute high-precision numerical values ofn-dimensional integrals such as these. But as it turn out, the Cn integrals canbe converted to one-dimensional integrals, which are amenable to evaluation
with the tanh-sinh scheme:
Z2n ∞pK0n (p) dp.
Cn =
n! 0
Here K0 is the modified Bessel function (see Section 7.4). 1000-digit values ofthese sufficed to identify the first few instances of Cn in terms of well-known
constants. For example, C4 = 7ζ(3)/12, where ζ denotes the Riemann zetafunction. For larger n, it quickly became clear that the Cn approach the limit
lim Cn = 0.630473503374386796122040192710 . . . .

n→∞This numerical value was quickly identified, using the Inverse Symbolic Calculator 2.0 (an online tool that unfortunately is no longer available), aslim Cn = 2e−2γ ,n→∞where γ is Euler’s constant. This identity was then proven in [11].

5Other specific results found in this study include the following:
D3
D4
E2
E3
E4
E5

8 + 4π 2 /3 − 27 L−3 (2)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 ζ is the Riemann zeta function and Lin (x) is the polylogarithm function(see Section 7.4).
E5 was computed by first reducing it to a 3-D integral of a 60-line integrand, which was evaluated using tanh-sinh quadrature to 250-digit arithmetic using over 1000 CPU-hours on a highly parallel system. The PSLQ
calculation required only seconds to produce the relation above. This formularemained a “numerical conjecture” for several years, but was proven in March2014 by Erik Panzer, who mentioned that he relied on these computational
results to guide his research [35].

1.4

Algebraic numbers in Poisson potential functions

The Poisson potential function arises in several contexts of mathematicalphysics and engineering, including the analysis of crystal structures and eventhe sharpening of photographic images. In two recent studies [13, 12], researchers explored a simple two-dimensional instance:
φ2 (x, y) =1 X cos(mπx) cos(nπy).
π 2 m,n odd
m2 + n2(3)

After detailed numerical experimentation, these researchers discovered andthen proved the remarkable fact that when x and y are rational,
φ2 (x, y) =1
log β(x, y),π(4)where β(x, y) is an algebraic number, namely the root of an algebraic equation
with integer coefficients of some degree m.
6

These experiments computed α = (β(x, y))8 = exp(8πφ2 (x, y)), using
some rapidly convergent formulas found in [13], for various simple rationals
x and y. Then the vector (1, α, α2 , · · · , αm ) was input to a program implementing the three-level multipair PSLQ program [15]. When successful, theprogram returned the vector of integer coefficients (a0 , a1 , a2 , · · · , am ) of apolynomial satisfied by α as output. Table 1 shows some examples [13].
s
5
6
7
8
9

10

Minimal polynomial corresponding to x = y = 1/s:
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

Table 1: Sample of polynomials produced in earlier study [13].

Using this and other data, Jason Kimberley of the University of Newcastle, Australia, conjectured a formula that gives the degree m as a functionof s [13]. These computations required prodigiously high precision (up to
10,000 digits) and very long run times.
In a subsequent study [14], the present author and three other researchers
extended these studies, confirming Kimberley’s formula for all 10 ≤ s ≤ 40
and also for s = 42, 44, 46, 48, 50, 52, 60 and 64. These runs, which required
precision levels up to 64,000 digits, were facilitated by an improved PSLQ
code, the faster MPFUN2015 package described below in Section 2, as wellas parallel execution facilitated by the thread-safe feature of this package.
Examination of these results led to additional discoveries, and ultimately toa full proof of Kimberley’s conjecture, at least for the special case x = y = 1/s
for integer s, as well as other facts about these curious polynomials, such as
the fact that when s is even, the polynomials for the case x = y = 1/s arepalindromic (left-to-right symmetric) [14].
In two more recent studies [2, 3], these results were further extended, using
the arbitrary precision package described in this paper, to include the vastly7

larger class of rationals x = p/s and y = q/s, for 1 ≤ p, q < s/2 ≤ 38, and also
for s = 40, 42 and 50. These computations revealed that Kimberley’s formula
does not give the correct degree for the larger class, but that a generalizationof Kimberley’s rule does hold. Two other new general features of thesepolynomials were also identified. See [3] for these results and additional
details.

2

Current high-precision software

By far the most common form of extra-precision arithmetic is roughly twice
the level of standard 64-bit IEEE floating-point arithmetic. One option is
the IEEE standard for 128-bit binary floating-point arithmetic, with 113
mantissa bits; sadly it is not yet widely implemented in hardware, although
it is supported, in software, by some compilers. Another option for thislevel of precision is “double-double” arithmetic (approximately 31 digits),
which consists of two 64-bit IEEE floats, or even quad-double arithmetic(approximately 62 digits), which consists of four IEEE 64-bit floats (see the
QD package in the list below) [30].
For higher-levels of precision, software packages typically represent a highprecision datum as a string of floats or integers, where the first few words
contain bookkeeping information and the binary exponent, and subsequentwords contain the binary mantissa.
Software for performing high-precision arithmetic has been available forquite some time, for example in the commercial packages Mathematica and
Maple, but until recently those using standard languages for high-performancecomputing, such as C, C++ or Fortran, often found necessary to rewrite theirapplications, replacing each arithmetic operation with a subroutine call to acustom-coded library. Today there are several freely available high-precision
software packages for high-precision computing, and even some high-level
language interfaces that make code conversions relatively painless:• ARPREC: [By the present author; superseded by MPFUN2015 and
MPFUN2020] Supports arbitrary precision integer, real, complex andtranscendental functions, with high-level C++ and Fortran-90 interfaces; available at https://www.davidhbailey.com/dhbsoftware/.
• CLN: Supports arbitrary precision integer, real, complex and transcendental functions in C++; available at http://www.ginac.de/CLN.
8• GMP: Supports low-level high-precision integer, rational and floatingpoint calculations; available at http://gmplib.org.
• Julia: High-level programming environment incorporating GMP and
MPFR; available at http://julialang.org.
• MPFR: Supports (at a low level) multiple-precision floating-point computations with correct rounding, based on GMP; includes numerous
algebraic and transcendental functions, and a thread-safe build option;achieves very fast timings; available at http://www.mpfr.org.
• MPFR++: High-level C++ interface to MPFR (although currently
available version is not up-to-date with MPFR); available at http://perso.ens-lyon.fr/nathalie.revol/software.html.
• MPFR C++: High-level C++ interface to MPFR with a thread-safe
option; available at http://www.holoborodko.com/pavel/mpfr.
• mpmath: Python library for arbitrary precision floating-point arithmetic, with transcendentals; available at https://code.google.com/
p/mpmath.
• NTL: C++ library for arbitrary precision integer and floating-point
arithmetic; available at http://www.shoup.net/ntl.
• Pari/GP: Computer algebra system including high-precision arithmetic
and transcendental functions; available at http://pari.math.u-bordeaux.
fr.
• QD: [By the present author] Supports “double-double” (approx. 31 digits) and “quad-double” (approx. 62 digits) arithmetic, with transcendental functions; includes high-level interfaces for C++ and Fortran-90;
available at https://www.davidhbailey.com/dhbsoftware/.
• Sage: Open-source symbolic computing system including high-precision
facilities; available at http://www.sagemath.org.
• MPFUN2015: [By the present author; superseded by MPFUN2020] A
thread-safe arbitrary precision package for Fortran application, withan all-Fortran version and a version based on MPFR for faster performance; available at https://www.davidhbailey.com/dhbsoftware/.
92.1

Thread safety and parallel implementations

The scientific computing world is moving rapidly into multicore and multinode parallel computing, so that future improvements in performance on
high-precision computations will only be obtained by aggressively exploiting
parallelism. On modern multicore processors, parallel computing is most efficiently performed using a shared memory, multithreaded environment suchas OpenMP [39] within a single node, even if a message passing system such
as MPI is employed for parallelism between nodes.
Computations that use a thread-parallel environment such as OpenMP
must be entirely “thread-safe.” One impediment to thread safety for multiprecision applications is the design of the operator overloading feature (i.e.,
extending +, −, ×, ÷ to multiprecision data) of modern computer languages,
which typically does not permit one to carry information such as the current working precision level. More importantly, most existing high-precision
packages generate a “context” of auxiliary data, such as the current working
precision level and data to support transcendental function evaluation, which
context typically ruins thread safety unless special care is taken.
Of the packages listed in Section 2, only one is both thread-safe andsupports arbitrary precision floating-point calculations with a high-level interface, namely the MPFR C++ package [33]. This package is built upon thelower-level MPFR package [27], which in turn is well-designed, features correct rounding to the last bit, includes numerous transcendental and special
functions, and achieves the the fastest overall timings of any floating-point
package in the above list [22]. (However, the MPFR package is only threadsafe when a special option is enabled in the installation step.)
There was, to this author’s knowledge, no high-level, thread-safe arbitrary
precision package to support Fortran applications, prior to MPFUN2015,
which was mentioned in the above list. However, the present author will be
the first to acknowledge that MPFUN2015 has some shortcomings, notablythat the run-time performance of the all-Fortran version is typically five times
or more slower than the version based on MPFR. Another weakness is thatin both versions of the package, significant memory is wasted for applications
where part of multiprecision arrays require only a medium level of precision.
Further, the MPFUN2015 package lacks a full package of special functions
and polynomial solutions routines.

10

3

MPFUN2020: A new thread-safe package

To that end, the present author has written a new thread-safe arbitrary
precision package. As with the MPFUFN2015 package, it is offered in two
versions: an easy-to-install all-Fortran version (MPFUN20-Fort), which runs
several times faster than the all-Fortran version of MPFUN2015, and alsoa version based on the MPFR package (MPFUN20-MPFR), which is evenfaster and more complete than the MPFR version of MPFUN2015.
Both versions feature: (a) a completely thread-safe design, so user codes
can be converted for parallel execution at the application level; (b) a fullfeatured high-level Fortran interface, so that most applications can be converted to multiprecision with only minor changes to source code; (c) full
support for both real and complex datatypes; (d) an extensive set of transcendental and special functions; (e) run-time checking to guard against problems with converting double precision data; (f) a medium precision datatypethat improves performance on large variable precision applications; and (g)interoperability, so that with a simple restriction, application codes written
for one version can be run using the other without any code changes.
The sharply improved performance of the all-Fortran MPFUN20-Fort version is due primarily to changing the underlying design of the package to be
based on 64-bit integer operations, rather than on 64-bit floating-point operations as in the earlier package. Additionally, the MPFUN20-Fort version
employs FFT-based arithmetic for significantly faster execution at very highprecision levels — see Sections 7.2. The new MPFUN20-MPFR version isfaster than the MPFUN20-Fort version on most applications, although it
is significantly more complicated to install, because the GMP and MPFR
packages must first be installed, usually requiring administrator privilege.

3.1

Data structure

For the MPFUN20-Fort version, the structure is a (N + 6)-long vector of
64-bit integers, where N is the number of mantissa words:• Word 0: Total space allocated for this array, in 64-bit integer words.
• Word 1: The working precision level (in words) associated with this
data.

11

• Word 2: The number of mantissa words N ; the sign of word 2 is thesign of the value.
• Word 3: The multiprecision exponent, base 260 .
• Word 4 to N + 3: N mantissa words (whole numbers between 0 and
260 − 1).
• Word N + 4 and N + 5: Scratch words for internal usage.
For the MPFUN20-MPFR version, the structure is a (N + 6)-long vectorof 64-bit integers, where N is the number of mantissa words:• Word 0: Total space allocated for this array, in 64-bit words.
• Word 1: The working precision level (in bits) associated with this data.
• Word 2: The sign of the value.
• Word 3: The exponent, base 2.
• Word 4: A pointer to the first word of the mantissa, which in MPFUN20MPFR always points to Word 5.
• Word 5 to N + 4: Mantissa words (unsigned integers between 0 and
264 − 1).
• N + 5: Not used at present.
Note that in the MPFUN20-MPFR version, words 1 through N + 4 correspond exactly to the data structure of the MPFR package.
For each version, a complex multiprecision datatype is a contiguous pairof real multiprecision data. The imaginary member of the real-imaginary
pair starts at an offset in the array equal to the value of word 0. Note that
this offset is independent of the working precision.

3.2

Modules

The MPFUN20-Fort and MPFUN20-MPFR versions both include the following separate modules, each in its own source file:121. MPFUNA: Contains compile-time global data, including the binary
values of log(2), π and γ (Euler’s constant), to 20,000-digit precision.
2. MPFUNB (MPFUN20-Fort only): Handles basic arithmetic functions,rounding, normalization, square roots and n-th roots.
3. MPFUNC (MPFUN20-Fort only): Handles binary-to-decimal conversion, decimal-to-binary conversion and input/output operations.
4. MPFUND: Includes routines for common transcendental constants andfunctions (see Table 3 and Section 7.3).
5. MPFUNE: Includes routines for special functions, including Bessel
functions, the error function, the gamma function, the zeta function
and numerous others (see Table 5 and Section 7.4).
6. MPFUNF: Defines the default standard precision and medium precisionlevels, in digits, and also the equivalent levels in words.
7. MPFUNG: A high-level user interface that provides support for the
standard high-precision datatype.
8. MPFUNH: A high-level user interface that provides support for the
medium precision datatype.
9. MPMODULE: The main module that references the others; in normal
usage it is the only module that is directly referenced by the user.

3.3

The MPFUN2020 solution to thread safety

All of the software modules above are 100% thread safe. There are no global
parameters or arrays, except for static, compile-time data, and no initialization is required unless extremely high precision is required (over 20,000
digits). The working precision level is always passed as a subroutine argument, ensuring thread safety. The MPFUN20-MPFR version is thread-safe
provided that the MPFR package is compiled with the thread-safe option.
Thread safety at the language interface or user level in both versions isachieved by assigning a working precision level to each multiprecision datum,
which then is passed through the multiprecision software. Note, in the data
structures given in Section 3.1 above, that word 1 (the second word of the13array) is the working precision level associated with that datum. This solves
the thread safety problem when precision is dynamically changed in the application, although it requires a slightly different programming style for some
applications, as we shall explain below (see Section 5.9).
All computations are performed to a fixed precision set by the user inmodule MPFUNF (default is 4,000 digits), unless the user specifies a lower
value for certain calculations (see Section 5.9 for details) or otherwise employsa variable precision design. The result of any operation involving multiprecision variables or array elements inherits the working precision level of the
input operands; if the operands have different working precision levels, thehigher precision level is chosen for the result. When assigning a double precision constant or variable to a multiprecision variable or array element, or
when reading multiprecision data from a file, the result is assigned the default
precision unless a lower precision level is specified.

4

Installation, compilation and linking

The two versions of the MPFUN2020 package, together with installation instructions, are available at https://www.davidhbhailey.com/dhbsoftware.
Installation, compilation and linking is relatively straightforward, provided that one has a Unix-based system, such as Linux or Apple OS X,
and a Fortran-2008-compliant compiler, such as GNU gfortran, Intel’s ifort
or NAG’s nagfor. Full details on installing the MPFUN2020 software (andalso for installing gfortran, if needed) are available in the README file forthe distribution, which is in the main directory of the distribution, or athttps://www.davidhbhailey.com/dhbsoftware.

4.1

The two variants

Each of the two versions of the software comes in two variants:• Variant 1: This is recommended for beginning users and for basic applications that do not dynamically change the working precision level(or do so only rarely).
• Variant 2: This is recommended for more sophisticated applications
that dynamically change the working precision level. It does not allow
some mixed-mode combinations, and requires one to explicitly specify
14

a working precision parameter for some functions. However, in thepresent author’s experience, these restrictions result in less overall effort
to produce a debugged, efficient application code. See Section 5.9.
The Fortran source files and scripts required for each of these variants arein the respective directories fortran-var1 and fortran-var2.
Compile/link scripts are available in the fortran-var1 and fortran-var2 directories for the gfortran and Intel ifort compilers, and, with MPFUN20-Fort,for the NAG nagfor compiler. For example, to compile variant 1 of either
the MPFUN20-Fort or MPFUN20-MPFR library using the GNU gfortrancompiler, go to the fortran-var1 directory and type./gnu-complib1.scr

To compile and link the application program tpslq1.f90 for variant 1, using
the GNU gfortran compiler, producing the executable file tpslq1, type./gnu-complink1.scr tpslq1

To execute the program, with output to tpslq1.txt, type./tpslq1 > tpslq1.txt

These scripts assume that the user program is in the same directory as thelibrary files; this can easily be changed by editing the script files.
Several sample test programs, together with reference output files, areincluded in the fortran-var1 and fortran-var2 directories — see Section 6.

5

Coding Fortran applications

A high-level Fortran interface is provided for each version of the package(MPFUN20-Fort and MPFUN20-MPFR).

5.1

Basic instructions for using the package

To use either MPFUN20-Fort or MPFUN20-MPFR in a user program, first
set the parameter mpipl, the default standard precision level in digits, whichis the maximum precision level to be used for subsequent computation; thisis used to specify the amount of storage required for multiprecision data.
mpipl is set in a parameter statement in file mpfunf.f90 in the fortran-var1
or fortran-var2 directory of the software. In the code as distributed, mpipl
is set to 4,000 digits (sufficient to run the six test programs of Section 6),
15

but it can be set to any level greater than 50 digits. mpipl is automaticallyconverted to mantissa words by the formulampwds = int (mpipl / mpdpw + 2)
where mpdpw is a double precision system parameter set in file mpfuna.f90 (inMPFUN20-Fort, mpdpw = log10 260 = 18.0617997 . . .; in MPFUN20-MPFR,
mpdpw = log10 264 = 19.2659197 . . .). The resulting parameter mpwds is theinternal default precision level, in words. All subsequent computations are
performed to mpwds words precision unless the user, in an application code,
specifies a lower precision.
After setting the value of mpipl, if needed, compile the appropriate version of the library, using one of the scripts mentioned in the previous section.
Next, place the following line in every subprogram of the user’s application code that contains a multiprecision variable or array, at the beginningof the declaration section, before any implicit or type statements:
use mpmodule
To designate a variable or array as multiprecision real in an application program, use a Fortran-90 type statement with the type mp real, as in thisexample:type (mp real) a, b(m), c(m,n)
Similarly, to designate a variable or array as multiprecision complex, use atype statement with the type mp complex. Thereafter when one of thesevariables or arrays appears, as in the coded = a + b(i) * sqrt(3.d0 - c(i,j))
the proper underlying multiprecision routines are automatically called.
Most common mixed-mode combinations (arithmetic operations, comparisons and assignments) involving multiprecision real (MPR), multiprecisioncomplex (MPC), double precision (DP), double complex (DC), and integer
operands are supported. A complete list of supported mixed-mode operations is given in Table 2. See Section 5.7 below about DP and DC constantsand expressions.
Input and output of MP variables or array elements are done by using
the special subroutines mpread and mpwrite. See Table 4 and Section 5.6.

5.2

Extra-high precision

The above instructions apply if the default precision level, namely mpipl, isless than or equal to 20,000 digits. If one requires a precision level greater
16

Operatora = b(assignment)

a**b
(ab )

<=, >=, <, >(≤, ≥, <, > tests)

Arg 1
MPR
DP
Int
MPR
MPC
MPC
DP
DC
MPR
MPR
MPR
MPC
MPC
MPR
MPR
MPC
MPC
MPR
MPC
MPR
DP
MPR
Int
MPR

Arg 2
MPR
MPR
MPR
MPC
MPR
MPC
MPC
MPC
DP {1}
Int {1}
Char {1}DP {1}
DC {1}
Int
MPR
Int
MPC
MPC
MPR
MPR
MPR
DP
MPR
Int

Operator+, -, *, /(+, −, ×, ÷)

Arg 1
MPR
DP
MPR
Int
MPR
MPC
DP
MPC
DC
MPC
MPR
MPC

Arg 2
MPR
MPR
DP
MPR
Int
MPC
MPC
DP
MPC
DC
MPC
MPR

==, /=(=, 6= tests)

MPR
DP
MPR
Int
MPR
MPC
DP
MPC
DC
MPC
MPR
MPC

MPR
MPR
DP
MPR
Int
MPC
MPC
DP
MPC
DC
MPC
MPR

Table 2: Supported mixed-mode operator combinations. MPR denotes multiprecision real, MPC denotes multiprecision complex, DP denotes double
precision, DC denotes double complex, Int denotes integer and Char denotes
arbitrary-length character string. Note:{1} These operations are not allowed in variant 2 — see Section 5.9.

17

than this, then in addition to changing the value of mpipl (and, thus, mpwds)
in file mpfunf.f90, one must include a call to mpinit at the start of execution,
in a single-threaded section of code, before any multiprecision operations are
performed:
call mpinit (nwds)
where nwds is the requested precision level in words (normally nwds = mpwds).
A call to mpinit initializes FFT tables (in MPFUN20-Fort) and calculateslog(2), π and γ. If the requested precision (mpwds or nwds) is less than
mpl2pi, a call to mpinit has no effect.
It should be noted, however, that some special functions internally requirelog(2), π and γ to higher precision. If one sees a run-time error message suchas “Pi must be precomputed to precision 1200 words,” insert a call to mpinit
at start of execution with the specified precision as the argument.

5.3

Functions and subroutines

As mentioned above, the MPFUN20-Fort and MPFUN20-MPFR packages
support most Fortran-2008 intrinsics, including all the well-known transcendentals (e.g., sin, exp, log, etc.), and, in addition, an extensive set of special
functions and a set of routines for polynomial solution. Further, the package
includes a set of I/O and conversion functions, such as functions to convert
between double precision and multiprecision real or between double complex and multiprecision complex. These functions and subroutines are listed
below in tables 3, 4, 5, 6 and 7 (some are listed in more than one table).
In these tables, “F” denotes function, “S” denotes subroutine, “MPR”
denotes multiprecision real, “MPC” denotes multiprecision complex, “DP”
denotes double precision, “DC” denotes double complex, “Int” denotes integer and “QP” denotes IEEE quad precision (if supported by the compiler).
The variable names r1,r2,r3 are MPR; z1 is MPC; d1 is DP; dc1 is DC;
q1 is QP; i1,i2,i3,n,nb,np,nq are integers; s1 is character(1); sn is character(n); rb is an MPR vector of length nb; ss is an MPR vector of lengthn; aa is an MPR vector of length np; and bb is an MPR vector of length nq.

5.4

Special functions

As briefly mentioned above, the package includes numerous special functions: the BesselI, BesselJ, BesselK and BesselY functions (each with integer18

or MPR order), the digamma function, the erf and erfc functions, the exponential integral (expint) function, the gamma function, the Hurwitz zetafunction, the hypergeometric pFq function, the incomplete gamma function,
the polygamma function, the polylogarithm function, the StruveH function
and the zeta function.
A complete listing of these functions and their arguments is presented inTables 5 and 6. The algorithms employed in these routines are described in
the Appendix (Section 7.4).

5.5

Polynomial evaluation and solution

As briefly mentioned above, the package includes several routines for polynomial evaluation and solution. Two routines solve an input MPR polynomial,
given an MPR or MPC initial approximation to the root, by means of Newton iterations with a dynamically varying level of precision. A third routine
returns the greatest common divisor (GCD) of two input MPR polynomials. A fourth routine divides two MPR polynomials, producing the quotient
and remainder polynomials. Two other routines evaluate a polynomial at an
MPR or MPC argument.
Polynomials with repeated roots often defeat Newton iteration schemes
for solution. Such polynomials can be handled as follows: first apply thepolynomial GCD routine to the input polynomial and its derivative, andthen divide the input polynomial by the GCD polynomial. This produces
a reduced polynomial that has all the roots of the original polynomial, but
no repeated roots. Then the Newton iteration-based solution routines canquickly produce high-precision roots.
Table 7 lists the polynomial evaluation and solution routines. Usageof these routines is demonstrated by the tpolysolve sample program; see
Section 6.

5.6

Input and output of multiprecision data

As briefly mentioned above in Table 4, the package includes routines forbinary-decimal conversion and input or output of multiprecision data. Theseroutines and their usage are described in Table 8.

19

Type
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPC
F:MPR
F:MPC
F:MPR
F:DP
F:DC
F:MPR
F:MPR
F:MPR
F:MPC
F:MPR
F:MPR
F:MPR
F:MPC
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPC
F:MPR
F:MPR
F:MPC
F:MPR
F:MPR

Name
abs(r1)abs(z1)acos(r1)
acosh(r1)aimag(z1)aint(r1)
anint(r1)asin(r1)
asinh(r1)atan(r1)
atan2(r1,r2)
atanh(r1)bessel j0(r1)bessel j1(r1)bessel jn(n,r1)
bessel y0(r1)bessel y1(r1)bessel yn(n,r1)
conjg(z1)cos(r1)cos(z1)cosh(r1)dble(r1)
dcmplx(z1)erf(r1)erfc(r1)
exp(r1)exp(z1)
gamma(r1)hypot(r1,r2)log(r1)log(z1)log10(r1)max(r1,r2)
min(r1,r2)
mod(r1,r2)
sign(r1,r2)
sin(r1)sin(z1)sinh(r1)sqrt(r1)sqrt(z1)
tan(r1)tanh(r1)

Description
Absolute value
Absolute value of complex arg
Inverse cosine
Inverse hyperbolic cosine
Imaginary part of complex arg
Truncates to integer
Rounds to closest integer
Inverse sine
Inverse hyperbolic sine
Inverse tangent
Arctangent with two args
Inverse hyperbolic tangent
BesselJ function, order 0
BesselJ function, order 1
BeselJ function, order nBesselY function, order 0
BesselY function, order 1
BeselY function, order nComplex conjugate
Cosine of real arg
Cosine of complex arg
Hyperbolic cosine
Converts MPR argument to DP
Converts MPC argument to DC
Error function
Complementary error function
Exponential function of real arg
Exponential function of complex arg
Gamma function
Hypotenuse of two args
Natural logarithm of real arg
Natural logarithm of complex arg
Base-10 logarithm of real arg
Maximum of two (or three) args
Minimum of two (or three) args
Mod function = r1 - r2*aint(r1/r2)
Transfers sign from r2 to r1
Sine function of real arg
Sine function of complex arg
Hyperbolic sine
Square root of real arg
Square root of complex arg
Tangent function
Hyperbolic tangent function

Table 3: Fortran-2008 intrinsic functions extended to multiprecision.
20

Type
F:MPC
F:MPC
F:MPC
F:MPC
S

Name
mpcmplx(r1,r2)
mpcmplx(dc1)
mpcmplx(z1)mpcmplxdc(dc1)
mpcssh(r1,r2,r3)

S

mpcssn(r1,r2,r3)

SS

mpdecmd(r1,d1,i1)
mpeform(r1,i1,i2,s1)

S

mpfform(r1,i1,i2,s1)

F:MPR
S
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR

mpegamma()
mpinit()mplog2()
mpnrt(r1,i1)
mppi()
mpprod(r1,d1)
mpquot(r1,d1)
mprand(r1)S

mpread(i1,r1)

S

mpread(i1,z1)

F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:Int
F:Int
S

mpreal(r1)mpreal(z1)mpreal(d1)mpreal(s1,i1)
mpreal(sn)mpreald(d1)mpwprec(r1)mpwprec(z1)mpwrite(i1,i2,i3,r1)

S

mpwrite(i1,i2,i3,z1)

Description
Converts (r1,r2) to MPC [1]Converts DC arg to MPC [1]Converts MPC arg to MPC [1]Converts DC to MPC, without checking [1][2]Returns both cosh and sinh of r1, in the sametime as calling just cosh or just sinh
Returns both cos and sin of r1, in the sametime as calling just cos or just sin
Converts r1 to the form d1*10^i1
Converts r1 to char(1) string in Ei1.i2
format, suitable for output (Sec. 5.6)
Converts r1 to char(1) string in Fi1.i2
format, suitable for output (Sec. 5.6)
Returns Euler’s γ constant [1]Initializes for extra-high precision (Sec. 5.2) [1]
Returns log(2) [1]
Returns the i1-th root of r1
Returns π [1]Returns r1*d1, without checking [2]Returns r1/d1, without checking [2]Returns pseudorandom number, based on r1
Start with an irrational, say r1 = mplog2()
Typical iterated usage: r1 = mprand(r1)Inputs r1 from Fortran unit i1; up to five
MPR args may be listed (Sec. 5.6) [1]
Inputs z1 from Fortran unit i1; up to five
MPC args may be listed (Sec. 5.6) [1]
Converts MPR arg to MPR [1]Converts MPC arg to MPR [1]Converts DP arg to MPR [1]Converts char(1) string of length i1 to MPR (Sec. 5.6) [1]
Converts char(n) string to MPR (Sec. 5.6) [1]
Converts DP to MPR, without checking [1][2]Returns precision in words assigned to r1
Returns precision in words assigned to z1
Outputs r1 in Ei2.i3 format to unit i1; up to
five MPR args may be listed (Sec. 5.6)
Outputs z1 in Ei2.i3 format to unit i1; up to
five MPC args may be listed (Sec. 5.6)

Table 4: Additional general routines (F: function, S: subroutine). Notes:[1]: In variant 1, an integer precision level argument (mantissa words) may optionally be
added as the final argument; this argument is required in variant 2. See Sec. 5.9.
[2]: These do not check DP or DC values. See Sec. 5.7.

21

Type
F:MPR
F:MPR
S

Name
agm(r1,r2)
airy(r1)
mpberne(nb,rb)

F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR

bessel i(r1,r2)
bessel in(n,r1)
bessel j(r1,r2)
bessel j0(r1)bessel j1(r1)bessel jn(n,r1)
bessel k(r1,r2)
bessel kn(n,r1)
bessel y(r1,r2)
bessel y0(r1)bessel y1(r1)bessel yn(n,r1)
digamma(r1)digamma be(nb,rb,r1)

F:MPR
F:MPR
F:MPR
F:MPR
F:MPR
F:MPR

erf(r1)erfc(r1)
expint(r1)gamma(r1)hurwitz zetan(n,r1)
hurwitz zetan be(nb,rb,n,r1)
hypergeom pfq(np,nq,aa,bb,r1)
incgamma(r1,r2)
polygamma(n,r1)
polygamma be(nb,rb,n,r1)

F:MPR
F:MPR
F:MPR
F:MPR

Description
Arithmetic-geometric mean
Airy function [1]Initialize array rb of length nb with first nb even
Bernoulli numbers [2][3]BesselI function, order r1, of r2
BesselI function, order n, of r1
BesselJ function, order r1, of r2
BesselJ function, order 0, of r1
BesselJ function, order 1, of r1
BesselJ function, order n, of r1
BesselK function, order r1, of r2
BesselK function, order n, of r1
BesselY function, order r1, of r2
BesselY function, order 0, of r1
BesselY function, order 1, of r1
BesselY function, order n, of r1
Digamma function of r1 [1]Digamma function of r1, using nb even
Bernoulli numbers in rb [3]Error function
Complementary error function
Exponential integral function
Gamma function
Hurwitz zeta function, order n ≥ 2, of r1 [4]Hurwitz zeta function, order n ≥ 2, of r1, using
nb even Bernoulli numbers in rb [3]Hypergeometric pFq function of aa, bb and r1;
dimensions are aa(np) and bb(nq) [5]
Incomplete gamma function [6]Polygamma function, order n ≥ 1, of r1 [4]Polygamma function, order n ≥ 1, of r1, using
nb even Bernoulli numbers in rb [3]Table 5: Special functions (F: function, S: subroutine) (continued in Table 6). Notes:[1]: Only available with MPFUN20-MPFR.
[2]: In variant 1, an integer precision level argument (mantissa words) may optionally be
added as the final argument; this argument is required in variant 2. See Sec. 5.9.
[3]: For most applications, set nb > 2X precision in decimal digits; see mpberne above.
[4]: For hurwitz zetan and polygamma, the argument r1 is limited to the range (0, 1].
[5]: For hypergeom pfq, the integers np and nq must not exceed 10.
[6]: For incgamma, r1 must not be zero, and must not be negative unless r2 = 0.

22

Spolylog ini(n,ss)

F:MPR

polylog neg(n,ss,r1)

F:MPR
F:MPR
F:MPR
F:MPR

polylog pos(n,r1)
struve hn(n,r1)zeta(r1)
zeta be(nb,rb,r1)

F:MPR

zeta int(n)Initialize array ss, of size |n|, for computing
polylogarithms of order n when n < 0 [2][6]Polylogarithm function of r1, for n < 0, using
precomputed data in ss [7]Polylogarithm function, order n ≥ 0, of r1 [8]StruveH function, order n ≥ 0, of r1 [9]Zeta function of r1
Zeta function of r1, using nb even Bernoulli
numbers in rb [3]Zeta function of integer argument n [2]Table 6: Special functions, continued (F: function, S: subroutine). Notes:[2]: In variant 1, an integer precision level argument (mantissa words) may optionally be
added as the final argument; this argument is required in variant 2. See Sec. 5.9.
[6]: For incgamma, r1 must not be zero, and must not be negative unless r2 = 0.
[7]: For polylog ini and polylog neg, the integer n is limited to the range [−1000, −1].
[8]: For polylog pos, the argument r1 is limited to the range (−1, 1).
[9]: For struve hn, the argument r1 is limited to the range [−1000, 1000].

5.7

Handling double precision values

Double precision constants and expressions are indispensable in high-precisionapplications. For one thing, the product, say, of a multiprecision value times
a double precision value is more rapidly computed with a routine dedicatedto this task than converting the double precision value to multiprecision andthen calling the full multiplication routine. Certainly the usage of double
precision constants such as modest-sized whole numbers and exact binary
fractions (e.g., 3., 12345., 2.5, 6.125), which are entirely safe in a multiprecision application, should be allowed.
Users should be aware, however, that there are some hazards in this typeof programming, inherent in conventions adopted by all Fortran compilers
and most other languages as well. For example, the code
r1 = 3.14159d0
where r1 is MPR, does NOT produce the true multiprecision equivalentof 3.14159. In fact, the above line of code will result in a run-time error
with either version of the MPFUN2020 software. To obtain the full MPR
converted value, write this as
r1 = ’3.14159’
or, if using variant 2, as
r1 = mpreal (’3.14159’, nwds)
23

S

S

S

S

S

S

mprpolysolve (nq, polyq, rootapp, root)
Given an input MPR polynomial with coefficients polyq(0:nq) (degree nq),
and an input MPR approximation in rootapp, this returns the full precision
MPR root in root; if no root is found, MPR zero is returned in root.
mpcpolysolve (nq, polyq, rootapp, root)
Given an input MPR polynomial with coefficients polyq(0:nq) (degree nq),
and an input MPC approximation in rootapp, this returns the full precision
MPC root in root; if no root is found, MPC zero is returned in root.
mprpolygcd (n1, poly1, n2, poly2, ng, polyg)
Given two input MPR polynomials poly1(0:n1) (degree n1) and poly2(0:n2)(degree n2), with n1 >= n2, this returns the MPR greatest common divisor(GCD) polynomial in polyg(0:ng) (degree ng). polyg must be dimensioned(0:n1) in calling program.
mprpolydiv (n1, poly1, n2, poly2, n3, poly3, n4, poly4)
Given two input MPR polynomials poly1(0:n1) (degree n1) and poly2(0:n2)(degree n2), with n1 >= n2, this returns the MPR polynomial quotient inpoly3(0:n3) (degree n3) and the MPR remainder in poly4(0:n4) (degree n4).
poly3 and poly4 must be dimensioned (0:n1) in calling program.
mprpolyeval (np, poly, arg, result)
Given an input MPR polynomial with coefficients poly(0:np) (degree np),
and an input MPR argument in arg, this returns the MPR result of polynomial
evaluation in result.
mpcpolyeval (np, poly, arg, result)
Given an input MPR polynomial with coefficients poly(0:np) (degree np),
and an input MPC argument in arg, this returns the MPC result of polynomial
evaluation in result.

Table 7: Polynomial evaluation and solution routines

24

S

S

S

S

F:MPR

F:MPR

S

S

Note

mpeform (r1,i1,i2,s1)
This converts the MPR number r1 into decimal character form in thecharacter(1) array s1. The argument i1 (input) is the lengthof the output string, and i2 (input) is the number of digits after the decimal
point. The format is analogous to Fortran E format. The result is left-justified
among the i1 cells of s1. The condition i1 >= i2 + 20 must hold.
mpfform (r1,i1,i2,s1)
This converts the MPR number r1 into decimal character form in thecharacter(1) array s1. The argument i1 (input) is the length ofthe output string, and i2 (input) is the number of digits after the decimal
point. The format is analogous to Fortran F format. The result is right-justified
among the i1 cells of s1. The condition i1 >= i2 + 20 must hold.
mpread (i1,r1)
This reads the MPR number r1, presumed in decimal format, from Fortran
logical unit i1. The digits of r1 may span more than one line, provided
that a backslash appears at the end of a line to be continued. Individual input
lines may not exceed 2048 characters in length. Format: The input number musthave a leading digit (possibly zero), must have a period somewhere; may include
an e, d, E or D followed by an integer exponent; but must not have embedded
blanks. Up to five MPR arguments may be included in argument list. See theNote below on an additional precision argument.
mpread (i1,z1)This is the same as the previous item, except that the input argument z1
is of type MPC (a pair of MPR). Up to five MPC arguments may be included inargument list. See the Note below on an additional precision argument.
mpreal (s1,i1)
This converts the string s1, which is of type character(1) and length i1,
to MPR. See mpread above for format. See the Note below on an additional
precision argument.
mpreal (sn)This converts the string sn, which may be of type character(n) for any
n, to MPR. See mpread above for format. On some systems, n may belimited, say to 2048; if this is a problem, use the previous item. See the Notebelow on an additional precision argument.
mpwrite (i1,i2,i3,r1)
This writes the MPR number r1 to Fortran logical unit i1, in a format
analogous to Fortran E format, left-justified in the field. The argument i2(input) is the length of the output field, and i3 (input) is the number ofdigits after the decimal point; i3 must exceed i2 by 20. Up to five
MPR arguments may be included in argument list.
mpwrite (i1,i2,i3,z1)This is the same as the previous item, except that the argument z1 is oftype MPC (a pair of MPR). Up to five MPC arguments may be included inargument list.
For the mpread and mpreal routines above, when using variant 1, an integer
working precision level argument (mantissa words) may optionally be added
as the final argument; this argument is required in variant 2. See Section 5.9.

25
Table 8: Input/output and conversion routines

where nwds is the working precision (in words) to be assigned to r1. Similarly,
the coder2 = r1 + 3.d0 * sqrt (2.d0)
where r1 and r2 are MPR, does NOT produce the true multiprecision value
one might expect, since the expression 3.d0 * sqrt (2.d0) will be performed in double precision, according to Fortran precedence rules. In fact,
the above line of code will result in a run-time error with either version of
the MPFUN2020 software. To obtain the fully accurate result, write this as
r2 = r1 + 3.d0 * sqrt (mpreal (2.d0))
or, if using variant 2, as
r2 = r1 + 3.d0 * sqrt (mpreal (2.d0, nwds))
where nwds is the working precision in words.
To help avoid such problems, both versions of the MPFUN2020 software
check every double precision value (constants, variables and expression values) in a multiprecision statement at execution time to see if it has morethan 40 significant bits. If so, it is flagged as an error, since very likely suchusage represents an unintended loss of precision in the application program.
This feature catches 99.99% of accuracy loss problems due to the usage of
inexact double precision values.
On the other hand, some applications (including several of the sample testcodes mentioned in Section 6) contain legitimate double precision constants
that are trapped by this test. In order to permit such usage, four special
functions have been provided: mpprod, mpquot, mpreald, mpcmplxdc (seeTable 4). The first and second return the product and quotient, respectively,of a MPR argument and a DP argument; the third converts a DP value
to MPR (with an optional precision level parameter — see Section 5.9);
and the fourth converts a DC value to MPC (with an optional precisionlevel parameter — see Section 5.9). These routines do not check the double
precision argument to see if it has more than 40 significant bits.

5.8

Support for quad precision

Several processor/compiler systems now support IEEE quad (128-bit) floatingpoint arithmetic. If one’s platform does support IEEE quad, limited support
is provided in the MPFUN2020 package. To determine whether or not IEEE
quad is supported on one’s platform, compile and execute the simple program

26

Type
F:MPR
F:MPR
F:MPR
F:MPR
F:Q

Name
mpprod(r1,q1)
mpquot(r1,q1)
mpreal(q1)mprealq(d1)qreal(r1)

Description
Returns r1*q1, without checking [2]Returns r1/q1, without checking [2]Converts quad real to MPR [1]Converts quad to MPR, without checking [1][2]Converts MPR to quad.

Table 9: Support for IEEE quad precision (if available) (F: function, S: subroutine).
Notes:[1]: In variant 1, an integer precision level argument (mantissa words) may optionally be
added as the final argument; this argument is required in variant 2. See Sec. 5.9.
[2]: These do not check quad values. See Sec. 5.8.

program testq
integer, parameter:: mprknd2 = selected real kind (33, 4931)
real (mprknd2) q1, q2
q1 = real (-1.d0, kind = mprknd2)
q2 = acos (q1)write (6, ’(1p,d42.34)’) q2
stop
end

If this code is successfully compiled and produces the output
3.1415926535897932384626433832795028D+00(except possibly for the final digit), then one’s system supports IEEE quad.
The MPFUN2020 function mpreal function accepts a quad argument for
conversion to multiprecision, and the qreal function accepts a multiprecision argument for conversion to quad. Other routines available are mpprod,
mpquot and mprealq. Note that the mpreal function checks the quad argument to see if it has more than 90 significant bits; if so, an error is flagged,
as mentioned above for double precision conversions. To convert without
checking, use mprealq. See Table 9.
When using variant 1, an integer working precision level argument (mantissa words) may optionally be added as the final argument for the functions
mpreal and mprealq; this argument is required in variant 2. See Section 5.9.

5.9

Dynamically changing the working precision

Some multiprecision applications run fine with a fixed precision level, butothers are more efficiently implemented with a working precision level that27

is changed frequently. Accordingly, for each version of the package, there are
two variants of the package, as mentioned above in Section 4:• Variant 1: This is recommended for beginning users and for basic applications that do not dynamically change the working precision level(or do so only rarely).
• Variant 2: This is recommended for more sophisticated applications
that dynamically change the working precision level. It does not allow
some mixed-mode combinations, and requires one to explicitly specify
a working precision parameter for some functions. However, in thepresent author’s experience, these restrictions result in less overall effort
to produce a debugged, efficient application code.
In particular, with variant 1:
1. Assignments of the form R = X, where R is MPR and X is DP, integer
or literal, are permitted. Assignments of the form Z = Y, where Z is
MPC and Y is DP or DC, are permitted.
2. The routines mpcmplx, mpcmplxdc, mpegamma, mpinit, mplog2,
mppi, mpread, mpreal, mpreald, mprealq and zeta int each havean (optional) integer argument as the final argument in the list. This
argument is the working precision level, in words, to be assigned to theresult(s). If this argument is not present, the default precision level(mpwds words, corresponding to mpipl digits) is assumed.
In contrast, with variant 2:
1. The assignments mentioned in item 1 above are not permitted. If anyof these appears in code, compile-time errors will result. Instead, onemust use mpreal and mpcmplx, as appropriate, with the precision level(in mantissa words) as the final argument, to perform these conversions.
2. The optional working precision level arguments mentioned in item 2
above are required in all cases. For example, if the full default precisionlevel (mpwds, corresponding to mpipl digits) is required, then a call to
one of the functions or subroutines in item 2 above must have mpwds
as the final argument.

28

Note that the mpreal function, with the precision level (in words) as thesecond argument, can be used to assign an MPR argument with one precisionlevel to an MPR variable or array element with a different working precisionlevel. The same is true of mpcmplx. The working precision currently assigned
to any MP variable or array element may be obtained by using the function
mpwprec — see Table 4.
Along this line, when one uses the precision level arguments, a precisionlevel of ndig digits can be converted to words by the formula nwds = int(ndig / mpdpw + 2). By using the global built-in variable mpdpw (which isdifferent between MPFUN20-Fort and MPFUN20-MPFR) in this way, theuser code remains completely portable between the two versions.
As it turns out, in most applications, even those that frequently require
the working precision to be changed, only a few changes need to be made to
the source code to implement variable precision. Consider, for example, thefollowing user code, where the default precision has beenn set to 2,500 digits:
integer k, nx
parameter (nx = 128)
type (mp real) x(nx)x(1) = 1.d0
do k = 2, nx
x(k) = 2.d0 * x(k-1) + 1.d0
enddoThis code, as written, is permissible with variant 1, but not with variant
2, because the assignment x(k) = 1.d0 is not allowed. Furthermore, alloperations are performed with the default precision level of 2,500 digits. So
with variant 2, where one wishes to perform this loop with a precision levelof, say, only 500 digits, this should be written as:
integer k, ndig, nwds, nx
parameter (nx = 128, ndig = 500, nwds = int (ndig / mpdpw + 2))
type (mp real) x(nx)x(1) = mpreal (1.d0, nwds)
do k = 2, nx
x(k) = 2.d0 * x(k-1) + 1.d0
enddo

Note that by changing x(1) = 1.d0 to x(1) = mpreal (1.d0, nwds),
29

the array element x(1) is assigned the value 1.0, with a working precision ofnwds words (i.e., 500 digits). In the loop, when k is 2, x(2) also inherits theworking precision level nwds words, since it is computed with an expressionthat involves x(1). By induction, all elements of the array x inherit theworking precision level nwds words (i.e., 500 digits).
This scenario is entirely typical of using variable precision: in most cases,
it is only necessary to make a few code changes, such as in assignments todouble precision values before a loop, to completely control dynamic precision. However, it is highly recommended that the user frequently employ thesystem function mpwprec, which returns the working precision level currently
assigned to an multiprecision variable or array element (see Table 4), to ensure that the working precision level the user thinks is assigned to a variable
is indeed the level being used by the program.
Using variant 2, with its stricter coding standards, requires a bit moreprogramming effort, but in the present author’s experience, when dealing
with applications that dynamically change the precision level, this additional
effort is more than repaid by fewer debugging and performance problems in
actual usage. A code written for variant 2 also works with variant 1, but notvice versa.
See the sample test codes mentioned in Section 6, all of which are written
to conform to the stricter standards of variant 2.

5.10

Medium precision datatype

In many high-precision applications, only part of the multiprecision variables
and arrays contain full precision data; others contain data with only modest
precision, typically only a fraction as high (although still higher than double
precision). Since all multiprecision data are allocated sufficient space toaccommodate full precision values, much of the storage and data movement
costs for modest precision data are wasted.
To reduce memory usage and improve performance in such applications,a medium precision real and a medium precision complex datatype have been
defined, and are available in both the MPFUN20-Fort and MPFUN20-MPFR
versions of the software. To use these datatypes, do the following:First, set the medium precision level mpiplm in file mpfunf.f90; by default
it is set to 400 digits. Then in each user subprogram that will include either
full precision or medium precision data, insert the following line:
use mpmodule30

as with standard full precision. To designate a variable or array as medium
precision real, use a Fortran-90 type statement with the type mp realm, as
in this example:
type (mp realm) a, b(m), c(m,n)
Similarly, use the type mp complexm for medium precision complex data.
Direct assignments between full precision and medium precision variables, as well as mixed-mode arithmetic operations between full precisionand medium precision variables, are not allowed. If one needs to converta regular full precision value to a medium precision value, use the function
mprealm, as in this example
type (mp real) a; type (mp realm) b
b = mprealm (a)Similarly, to convert a medium precision value to a full precision value, use
the function mpreal, as in
type (mp real) a; type (mp realm) b
a = mpreal (b)Note however, that with either mpreal or mprealm, in variant 1, an integer
working precision level argument (mantissa words) may optionally be added
as the final argument; this argument is required in variant 2. See Section 5.9
for details.
Each of the intrinsic functions and subroutines listed in tables 2, 3, 4, 5, 6,
7, 8 and 9 is also supported with medium precision arguments, except for thefollowing: mpreald, mprealq, mpcmplx, mpcmplxd, mpegamma, mplog2,
mppi and zeta int. The medium precision equivalents of these functions are:mprealdm, mprealqm, mpcmplxm, mpcmplxdm, mpegammam, mplog2m,
mppim and zeta intm, respectively. See Table 10. Note that the comment
regarding the working precision argument in the previous paragraph applies
to each of these functions.
Usage of the medium precision datatype is illustrated in the programs
tpslqm3.f90 and tpphix3.f90 in the application programs in Section 6.

6

Sample applications and performance

The current release of the software includes a set of sample application programs in the fortran-var1 and fortran-var2 directories:
1. testmpfun.f90: An overall test suite: This briefly tests most individual MPFUN2020 operations and functions (including mixed mode
31

Type
F:MPRM
F:MPRM
F:MPRM
F:MPCM
F:MPCM
F:MPCM
F:MPRM
F:MPRM
F:MPRM
F:MPRMName
mprealm(r1)mprealm(d1)mprealdm(d1)mprealqm(q1)mpcmplxm(dc1)
mpcmplxdm(dc1)mpegammam()
mplog2m()
mppim()
zeta intm(n)Description
Converts MPR to med prec real [1]Converts DP to med prec real [1]Converts DP to med prec real, without checking [1][2]Converts quad to med prec real , without checking [1][2]Converts DC to med prec complex [1]Converts DC to med prec complex, without checking [1][2]Returns Euler’s gamma constant to med prec real [1]Returns log(2) to med prec real [1]Returns π to med prec real [1]Returns ζ(n) to med precision real [1]Table 10: Routines to support medium precision datatytpe (F: function, S: subroutine),
in addition to the functions listed in Tables 2, 3, 4, 5 and 9. Notes:[1]: In variant 1, an integer precision level argument (mantissa words) may optionally be
added as the final argument; this argument is required in variant 2. See Sec. 5.9.
[2]: These do not check double or quad values. See Sec. 5.7 and Sec. 5.8.
arithmetic, comparison operations, transcendental functions, special
functions and polynomial routines), by comparing each result with reference results in the file testmpfun.ref.txt, which must be present in thesame directory. This is not an exhaustive test of all possible scenarios,
but it often detects bugs and compiler issues. Size of code: 1,093 lines.
Precision: 500 digits.
2. tpolysolve.f90: A polynomial solution test program: This computes
each of the 12 real and 26 complex roots of an input degree-40 realpolynomial with large integer coefficients and two repeated roots. The
computation includes a polynomial greatest common divisor (GCD)
calculation to produce a reduced polynomial with no repeated roots;
this reduced polynomial is then solved using a Newton iteration schemewith dynamically varying precision. This code employs the polynomial
evaluation and solution routines described above in Table 7. Size ofcode: 274 lines. Precision level: 4,000 digits.
3. tpslq1.f90: A one-level standard PSLQ program: This finds the coefficients of the degree-30 polynomial satisfied by 31/5 − 21/6 . Size ofcode: 866 lines. Precision level: 250 digits.
4. tpslqm1.f90: A one-level multipair PSLQ program: This finds thecoefficients of the degree-30 polynomial satisfied by 31/5 − 21/6 . Size ofcode: 946 lines. Precision level: 250 digits.
32

5. tpslqm2.f90: A two-level multipair PSLQ program: This finds thecoefficients of the degree-56 polynomial satisfied by 31/7 − 21/8 . Sizeof code: 1,796 lines. Precision level: 750 digits; switches frequently
between multiprecision and double precision.
6. tpslqm3.f90: A three-level multipair PSLQ program: This finds thecoefficients of the degree-72 polynomial satisfied by 31/8 − 21/9 . Sizeof code: 2,365 lines. Precision level: 1200 digits; switches frequently
between full precision (1200 digits), medium precision (varies from 50
to 110 digits) and double precision.
7. tpphix3.f90: A Poisson phi program: This computes the value ofφ2 (x, y) and then employs a three-level multipair PSLQ algorithm tofind the minimal polynomial of degree m satisfied by exp(8πφ2 (1/k, 1/k))
for a given k (see Section 1.4). In the code as distributed, k = 28, m =96, and a palindromic option is employed so that the multipair PSLQ
routines (which are part of this application) search for a relation of size
49 instead of 97. This computation involves transcendental functions
and both real and complex multiprecision arithmetic. Size of code:
2,728 lines. Precision level: 2500 digits; switches frequently between
full precision (2500 digits), medium precision (varies from 50 to 250
digits) and double precision.
8. tquad.f90: A quadrature (numerical integration) program: This performs the tanh-sinh, the exp-sinh or the sinh-sinh quadrature algorithm, as appropriate, on a suite of 18 problems involving numerous
transcendental function references, producing results correct to 500digit accuracy. Size of code: 1,565 lines. Precision level: 1000 digits,
but most computation is done to 500 digits; switches frequently between 500 and 1000 digits.
9. tquadgs.f90: A quadrature (numerical integration) program: This
performs the Gaussian quadrature algorithm on many of the same 18
problems as in tquad, producing results correct to 500-digit accuracy.
This code runs much longer than tquad, due to the much greater expense of computing weights and abscissas. Once the weights and abscissas are computed, they are written to a file, so they can be reusedin future runs. Size of code: 759 lines. Precision level: 1000 digits, but

33most computation is done to 500 digits; switches frequently between
500 and 1000 digits.
In addition, the fortran-var1 and fortran-var2 directories include test
scripts that compile the library and run each of the above sample programs
above (except tquadgs.f90): In directory fortran-var1, these scripts are:• gnu-mpfun-tests1.scr• intel-mpfun-tests1.scr• nag-mpfun-tests1.scr (MPFUN20-Fort only)
while in directory fortran-var2 they are:• gnu-mpfun-tests2.scr• intel-mpfun-tests2.scr• nag-mpfun-tests2.scr (MPFUN20-Fort only)

For each of the eight test programs, the script outputs either TEST PASSEDor TEST FAILED. If all tests pass, then one can be fairly confident that the
MPFUN2020 software and underlying compilers are working properly.
These programs are provided, in part, as examples of programming techniques when using the MPFUN2020 package. Users may feel free to adapt
these codes, although the present author asks to be notified and credited
when this is done. All application programs and library codes are publicly
available but are subject to copyright and other legal conditions. For details,
see the file disclaimer.txt in the distribution package.

6.1

Timings

Table 11 presents some performance timings comparing the two versions of
the package for the first seven test programs listed above, plus an additional
run using the tpphix3.f90 code, with different parameters k and m, andwithout the palindromic option, which is not available when k is odd. Allruns except the last were performed with default precision level mpipl set to
4,000 digits, even though most of these runs required less.
These runs were performed using the GNU gcc/gfortran compiler (version
13.2.0). For the MPFUN20-MPFR version, the GNU gcc compiler (version
13.2.0) was used to build the GMP and MPFR libraries. These runs were34

Code nametestmpfun
tpolysolve
tpslq1
tpslqm1
tpslqm2
tpslqm3
tquad
tpphix3
tpphix3

Precision(digits)
500
4000
250
250
750
1200
1000
2500
5200

k

m

28
25

96
100

MPFUN20Fort
11.57
1.51
4.38
2.31
3.20
11.30
18.68
10.51
288.99MPFUNMPFR
3.16
0.28
2.52
2.12
3.08
11.09
3.56
8.13
219.29

Table 11: Timings on a suite of test programs (seconds).

executed on an Apple Mac Studio system, with OS version 14.4.1, an AppleM2 Max processor and 32 Gbyte of main memory. For uniformity, the timings
are listed in the Table 11 to two decimal place accuracy, but as with allcomputer run time measurements, they should not be considered repeatablebeyond about two or three significant digits.

35

7

Appendix: Numerical algorithms

This section presents a brief overview of the algorithms used in MPFUN20Fort, and, in some cases, also in MPFUN20-MPFR; other algorithms usedin MPFUN20-MPFR are described in [27].

7.1

Algorithms for basic arithmetic

Multiplication. [MPFUN20-Fort.] For modest levels of precision, MPFUN20Fort employs adaptations of the usual schemes we all learned in grade school,
where the number base is 260 = 1152921504606846976. Note that if two nword arguments are multiplied, and the working precision is also n words,
then since only an n-word result is returned, only slightly more than half of
the “multiplication pyramid” need be calculated.
Division. [MPFUN20-Fort.] Division of two extra-high-precision argumentsa and b can be performed by the following scheme. First, note that this
Newton-Raphson algorithm iteration converges to 1/b [17, pg. 226]:
xk+1 = xk + (1 − xk · b) · xk ,(5)

where the multiplication () · xk is performed with only half of the normal
level of precision. These iterations are performed with a working precisionlevel that is approximately doubles with each iteration, except that at three
iterations before the final iteration, the iteration is repeated without doubling
the precision, in order to enhance accuracy. The final iteration is performed
as follows (due to A. Karp):
a/b ≈ (a · xn ) + [a − (a · xn ) · b] · xn ,(6)

where the multiplications a · xn and [] · xn are performed with only half of the
final level of precision. The total cost of this procedure is only about threetimes the cost of a single full-precision multiplication.
Square roots. [MPFUN20-Fort.] Square roots are calculatedby the following√
Newton-Raphson iteration, which converges to 1/ a [17, pg. 227]:
xk+1 = xk + 1/2 · (1 − x2k · a) · xk ,(7)

where the multiplication () · xk is performed with only half of the normal
level of precision. These iterations are performed with a working precision36

level that approximately doubles with each iteration, except that at three
iterations before the final iteration, the iteration is repeated without doubling
the precision, in order to enhance accuracy. The final iteration is performed
as follows (due to A. Karp):√a ≈ (a · xn ) + 1/2 · [a − (a · xn )2 ] · xn ,(8)where the multiplications (a · xn ) and [] · xn are performed with only half the
final level of precision. If this is done properly, the total cost of the calculationis only about three times the cost of a single full-precision multiplication.
n-th roots. [MPFUN20-Fort.] A similar scheme is used to compute n-th roots
for any integer n. Computing xnk , which is required here, can be efficiently
performed using the binary algorithm for exponentiation. This is merely
the observation that exponentiations can be accelerated based on the binary
expansion of the exponent: for example, 317 can be computed as ((((3)2 )2 )2 )2 ·3 = 129140163.
Pseudorandom number generation. [MPFUN20-Fort and MPFUN-MPFR.]
The pseudorandom number included in the package versions is the following:
xn+1 = frac (5501758857861179 · xn ),(9)where the calculation is performed with slightly more than the standard precision for full accuracy. Here the large prime 5501758857861179 is specified
as the double precision constant mprandx, set in module MPFUNA in file
mpfuna.f90. Its value may be changed if desired, but must not exceed 253 .
To use the generator, firstpset a MPR variable r1 to some irrational value
between 0 and 1, such as 1/2, log(2) or π/4. Then successive iterates canbe generated, say, by typing
r1 = mprand (r1)
The same scheme is used for both MPFUN20-Fort and MPFUN20-MPFR.
However, the specific sequence is not the same between the two versions, nor
is it same if the precision level is changed. Note that the above algorithms
are trivially thread-safe, since no auxiliary data or saved context is involved.

7.2

FFT-based multiplication[MPFUN20-Fort only.] Although the multiplication algorithm describedabove in Section 7.1 is very efficient, for higher levels of precision (above37

approximately 20,000 digits, based on the present author’s implementation),significantly faster performance can be achieved by employing a convolution
scheme based on fast Fourier transforms [20][17, pg. 223–224].
Suppose one wishes to multiply two n-precision values whose mantissa
words are given by a = (a0 , a1 , a2 , · · · , an−1 ) and b = (b0 , b1 , b2 , · · · , bn−1 ). It
is easy to see that the desired result, except for releasing carries, is an acyclic
convolution. In particular, assume that a and b are extended to 2n words
each by padding with zeroes. Then the product c = (ck ) is given by
ck =

2n−1
X

aj bk−j ,0 ≤ k < 2n,

(10)

j=0

where bk−j is read as bk−j+2n when k − j is negative. This convolution canbe calculated as(c) = F −1 [F (a) · F (b)],(11)

where F (a) and F (b) denote a real-to-complex discrete Fourier transform(computed using an FFT algorithm), the dot means element-by-element complex multiplication, and F −1 [] means an inverse complex-to-real FFT. The ck
results from this process are floating-point numbers. Rounding these values
to the nearest integer, and then releasing carries beginning at c2n−1 gives thedesired multiplication result.
One important detail was omitted above: to avoid loss of numerical significance, the 60-bit mantissa words of the input multiprecision values are firstdivided into 15-bit sections. Then the FFT-convolution procedure is performed, and the result is suitably reconstituted to an output multiprecision
value of 60-bit words at the end.
In contrast to the basic arithmetic algorithms, FFT-based multiplication
requires precomputed FFT root-of-unity data. Thus, if one requires a precision level greater than 20,000 digits, one must call mpinit at the start of theuser’s program, in a single-threaded initialization section, before any parallel
execution — see Table 4 and Section 5.2 for details.

7.3

Transcendental constants and functions

Most arbitrary precision packages require a significant “context” of global
data to support transcendental function evaluation at a particular precision38

level, and this data is often problematic for thread safety. Further, if thiscontext data must be created and freed within each running thread, thislimits the efficiency in a multithreaded environment. With this in mind,
the transcendental functions in MPFUN2020 were designed to require only aminimum of global data, namely the three constants log(2), π and γ, which
are provided in compile-time global data statements unless extremely highprecision is required. If a precision over 20,000 digits is required, mpinit
must be called at the start of execution, as noted above in Section 5.2, tocompute these three constants.
π is computed in MPFUN20-Fort using a quadratically convergent scheme
discovered by Richard
Brent and Eugene Salamin in 1976 [8, 109–110]. Set√
a0 = 1, b0 = 1/ 2 and s0 = 1/2. Iterate, beginning with k = 0,p
ak + b k, bk+1 = ak bk ,2
ck+1 = a2k+1 − b2k+1 , sk+1 = sk − 2k+1 ck+1 ,2a2
pk+1 = k+1 .
sk+1

ak+1 =(12)

Then pk converges quadratically to π: each iteration approximately doubles
the number of correct digits.
log(2) is computed in MPFUN20-Fort using the following quadraticallyconvergent algorithm, also due to Salamin and Brent [8, 228]: Select N > 2B ,where B is the number of bits of desired precision. Thenlog(2) =π,
2 AGM (1, 4/2N )(13)

where AGM (a, b) denotes the arithmetic-geometric mean: Set a0 = a andb0 = b, then iterate
ak+1 = (ak + bk )/2
p
bk+1 = ak · bk

(14)

until convergence, i.e., until ak = bk to available precision.
Euler’s constant γ = 0.57721566490153286 . . . is calculated in MPFUN20Fort using the following formula, which is an improvement of a techniquepreviously used by Sweeney [37]. If a result accurate to at least B bits is39

desired, first select the integer N = dlog2 (B log 2)e. Then∞m
2N X 2mN X 1− N log 2.
γ ≈ 2Ne m=0 (m + 1)! t=0 t + 1

(15)

N

The error in this approximation is less than 1/(2N e2 ).
The binary values of log(2), π and γ are stored to 20,000 digit precision indata statements in module MPFUNA, in both versions. If higher precision isrequired, mpinit must be called at the start of execution in a single-threaded
section of code, as explained in Section 5.2.
Exponential and logarithm. [MPFUN20-Fort only; MPFUN20-MPFR uses
the MPFR routines.] The exponential function routine first reduces the input
argument to the interval (− log(2)/2, log(2)/2]. Then it divides this value by
2q , producing a very small value, which is then input to the Taylor series forexp(x):x2 x3++ ··· ,
(16)
2!
3!
where the terms are computed with a precision level that decreases as theterms decrease in size, thus saving approximately one-half of the total run
time. When complete, the result is squared q times, and then corrected forthe initial reduction. In the current implementation, q is set to the nearest
integer to B 2/5 , where B is the number of bits of precision.
Since the Taylor series for the logarithm function converges much more
slowly than that of the exponential function, the Taylor series is not used for
logarithms unless the argument is extremely close to one. Instead, logarithms
are computed based on the exponential function, by employing the following
Newton iteration with a level of precision that approximately doubles witheach iteration:exk − a.
(17)
xk+1 = xk −e xkexp(x) = x +Trigonometric functions. [MPFUN20-Fort only; MPFUN20-MPFR uses the
MPFR routines.] The sine/cosine routine first reduces the input argument towithin the interval (−π, π]. This value is then divided by 2q and then input
to the Taylor series for sin(x):x3 x5
sin(x) = x −++ ··· ,
3!
5!
40

(18)

where, as with exp(x), the terms are computed with a precision level that
decreases as the terms decrease in size. Then the double-angle formulascos(2x) = 1 − 2 sin2 (x)cos(2x) = 2 cos2 (x) − 1,

(19)
(20)

are applied q times (formula (19) is used once, and (20) thereafter).
In the√current implementation, q is set to the greatest integer in 2B, where B
is the precision in bits, unless the reduced argument is very pclose to zero,in which case q = 0. When complete, sin(x) is computed as 1 − cos2 (x),withsign, except for the case q = 0, when cos(x) is computed asp corrected
1 − sin2 (x).
The inverse cos/sin function is based on the sine routine, by employing a
Newton iteration with a level of numeric precision that roughly doubles witheach iteration.
Power function. [MPFUN20-Fort only; MPFUN20-MPFR uses the MPFR
routine.] The power function, namely ab for real a > 0 and b, can be computed as eb log a . To further accelerate this operation, the software first examines the value of b to see if it is a rational number with numerator anddenominator up to 107 size, using the extended Euclidean algorithm performed in double precision. If it is, ab is performed using a combination of
the binary algorithm for exponentiation for the numerator, and the n-th root
function for the denominator.

7.4

Special functions

Modern mathematical and scientific computing frequently often involves moresophisticated functions, collectively termed “special functions” [23]. Many ofthese functions have been implemented in the MPFUN20-Fort and MPFUN20MPFR packages, and others will be added as they are developed. Here is abrief description of the functions that have been implemented and the algorithms employed. In each case, care is taken to preserve thread safety, andto avoid, as far as possible, any need to precalculate auxiliary data.
Even Bernoulli coefficients. [MPFUN20-Fort and MPFUN20-MPFR.] The
subroutine mpberne returns an array containing Bernoulli numbers B2k for
k = 1 to n. It employs an advanced algorithm, due to the present author[10]: First, compute {ζ(2k), 0 ≤ k ≤ n} to P words precision based on the41

following well-known formulas:coth(πx) = cosh(πx)/ sinh(πx)
1 1 + (πx)2 /2! + (πx)4 /4! + · · ··=πx 1 + (πx)2 /3! + (πx)4 /5! + · · ·
1=· 1 + (πx)2 /3 − (πx)4 /45 + 2(πx)6 /945 − · · ·
πx
∞
2 X
=
(−1)k+1 ζ(2k)x2k .
πx k=1(21)

The strategy is to calculate the coefficients of the series by polynomial operations. Polynomial division is performed by computing the reciprocal ofthe denominator polynomial, by a polynomial Newton iteration, as follows.
Let N (x) be the polynomial approximation to the numerator series; let D(x)be a polynomial approximation to the numerator numerator series; and let
Qk (x) be polynomial approximations to R(x) = 1/D(x). Then iterate:
Qk+1 = Qk (x) + [1 − D(x)Qk (x)]Qk (x).

(22)

In these iterations, both the degree of the polynomial Qk (x) and the precisionlevel in words are initially set to 4. When convergence is achieved at this
precision level, the degree is doubled, and iterations are continued, etc., until
the final desired degree is achieved. Then the precision level is doubled anditerations are performed in a similar way, until the final desired precision level
P is achieved. The reciprocal polynomial R(x) produced by this process isthen multiplied by the numerator polynomial N (x) to yield an approximation
to the quotient series. The even zeta values are then the coefficients of thisseries, scaled according to the formula (21)above.
Once the even integer zeta values have been computed, the even Bernoulli
coefficients are computed via this well-known formula (for n > 0):B2n =(−1)n−1 2(2n)!ζ(2n).
(2π)2n(23)

Aside from requiring π, which is a compile-time constant, the code implementing this function does not rely on any stored data, and so is thread-safe.
The library subroutine name is mpberne (see Table 5).
BesselI function. [MPFUN20-Fort and MPFUN20-MPFR.] The BesselI function, or, more formally, the modified Bessel function of the first kind, is42defined as [23, 10.25.2]:
Iν (z) = 2 kz4∞
 z ν X

2

k=0

k!Γ(ν + k + 1)

(24)

For modest-sized z, this function can be evaluated directly according to thisdefinition, terminating the summation when terms are smaller than epsilon
times the current sum. For larger z (the author uses the condition |z| > 0.4B,
where B is the precision in bits), the following asymptotic formula is faster[23, 10.40.1]:∞

ez X (−1)k ak (ν),
Iν (z) ≈ √zk
2πz k=0

(25)

where
(4ν 2 − 12 )(4ν 2 − 32 ) · · · (4ν 2 − (2k − 1)2 )
.
(26)
8k k!
Aside from requiring π, which is a compile-time constant, the code implementing this algorithm does not rely on any stored data, and so is thread-safe.
The library function names are bessel i and bessel in (see Table 5).
ak (ν) =BesselJ function. [MPFUN20-Fort only; MPFUN20-MPFR uses the MPFR
routine.] The BesselJ function, or, more formally, the Bessel function of thefirst kind, is defined as [23, 10.2.2]:
 2 k∞ z ν X(−1)k z4.
(27)
Jν (z) =2 k=0 k!Γ(ν + k + 1)
For modest-sized z, this function can be evaluated directly according to thisdefinition, terminating the summation when terms are smaller than epsilon
times the current sum. However, large amounts of cancellation occurs in thisformula, so the present author’s implementation employs a working precision
P = B + 1.4|z| bits, where B is desired result precision in bits, up to 3X therequested precision.
For larger z (the author uses the condition |z| > 0.4B, where B is theprecision in bits), the following asymptotic formula [23, 10.17.3] is faster:
!
r
∞
∞
X
X
2
(−1)k a2k (ν)(−1)k a2k+1 (ν)
Jν (z) ≈cos(ω)− sin(ω), (28)πzz 2k
z 2k+1k=0
k=043

where ω = z − πν/2 − π/4, and where ak (ν) is given by (26). Aside fromrequiring π, which is a compile-time constant, the code implementing thesealgorithms does not rely on any stored data, and so is thread-safe. The
library function names are bessel j, bessel jn, bessel j0 and bessel j1(see Tables 3 and 5).
BesselK function. [MPFUN20-Fort only; MPFUN20-MPFR uses the MPFR
routine.] The BesselK function, or, more formally, the modified Bessel function of the second kind, is defined as [23, 10.31.1]:
n−1

1  z −n X (n − k − 1)!
Kn (z) =2 2k!
k=0+ (−1)n1∞
 z n X

2 2

 2 kz z+ (−1)n+1 log
In (z)−
42
 2 k(ψ(k + 1) + ψ(n + k + 1))

k=0

z
4

k!(n + k)!

,

(29)

where the ψ(n) (digamma) function for positive
integer arguments can bePn−1
evaluated from the simple formula ψ(n) = −γ + k=1 1/k, where γ is Euler’s
constant. For modest-sized z, this function can be evaluated directly according to this definition, terminating the summation when terms are smaller
than epsilon times the current sum. However, large amounts of cancellation
occurs in this formula, so the present author’s implementation employs aworking precision B + 2.8|z| bits, where B is desired result precision in bits,
up to 3X the requested precision.
For larger z (the author uses the condition |z| > 0.4B, where B is theprecision in bits), the following asymptotic formula [23, 10.40.1] is faster:rKν (z) ≈∞π −z X ak (ν)
e
,
2z
zk
k=0

(30)

where ak (ν) is given by (26).
In the case when the order ν is not an integer, this formula is used instead:Kν (z) =π(I−ν − Iν ).
2 sin(νπ)(31)

Aside from requiring γ and π, which are compile-time constants, the code
implementing these algorithms does not rely on any stored data, and so is44thread-safe. The library function names are bessel k and bessel kn (seeTable 5).
BesselY function. [MPFUN20-Fort and MPFUN20-MPFR.] The BesselY
function, or, more formally, the Bessel function of the second kind, is definedas [23, 10.8.1]:
"
 k
n−1
z 
1  z −n X (n − k − 1)! z 2+ 2 logJn (z)Yn (z) =
π
2
k!
4
2
k=0
 2 k ∞ z n X
− z4

+
(32)
(ψ(k + 1) + ψ(n + k + 1)),2 k=0
k!(n + k)!
where the ψ(n) (digamma) function for positive
integer arguments can bePn−1
evaluated from the simple formula ψ(n) = −γ + k=1 1/k, where γ is Euler’s
constant. For modest-sized z, this function can be evaluated directly according to this definition, terminating the summation when terms are smaller
than epsilon times the current sum. However, large amounts of cancellation
occurs in this formula, so the present author’s implementation employs aworking precision P = B + 1.4|z| bits, where B is desired result precision inbit, up to 3X the requested precision.
In the case when the order ν is not an integer, this formula is used instead:Yν (z) =Jν cos(νπ) − J−ν.
sin(νπ)(33)

For larger z (the author uses the condition |z| > 0.4B, where B is theprecision in bits), the following asymptotic formula [23, 10.17.4] is faster:
!
r
∞
∞
X
X
2
(−1)k a2k (ν)(−1)k a2k+1 (ν)
Jν (z) ≈sin(ω)+ cos(ω), (34)
2k
2k+1
πz
z
z
k=0
k=0
where ω = z − πν/2 − π/4, and where ak (ν) is given by (26).
Aside for requiring γ and π, which are compile-time constants, the code
implementing these algorithms does not rely on any stored data, and so isthread-safe. The library function names are bessel y, bessel y0, bessel y1
and bessel yn (see Tables 3 and 5).
45

Digamma function. [MPFUN20-Fort and MPFUN20-MPFR; MPFUN20MPFR also has a separate routine from MPFR.] This evaluates the digammafunction, using the following asymptotic formula [23, 5.11.2]:
∞
X
B2k
1
−
ψ(x) ≈ log(x) −,2x k=1 2kx2k(35)

where B2k are the even Bernoulli coefficients (see Even Bernoulli coefficientsabove). The code implementing this algorithm is entirely thread-safe, provided the requisite even Bernoulli numbers are precomputed.
The library function name is digamma be. Before calling digamma be,subroutine mpberne must be called to compute the even Bernoulli numbers,
where the parameter nb is at least twice the desired precision in digits, and theresulting array and size must be included in subsequent calls to digamma be.
See Table 5.
Error functions. [MPFUN20-Fort only; MPFUN20-MPFR uses the MPFR
routines.] The error function is defined as [23, 7.6.2]:∞
2n z 2n+1
2 X.
erf (z) = √ z2
πe n=0 1 · 3 · · · (2n + 1)

(36)

This formula is satisfactory for computation unless the argument z is large,
in which case the following asymptotic expression for the error function complement works better:
erfc (z) ≈ √∞1 X (−1)n 1 · 3 · · · (2n − 1).
2n z 2n+1
πez2 n=0(37)where erfc (z) = 1 − erf (z), and where the product in the numerator is 1
when n = 0. In the author’s code, (36) is used when |z| < B/100, and (37)is used when |z| ≥ B/100, where B is the number of bits of precision. Asidefrom requiring π, which is a compile-time constant, the code implementing
this algorithm does not rely on any stored data, and so is thread-safe. The
library function names are erf and erfc (see Tables 3 and 5).
Exponential integral function. [MPFUN20-Fort only; MPFUN20-MPFR uses
the MPFR routine.] The exponential integral function is evaluated using the46

following formula due to Ramanujan [24]:Ein(z) = γ + log(|z|) + exp(z/2)

b(k−1)/2c∞
X(−1)k−1 z k X
k=1

2k−1 k!

j=0

1
2j + 1= −Γ(0, −z)

(38)
(39)

See incomplete gamma function below for Γ(0, −z)).
Gamma function. [MPFUN20-Fort only; MPFUN20-MPFR uses the MPFR
routine.] The gamma function employs a very efficient but little-known formula due to Ronald W. Potter [36] for most values, combined with a separate formula for values close to one. If the input t is a positive integer, thenΓ(t) = (t − 1)!. If not, use the recursion Γ(t + 1) = tΓ(t) to reduce theargument (positive or negative) to the interval (0, 1). If 1 − 1/1010 < t < 1,then employ this formula from the DLMF [23, 5.7.3]:∞
Xxk(−1) (ζ (k) − 1) ,log(Γ (1 + x)) = − log (1 + x) + x(1 − γ) +
k
k=2
k

(40)

where x = t − 1 and γ is Euler’s constant, and where ζ(k) is computed using
Peter Borwein’s algorithm, as described below under Zeta. For other valuesof t, define α = 2 int (b/4 · log 2 + 1), where b is the number of bits of precisionand int is greatest integer. Define the Pochhammer symbol as(t)k = t(t + 1)(t + 2) · · · (t + k − 1).

(41)

Then define the functions∞
 α t X(α2 /4)kt
A(t, α) =2
k!(t)k+1
k=0∞ α −tX(α2 /4)k
B(t, α) =(−t).
2
k!(−t)k+1
k=0With these definitions, the gamma function can then be computed assπA(t, α)Γ(t) =.
t sin(πt)B(t, α)
47

(42)

(43)

Aside from requiring π and γ, which are compile-time constants, the code
implementing this function does not rely on any stored data, and so is threadsafe. The library function name is gamma (see Tables 3 and 5).
Hurwitz zeta function. [MPFUN20-Fort and MPFUN20-MPFR.] The Hurwitz zeta function is computed by the following formula, which is a modification of a scheme originally due to Richard Crandall [6]: Let the notation
σ(x) denote sign(x). Then for 0 < a < 1 and for s 6= 1, 0, −1, −2, · · · ,√ (s−1)/2ππ
ζ (s, a) =(s − 1) Γ (s/2)
∞ 
1 X Γ( 2s ,π (n+a)2 )1
Γ( s+1,π (n+a)2 )2++ σ(n + a) Γ((s+1)/2)
Γ(s/2)
2 n=−∞|n + a|s  1−s m2 π2 ∞s m2 π 2
X
Γ
,
Γ
1−
,
1
2
π
2
π
s−1/2
+π
cos (2 πma) + Γ((s+1)/2) sin (2 π ma) .
Γ(s/2)
m1−s
m=1
(44)
Although each of the individual terms are relatively expensive to computedue to the incomplete gamma functions, convergence is quite rapid: 30 terms
produce over 1000 digits. Note that when a = 1, ζ (s, a) = ζ(s). Aside fromrequiring π, which is a compile-time constant, the code implementing thisfunction does not rely on any stored data, and so is thread-safe. The libraryfunction name is hurwitz zetan (see Table 5).
[MPFUN20-Fort and MPFUN20-MPFR.] An even faster formula, which
relies on precomputed even Bernoulli numbers (see Even Bernoulli numbers
above), is as follows. Select an integer q > 0 (the present author uses q =nint (0.6d), where d is the precision in digits and nint is nearest integer).
Thenζ(s, a) =
+

q−1
X

1
1
1
+
+
s
s−1
(a + k)(s − 1)(a + q)2(a + q)s
k=0∞
XB2k s(s + 1) · · · (s + 2k − 2)
k=1(2k)!(a + q)s+2k−1

.

(45)

The code implementing this function (hurwitz zetan be) relies only on precomputed even Bernoulli numbers, and thus is thread-safe. Before callinghurwitz zetan be, subroutine mpberne must be called to compute the even48

Bernoulli numbers, where the parameter nb is at least twice the desired precision in digits, and the resulting array and size must be included in subsequentcalls to hurwitz zetan be (see Table 5).
Hypergeometric pFq function. [MPFUN20-Fort and MPFUN20-MPFR.] The
hypergeometric pFq function is defined as: X∞(a1 )k (a2 )k · · · (ap )k xk
a1 , a2 , · · · , ap ; x
,
(46)
=
p Fq
b1 , b2 , · · · , bq(b1 )k (b2 )k · · · (bq )k k!
k=0
where (s)n = s(s+1)(s+2) · · · (s+n−1) is the Pochhammer symbol as before.
This formula is evaluated directly as stated, terminating the summation whenterms are smaller than epsilon times the current sum. Note that the terms
satisfy a relatively simple recursion. The code implementing this function
does not rely on any stored data, and so is thread-safe. The library function
name is hypergeom pfq (see Table 5).
Incomplete gamma function. [MPFUN20-Fort only; MPFUN20-MPFR uses
the MPFR routine.] For modest-sized positive arguments (the author uses
the condition |z| < 0.833B, where B is the precision level in bits), the
MPFUN20-Fort incomplete gamma function is evaluated using the following formula [23, 8.7.3]:
!
∞
za X
zk
.
(47)
Γ(a, z) = Γ(a) 1 − z
e k=0 Γ(a + k + 1)
Note, as with the BesselJ function, that although formula (47) involves the
gamma function, this is only called once to compute Γ(a + 1), after which
the recursion Γ(t + 1) = tΓ(t) can be applied for all other terms.
For large values of z, the following asymptotic formula is used [23, 8.11.2]:∞

Γ(a, z) ≈z a−1 X (−1)k (1 − a)k.
ez k=0zk(48)

For the case a = 0, the following formula is employed, which is a modificationof a formula due to Ramanujan for the exponential integral function [24]:Γ(0, z) = −γ − log(|z|) + exp(−z/2)

49

∞
X

zk

2k−1 k!
k=1b(k−1)/2c

X
j=0

1
.
2j + 1(49)

For z < 0, formula (49) returns the correct real part of the complex value.
For a = −1, −2, · · · , the following formula is employed [31]:
!
n−1
X
1
1
(−1)n Γ(0, z) + z nΓ(−n, z) =(−1)k (n − k − 1)!z k .
(50)n!
e z k=0
One difficulty is that cancellation may occur in formulas (49) and (50), so forcertain arguments the working precision is increased, up to 3X the requested
precision. Aside from requiring Euler’s constant γ in the formula for Γ(0, z),
the code implementing this function does not require any saved data andthus is thread-safe. The library function name is incgamma (see Table 5).
Polygamma function. [MPFUN20-Fort and MPFUN20-MPFR.] The polygammafunction is evaluated using the Hurwitz zeta function via the formula [29]
ψ (m) = (−1)m+1 m!ζ(m + 1, z).

(51)

See Hurwitz zeta function above for details. The library function name ispolygamma; an alternate function name is polygamma be, requiring precomputed even Bernoulli numbers (see Table 5).
Polylogarithm function. [MPFUN20-Fort and MPFUN20-MPFR.] The polylogarithm function for integer order n is defined as [23, 25.12.10]:
Lin (z) =

∞
X
zk
k=1

kn

(52)

This is satisfactory for computation for n ≥ 1 and |z| < 1. When n = 0, one
uses the formula Li0 (z) = z/(1 − z). For negative order n, one can use theformula [40]:nX
1an,k z k ,Li−n (z) =n+1(1 − z)k=1(53)where the coefficients an,k are given by the recurrencean,k = (n + 1 − k)an−1,k−1 + kan−1,k ,(54)

with an,1 = 1. One difficulty is that cancellation may occur in formula (53),so for certain arguments the working precision is increased, up to 3X the50

requested precision. Since the requisite auxiliary data is precomputed andpassed via argument, the code implementing these algorithms is thread-safe.
Because of the different algorithms and data requirements, MPFUN2020
includes three routines for polylogarithms, namely polylog ini, polylog neg
and polylog pos. Before calling polylog neg, subroutine polylog ini must
be called to compute the coefficients an,k for the given value of n, and theresulting array and size must be included in subsequent calls to polylog neg(see Table 6).
StruveH fuction. [MPFUN20-Fort and MPFUN20-MPFR.] The StruveH
function is defined for positive integer order n as [23, 11.2.1]:
Hn (z) =∞ z n+1 X2(−1)k (z/2)2k.
Γ(k+3/2)Γ(n
+
k
+
3/2)
k=0

(55)

This formula is evaluated as stated. As with the BesselJ function, the gammafunction references√ here can be computed iteratively, starting with the factthat Γ(3/2) = π/2. One difficulty is that due to the alternating sign in the
numerator, precision is lost due to cancellations. In the author’s code, thisis overcome by employing B + B|z|/1000 bits of precision, where B is thenumber of bits desired for the result, up to 3X the requested precision. Asidefrom requiring π, which is a compile-time constant, the code implementing
this function does not rely on any stored data, and so is thread-safe. The
library function name is struve hn (see Table 6).
Zeta function. [MPFUN20-Fort only; MPFUN20-MPFR uses the MPFR
routine.] For large positive arguments s (the present author uses the condition s > B log(2)/ log(2B/3)), where B is the precision in digits), it suffices
to use the definition of zeta, namely
ζ(s) =

∞
X
1
n=1

ns

.

(56)

For modest-sized arguments, the zeta function can be evaluated by means of
this formula, due to Peter Borwein [18]. Select n to be the number of digitsof precision required for the result. Define
!
j−n
X
n!
ej = (−1)j− 2n ,
(57)
k!(n
−
k)!
k=0
51

where the empty sum is zero. Then
2n−1
X
ej
−1
ζ(s) ≈ n.
1−s
2 (1 − 2 ) j=0 (j + 1)s(58)

There is some loss of precision due to cancelation with these formulas, so theauthor’s code employs a working precision 12.5% higher than the requested
precision.
The above formulas are used for positive real arguments (except s = 1, forwhich the zeta function is undefined). For negative s, Riemann’s reflection
formula is used to convert the calculation to a positive argument:ζ(s) =

2 cos(π(1 − s)/2)Γ(1 − s)ζ(1 − s).
(2π)1−s(59)

Aside from the value of π, no auxiliary data for this algorithm required, so
the code implementing these algorithms is thread-safe. The library formula
names are zeta and zeta int (see Table 6).
[MPFUN20-Fort and MPFUN20-MPFR.] An even faster algorithm forthe zeta function, based on the Euler-Maclaurin summation formula, can bederived from the following [23, 25.2.9]: Select an integer parameter N > 0(the present author uses N = nint (0.6d), where d is the precision in digits
and nint is nearest integer). Then

∞ 
X
1
1
s + 2k − 2
B2k
+
−
+
, (60)
ζ(s) =ss−1sk(s − 1)N2N
2k − 1
2kN s−1+2k
k=1
k=1
N
X
1

where B2k are the even Bernoulli numbers (see Even Bernoulli numbers
above). The library function name is zeta be. Before calling zeta be, thesubroutine mpberne must be called to compute the even Bernoulli numbers,
where the parameter nb is greater than twice the desired precision in digits, and the resulting array and size must be included in subsequent calls tozeta be (see Table 6).
Polynomial solution routines. As briefly mentioned above, the package includes several routines for polynomial evaluation and solution. mprpolysolve
and mpcpolysolve, respectively, solve an input MPR polynomial, given an
MPR or MPC initial approximation. Both routines employ Newton iterations
with a working precision that approximately doubles with each iteration.
52

Polynomials with repeated roots often defeat Newton iteration schemes
for solution. Such polynomials may be solved by first applying the polynomial greatest common divisor (GCD) routine, provided in the package, to the
input polynomial and its formal derivative, and then divide the input polynomial by the GCD polynomial, using the polynomial division routine. This
produces a reduced polynomial that has all the roots of the original polynomial, but no repeated roots. Then the Newton iteration-based solution
routines can quickly produce high-precision roots.
The polynomial GCD routine is mprpolygcd, and the polynomial division routine, which produces the quotient and remainder polynomials, ismprpolydiv. These routines employ straightforward implementations of thepolynomial GCD and division algorithms. Because the polynomial GCD process is highly sensitive numerically, the author’s code for this routine employsa working precision 12.5% higher than the requested precision.
Two other routines (mprpolyeval) and (mpcpolyeval) evaluate a polynomial at an MPR or MPC argument, respectively. See Table 7 for details.

53

References[1] D. H. Bailey, “MPFUN2015: A thread-safe arbitrary precision package(full documentation),” manuscript, 18 Nov 2020, preprint at
https://www.davidhbailey.com/dhbpapers/mpfun2015.pdf[2] D. H. Bailey, “Finding large Poisson polynomials using four-levelvariable precision,” Correctness 2021: 5th International Workshop on
Software Correctness for HPC Applications, 20 Nov 2021, preprint at
https://www.davidhbailey.com/dhbpapers/correct-2021.pdf.
[3] D. H. Bailey, “Large Poisson polynomials: Computation, results andanalysis,” manuscript, 5 Mar 2023, available athttps://www.davidhbailey.com/dhbpapers/poisson-2023.pdf.
[4] D. H. Bailey, X. S. Li and K. Jeyabalan, “A comparison of threehigh-precision quadrature schemes,” Experimental Mathematics, vol.
14 (2005), no. 3, 317–329.
[5] D. H. Bailey, R. Barrio and J. M. Borwein, “High precisioncomputation: Mathematical physics and dynamics,” Applied
Mathematics and Computation, vol. 218 (2012), 10106-10121.
[6] D. H. Bailey and J. M. Borwein, “Crandall’s computation of theincomplete gamma function and the Hurwitz zeta function withapplications to Dirichlet L-series,” Applied Mathematics andComputation, vol. 268C (Oct 2015), 462–477, preprint at
https://www.davidhbailey.com/dhbpapers/lerch.pdf.
[7] D. H. Bailey and J. M. Borwein, “High-precision arithmetic inmathematical physics,” Mathematics, vol. 3 (2015), 337–367,
http://www.mdpi.com/2227-7390/3/2/337/pdf.
[8] Jonathan M. Borwein and David H. Bailey, Mathematics byExperiment: Plausible Reasoning in the 21st Century, AK Peters,
Natick, MA, 2008.
[9] D. H. Bailey and J. M. Borwein, “Hand-to-hand combat withthousand-digit integrals,” Journal of Computational Science, vol. 3(2012), 77–86.
54

[10] D. H. Bailey, J. M. Borwein and R. E. Crandall, “On the Khintchine
constant,” Mathematics of Computation, vol. 66 (1997), 417–431,
preprint at
https://www.davidhbailey.com/dhbpapers/khinchine.pdf.
[11] 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), 12271–12302.
[12] D. H. Bailey and J. M. Borwein, “Compressed lattice sums arisingfrom the Poisson equation: Dedicated to Professor Hari Sirvastava,”
Boundary Value Problems, vol. 75 (2013), DOI:
10.1186/1687-2770-2013-75,
http://www.boundaryvalueproblems.com/content/2013/1/75.
[13] D. H. Bailey, J. M. Borwein, R. E. Crandall and J. Zucker, “Latticesums arising from the Poisson equation,” Journal of Physics A:
Mathematical and Theoretical, vol. 46 (2013), 115201, preprint at
http://www.davidhbailey.com/dhbpapers/PoissonLattice.pdf.
[14] D. H. Bailey, J. M. Borwein, J. Kimberley and W. Ladd, “Computer
discovery and analysis of large Poisson polynomials,” Experimental
Mathematics, 27 Aug 2016, vol. 26 (2016), 349–363, preprint at
https://www.davidhbailey.com/dhbpapers/poisson-res.pdf.
[15] D. H. Bailey and D. J. Broadhurst, “Parallel integer relation detection:
Techniques and applications,” Mathematics of Computation, vol. 70,
no. 236 (Oct 2000), 1719–1736.
[16] D. H. Bailey, X. S. Li and B. Thompson, “ARPREC: An arbitrary
precision computation package,” Sep 2002, preprint at
http://crd.lbl.gov/~dhbailey/dhbpapers/arprec.pdf.
[17] J. M. Borwein and D. H. Bailey, Mathematics by Experiment: Plausible
Reasoning in the 21st Century, 2nd ed., A. K. Peters, Natick, MA,
2008.
[18] P. Borwein, “An efficient algorithm for the Riemann zeta function,”
1995, http://www.cecm.sfu.ca/~pborwein/PAPERS/P155.pdf.

55

[19] R. P. Brent, “Fast multiple-precision evaluation of elementaryfunctions,” J. of the ACM, vol. 23 (1976), 242–251.
[20] R. P. Brent and P. Zimmermann, Modern Computer Arithmetic,
Cambridge Univ. Press, 2010.
[21] Jingwei Chen, Damien Stehle and Gilles Villard, “A new view on
HJLS and PSLQ: Sums and projections of lattices,” Proc. of
ISSAC’13, 149–156.
[22] “Comparison of multiple-precision floating-point software,”http://www.mpfr.org/mpfr-3.1.0/timings.html.
[23] “Digital library of mathematical functions,” National Institute ofStandards and Technology, 2015, http://dlmf.nist.gov.
[24] “Exponential integral,” Wikipedia article, viewed 3 Feb 2022,
https://en.wikipedia.org/wiki/Exponential_integral.
[25] 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), 351–369.
[26] “Fortran 2008,”
http://fortranwiki.org/fortran/show/Fortran+2008.
[27] “The GNU MPFR library,” http://www.mpfr.org.
[28] “GNU MPFR library: Comparison of multiple-precision floating-point
software,” http://www.mpfr.org/mpfr-current/timings.html.
[29] “Polygamma function,” Wikipedia page, viewed 28 Jan 2022,
https://en.wikipedia.org/wiki/Polygamma_function.
[30] Y. Hida, X. S. Li and D. H. Bailey, “Algorithms for Quad-Double
Precision Floating Point Arithmetic,” Proc. of the 15th IEEE
Symposium on Computer Arithmetic (ARITH-15), 2001.
[31] “Incomplete gamma function,” Wikipedia article, viewed 28 Jan 2022,
https://en.wikipedia.org/wiki/Incomplete_gamma_function.

56

[32] A. K. Lenstra, H. W. Lenstra, and L. Lovasz, “Factoring polynomials
with rational coefficients,” Mathematische Annalen, vol. 261 (1982),
515–534.
[33] “MPFR C++,” http://www.holoborodko.com/pavel/mpfr.
[34] NIST Digital Library of Mathematical Functions, version 1.0.6 (May2013), http://dlmf.nist.gov.
[35] E. Panzer, “Algorithms for the symbolic integration of
hyperlogarithms with applications to Feynman integrals,” manuscript,
2014, https://arxiv.org/abs/1403.3385.
[36] R. W. Potter, Arbitrary Precision Calculation of Selected Higher
Functions, Lulu.com, San Bernardino, CA, 2014.
[37] D. W. Sweeney, “On the computation of Euler’s constant”,
Mathematics of Computation, vol. 17 (1963), 170–178.
[38] H. Takahasi and M. Mori, “Double exponential formulas for numericalintegration,” Publications of RIMS, Kyoto University, vol. 9 (1974),
721–741.
[39] S. W. Williams and D. H. Bailey, “Parallel computer architecture,” inDavid H. Bailey, Robert F. Lucas and Samuel W. Williams, ed.,
Performance Tuning of Scientific Applications, CRC Press, Boca
Raton, FL, 2011, 11–33.
[40] D. C. Wood, “The computation of polylogarithms,” undated
manuscript,
https://www.cs.kent.ac.uk/pubs/1992/110/content.pdf.

57