Reproducibility in scientific computing
A December 2012 workshop on reproducibility in computing, held at Brown University
in Rhode Island, USA, noted that
Science is built upon the foundations of theory and experiment validated andimproved through open, transparent communication. With the increasingly
central role of computation in scientific discovery, this means communicating
all details of the computations needed for others to replicate the experiment.
... The “reproducible research” movement recognizes that traditional
scientific research and publication practices now fall short of this ideal, andencourages all those involved in the production of computational science ...
to facilitate and practice really reproducible research.

I V. Stodden, D. H. Bailey, J. Borwein, R. J. LeVeque, W. Rider and W. Stein, “Setting the default toreproducible: Reproducibility in computational and experimental mathematics,”
http://www.davidhbailey.com/dhbpapers/icerm-report.pdf.

Reproducibility in scientific computing, continued
Issues identified in the ICERM report and other studies include:
I The need to carefully document the full context of computational

experiments—system environment, input data, code used, computed results, etc.
I The need to save the code and output data in a permanent repository.
I The need for reviewers, research institutions and funding agencies to recognize theimportance of computing and computing professionals, and to allocate funding forafter-the-grant support and repositories.
I The increasing importance of numerical reproducibility, and the need for tools toensure and enhance numerical reliability.
I The need to encourage publication of negative results—other researchers canoften learn from them.
I The re-emergence of the need to ensure responsible reporting of performance.

Reproducibility in biomedicine
The biomedical field has been stung by numerous cases where pharma products or lab
discoveries look great, but then fall flat or cannot be reproduced. Examples:
I In 2004, GlaxoSmithKline acknowledged that while some trials of Paxil found it
effective for depression in children, other unpublished studies showed no benefit.
I In 2011, Bayer researchers reported that they were able to reproduce the results of
only 17 of 67 published studies they examined.
I In 2012, Amgen researchers reported that they were able to reproduce the resultsof only 6 of 53 published cancer studies.
I In 2014, a review of Tamiflu found that while it made flu symptoms disappear a
bit sooner, it did not stop serious complications or keep people out of the hospital.
I In June 2014, Japanese retracted an earlier, highly publicized claim that they had
produced pluripotent stem cells from ordinary skin or blood cells.
These experiences have exposed a fundamental flaw in methodology:
Only publicizing the results of successful trials introduces a bias into the results.
The AllTrials movement would require all results to be public: http://www.alltrials.net

Reproducibility in physics
In March 2014, a team of researchers from Harvard
University made the dramatic announcement thatthey had discovered an interesting “twisting”
pattern in cosmic microwave background data,measured using their BICEP2 experimental system.
This pattern fit very well with the hypothesizedpattern of the most commonly assumed model of“inflation” in the first tiny fraction of a second after
the big bang, and thus has been trumpeted as thefirst experimental evidence of the inflationary
cosmology.
But other researchers had difficulty reconstructing
the claimed results. Finally, two teams challenged
the BICEP2 findings, saying that the results could
more readily be explained by dust in the Milky Way.

I Ron Cowen, “Doubt grows about
gravitational waves detection,”
Scientific American, 2 Jun 2014.

Reproducibility in social science
The “blank slate” paradigm (1920–1990):
I The human mind at birth is a tabula rasa (“blank slate”).
I Heredity and biology play no significant role in human psychology; all personality

and behavioral traits are socially constructed.
I Pre-modern societies were serene, devoid of most psychological and social

problems that afflict modern societies.
As Ashley Montagu wrote in 1973,
Man is man because he has no instincts, because everything he is and has
become he has learned, acquired, from his culture, from the man-made partof the environment, from other human beings.

The fall of the blank slate
The current consensus in the field, based on latest research:
I Humans at birth possess sophisticated facilities for social interaction, language

acquisition, pattern recognition, navigation and counting.
I Heredity, evolution and biology are major factors in human personality.
I Some personality traits are more than 50% heritable.
I Pre-modern societies had more crime, war and social problems than modern

