Code supplement to "Computing Science: Pencil, Paper, and Pi," American Scientist 102(5):342-345.
Brian Hayes. brian@bit-player.org
Python 2.7.6
# get floating-point division, decimal fractions, and regular expressions
from __future__ import division
import decimal
from decimal import Decimal as D # thus we have D('3.14') -> Decimal(3.14)
import re
# set the decimal precision globally ; can override locally in a 'with' statement
decimal.getcontext().prec = 709
Entered as strings of decimal digits.
# the correct value of arctan(1/5), to 709 decimal places,
# calculated with this Sage command: n(atan(1/5), digits=709)
atan_5_true = ' 0.1973955598498807583700497651947902934475851037878521015176889402410339\
699782437857326978280372880441126281180736913601044564798867942393557475\
654952163032700522107470015645015560061286185526633257318692806643896806\
189528405825931124251613297313993397113233537821796084176648310525473039\
665725650488878155309384290579311695934192851806364919697519401708560949\
527368673738508400812367856158009329822514023246675549211026704574378815\
474839079978985020075223696837961392278354193255722328413846477441352909\
705465122438302697560518377574220877835853152464749330914587633823112490\
332030126805100670223312575050942448460267162254894079226140467995062365\
969282873058287872053603034570766066681374312566267431408992606'
# Shanks's value of arctan(1/5).
# transcribed by hand from ProcRoySoc (1873) Vol. 21 p. 319, with the 75th decimal
# place corrected from 7 to 8, per ProcRoySoc (1873) Vol. 22 p. 45
atan_5_shanks = ' 0.1973955598498807583700497651947902934475851037878521015176889402410339\
699782437857326978280372880441126281180736913601044564798867942393557475\
654952163032700522107470015645015560061286185526633257318692806643896806\
189528405825931124251613297313993397113233537821796084176648310525473039\
665725650488878155309384290579311695934192851806364919697519401708560949\
527368673738508400812367856158009329822514023246675549211026704574378815\
474839079978985020075223696837961392278354193255722328413846477441352909\
705465122438302697560518377617781642423378303370181926488028277686119150\
985606759012135985563630343210056649978267636088711523275661084900937733\
802319504706876572938513592431975937947360575063620935078532833'
# Shanks's June 1853 value of arctan(1/5), to 609 places,
# transcribed by hand from _Rectification of the Circle_.
atan_5_shanks_609 = ' 0.\
1973955598498807583700497651947902934475851037878\
5210151768894024103396997824378573269782803728804\
4112628118073691360104456479886794239355747565495\
2163032700522107470015645015560061286185526633257\
3186928066438968061895284058259311242516132973139\
9339711323353782179608417664831052547303966572565\
0488878155309384290579311695934192851806364919697\
5194017085609495273686737385084008123678561580093\
2982251402324667554921102670457437881547483907997\
8985020075223696837961392278354193255722328413846\
4774413529097054651224383026975605183776177816424\
2337830337018192648802827768611915098560675901213\
598556363034347839926'
# the correct value of arctan(1/239), to 709 decimal places
# from the Sage command: n(atan(1/239), digits=707)
# (Note that the 'digits' parameter in Sage is number of significant digits;
# I have made an adjustment to get the right number of decimal places.)
atan_239_true = ' 0.0041840760020747238645382149592854527410480653076319508270196128871817\
783414228932737826058136229094549754506664448637560524583947893118650589\
221288330928008462719623307733759476346033184734145703319860154548148059\
924498302114603912539495276077968815588812733978533465180457425481358674\
644751979102328309770020646528276346532969104818386543560789195914512322\
209446368627661552083167964264657465511032510343526282445126935567049968\
444524790433177283930708631401938695195037058641077085585540452235538814\
237677083651569182527020022930895449500435854409344964401424187249509228\
386239545533356511719737470202349475977909746950111888547667397957315370\
930327821130898425830836771909100839098516551041922416780920533'
# Shanks's value of arctan(1/239)
# transcribed by hand from ProcRoySoc (1873) Vol. 21 p. 319,
atan_239_shanks = ' 0.0041840760020747238645382149592854527410480653076319508270196128871817\
783414228932737826058136229094549754506664448637560524583947893118650589\
221288330928008462719623307733759476346033184734145703319860154548148059\
924498302114603912539495276077968815588812733978533465180457425481358674\
644751979102328309770020646528276346532969104818386543560789195914512322\
209446368627661552083167964264657465511032510343526282445126935567049968\
444524790433177283930708631401938695195037058641077085585540452235538814\
237677083651569182527020022930895449500435854409344964401424187249509228\
386239545533356516494212200685238821940064584929313238867346764889187318\
168283021211013789711546961918469218237339030468204140799856684'
# Shanks's June 1853 value of arctan(1/239), 609 digits,
# transcribed by hand from _Rectification of the Circle_.
atan_239_shanks_609 = ' 0.\
0041840760020747238645382149592854527410480653076\
3195082701961288718177834142289327378260581362290\
9454975450666444863756052458394789311865058922128\
8330928008462719623307733759476346033184734145703\
3198601545481480599244983021146039125394952760779\
6881558881273397853346518045742548135867464475197\
9102328309770020646528276346532969104818386543560\
7891959145123222094463686276615520831679642646574\
6551103251034352628244512693556704996844452479043\
3177283930708631401938695195037058641077085585540\
4522355388142376770836515691825270200229308954495\
0043585440934496440142418724950922838623954553335\
651649421220068523882'
# the correct value of pi, with 707 decimal places
# from Sage: n(pi, digits=708)
pi_true = ' 3.1415926535897932384626433832795028841971693993751058209749445923078164\
062862089986280348253421170679821480865132823066470938446095505822317253\
594081284811174502841027019385211055596446229489549303819644288109756659\
334461284756482337867831652712019091456485669234603486104543266482133936\
072602491412737245870066063155881748815209209628292540917153643678925903\
600113305305488204665213841469519415116094330572703657595919530921861173\
819326117931051185480744623799627495673518857527248912279381830119491298\
336733624406566430860213949463952247371907021798609437027705392171762931\
767523846748184676694051320005681271452635608277857713427577896091736371\
7872146844090122495343014654958537105079227968925892354201996'
# Shanks's value of pi from his April 1873 paper, the first to proclaim 707 digits.
# Transcribed by hand from ProcRoySoc (1873) Vol. 21 p. 319.
pi_shanks_1873a = ' 3.1415926535897932384626433832795028841971693993751058209749445923078164\
062862089986280348253421170679821480865132823066470938446095505822317253\
594081284811174502841027019385211055596446229489549303819644288109756659\
334461284756482337867831652712019091456485669234603486104543266482133936\
072602491412737245870066063155881748815209209628292540917153643678925903\
600113305305488204665213841469519415116094330572703657595919530921861173\
819326117931051185480744623798347495673518857527248912279381830119491298\
336733624419366430860213950160924480772309436285530966202755693979869502\
224749962060749703041236688619951100892023837702131416941190298858254468\
1639799904659700081700296312377381342084130791451183980570985'
# Shanks's final value of pi, eight months later, correcting minor errors in the April
# value but introducing one further typo.
# transcribed by hand from ProcRoySoc (1873) Vol. 22 p. 45
pi_shanks_1873b = ' 3.1415926535897932384626433832795028841971693993751058209749445923078164\
062862089986280348253421170679821480865132823066470938446095505822317253\
594081284811174502841027019385211055596446229489549303819644288109756659\
334461284756482337867831652712019091456485669234603486104543266482133936\
072602491412737245870066063155881748815309209628292540917153643678925903\
600113305305488204665213841469519415116094330572703657595919530921861173\
819326117931051185480744623799627495673518857527248912279381830119491298\
336733624406566430860213950160924480772309436285530966202755693979869502\
224749962060749703041236688619951100892023837702131416941190298858254468\
1639799904659700081700296312377387342084130791451183980570985'
# Shanks's "best" value of pi. He never actually published this number.
# It's the version of 1873b with the typo at position 327 corrected
# (309209 --> 209209). This is the value I use in most of what follows.
pi_shanks = ' 3.1415926535897932384626433832795028841971693993751058209749445923078164\
062862089986280348253421170679821480865132823066470938446095505822317253\
594081284811174502841027019385211055596446229489549303819644288109756659\
334461284756482337867831652712019091456485669234603486104543266482133936\
072602491412737245870066063155881748815209209628292540917153643678925903\
600113305305488204665213841469519415116094330572703657595919530921861173\
819326117931051185480744623799627495673518857527248912279381830119491298\
336733624406566430860213950160924480772309436285530966202755693979869502\
224749962060749703041236688619951100892023837702131416941190298858254468\
1639799904659700081700296312377387342084130791451183980570985'
The aim here is to visually present pairs of large numbers in a way that makes discrepancies conspicuous. The numbers are easier to read if we break them into groups of digits. Then digits that differ are highlighted by color.
First off, a bit of a false start. I thought I would do the comparison of two numeric strings by constructing a third string that records digit-by-digit differences. Later I decided to take another approach, but I'm leaving the code here since the difference strings might be useful for other purposes.
def digit_diff(a, b):
"""For digit a and b, return the digit that when
added to a modulo 10 produces b."""
if a.isdigit() and b.isdigit():
if b < a:
b = '1' + b # concatenation not addition!
return str(int(b) - int(a))
elif a == b == '.':
return '.'
else:
return '?'
def numstr_diffs(astr, bstr):
"""Return a string of digits representing the digit-by-digit
differences between digit strings bstr and astr."""
return(''.join([digit_diff(a, b) for a, b in zip(astr, bstr)]))
def flag_diffs(diffstr):
"""Return a list of tuples giving start and end indices for
regions where the diffstr is nonzero."""
return [tpl.span() for tpl in re.finditer('[1-9]+', diffstr)]
Present two long strings of digits in a tidily formatted table that makes it easy to see where they differ.
# We're going to re-encode discrepant digits in a different alphabet
# Seems a little cheesy, but it's a way of flagging erroneous digits
# in a pure-ascii string without disturbing character counts or indexing.
from IPython.display import HTML
from string import maketrans
n_to_a = maketrans('0123456789', 'abcdefghij')
a_to_n = maketrans('abcdefghij', '0123456789')
# Here we do the re-encoding...
def flag_discrepancies(astr, bstr):
"""Return a new version of bstr in which every digit that differs
from the corresponding digit in astr has been re-encoded by mapping
the alphabet [0-9] onto [a-j]."""
diffstr = ''
for a, b in zip(astr, bstr):
if a == b:
diffstr += b
else:
diffstr += b.translate(n_to_a)
return(diffstr)
# an example of what the re-encoding routine produces
flag_discrepancies(pi_true, pi_shanks)
# ... and here we undo it, restoring the normal digits but
# surrounding each block with a styled span tag
def color_error_digits(flagstr):
"""Take a string of digits with discrepancies marked by alphabetic encoding;
add HTML background colors to any such alphabetic spans; finally restore
those digits to normal decimal encoding."""
span_tag = '<span style="background-color:orange;">'
parts = re.split('([a-j]+)', flagstr)
html = ''
for item in parts:
if item.isalpha():
html += span_tag + item.translate(a_to_n) + '</span>'
else:
html += item
return(html)
# Construct and print the styled comparison table, and display
# in the cell as HTML.
def print_diff_table(astr, bstr, digits_per_line=60, digits_per_group=5):
"""Print the two digit strings in alternating lines, broken
into groups, with discrepancies in the second string highlighted
in color."""
astr = astr.partition('.')[2] # discard integer part and decimal point
bstr = bstr.partition('.')[2]
bstr = flag_discrepancies(astr, bstr) # do the re-encoding of error digits
alst = re.findall('.{1,' + str(digits_per_line) + '}', astr) # break up into lines
blst = re.findall('.{1,' + str(digits_per_line) + '}', bstr)
alst = map(lambda(s): re.findall('.{1,' + str(digits_per_group) + '}', s), alst) # groups of digits
blst = map(lambda(s): re.findall('.{1,' + str(digits_per_group) + '}', s), blst)
alst = map(lambda(lst): ' '.join(lst), alst) # convert lists of groups in string with spaces
blst = map(lambda(lst): ' '.join(lst), blst)
blst = map(color_error_digits, blst) # remove alpha encoding, replace with styled spans
tbl = """"""
tbl += '<pre style="font-size:10pt; line-height:10pt">'
for aline, bline in zip(alst, blst):
tbl += aline + '\n'
tbl += bline + '\n'
tbl += '\n'
tbl += '</pre>'
# print tbl # uncomment if you want the raw HTML, perhaps to export
return(HTML(tbl))
Here we have digit-by-digit comparisons for the true values and Shanks's values of atan(1/5), atan(1/239) and pi. The upper number in each row is the correct value; the lower one is Shanks's value. Digits that differ are flagged in orange.
print_diff_table(atan_5_true, atan_5_shanks)
print_diff_table(atan_239_true, atan_239_shanks)
print_diff_table(pi_true, pi_shanks_1873a)
# The first two isolated errors (962 --> 834 and 065 --> 193) did not appear in Shanks's publications of the 1850s.
print_diff_table(pi_true, pi_shanks_1873b)
# In his followup publication of December 1873, Shanks corrected the two three-digit error blocks from the April paper.
# Unfortunately, he (or the typesetter) introduced a new single-digit error.
print_diff_table(pi_true, pi_shanks)
# Finally we have the best version Shanks might have achieved, cherry-picking all the correct digits he ever published.
# The erroneous tail of the value begins with decimal place 528.
# Machin's formula: given atan 1/5 and atan 1/239, construct pi
# args can be either Decimals or strings
def make_pi(a5, a239):
a5 = D(a5)
a239 = D(a239)
return(16 * a5 - 4 * a239)
make_pi(atan_5_true, atan_239_true)
In his calculation of $\pi$ Shanks used the Machin series:
$$\frac{\pi}{4} = 4 \arctan \left( \frac{1}{5} \right) - \arctan \left( \frac{1}{239} \right)$$
And he computed the arctangents by summing terms of this series:
$$\arctan(x) = \frac{x^{1}}{1} - \frac{x^{3}}{3} + \frac{x^{5}}{5} - \frac{x^{7}}{7} + \cdots = \sum_{k=0}^\infty \frac{(-1)^{k}}{2k+1}x^{2k+1}$$
With a modern computer, the series sum for arctangent is just a matter of directly translating the mathematical formula into a suitable programming language.
# This one loops through the terms of the series, starting from scratch each time
# to compute the sign as (-1)^k, the exponent and denominator m as 2k+1, the power
# x^m and finally the quotient (x^m)/m. The result is added to a running sum s, and
# the loop continues until the terms are too small to contibute to the sum at the
# specified precision.
# Note the use of the Python decimal module here: D(a) is decimal.Decimal(a), and
# ctx.prec = p is locally binding the decimal precision to p places.
def arctan_silicon(a, b, p):
"""Return p decimal digits of arctan of a/b."""
with decimal.localcontext() as ctx:
ctx.prec = p
x = D(a) / D(b)
assert abs(x) <= 1
min_term = D(10) ** -p
s = D('0')
k = 0
term = x
while abs(term) > min_term:
m = D((2 * k) + 1)
term = ((-1)**k / m) * (x**m)
s += term
k += 1
return(s)
arctan_silicon(1, 5, 709)
# The value calculated by arctan_silicon differs from the canonical Sage value
# in the last decimal place, because of truncation error. But the 707-place prefix
# matches correctly:
'{:.709f}'.format(arctan_silicon(1, 5, 709))[:709] == atan_5_true[:709]
# It's also handy to have a procedure
# for generating a single term of an arctan series.
def arctan_term(a, b, p, k):
"""Return p digits of term k of the series expansion for arctan(a/b)."""
with decimal.localcontext() as ctx:
ctx.prec = p
x = D(a) / D(b)
m = D(k + k + 1)
return((-1)**k * (x**m) / m)
'{: .709f}'.format(arctan_term(1, 239, 709, 2))
No human computer would use a method anything like the one above for a paper-and-pencil calculation. We do not keep track of alternating signs by raising -1 to the kth power. And if you have just computed $x^{25}$ and you now need $x^{27}$, you don't start all over; you multiply $x^{25} \times x^{2}$. Here's a first attempt at a people-friendly arctan algorithm.
# A naive algorithm for computing arctan on paper. Each time through the loop we
# update a running value of m, x^m and (x^m)/m; we also toggle the sign.
def arctan_pencil(a, b, p):
"""Return p decimal digits of arctan(a/b)."""
with decimal.localcontext() as ctx:
ctx.prec = p
x = D(a) / D(b)
assert abs(x) <= 1
min_term = D(10) ** -p
xsq = x * x # x^2 multiplier for forming odd powers of x
s = D('0') # above here are constants; below are loop variables
k = 0 # not actually needed here, but convenient term counter
sign = 1 # toggled plus and minus
m = 1 # = 2k+1; each term is (x^m)/m
x_to_m = x # current power of x
term = x # current term of series
while abs(term) > min_term:
s += term
k += 1
sign *= -1
m += 2
x_to_m *= xsq
term = (sign * x_to_m) / m
return(s)
# The silicon and pencil-lead procedures should give the same result.
arctan_pencil(1, 5, 709) == arctan_silicon(1, 5, 709)
There's still something inhuman about the algorithm above. Each time through the loop we overwrite all the variables, keeping only the latest value of m, x^m, and so on. Overwriting is not easy to do on paper, and it would be ill-advised anyway. A human computer is more likely to keep all the intermediate results in a table.
# A version of the arctan algorithm that prints out a table of intermediate
# results.
def arctan_pencil_table(a, b, p):
"""Return p decimal digits of arctan(a/b), printing out
a table of intermediate values."""
with decimal.localcontext() as ctx:
ctx.prec = p
x = D(a) / D(b)
assert abs(x) <= 1
min_term = D('10') ** -p
xsq = x * x
s = D('0')
number_fmt = '{: >3}{: >4}{: >4} {: <.' + `p` + 'f} {: < .' + `p` + 'f}'
header_fmt = '{: >3}{: >4}{: >4} {: ^' + `p` + '} {: ^' + `p` + '}'
k = 0
sign = 1
m = 1
x_to_m = x
term = x
print header_fmt.format('k', 'sgn', 'm', 'x^m', '(x^m)/m')
while abs(term) > min_term:
print number_fmt.format(k, sign, m, x_to_m, term)
s += term
k += 1
sign *= -1
m += 2
x_to_m *= xsq
term = (sign * x_to_m) / m
return(s)
# If we were to continue this to 700 digits, the table would get somewhat unruly!
arctan_pencil_table(1, 5, 40)
arctan_pencil_table(1, 239, 40)
In the series calculation for $\arctan(x)$, we're guaranteed that $-1 \le x \le 1$; in all the examples considered here, $x$ can be expressed as a fraction with the form $1/z$. So far, we've been evaluating these terms as
$$\frac{\left( \frac{1}{z} \right)^{m}}{m}.$$
But instead we could write the expression as
$$\frac{1}{m \cdot z^{m}}.$$
Mathematically, the two expressions are identical, but computationally they're quite different. Would you rather divide integers into tiny numbers, making them even tinier, or would you prefer to multiply large numbers, making them still larger? In the latter case, you've still got to do a division at the end, but it's a reciprocal: You're dividing into 1.0.
# A version of the arctan algorithm that does more multiplying
# and less dividing; optionally prints out a table of intermediate
# results.
def arctan_pencil_recip(b, p, table=False):
"""Return p decimal digits of arctan(1/b). If table=True, print out
a table of intermediate values."""
with decimal.localcontext() as ctx:
ctx.prec = p
min_term = D('10') ** -p
bsq = b * b
s = D('0')
number_fmt = '{: >3}{: >4}{: >4} {:< ' + `p` + 'd} {:< ' + `p` + 'd} {:< .' + `p` + 'f}'
header_fmt = '{: >3}{: >4}{: >4} {: ^' + `p` + '} {: ^' + `p` + '} {: ^' + `p` + '}'
k = 0
sign = 1
m = 1
b_to_m = b
m_times_b_to_m = b
term = D(1) / b
if table:
print header_fmt.format('k', 'sgn', 'm', 'b^m', 'm(b^m)', '1/m(b^m)')
while abs(term) > min_term:
if table:
print number_fmt.format(k, sign, m, b_to_m, m_times_b_to_m, term)
s += term
k += 1
sign *= -1
m += 2
b_to_m *= bsq
m_times_b_to_m = m * b_to_m
term = D(1) / (sign * m_times_b_to_m)
return(s)
arctan_pencil_recip(5, 25, table=True)
# Just checking...
arctan_pencil_recip(5, 709) == arctan_silicon(1, 5, 709)
How and where did Shanks stray into error?
# Rumor says that Shanks might have omitted a term while adding up
# the hundreds of terms in an arctan series. Let's also check for
# including a term twice.
def search_omitted__or_duped_terms(a, b, p, bogus_atan, match_range_tuple):
"""
In the series for arctan(a/b) at precision p decimal places, try
deleting each term and duplicating each term; if the resulting digits
in the range specified as match_range_tuple matches the corresponding
digits in bogus_atan, report success.
"""
with decimal.localcontext() as ctx:
ctx.prec = p
fmt = '{:.' + str(p) + 'f}'
true_atan = arctan_silicon(a, b, p)
left, right = match_range_tuple
frame_left = max(0, left-20)
frame_right = min(p, right+20)
match_str = bogus_atan[left:right]
min_term = D(10) ** -p
k = 0
term = arctan_term(a, b, p, k)
print 'match string = {}'.format(match_str)
while abs(term) > min_term:
adj_atan_minus = true_atan - term
adj_atan_plus = true_atan + term
if fmt.format(adj_atan_minus)[left:right] == match_str:
print 'k = {}, possible omitted term'.format(k)
print 'term = ' + fmt.format(term)[frame_left:frame_right]
print 'bogus = ' + fmt.format(D(bogus_atan))[frame_left:frame_right]
print 'true = ' + fmt.format(true_atan)[frame_left:frame_right]
print 'adj = ' + fmt.format(adj_atan_minus)[frame_left:frame_right]
if fmt.format(adj_atan_plus)[left:right] == match_str:
print 'k = {}, possible duplicated term'.format(k)
print 'term = ' + fmt.format(term)[frame_left:frame_right]
print 'bogus = ' + fmt.format(D(bogus_atan))[frame_left:frame_right]
print 'true = ' + fmt.format(true_atan)[frame_left:frame_right]
print 'adj = ' + fmt.format(adj_atan_plus)[frame_left:frame_right]
k += 1
term = arctan_term(a, b, p, k)
print 'Searched {} omissions and duplicates.'.format(k+1)
Go fishing for omitted or duplicated terms. None to be found in the example below.
search_omitted__or_duped_terms(1, 5, 709, atan_5_shanks, (532,541))
search_omitted__or_duped_terms(1, 5, 709, atan_5_shanks_609, (603,608))
search_omitted__or_duped_terms(1, 239, 709, atan_239_shanks, (594,600))
# One plausible error mechanism is omitting a digit or a group of digits
# while copying numbers from one sheet to another. Our digits might be
# mistakenly duplicated while copying. Errors of this kind produce a
# characteristic signature in the affected term: Beyond the point of error,
# all digits are shifted left or right.
def search_shifts(a, b, p, bogus_atan, match_range_tuple, display_range_tuple):
"""
Compute arctan(a/b), then subtract bogus_atan to get a difference delta;
then, for each term in the arctan series, subtract delta,
extract a pattern from resulting adjusted term, and look for a matching
pattern at any position within the true term.
"""
with decimal.localcontext() as ctx:
ctx.prec = p
fmt = '{: .' + `p` + 'f}'
true_atan_d = arctan_silicon(a, b, p)
bogus_atan_d = D(bogus_atan)
atan_diff_d = true_atan_d - bogus_atan_d
start, end = match_range_tuple
left, right = display_range_tuple
min_term = D(10) ** -p
k = 0
term = arctan_term(a, b, p, k)
while abs(term) > min_term:
adj_term = term - atan_diff_d
tru_term_str = fmt.format(term)
diff_term_str = fmt.format(atan_diff_d)
adj_term_str = fmt.format(adj_term)
pattern = adj_term_str[start:end]
position = tru_term_str.find(pattern)
if position != -1:
print 'k = {}, m = {}, pattern = {}'.format(k, k+k+1, pattern)
print 'match at position {}, shifted by {}'.format(position, start - position)
print 'true = ', tru_term_str[left:right]
print 'adj = ', adj_term_str[left:right]
print 'diff = ', diff_term_str[left:right]
k += 1
term = arctan_term(a, b, p, k)
print 'Searched for shifts in {} terms.'.format(k+1)
search_shifts(1, 5, 709, atan_5_shanks, (537, 543), (521, 601))
# finds a shift by -1 in term 248
Now we have a candidate error in the computation of arctan(1/5). In term 248, at decimal place 530 (or string index 533, because of the " 0." prefix), it looks like a zero digit was omitted. There's also a substitution 2 --> 3 two digits later. If this really is the source of Shanks's error, then by injecting the change into the computation of arctan(1/5), we should get the erroneous value that Shanks reported.
# Inject the "uncorrection" into term 248
def uncorrect_term_248_at_533_and_536():
fmt = '{: .709f}'
t248 = fmt.format(arctan_term(1, 5, 709, 248))
t248_ucx = t248[:533] + '483' + t248[537:] + '0'
return(t248_ucx)
uncorrect_term_248_at_533_and_536()
# store the uncorrection in a dictionary
# k = 248
# 4446680080482897 --> 444668008483897
# at end, append 3
atan_5_error_dict_248 = {248: ' 0.000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000823287376231424011572733827169919303\
87742506188205854866279963986321272545389199066480636124896985843835178172222335789\
9702529312345420788060494866543259557344064386317907444668008483897384305835010060\
36217303822937625754527162977867203219315895372233400402414486921529175050301810865\
191146881287726358148893360160965794768611670020120724346076458752515090543'}
# Now create a new version of the arctan function that will inject
# any uncorrections in the appropriate places.
def inject_error_terms(a, b, p, error_dict):
"""Return arctan(a/b) with termsubstitutions specified in error_dict."""
x = D(a) / D(b)
fmt = '{: .' + `p` + 'f}'
with decimal.localcontext() as ctx:
ctx.prec = p
min_term = D('10') ** -p
xsq = x * x
s = D('0')
k = 0
sign = 1
m = 1
x_to_m = x
term = x
while abs(term) > min_term:
substitute_term = error_dict.get(k)
if substitute_term:
# print k
# print fmt.format(term)[500:550]
# print substitute_term[500:550]
s += D(substitute_term)
else:
s += term
k += 1
sign *= -1
m += 2
x_to_m *= xsq
term = (sign * x_to_m) / m
return(fmt.format(s))
# give it a try....
inject_error_terms(1, 5, 709, atan_5_error_dict_248)
# And now let's see if the 'uncorrection' transforms the true
# arctan value into the Shanks value
# First, for comparison, here again is the difference table comparing
# the true arctan(1/5) with the Shanks result
print_diff_table(atan_5_true, atan_5_shanks)
# And now compare the 'uncorrected' true value with Shanks
print_diff_table(inject_error_terms(1, 5, 709, atan_5_error_dict_248), atan_5_shanks)
Progress, sort of: We have uncorrected all the digits from decimal place 530 through decimal place 568. But starting at place 569, Shanks departs from the uncorrected version. Evidently there's another error.
In 1946 D. F. Ferguson suggested an explanation for the error at decimal place 569: Shanks had truncated term 72, effectively zeroing out all its digits starting at decimal place 569. We can test this idea.
# Set all digits starting at decimal place 569 (string index 572) to 0
def uncorrect_term_72_zeros_569():
fmt = '{: .709f}'
t72 = fmt.format(arctan_term(1, 5, 709, 72))
t72_ucx = t72[:571] + '0'*(709-571)
return(t72_ucx)
uncorrect_term_72_zeros_569()
# add the uncorrection to our error dictionary
# k = 248
# 4446680080482897 --> 444668008483897
# at end, append 3
# k = 72
# [:571] --> 0
atan_5_error_dict_248_and_72_zeroes = {248: ' 0.000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000823287376231424011572733827169919303\
87742506188205854866279963986321272545389199066480636124896985843835178172222335789\
9702529312345420788060494866543259557344064386317907444668008483897384305835010060\
36217303822937625754527162977867203219315895372233400402414486921529175050301810865\
191146881287726358148893360160965794768611670020120724346076458752515090543',
72: ' 0.000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000307596485496974112297044389967563\
6069790402206896551724137931034482758620689655172413793103448275862068965517241379310\
3448275862068965517241379310344827586206896551724137931034482758620689655172413793103\
4482758620689655172413793103448275862068965517241379310344827586206896551724137931034\
4827586206896551724137931034482758620689655172413793103448275862068965517241379310344\
8275862068965517241379310344827586206896551724137931034482758620689655172413793103448\
2758620000000000000000000000000000000000000000000000000000000000000000000000000000000\
000000000000000000000000000000000000000000000000000000000000'}
# try applying the two uncorrections
inject_error_terms(1, 5, 709, atan_5_error_dict_248_and_72_zeroes)
# And now see if we've gotten closer to Shanks
print_diff_table(inject_error_terms(1, 5, 709, atan_5_error_dict_248_and_72_zeroes), atan_5_shanks)
Looks no better, does it? We are no closer to reproducing Shanks's value. As a matter of fact, we're one digit farther. So now the forensic mathematician has two mysteries to solve: What Shanks did wrong, and how Ferguson came up with a fix that doesn't fix anything. I think I have a clue about the first question, but the second stumps me.
For starters, let's look at the various pairs of numbers involved, and the differences between them. I'm displaying the sign of each number, the first few digits, and then the range of digits where all the action is, from decimal place 520 to 610. (The code for producing these numbers is in the cell below.)
pi true 3.141 ... 430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681 pi shanks 3.141 ... 430860213950160924480772309436285530966202755693979869502224749962060749703041236688619951 pi diff -0.000 ... 000000000000696972233400402414486921529175050301808106570457226115312565026347185368614269
a5 true 0.197 ... 697560518377574220877835853152464749330914587633823112490332030126805100670223312575050942 a5 shanks 0.197 ... 697560518377617781642423378303370181926488028277686119150985606759012135985563630343210056 diff0 -0.000 ... 000000000000043560764587525150905432595573440643863006660653576632207035315340317768159114
The numbers above are what we know at the outset, because Shanks published his values for pi and for arctan 1/5, and because we can calculate the true values. In the case of arctan 1/5, the difference diff0 tells us how we need to "uncorrect" the true value in order to arrive at Shanks's value. And, as we've already seen, there's a simple and fairly plausible change to term 248 in the arctan 1/5 summation that gets us partway there. After making that uncorrection, we're left with this residual unerror in the arctan 1/5 value:
a5 248ucx 0.197 ... 697560518377617781642423378303370181926488028277686291564778710207287998054529147585111304 a5 shanks 0.197 ... 697560518377617781642423378303370181926488028277686119150985606759012135985563630343210056 diff1 0.000 ... 000000000000000000000000000000000000000000000000000172413793103448275862068965517241901247
So now the task is to find some plausible change that will have the effect of subtracting diff1 from one of the terms of the arctan 1/5 expansion. Let's look at term 72, the one that Ferguson mentioned. Here I show the correct value of term 72, the "uncorrected" version after subtracting diff1, and, for convenience I have repeated the value of diff1 itself. The orange bars show one cycle of 28 digits in the repeating decimal expansion of term 72, beginning at decimal place 569, and the corresponding ranges in the other two numbers.
t72 0.000 ... 517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241 t72-diff1 0.000 ... 517241379310344827586206896551724137931034482758620517241379310344827586206896551723615993 diff1 0.000 ... 000000000000000000000000000000000000000000000000000172413793103448275862068965517241901247
What catches the eye here is the repetition of the same sequence of digits, with different offsets or phases, in all three orange bands. What's that all about? The answer turns out to be interesting and important, and I'll return to it below; but for now I want to focus solely on what the pattern implies in our attempt to understand the errors in Shanks's computation. We have another hypothetical "uncorrection" here that transforms the true arctan value into something closer to Shanks's result. In the earlier uncorrection of term 248, we assumed that Shanks skipped over a single digit while recording or transcribing the term. In this case, the hypothesis requires that he skip five digits. In term 72, omit the first five digits inside the orange band (68965), and t72 will exactly match t72-diff1 throughout the orange region and for another four digits.
Skip over the next cell (which generates the numbers displayed above) to see how this works out.
def compare_atan_terms(n, left, right):
fmt = '{: .709f}'
diff_pi = fmt.format(D(pi_true) - D(pi_shanks))
diff0 = fmt.format(D(atan_5_true) - D(atan_5_shanks))
atan_5_248ucx = inject_error_terms(1, 5, 709, atan_5_error_dict_248)
diff1 = fmt.format(D(atan_5_248ucx) - D(atan_5_shanks))
tn = fmt.format(arctan_term(1, 5, 709, n))
tn_minus_diff1 = fmt.format(D(tn) - D(diff1))
print
print 'pi true {} ... {}'.format(pi_true[0:6], pi_true[left:right])
print 'pi shanks {} ... {}'.format(pi_shanks[0:6], pi_shanks[left:right])
print 'pi diff {} ... {}'.format(diff_pi[0:6], diff_pi[left:right])
print
print 'a5 true {} ... {}'.format(atan_5_true[0:6], atan_5_true[left:right])
print 'a5 shanks {} ... {}'.format(atan_5_shanks[0:6], atan_5_shanks[left:right])
print 'diff0 {} ... {}'.format(diff0[0:6], diff0[left:right])
print
print 'a5 248ucx {} ... {}'.format(atan_5_248ucx[0:6], atan_5_248ucx[left:right])
print 'a5 shanks {} ... {}'.format(atan_5_shanks[0:6], atan_5_shanks[left:right])
print 'diff1 {} ... {}'.format(diff1[0:6], diff1[left:right])
print
print 'tn {} ... {}'.format(tn[0:6], tn[left:right])
print 'tn-diff1 {} ... {}'.format(tn_minus_diff1[0:6], tn_minus_diff1[left:right])
print 'diff1 {} ... {}'.format(diff1[0:6], diff1[left:right])
print
# try n = 14, 72; 362 is part of the same family, but the cyclic part of the decimal fraction
# doesn't begin until after digit 709
compare_atan_terms(14, 520, 612)
# Hypothesis: Shanks's second error was omitting five digits in term 72
# while extending the value beyond decimal place 569 (string index 572)
def uncorrect_term_72_skip_5():
fmt = '{: .709f}'
t72 = fmt.format(arctan_term(1, 5, 709, 72))
t72_ucx = t72[:571] + t72[576:-1] + '689655'
return(t72_ucx)
uncorrect_term_72_skip_5()
# prepare an error dictionary with the five-digit skip uncorrection
# k = 248
# 4446680080482897 --> 444668008483897
# at end, append 3
# k = 72
# remove [571:576]
# at end, complete the 28-digit period
atan_5_error_dict_248_and_72_skip = {248: ' 0.000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000823287376231424011572733827169919303\
87742506188205854866279963986321272545389199066480636124896985843835178172222335789\
9702529312345420788060494866543259557344064386317907444668008483897384305835010060\
36217303822937625754527162977867203219315895372233400402414486921529175050301810865\
191146881287726358148893360160965794768611670020120724346076458752515090543',
72: ' 0.000000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000003075964854969741122970443899675636069790402206896551724\
13793103448275862068965517241379310344827586206896551724137931034482758620689655172\
41379310344827586206896551724137931034482758620689655172413793103448275862068965517\
24137931034482758620689655172413793103448275862068965517241379310344827586206896551\
724137931034482758620689655172413793103448275862068965517241379310344827586206896551\
724137931034482758620689655172413793103448275862068965517241379310344827586205172413\
793103448275862068965517241379310344827586206896551724137931034482758620689655172413\
79310344827586206896551724137931034482758620689655'}
inject_error_terms(1, 5, 709, atan_5_error_dict_248_and_72_skip)
# Now test our "uncorrection" to see how close we've gotten to the Shanks
# value of arctan 1/5
# We're good now through decimal place 601.
print_diff_table(inject_error_terms(1, 5, 709, atan_5_error_dict_248_and_72_skip), atan_5_shanks)
We have found two potential errors, in terms 248 and 72, that might explain some of the problems Shanks encountered in extending his calculations beyond 530 decimal places. Starting with the correct values for the terms in the series for arctan 1/5, skipping a single digit in term 248 uncorrects decimal places 530 through 568, making them match the digits given by Shanks. Then a five-digit omission in term 72 uncorrects decimal places 569 through 601. Beyond that point, there appear to be further unidentified errors in the arctan 1/5 series; there's also a mistake at decimal place 592 in the Shanks value for arctan 1/239.
The narrative above reflects my state of knowledge when the American Scientist column went to press in mid-July 2014. In the next few weeks, some loose ends in the story continued to nag at me.
First, skipping one digit while transcribing a huge sequence of digits seems like a typical human failure mode; skipping a block of five digits seems less likely. Could there be some other explanation of the error at decimal place 569?
Second, how do we explain Ferguson's faulty diagnosis in term 72? Here is what he writes on the subject at the end of his 1946 paper in The Mathematical Gazette:
In calculations made since the preceding part of this Note was written, from the 569th decimal place onwards I have a further disagreement with Shanks' figures, namely 6896551724137931034482758620 (recurring). These are precisely the figures which should occur at this exact spot (i.e. from the 569th d.p. onwards) in the term 145 which comes in the series for $\tan^{-1}(1/5)$. This surely cannot be a mere coincidence. It looks very much as if Shanks omitted the whole of this term from the 569th place onwards.
In this passage Ferguson refers to term 145, but that's the same term I've been calling 72. (He counts by $m$, I count by $k$, where $m=2k+1$.) The sequence of digits he cites does indeed occur exactly where he describes it. But removing those digits (or setting them all to 0) does not recover the Shanks value of arctan 1/5. Omitting the first five digits and shifting the rest to the left does accomplish what we want. Ferguson was a skillful and thorough worker, so this proposal of the wrong fix at the right place is puzzling.
Finally, the very digit sequence mentioned by Ferguson has already been mentioned above as an object of curiosity:
t72 0.000 ... 517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241 t72-diff1 0.000 ... 517241379310344827586206896551724137931034482758620517241379310344827586206896551723615993 diff1 0.000 ... 000000000000000000000000000000000000000000000000000172413793103448275862068965517241901247
In the orange bands, we take that digit sequence, subtract a cyclic permutation of the same sequence, and get a third cyclic permutation as the result. That's not something that happens with just any string of numbers!
# See for yourself. Try a random sequence of digits, subtract
# all possible cyclic permutations of the same sequence, and see
# if you're left with another cyclic permutation.
def rotate(s, n):
return(s[n:] + s[0:n])
def cyclic_diffs(s):
sss = s + s + s
for n in range(1, len(s)):
print '{: }'.format(int(sss) - int(rotate(sss, n)))[len(s)+1:2*len(s)+1]
# Nope, this one doesn't work.
cyclic_diffs('123456')
# Neither does this.
cyclic_diffs('910672')
# But this one is perfect.
cyclic_diffs('142857')
# And so does the 28-digit sequence plucked out of term 72.
cyclic_diffs('6896551724137931034482758620')
Obviously there's a trick here. The six-digit sequence 142857 is the recurrent pattern in the decimal fraction for 1/7. The 28-digit sequence 0344827586206896551724137931 (a cyclic permutation of 6896551724137931034482758620) is the recurrent decimal fraction for 1/29. Both 7 and 29 are full-repetend primes: primes $p$ with the property that the decimal expansion of $1/p$ has a period of $p-1$ digits; multiplying those digits by any number in the range $1, 2, 3, \ldots n-1$ produces a cyclic permutation of the same digits. In other words, it is a cyclic number,
This is a lovely bit of number theory, which goes back to Gauss (who else?) and still has open questions (see Artin's conjecture). But what does it have to do with Shanks, pi, and arctangents?
Term 72 in the series for arctangent 1/5 is $1/(145 \cdot 5^{145})$. The denominator has the prime factorization $29 \cdot 5^{146}$. In the decimal expansion of this fraction, the factors of 5 contribute a long nonrecurrent prefix; the repetend—the repeating part of the expansion—is determined entirely by the prime 29. This explains the spooky behavior of term 72—why shifting the sequence and subtracting yields the same sequence again.
Making this connection led me to look at another term in the arctan 1/5 series. Term 14 is $1/29 \cdot 5^{29}$, and its repetend must have the same sequence of digits as term 72 (in some cyclic permutation). Applying the same procedure to term 14 yields this result:
tn 0.000 ... 379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310 tn-diff1 0.000 ... 379310344827586206896551724137931034482758620689654999999999999999999999999999999999478062 diff1 0.000 ... 000000000000000000000000000000000000000000000000000172413793103448275862068965517241901247
This was an aha moment, suggesting a different explanation of the error at position 569, more along the lines of what Ferguson reported seeing.
# In term 14 set all digits starting at decimal place 569 (string index 572) to 0
def uncorrect_term_14_zeros_569():
fmt = '{: .709f}'
t14 = fmt.format(arctan_term(1, 5, 709, 14))
t14_ucx = t14[:571] + '0'*(709-571)
return(t14_ucx)
uncorrect_term_14_zeros_569()
# prepare an error dictionary with the term 14 omission uncorrection
# k = 248
# 4446680080482897 --> 444668008483897
# at end, append 3
# k = 14
# from 571 on, set digits to 0
atan_5_error_dict_248_and_14_zeroes = {248: ' 0.000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000823287376231424011572733827169919303\
87742506188205854866279963986321272545389199066480636124896985843835178172222335789\
9702529312345420788060494866543259557344064386317907444668008483897384305835010060\
36217303822937625754527162977867203219315895372233400402414486921529175050301810865\
191146881287726358148893360160965794768611670020120724346076458752515090543',
14: ' 0.000000000000000000000185127900689655172413793103448275862068965517241379310\
34482758620689655172413793103448275862068965517241379310344827586206896551724137931\
03448275862068965517241379310344827586206896551724137931034482758620689655172413793\
10344827586206896551724137931034482758620689655172413793103448275862068965517241379\
31034482758620689655172413793103448275862068965517241379310344827586206896551724137\
93103448275862068965517241379310344827586206896551724137931034482758620689655172413\
79310344827586206896551724137931034482758620689655172413793103448275862068965500000\
00000000000000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000'}
inject_error_terms(1, 5, 709, atan_5_error_dict_248_and_14_zeroes)
print_diff_table(inject_error_terms(1, 5, 709, atan_5_error_dict_248_and_14_zeroes), atan_5_shanks)
This uncorrection, blanking out term 14 beyond position 569, produces the same effect as a five-digit shift in term 72. Knowing something about the cyclic structure of the two terms, we can understand why those two quite different actions have the same result.
Could there be other terms with the same cyclic repetend? Yes, term 362 has a denominator that factors as $29 \cdot 5^{729}$, but that number is so large that the repetent does not appear until beyond the Shanks limit of 709 decimal places.
Given two alternative uncorrections with the same effect, which one is more likely to explain Shanks's misstep? Is it more plausible that he skipped over a few digits while calculating or copying a single term, or that he dropped all the trail digits of one term while summing up the whole series? In the absence of documentary evidence—unlikely to turn up 140 years later—the answer must remain a matter of conjecture. Of course there may well be additional candidates I haven't spotted at all. And there remain at least two more errors—one in the arctan 1/5 series, the other in arctan 1/239—that remain untirely unexplained. I also remain utterly muddled about Ferguson's suggested repair, which seems to work successfully, but when applied to a different term than the one he identified.