western societies.
How did the early 20th century social scientists err so badly?
I Sloppy experimental methodology and analysis.
I Pervasive wishful thinking and politically correct biases.
I Ignoring or dismissing data that runs counter to predisposition.
I S. Pinker, The Blank Slate: The Modern Denial of Human Nature, Penguin Books, 2003.

Science adds statistical checks to peer-review
The journal Science is adding an extra round of statistical checks to its
peer-review process, editor-in-chief Marcia McNutt announced on July 3.
The policy follows similar efforts from other journals, after widespread
concern that basic mistakes in data analysis are contributing to theirreproducibility of many published research findings.
“Readers must have confidence in the conclusions published in our journal,”
writes McNutt in an editorial today. Working with the American Statistical
Association, the journal has appointed seven experts to a statistics board ofreviewing editors (SBoRE). Manuscript will be flagged up for additional
scrutiny by the journal’s internal editors, or by its existing Board of Reviewing
Editors (more than 100 scientists whom the journal regularly consults on
papers) or by outside peer reviewers. The SBoRE panel will then findexternal statisticians to review these manuscripts.
I Richard Van Noorden, “Major scientific journal joins push to screen statistics in papers it publishes,”
Nature, 6 Jul 2014.

Reproducibility in finance
Finance, like the pharmaceutical world, has been stung with numerous instances ofinvestment strategies that look great on paper, but fall flat in practice. A primarycause is statistical overfitting of backtest (historical market) data.
When a computer can analyze thousands or millions of variations of a given strategy, it
is almost certain that the best such strategy, measured by backtests, will be overfit andthus of dubious value.
In two 2014 papers by myself and three co-authors, we derive (a) a formula relating
the number of trials to the minimum backtest length, and (b) a formula for the
probability of backtest overfitting. We also show that under the assumption of memory
in markets, overfit strategies are actually somewhat prone to lose money.
I D. H. Bailey, J. M. Borwein, M. Lopez de Prado and Q. J. Zhu, “Pseudo-mathematics and financial
charlatanism: The effects of backtest over fitting on out-of-sample performance,” Notices of theAmerican Mathematical Society, May 2014, pg. 458–471.
I D. H. Bailey, J. M. Borwein, M. Lopez de Prado and Q. J. Zhu, “The probability of backtest overfitting,”
manuscript, 12 Feb 2014, http://ssrn.com/abstract=2326253.

How easy is it to overfit a backtest? Answer: Very easy!
I If only 2 years of daily backtest data are available, then no more than 7 strategy

variations should be tried.
I If only 5 years of daily backtest data are available, then no more than 45 strategy

variations should be tried.
I A backtest that does not report the number of trials N makes it impossible toassess the risk of overfitting.
I Given any desired performance level, a financial researcher just needs to keep

trying alternative parameters for that strategy!

MinBTL ≈




 !2
(1 − γ)Z −1 1 − N1 + γZ −1 1 − N1 e −1
E [maxN ]An absurd investment strategy
I A financial advisor sends letters to 10, 240 = 10 × 210 potential clients, with 5120

letters predicting a certain security will go up, and the other half predicting it will
go down.
I One month later, the advisor sends letters only to the 5120 investors who werepreviously sent the correct prediction, with 2560 letters predicting a certain
security will go up, and the other half predicting it will go down.
I The advisor continues this process for 10 months.
I The remaining ten investors, so impressed by the advisor’s ten consecutive spot-on

predictions, will entrust to him/her all of their assets!
This strategy is absurd; even fraudulent.
But why is marketing a statistically overfit strategy, where potential investors are notinformed of the number of trials behind the strategy, any different?

Why the silence in the mathematical finance community?
I Historically scientists have led the way in exposing those who utilize pseudoscience

to extract a commercial benefit: e.g., in the 18th century, physicists exposed thenonsense of astrologers.
I Yet financial mathematicians in the 21st century have remained disappointingly
silent with the regards to those in the community who, knowingly or not:1. Fail to disclose the number of models or variations that were used to develop an
investment strategy.
2. Make vague predictions that do not permit rigorous testing and falsification.
3. Misuse probability theory, statistics and stochastic calculus.
4. Use dubious technical jargon: “stochastic oscillators,” “Fibonacci ratios,” “cycles,”“Elliot wave,” “Golden ratio,” “parabolic SAR,” “pivot point,” “momentum,” etc.

As we recently wrote in our paper “Pseudo-Mathematics and Financial Charlatanism”:“Our silence is consent, making us accomplices in these abuses.”
I D. H. Bailey, J. M. Borwein, M. Lopez de Prado and Q. J. Zhu, “Pseudo-mathematics and financial
charlatanism: The effects of backtest over fitting on out-of-sample performance,” Notices of theAmerican Mathematical Society, May 2014, pg. 458–471.

Reproducibility in highly parallel computing: History from 1990-1994
Background:
I Many new parallel systems had been introduced; each claimed theirs was best.
I Many researchers were excited about the potential of highly parallel systems.
I Few standard benchmarks and testing methodologies had been established.
I It was hard to reproduce published performance results; much confusion reigned.
I Overall, the level of rigor and peer review in the field was rather low.

In response, in 1991 I published a humorous essay ”Twelve Ways to Fool the Masses,”
poking fun at some of the abuses.
Since abuses continued, I presented a talk at Supercomputing 1992 and published apaper with specific examples.
I D. H. Bailey, “Misleading performance reporting in the supercomputing field,” Scientific Programming,
vol. 1., no. 2 (Winter 1992), pg. 141–151.

1991 paper: “Twelve ways to fool the masses in highly parallel computing”
1. Quote 32-bit performance results, not 64-bit results, but don’t mention this in paper.
2. Present performance figures for an inner kernel, then represent these figures as theperformance of the entire application.
3. Quietly employ assembly code and other low-level language constructs.
4. Scale up the problem size with the number of processors, but omit any mention of this.
5. Quote performance results projected to a full system.
6. Compare your results against scalar, unoptimized code on conventional systems.
7. When run times are compared, compare with an old code on an obsolete system.
8. Base Mflop/s rates on the operation count of the parallel implementation, instead of thebest practical serial algorithm.
9. Quote performance as processor utilization, parallel speedups or Mflop/s per dollar.
10. Mutilate the algorithm used in the parallel implementation to match the architecture.
11. Measure parallel run times on a dedicated system, but measure conventional run times ina busy environment.
12. If all else fails, show pretty pictures and animated videos, and don’t discuss performance.

1992 paper: Scaling performance results to full-sized system
In some published papers and conference presentations, performance results on
small-sized parallel systems were linearly scaled to full-sized systems, without evenclearly disclosing this fact.
Example: 8,192-CPU performance results were linearly scaled to 65,536-CPU results,
simply by multiplying by eight.
Excuse: “We can’t afford a full-sized system.”This and the other examples mentioned in the next few viewgraphs are presented in:
I D. H. Bailey, “Misleading performance reporting in the supercomputing field,” Scientific Programming,
vol. 1., no. 2 (Winter 1992), pg. 141–151.

1992 paper: Using inefficient algorithms on highly parallel systems
In many cases, inefficient algorithms were employed for the highly parallel
implementation, requiring many more operations, thus producing artificially highMflop/s rates:
I Numerous researchers cited parallel PDE performance based explicit schemes,

where implicit schemes were known to be much better.
Excuse: Explicit schemes “run better” on the researchers’ parallel system.
I One paper cited performance for computing a 3D discrete Fourier transform bydirect evaluation of the defining formula (8n2 operations), rather than by using afast Fourier transform (5n log2 n).
Excuse: Direct computation was “more appropriate” for the architecture being
analyzed.
Both examples violate a rule of professional performance reporting, namely to base theoperation count (when computing Mflop/s or Gflop/s rates) on the best practical serial
algorithm, no matter what scheme was actually used on the parallel system.

1992 paper: Not actually performing a claimed computation
Abstract of published paper: “The current Connection Machine implementation runs
at 300-800 Mflop/s on a full [64K] CM-2, or at the speed of a single processor of a
Cray-2 on 1/4 of a CM-2.”
I Excerpt from text: “This computation requires 568 iterations (taking 272

seconds) on a 16K Connection Machine.”
In other words, the computation was run on a 16K system, not on a 64K system;
the figures cited in the Abstract were merely multiplied by four.
I Excerpt from text: “In contrast, a Convex C210 requires 909 seconds to compute

this example. Experience indicates that for a wide range of problems, a C210 isabout 1/4 the speed of a single processor Cray-2.”
In other words, the computation mentioned in the Abstract was not actually run
on a Cray-2; instead, it was run on a Convex system, and a questionable
rule-of-thumb scaling factor was used to produce the Cray-2 rate.

1992 paper: Performance plot [parallel (lower) vs vector (upper)]3.5

3Time (Hours)

2.5

2

1.5

1

0.5

0

0

1000

2000

3000

4000

5000

6000

Number of Objects

7000

8000

9000

10000

1992 paper: Data for performance plot
Problem size(x axis)
20
40
80
160
990
9600

Parallel system
run time
8:18
9:11
11:59
15:07
21:32
31:36

Vector system
run time
0:16
0:26
0:57
2:11
19:00
3:11:50*

Details in text of paper:
I In last entry, the 3:11:50 figure is an “estimate.”
I The vector system code is “not optimized.”

Note that the parallel system is actually slower than the vector system for all cases,
except for the last (estimated) entry. Also, except for the last entry, all real data in thegraph is in the lower left corner.

Fast forward to 2014: New ways to fool the masses
I Cite performance rates for a run with only one processor core active in a shared-

memory multi-core node, producing artificially inflated performance (since there is
no shared memory interference) and wasting resources (since most cores are idle).
I

Example: Cite performance on “1024 cores,” even though the code was run on 1024
multicore nodes, one core per node, with 15 out of 16 cores idle on each node.

I Claim that since one is using a graphics processing unit (GPU) system, that

efficient parallel algorithms must be discarded in favor of more basic algorithms.
I Cite performance rates only for a core algorithm (such as FFT or linear systemsolution), even though full-scale applications have been run on the system.
I List only the best performance figure in the paper, even though the run was madenumerous times (recall the experience of pharmaceutical tests).
I Employ special hardware, operating system or compiler settings that are notappropriate for real-world production usage.
I Define “scalability” as successful execution on a large number of CPUs, regardlessof performance.

Numerical reproducibility in high-performance computing
The report mentioned above on reproducibility in high-performance computing noted:
Numerical round-off error and numerical differences are greatly magnified ascomputational simulations are scaled up to run on highly parallel systems. As
a result, it is increasingly difficult to determine whether a code has beencorrectly ported to a new system, because computational results quickly
diverge from standard benchmark cases. And it is doubly difficult for otherresearchers, using independently written codes and distinct computer
systems, to reproduce published results.

I V. Stodden, D. H. Bailey, J. Borwein, R. J. LeVeque, W. Rider and W. Stein, “Setting the default toreproducible: Reproducibility in computational and experimental mathematics,” Jan 2013, available at
http://www.davidhbailey.com/dhbpapers/icerm-report.pdf.

Numerical reliability
Many applications routinely use either 32-bit or 64-bit IEEE arithmetic, and employ
fairly simple algorithms, assuming that all is well. But problems can arise:
1. Large-scale, highly parallel simulations, running on systems with hundreds ofthousands or millions of processors — numerical sensitivities are greatly magnified.
2. Certain applications with highly ill-conditioned linear systems.
3. Large summations, especially those involving +/- terms and cancellations.
4. Long-time, iterative simulations (such as molecular dynamics or climate models).
5. Computations to resolve small-scale phenomena.
6. Studies in computational physics or experimental mathematics often require hugeprecision levels.
I D. H. Bailey, R. Barrio, and J. M. Borwein, “High precision computation: Mathematical physics anddynamics,” Applied Mathematics and Computation, vol. 218 (2012), pg. 10106–10121.

Analysis of collisions at the Large Hadron Collider
I The 2012 discovery of the Higgs boson at the ATLAS experiment in the LHC

relied crucially on the ability to track charged particles with exquisite precision (10microns over a 10m length) and high reliability (over 99% of roughly 1000
charged particles per collision correctly identified).
I Software: 5 millions line of C++ and python code, developed by roughly 2000

physicists and engineers over 15 years.
I Recently, in an attempt to speed up the calculation, researchers found that merely

changing the underlying math library resulted in some collisions being missed or
misidentified.
Questions:
I How serious are these numerical difficulties?
I How can they be tracked down?
I How can the library be maintained, producing numerically reliable results?

Enhancing reproducibility with selective usage of high-precision arithmetic
Problem: Findof the irregular function
P the arc length−kg (x) = x + 0≤k≤10 2 sin(2k x), over the interval (0, π)(using 107 abscissa points).
3.0

I If this computation is done with ordinary double

precision arithmetic, the calculation takes 2.59 seconds
and yields the result 7.073157029008510.
I If it is done using all double-double arithmetic (31-digit

2.5

2.0

1.5

1.0

0.5

accuracy), it takes 47.39 seconds seconds and yields theresult 7.073157029007832.
I But if only the summation is changed to double-double,

the result is identical to the double-double result (to 15
digits), yet the computation only takes 3.47 seconds.

0.5

1.0

1.5

2.0

2.5

3.0

Graph
P of g (x) =x+ 0≤k≤10 2−k sin(2k x),
over (0, π).

U.C. Berkeley’s CORVETTE project and the “Precimonious” toolObjective: Develop software facilities to find and ameliorate numerical anomalies inlarge-scale computations:
I Facilities to test the level of numerical accuracy required for an application.
I Facilities to delimit the portions of code that are inaccurate.
I Facilities to search the space of possible code modifications.
I Facilities to repair numerical difficulties, including usage of high-precision

arithmetic.
I Facilities to navigate through a hierarchy of precision levels (32-bit, 64-bit, 80-bit

or higher as needed).
The current version of this tool is known as “Precimonious.”

I C. Rubio-Gonzalez, C. Nguyen, H. D. Nguyen, J. Demmel, W. Kahan, K. Sen, D. H. Bailey and C. Iancu,“Precimonious: Tuning assistant for floating-point precision,” manuscript, May 2013.

Innocuous example where high precision is required for reproducible results
Problem: Find a polynomial to fit the data (1, 1048579, 16777489,
84941299, 268501249, 655751251, 1360635409, 2523398179, 4311748609) forarguments 0, 1, · · · , 8. The usual approach is to solve the linear system:






Pn n

k=1 xk..
P .

nn
k=1 xkPn
Pnk=1 xk2
k=1 xk..
P . n+1
k=1 xk

···
···
..
.
···

Pn

  Pnna0k
Pn k=1 xn+1
Pn k=1 yk  a1  
k=1 xkk=1 xk yk    ..  = ..
..
 .  

.
.
P
Pn
n
n
2n
a
x
y
x
n
k=1 k kk=1 k

A 64-bit computation (e.g., using Matlab, Linpack or LAPACK) fails to find thecorrect polynomial in this instance, even if one rounds results to nearest integer.
However, if Linpack routines are converted to use double-double arithmetic (31-digit
accuracy), the above computation quickly produces the correct polynomial:f (x) = 1 + 1048577x 4 + x 8 = 1 + (220 + 1)x 4 + x 8

Innocuous example, continued
The result on the previous page can be obtained with double precision using Lagrange
interpolation or the Demmel-Koev algorithm. But few scientists, outside of expert
numerical analysts, are aware of these schemes.
Besides, even these schemes fail for higher-degree problems. For example:(1, 134217731, 8589938753, 97845255883, 549772595201,
2097396156251, 6264239146561, 15804422886323, 35253091827713,
71611233653971, 135217729000001, 240913322581691, 409688091758593)
is generated by:
f (x) = 1 + 134217729x 6 + x 12 = 1 + (227 + 1)x 6 + x 12
Neither the Lagrange or Demmel-Koev algorithms get the right answer with double
precision. In contrast, a straightforward Linpack scheme, implemented withdouble-double arithmetic, works fine for this and a wide range of similar problems.

Other applications where high-precision arithmetic is useful or essential
1. Planetary orbit calculations (32 digits).
2. Supernova simulations (32–64 digits).
3. Certain components of climate modeling (32 digits).
4. Coulomb n-body atomic system simulations (32–120 digits).
5. Schrodinger solutions for lithium and helium atoms (32 digits).
6. Electromagnetic scattering theory (32–100 digits).
7. Scattering amplitudes of fundamental particles (32 digits).
8. Discrete dynamical systems (32 digits).
9. Theory of nonlinear oscillators (64 digits).
10. The Taylor algorithm for ODEs (100–600 digits).
11. Ising integrals from mathematical physics (100–1000 digits).
12. Problems in experimental mathematics (100–50,000 digits).
I D. H. Bailey, R. Barrio, and J. M. Borwein, “High precision computation: Mathematical physics anddynamics,” Applied Mathematics and Computation, vol. 218 (2012), pg. 10106–10121.

Climate modeling
I Climate and weather simulations are fundamentallychaotic: if microscopic changes are made to the current
state, soon the future state is quite different.
I In practice, computational results are altered even if

minor changes are made to the code or the system.
I This numerical variation is a major nuisance for code

maintenance.
I He and Ding found that by using double-double

arithmetic in two key inner loops of an atmosphericmodeling code, most of this numerical variation
disappeared.
I Y. He and C. Ding, “Using accurate arithmetics to improve numerical reproducibility and stability inparallel applications,” Journal of Supercomputing, vol. 18, no. 3 (Mar 2001), pg. 259–277.

Using high-precision arithmetic to discover new mathematical relations
Very high-precision arithmetic is essential to obtain reproducible results in experimental
mathematics.
Methodology:
1. Compute various mathematical entities (limits, infinite series sums, definite
integrals, etc.) to high precision, typically 100–10,000 digits.
2. Use an algorithm such as “PSLQ” to recognize these numerical values in terms ofwell-known mathematical constants.
3. When results are found experimentally, seek formal mathematical proofs of the
discovered relations.
Many results have recently been found using this methodology, both in puremathematics and in mathematical physics.
“If mathematics describes an objective world just like physics, there is no reason why
inductive methods should not be applied in mathematics just the same as inphysics.” – Kurt GodelIsing integrals from mathematical physics
We recently applied this methodology to study three classes of integrals, two of whicharise in mathematical physics:
Z
Z ∞
du1
4 ∞
1
dun
Cn :=
···
···
P
2
n! 0
u1
un
n
0
j=1 (uj + 1/uj )Q  ui −uj 2
Z ∞Z ∞i<j ui +uj
dun
4
du1
···
···
Dn :=
P
2
n! 0
u1
un
n
0
j=1 (uj + 1/uj )

2
Z 1
Z 1
Y uk − uj dt2 dt3 · · · dtn
En = 2···
uk + uj
0
0
1≤j<k≤nwhere in the last line uk = t1 t2 · · · tk .
I D. H. Bailey, J. M. Borwein and R. E. Crandall, “Integrals of the Ising class,” Journal of Physics A:
Mathematical and General, vol. 39 (2006), pg. 12271–12302.

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

n→∞where γ = 0.57714 . . . is Euler’s constant. This result and numerous related results
were subsequently proven.

Other Ising integral evaluations found using PSLQ and 500-digit arithmetic
D2 = 1/3
D3 = 8 + 4π 2 /3 − 27 Li−3 (2)D4 = 4π 2 /9 − 1/6 − 7ζ(3)/2E2 = 6 − 8 log 2E3 = 10 − 2π 2 − 8 log 2 + 32 log2 2
E4 = 22 − 82ζ(3) − 24 log 2 + 176 log2 2 − 256(log3 2)/3+16π 2 log 2 − 22π 2 /3
E5 = 42 − 1984 Li4 (1/2) + 189π 4 /10 − 74ζ(3) − 1272ζ(3) log 2+40π 2 log2 2 − 62π 2 /3 + 40(π 2 log 2)/3 + 88 log4 2+464 log2 2 − 40 log 2
where ζ(x) is the Riemann zeta function and Lin (x) is the polylogarithm function.

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

where α is an algebraic number, i.e., the root of an integer polynomial
0 = a0 + a1 α + a2 α 2 + · · · + an α n
The minimal polynomials for these α were found by PSLQ calculations, with the(n + 1)-long vector (1, α, α2 , · · · , αn ) as input, where α = exp(8πφ2 (x, y )). PSLQ
returned the vector of integer coefficients (a0 , a1 , a2 , ·, an ) as output.
I D. H. Bailey, J. M. Borwein, R. E. Crandall and J. Zucker, “Lattice sums arising from the Poisson
equation,” Journal of Physics A: Mathematical and Theoretical, vol. 46 (2013), pg. 115201.

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

10

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

The minimal polynomial corresponding to x = y = 1/32 has degree 128, withindividual coefficients ranging from 1 to over 1056 . This PSLQ computation required
10,000-digit precision.
Other polynomials required up to 50,000-digit precision.

Reproducibility in symbolic computing
Mathematicians, physicists and many others utilize symbolic computing packages such
as Maple and Mathematica to perform symbolic computing. But, like all software, they
have bugs.
For example, consider
Z 1Z 1W2 =0e 2πix + e 2πiy dx dy ,(1)0

The currently available versions of Maple and Mathematica both report that thisintegral is zero, despite the obvious fact that it must be strictly positive. The correct
result is 4/π = 1.2732395447351626862 . . .
Numerous other examples could be cited. Sometimes the bugs are not even consistent— correct results in some runs, incorrect results in others.

Potential sources of computational errors

I User programming errors.
I Bugs in symbolic or numerical software packages.
I Bugs in compilers.
I Bugs in operating systems or other system software.
I Hardware errors — memory errors, processor errors, data communication errors. In

current exascale and petascale systems, such errors must be expected, not ignored.
I Errors in data storage media, or in communication to/from such media.

So why should anyone believe the results of any computation?

Methods to enhance reliability of computed results
1. Perform a computation using two different algorithms, then compare final results.
I

This is widely employed, say in computing digits of π: if the results of two suchcalculation are the same, except for a few trailing digits, then this is strong evidence
that both calculations are almost certainly correct.

2. Perform a computation using two different software packages (e.g., Maple andMathematica).
3. Perform a computation with two teams of programmers on different systems —climate models are sometimes validated in this way.
4. Perform a computation using two different levels of numeric precision, to see how
many digits agree.
5. Perform internal validity checks:
I

For example, if a matrix and its inverse are both advanced in the computation,
periodically multiply them together to ensure that the result is still the identity, towithin acceptable numerical error.

Summary
I The issue of reproducibility has recently emerged in scientific computing, just as it

has in the pharmaceutical industry, finance and others.
I The field of scientific computing is woefully behind other fields in reproducible

practices — making careful records of runs, saving source code, analyzing results,
etc. New tools may help.
I The issue of questionable performance reporting practices is again emerging,particularly with the deployment of new advanced architectures.
I Numerical reproducibility is of particular concern, in light of huge petascale andexadscale calculations being attempted.
I High-precision arithmetic may alleviate some numerical problems. Other
applications require tens, hundreds or even thousands of digits.
I All studies involving computations need to present evidence of internal checks or
replications to ensure that the results are reliable.
The credibility of high-performance scientific computing is at stake. Have welearned the lessons of history?