# William Shanks's computations of pi, 1850-1873¶

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

## Preliminaries¶

In [1]:
# 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


### The Numbers

Entered as strings of decimal digits.

In [2]:
# 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'

In [3]:
# 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'

In [4]:
# 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'

In [5]:
# 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'

In [6]:
# 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'

In [7]:
# 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'

In [8]:
# 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'

In [9]:
# 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'

In [10]:
# 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'

In [11]:
# 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'


### Comparison and display of numbers

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.

In [12]:
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 '?'

In [13]:
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)]))

In [14]:
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)]


#### Display a table of digit discrepancies

Present two long strings of digits in a tidily formatted table that makes it easy to see where they differ.

In [15]:
# 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')

In [16]:
# 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)

In [17]:
# an example of what the re-encoding routine produces
flag_discrepancies(pi_true, pi_shanks)

Out[17]:
' 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139fab6a9ceeiah7cd0jedgciffdajggcac7f5g9dj7ji6jfacccehejjgcagahejhadaebcdggiigbjjf1baaij2acdidhhacbdbe1gjebbjac9iificfeegibgdjhjjja4gfjhaaaibhaacjgdbcdhhdi7dec0iebdahjbefbbidjiafha9if'

In [18]:
# ... 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)

In [19]:
# 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))


#### Output of the table-building routine

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.

In [20]:
print_diff_table(atan_5_true, atan_5_shanks)

Out[20]:
19739 55598 49880 75837 00497 65194 79029 34475 85103 78785 21015 17688
19739 55598 49880 75837 00497 65194 79029 34475 85103 78785 21015 17688

94024 10339 69978 24378 57326 97828 03728 80441 12628 11807 36913 60104
94024 10339 69978 24378 57326 97828 03728 80441 12628 11807 36913 60104

45647 98867 94239 35574 75654 95216 30327 00522 10747 00156 45015 56006
45647 98867 94239 35574 75654 95216 30327 00522 10747 00156 45015 56006

12861 85526 63325 73186 92806 64389 68061 89528 40582 59311 24251 61329
12861 85526 63325 73186 92806 64389 68061 89528 40582 59311 24251 61329

73139 93397 11323 35378 21796 08417 66483 10525 47303 96657 25650 48887
73139 93397 11323 35378 21796 08417 66483 10525 47303 96657 25650 48887

81553 09384 29057 93116 95934 19285 18063 64919 69751 94017 08560 94952
81553 09384 29057 93116 95934 19285 18063 64919 69751 94017 08560 94952

73686 73738 50840 08123 67856 15800 93298 22514 02324 66755 49211 02670
73686 73738 50840 08123 67856 15800 93298 22514 02324 66755 49211 02670

45743 78815 47483 90799 78985 02007 52236 96837 96139 22783 54193 25572
45743 78815 47483 90799 78985 02007 52236 96837 96139 22783 54193 25572

23284 13846 47744 13529 09705 46512 24383 02697 56051 83775 74220 87783
23284 13846 47744 13529 09705 46512 24383 02697 56051 83776 17781 64242

58531 52464 74933 09145 87633 82311 24903 32030 12680 51006 70223 31257
33783 03370 18192 64880 28277 68611 91509 85606 75901 21359 85563 63034

50509 42448 46026 71622 54894 07922 61404 67995 06236 59692 82873 05828
32100 56649 97826 76360 88711 52327 56610 84900 93773 38023 19504 70687

78720 53603 03457 07660 66681 37431 25662 67431 40899 2606
65729 38513 59243 19759 37947 36057 50636 20935 07853 2833


In [21]:
print_diff_table(atan_239_true, atan_239_shanks)

Out[21]:
00418 40760 02074 72386 45382 14959 28545 27410 48065 30763 19508 27019
00418 40760 02074 72386 45382 14959 28545 27410 48065 30763 19508 27019

61288 71817 78341 42289 32737 82605 81362 29094 54975 45066 64448 63756
61288 71817 78341 42289 32737 82605 81362 29094 54975 45066 64448 63756

05245 83947 89311 86505 89221 28833 09280 08462 71962 33077 33759 47634
05245 83947 89311 86505 89221 28833 09280 08462 71962 33077 33759 47634

60331 84734 14570 33198 60154 54814 80599 24498 30211 46039 12539 49527
60331 84734 14570 33198 60154 54814 80599 24498 30211 46039 12539 49527

60779 68815 58881 27339 78533 46518 04574 25481 35867 46447 51979 10232
60779 68815 58881 27339 78533 46518 04574 25481 35867 46447 51979 10232

83097 70020 64652 82763 46532 96910 48183 86543 56078 91959 14512 32220
83097 70020 64652 82763 46532 96910 48183 86543 56078 91959 14512 32220

94463 68627 66155 20831 67964 26465 74655 11032 51034 35262 82445 12693
94463 68627 66155 20831 67964 26465 74655 11032 51034 35262 82445 12693

55670 49968 44452 47904 33177 28393 07086 31401 93869 51950 37058 64107
55670 49968 44452 47904 33177 28393 07086 31401 93869 51950 37058 64107

70855 85540 45223 55388 14237 67708 36515 69182 52702 00229 30895 44950
70855 85540 45223 55388 14237 67708 36515 69182 52702 00229 30895 44950

04358 54409 34496 44014 24187 24950 92283 86239 54553 33565 11719 73747
04358 54409 34496 44014 24187 24950 92283 86239 54553 33565 16494 21220

02023 49475 97790 97469 50111 88854 76673 97957 31537 09303 27821 13089
06852 38821 94006 45849 29313 23886 73467 64889 18731 81682 83021 21101

84258 30836 77190 91008 39098 51655 10419 22416 78092 0533
37897 11546 96191 84692 18237 33903 04682 04140 79985 6684


In [22]:
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.

Out[22]:
14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944
14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944

59230 78164 06286 20899 86280 34825 34211 70679 82148 08651 32823 06647
59230 78164 06286 20899 86280 34825 34211 70679 82148 08651 32823 06647

09384 46095 50582 23172 53594 08128 48111 74502 84102 70193 85211 05559
09384 46095 50582 23172 53594 08128 48111 74502 84102 70193 85211 05559

64462 29489 54930 38196 44288 10975 66593 34461 28475 64823 37867 83165
64462 29489 54930 38196 44288 10975 66593 34461 28475 64823 37867 83165

27120 19091 45648 56692 34603 48610 45432 66482 13393 60726 02491 41273
27120 19091 45648 56692 34603 48610 45432 66482 13393 60726 02491 41273

72458 70066 06315 58817 48815 20920 96282 92540 91715 36436 78925 90360
72458 70066 06315 58817 48815 20920 96282 92540 91715 36436 78925 90360

01133 05305 48820 46652 13841 46951 94151 16094 33057 27036 57595 91953
01133 05305 48820 46652 13841 46951 94151 16094 33057 27036 57595 91953

09218 61173 81932 61179 31051 18548 07446 23799 62749 56735 18857 52724
09218 61173 81932 61179 31051 18548 07446 23798 34749 56735 18857 52724

89122 79381 83011 94912 98336 73362 44065 66430 86021 39494 63952 24737
89122 79381 83011 94912 98336 73362 44193 66430 86021 39501 60924 48077

19070 21798 60943 70277 05392 17176 29317 67523 84674 81846 76694 05132
23094 36285 53096 62027 55693 97986 95022 24749 96206 07497 03041 23668

00056 81271 45263 56082 77857 71342 75778 96091 73637 17872 14684 40901
86199 51100 89202 38377 02131 41694 11902 98858 25446 81639 79990 46597

22495 34301 46549 58537 10507 92279 68925 89235 42019 96
00081 70029 63123 77381 34208 41307 91451 18398 05709 85


In [23]:
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.

Out[23]:
14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944
14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944

59230 78164 06286 20899 86280 34825 34211 70679 82148 08651 32823 06647
59230 78164 06286 20899 86280 34825 34211 70679 82148 08651 32823 06647

09384 46095 50582 23172 53594 08128 48111 74502 84102 70193 85211 05559
09384 46095 50582 23172 53594 08128 48111 74502 84102 70193 85211 05559

64462 29489 54930 38196 44288 10975 66593 34461 28475 64823 37867 83165
64462 29489 54930 38196 44288 10975 66593 34461 28475 64823 37867 83165

27120 19091 45648 56692 34603 48610 45432 66482 13393 60726 02491 41273
27120 19091 45648 56692 34603 48610 45432 66482 13393 60726 02491 41273

72458 70066 06315 58817 48815 20920 96282 92540 91715 36436 78925 90360
72458 70066 06315 58817 48815 30920 96282 92540 91715 36436 78925 90360

01133 05305 48820 46652 13841 46951 94151 16094 33057 27036 57595 91953
01133 05305 48820 46652 13841 46951 94151 16094 33057 27036 57595 91953

09218 61173 81932 61179 31051 18548 07446 23799 62749 56735 18857 52724
09218 61173 81932 61179 31051 18548 07446 23799 62749 56735 18857 52724

89122 79381 83011 94912 98336 73362 44065 66430 86021 39494 63952 24737
89122 79381 83011 94912 98336 73362 44065 66430 86021 39501 60924 48077

19070 21798 60943 70277 05392 17176 29317 67523 84674 81846 76694 05132
23094 36285 53096 62027 55693 97986 95022 24749 96206 07497 03041 23668

00056 81271 45263 56082 77857 71342 75778 96091 73637 17872 14684 40901
86199 51100 89202 38377 02131 41694 11902 98858 25446 81639 79990 46597

22495 34301 46549 58537 10507 92279 68925 89235 42019 96
00081 70029 63123 77387 34208 41307 91451 18398 05709 85


In [24]:
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.

Out[24]:
14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944
14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944

59230 78164 06286 20899 86280 34825 34211 70679 82148 08651 32823 06647
59230 78164 06286 20899 86280 34825 34211 70679 82148 08651 32823 06647

09384 46095 50582 23172 53594 08128 48111 74502 84102 70193 85211 05559
09384 46095 50582 23172 53594 08128 48111 74502 84102 70193 85211 05559

64462 29489 54930 38196 44288 10975 66593 34461 28475 64823 37867 83165
64462 29489 54930 38196 44288 10975 66593 34461 28475 64823 37867 83165

27120 19091 45648 56692 34603 48610 45432 66482 13393 60726 02491 41273
27120 19091 45648 56692 34603 48610 45432 66482 13393 60726 02491 41273

72458 70066 06315 58817 48815 20920 96282 92540 91715 36436 78925 90360
72458 70066 06315 58817 48815 20920 96282 92540 91715 36436 78925 90360

01133 05305 48820 46652 13841 46951 94151 16094 33057 27036 57595 91953
01133 05305 48820 46652 13841 46951 94151 16094 33057 27036 57595 91953

09218 61173 81932 61179 31051 18548 07446 23799 62749 56735 18857 52724
09218 61173 81932 61179 31051 18548 07446 23799 62749 56735 18857 52724

89122 79381 83011 94912 98336 73362 44065 66430 86021 39494 63952 24737
89122 79381 83011 94912 98336 73362 44065 66430 86021 39501 60924 48077

19070 21798 60943 70277 05392 17176 29317 67523 84674 81846 76694 05132
23094 36285 53096 62027 55693 97986 95022 24749 96206 07497 03041 23668

00056 81271 45263 56082 77857 71342 75778 96091 73637 17872 14684 40901
86199 51100 89202 38377 02131 41694 11902 98858 25446 81639 79990 46597

22495 34301 46549 58537 10507 92279 68925 89235 42019 96
00081 70029 63123 77387 34208 41307 91451 18398 05709 85


In [25]:
# 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)

Out[25]:
Decimal('3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019957')


## Arctangent computations

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.

In [26]:
# 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)

In [27]:
arctan_silicon(1, 5, 709)

Out[27]:
Decimal('0.1973955598498807583700497651947902934475851037878521015176889402410339699782437857326978280372880441126281180736913601044564798867942393557475654952163032700522107470015645015560061286185526633257318692806643896806189528405825931124251613297313993397113233537821796084176648310525473039665725650488878155309384290579311695934192851806364919697519401708560949527368673738508400812367856158009329822514023246675549211026704574378815474839079978985020075223696837961392278354193255722328413846477441352909705465122438302697560518377574220877835853152464749330914587633823112490332030126805100670223312575050942448460267162254894079226140467995062365969282873058287872053603034570766066681374312566267431408992603')

In [28]:
# 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]

Out[28]:
False

In [29]:
# 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)


In [30]:
'{: .709f}'.format(arctan_term(1, 239, 709, 2))

Out[30]:
' 0.0000000000002564723144246473657052071108829952397959831309702859408384067108509352766796256526773237818220189053275836401771021551025640430246795106737638335808357386545495569466345656953541462710525538407816432168166396004211207689032255349482007271917366509046832075065856373672538036032558379147704262196836630492673776047538783075337828273111407506598189378794060269980816455984150889714518259680194368984450770369538450687118309727499841085557454243718916532322645544326412797690276830329752816233526954204140839275226388252072697023474690753958955106864438453237754583868042529252587209904867359448739683787077380738237477542839474073223274941313611680998980279518279239497601296308507856749936712579758'


#### Paper-and-pencil algorithms

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.

In [31]:
# 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)

In [32]:
# The silicon and pencil-lead procedures should give the same result.

arctan_pencil(1, 5, 709) == arctan_silicon(1, 5, 709)

Out[32]:
True


#### Tabulation

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.

In [33]:
# 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)

In [34]:
# If we were to continue this to 700 digits, the table would get somewhat unruly!

arctan_pencil_table(1, 5, 40)

  k sgn   m                    x^m                                     (x^m)/m
0   1   1  0.2000000000000000000000000000000000000000   0.2000000000000000000000000000000000000000
1  -1   3  0.0080000000000000000000000000000000000000  -0.0026666666666666666666666666666666666667
2   1   5  0.0003200000000000000000000000000000000000   0.0000640000000000000000000000000000000000
3  -1   7  0.0000128000000000000000000000000000000000  -0.0000018285714285714285714285714285714286
4   1   9  0.0000005120000000000000000000000000000000   0.0000000568888888888888888888888888888889
5  -1  11  0.0000000204800000000000000000000000000000  -0.0000000018618181818181818181818181818182
6   1  13  0.0000000008192000000000000000000000000000   0.0000000000630153846153846153846153846154
7  -1  15  0.0000000000327680000000000000000000000000  -0.0000000000021845333333333333333333333333
8   1  17  0.0000000000013107200000000000000000000000   0.0000000000000771011764705882352941176471
9  -1  19  0.0000000000000524288000000000000000000000  -0.0000000000000027594105263157894736842105
10   1  21  0.0000000000000020971520000000000000000000   0.0000000000000000998643809523809523809524
11  -1  23  0.0000000000000000838860800000000000000000  -0.0000000000000000036472208695652173913043
12   1  25  0.0000000000000000033554432000000000000000   0.0000000000000000001342177280000000000000
13  -1  27  0.0000000000000000001342177280000000000000  -0.0000000000000000000049710269629629629630
14   1  29  0.0000000000000000000053687091200000000000   0.0000000000000000000001851279006896551724
15  -1  31  0.0000000000000000000002147483648000000000  -0.0000000000000000000000069273666064516129
16   1  33  0.0000000000000000000000085899345920000000   0.0000000000000000000000002603010482424242
17  -1  35  0.0000000000000000000000003435973836800000  -0.0000000000000000000000000098170681051429
18   1  37  0.0000000000000000000000000137438953472000   0.0000000000000000000000000003714566310054
19  -1  39  0.0000000000000000000000000005497558138880  -0.0000000000000000000000000000140963029202
20   1  41  0.0000000000000000000000000000219902325555   0.0000000000000000000000000000005363471355
21  -1  43  0.0000000000000000000000000000008796093022  -0.0000000000000000000000000000000204560303
22   1  45  0.0000000000000000000000000000000351843721   0.0000000000000000000000000000000007818749
23  -1  47  0.0000000000000000000000000000000014073749  -0.0000000000000000000000000000000000299441
24   1  49  0.0000000000000000000000000000000000562950   0.0000000000000000000000000000000000011489
25  -1  51  0.0000000000000000000000000000000000022518  -0.0000000000000000000000000000000000000442
26   1  53  0.0000000000000000000000000000000000000901   0.0000000000000000000000000000000000000017


Out[34]:
Decimal('0.1973955598498807583700497651947902934476')

In [35]:
arctan_pencil_table(1, 239, 40)

  k sgn   m                    x^m                                     (x^m)/m
0   1   1  0.0041841004184100418410041841004184100418   0.0041841004184100418410041841004184100418
1  -1   3  0.0000000732497753612514108822356769037379  -0.0000000244165917870838036274118923012460
2   1   5  0.0000000000012823615721232368285260355544   0.0000000000002564723144246473657052071109
3  -1   7  0.0000000000000000224499146044928630193105  -0.0000000000000000032071306577846947170444
4   1   9  0.0000000000000000000003930238371963527077   0.0000000000000000000000436693152440391897
5  -1  11  0.0000000000000000000000000068805489609137  -0.0000000000000000000000000006255044509922
6   1  13  0.0000000000000000000000000000001204556811   0.0000000000000000000000000000000092658216
7  -1  15  0.0000000000000000000000000000000000021088  -0.0000000000000000000000000000000000001406


Out[35]:
Decimal('0.004184076002074723864538214959285452741046')


#### Reciprocals

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.

In [36]:
# 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)

In [37]:
arctan_pencil_recip(5, 25, table=True)

  k sgn   m             b^m                       m(b^m)                     1/m(b^m)
0   1   1   5                          5                          0.2000000000000000000000000
1  -1   3   125                        375                       -0.0026666666666666666666667
2   1   5   3125                       15625                      0.0000640000000000000000000
3  -1   7   78125                      546875                    -0.0000018285714285714285714
4   1   9   1953125                    17578125                   0.0000000568888888888888889
5  -1  11   48828125                   537109375                 -0.0000000018618181818181818
6   1  13   1220703125                 15869140625                0.0000000000630153846153846
7  -1  15   30517578125                457763671875              -0.0000000000021845333333333
8   1  17   762939453125               12969970703125             0.0000000000000771011764706
9  -1  19   19073486328125             362396240234375           -0.0000000000000027594105263
10   1  21   476837158203125            10013580322265625          0.0000000000000000998643810
11  -1  23   11920928955078125          274181365966796875        -0.0000000000000000036472209
12   1  25   298023223876953125         7450580596923828125        0.0000000000000000001342177
13  -1  27   7450580596923828125        201165676116943359375     -0.0000000000000000000049710
14   1  29   186264514923095703125      5401670932769775390625     0.0000000000000000000001851
15  -1  31   4656612873077392578125     144354999065399169921875  -0.0000000000000000000000069
16   1  33   116415321826934814453125   3841705620288848876953125   0.0000000000000000000000003


Out[37]:
Decimal('0.1973955598498807583700499')

In [38]:
# Just checking...

arctan_pencil_recip(5, 709) == arctan_silicon(1, 5, 709)

Out[38]:
True


## Forensic mathematics¶

How and where did Shanks stray into error?

### Search for omitted terms

In [39]:
# 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:
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 '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]
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.

In [40]:
search_omitted__or_duped_terms(1, 5, 709, atan_5_shanks, (532,541))

match string = 617781642
Searched 506 omissions and duplicates.


In [41]:
search_omitted__or_duped_terms(1, 5, 709, atan_5_shanks_609, (603,608))

match string = 34783
Searched 506 omissions and duplicates.


In [42]:
search_omitted__or_duped_terms(1, 239, 709, atan_239_shanks, (594,600))

match string = 649421
Searched 150 omissions and duplicates.



### Search for shifts (omitted or inserted digits)

In [43]:
# 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:
tru_term_str = fmt.format(term)
diff_term_str = fmt.format(atan_diff_d)
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 'diff  = ', diff_term_str[left:right]
k += 1
term = arctan_term(a, b, p, k)
print 'Searched for shifts in {} terms.'.format(k+1)

In [44]:
search_shifts(1, 5, 709, atan_5_shanks, (537, 543), (521, 601))

# finds a shift by -1 in term 248

k = 248, m = 497, pattern = 973843
match at position 538, shifted by -1
true  =  90744466800804828973843058350100603621730382293762575452716297786720321931589537
diff  =  00000000000043560764587525150905432595573440643863006660653576632207035315340317
Searched for shifts in 506 terms.



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.

In [45]:
# 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()

Out[45]:
' 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000082328737623142401157273382716991930387742506188205854866279963986321272545389199066480636124896985843835178172222335789970252931234542078806049486654325955734406438631790744466800848389738430583501006036217303822937625754527162977867203219315895372233400402414486921529175050301810865191146881287726358148893360160965794768611670020120724346076458752515090540'

In [46]:
# 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'}

In [47]:
# 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))

In [48]:
# give it a try....

inject_error_terms(1, 5, 709, atan_5_error_dict_248)

Out[48]:
' 0.1973955598498807583700497651947902934475851037878521015176889402410339699782437857326978280372880441126281180736913601044564798867942393557475654952163032700522107470015645015560061286185526633257318692806643896806189528405825931124251613297313993397113233537821796084176648310525473039665725650488878155309384290579311695934192851806364919697519401708560949527368673738508400812367856158009329822514023246675549211026704574378815474839079978985020075223696837961392278354193255722328413846477441352909705465122438302697560518377617781642423378303370181926488028277686291564778710207287998054529147585111304621498496538512439350855919140027255524923005207062312016922818326321269084790026224035080308672574092'

In [49]:
# 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)

Out[49]:
19739 55598 49880 75837 00497 65194 79029 34475 85103 78785 21015 17688
19739 55598 49880 75837 00497 65194 79029 34475 85103 78785 21015 17688

94024 10339 69978 24378 57326 97828 03728 80441 12628 11807 36913 60104
94024 10339 69978 24378 57326 97828 03728 80441 12628 11807 36913 60104

45647 98867 94239 35574 75654 95216 30327 00522 10747 00156 45015 56006
45647 98867 94239 35574 75654 95216 30327 00522 10747 00156 45015 56006

12861 85526 63325 73186 92806 64389 68061 89528 40582 59311 24251 61329
12861 85526 63325 73186 92806 64389 68061 89528 40582 59311 24251 61329

73139 93397 11323 35378 21796 08417 66483 10525 47303 96657 25650 48887
73139 93397 11323 35378 21796 08417 66483 10525 47303 96657 25650 48887

81553 09384 29057 93116 95934 19285 18063 64919 69751 94017 08560 94952
81553 09384 29057 93116 95934 19285 18063 64919 69751 94017 08560 94952

73686 73738 50840 08123 67856 15800 93298 22514 02324 66755 49211 02670
73686 73738 50840 08123 67856 15800 93298 22514 02324 66755 49211 02670

45743 78815 47483 90799 78985 02007 52236 96837 96139 22783 54193 25572
45743 78815 47483 90799 78985 02007 52236 96837 96139 22783 54193 25572

23284 13846 47744 13529 09705 46512 24383 02697 56051 83775 74220 87783
23284 13846 47744 13529 09705 46512 24383 02697 56051 83776 17781 64242

58531 52464 74933 09145 87633 82311 24903 32030 12680 51006 70223 31257
33783 03370 18192 64880 28277 68611 91509 85606 75901 21359 85563 63034

50509 42448 46026 71622 54894 07922 61404 67995 06236 59692 82873 05828
32100 56649 97826 76360 88711 52327 56610 84900 93773 38023 19504 70687

78720 53603 03457 07660 66681 37431 25662 67431 40899 2606
65729 38513 59243 19759 37947 36057 50636 20935 07853 2833


In [50]:
# 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)

Out[50]:
19739 55598 49880 75837 00497 65194 79029 34475 85103 78785 21015 17688
19739 55598 49880 75837 00497 65194 79029 34475 85103 78785 21015 17688

94024 10339 69978 24378 57326 97828 03728 80441 12628 11807 36913 60104
94024 10339 69978 24378 57326 97828 03728 80441 12628 11807 36913 60104

45647 98867 94239 35574 75654 95216 30327 00522 10747 00156 45015 56006
45647 98867 94239 35574 75654 95216 30327 00522 10747 00156 45015 56006

12861 85526 63325 73186 92806 64389 68061 89528 40582 59311 24251 61329
12861 85526 63325 73186 92806 64389 68061 89528 40582 59311 24251 61329

73139 93397 11323 35378 21796 08417 66483 10525 47303 96657 25650 48887
73139 93397 11323 35378 21796 08417 66483 10525 47303 96657 25650 48887

81553 09384 29057 93116 95934 19285 18063 64919 69751 94017 08560 94952
81553 09384 29057 93116 95934 19285 18063 64919 69751 94017 08560 94952

73686 73738 50840 08123 67856 15800 93298 22514 02324 66755 49211 02670
73686 73738 50840 08123 67856 15800 93298 22514 02324 66755 49211 02670

45743 78815 47483 90799 78985 02007 52236 96837 96139 22783 54193 25572
45743 78815 47483 90799 78985 02007 52236 96837 96139 22783 54193 25572

23284 13846 47744 13529 09705 46512 24383 02697 56051 83776 17781 64242
23284 13846 47744 13529 09705 46512 24383 02697 56051 83776 17781 64242

33783 03370 18192 64880 28277 68629 15647 78710 20728 79980 54529 14758
33783 03370 18192 64880 28277 68611 91509 85606 75901 21359 85563 63034

51113 04621 49849 65385 12439 35085 59191 40027 25552 49230 05207 06231
32100 56649 97826 76360 88711 52327 56610 84900 93773 38023 19504 70687

20169 22818 32632 12690 84790 02622 40350 80308 67257 4092
65729 38513 59243 19759 37947 36057 50636 20935 07853 2833



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.

In [51]:
# 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()

Out[51]:
' 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000307596485496974112297044389967563606979040220689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000'

In [52]:
# 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'}

In [53]:
# try applying the two uncorrections

inject_error_terms(1, 5, 709, atan_5_error_dict_248_and_72_zeroes)

Out[53]:
' 0.1973955598498807583700497651947902934475851037878521015176889402410339699782437857326978280372880441126281180736913601044564798867942393557475654952163032700522107470015645015560061286185526633257318692806643896806189528405825931124251613297313993397113233537821796084176648310525473039665725650488878155309384290579311695934192851806364919697519401708560949527368673738508400812367856158009329822514023246675549211026704574378815474839079978985020075223696837961392278354193255722328413846477441352909705465122438302697560518377617781642423378303370181926488028277685601909606296414184549778667078619594063242188151710926232454304195002096221042164384517407139603129714878045407015824508982655769963844987885'

In [54]:
# 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)

Out[54]:
19739 55598 49880 75837 00497 65194 79029 34475 85103 78785 21015 17688
19739 55598 49880 75837 00497 65194 79029 34475 85103 78785 21015 17688

94024 10339 69978 24378 57326 97828 03728 80441 12628 11807 36913 60104
94024 10339 69978 24378 57326 97828 03728 80441 12628 11807 36913 60104

45647 98867 94239 35574 75654 95216 30327 00522 10747 00156 45015 56006
45647 98867 94239 35574 75654 95216 30327 00522 10747 00156 45015 56006

12861 85526 63325 73186 92806 64389 68061 89528 40582 59311 24251 61329
12861 85526 63325 73186 92806 64389 68061 89528 40582 59311 24251 61329

73139 93397 11323 35378 21796 08417 66483 10525 47303 96657 25650 48887
73139 93397 11323 35378 21796 08417 66483 10525 47303 96657 25650 48887

81553 09384 29057 93116 95934 19285 18063 64919 69751 94017 08560 94952
81553 09384 29057 93116 95934 19285 18063 64919 69751 94017 08560 94952

73686 73738 50840 08123 67856 15800 93298 22514 02324 66755 49211 02670
73686 73738 50840 08123 67856 15800 93298 22514 02324 66755 49211 02670

45743 78815 47483 90799 78985 02007 52236 96837 96139 22783 54193 25572
45743 78815 47483 90799 78985 02007 52236 96837 96139 22783 54193 25572

23284 13846 47744 13529 09705 46512 24383 02697 56051 83776 17781 64242
23284 13846 47744 13529 09705 46512 24383 02697 56051 83776 17781 64242

33783 03370 18192 64880 28277 68560 19096 06296 41418 45497 78667 07861
33783 03370 18192 64880 28277 68611 91509 85606 75901 21359 85563 63034

95940 63242 18815 17109 26232 45430 41950 02096 22104 21643 84517 40713
32100 56649 97826 76360 88711 52327 56610 84900 93773 38023 19504 70687

96031 29714 87804 54070 15824 50898 26557 69963 84498 7885
65729 38513 59243 19759 37947 36057 50636 20935 07853 2833



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.

In [73]:
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)


pi true     3.141 ... 43086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127
pi shanks   3.141 ... 43086021395016092448077230943628553096620275569397986950222474996206074970304123668861995110
pi diff    -0.000 ... 00000000000069697223340040241448692152917505030180810657045722611531256502634718536861426982

a5 true     0.197 ... 69756051837757422087783585315246474933091458763382311249033203012680510067022331257505094244
a5 shanks   0.197 ... 69756051837761778164242337830337018192648802827768611915098560675901213598556363034321005664
diff0      -0.000 ... 00000000000004356076458752515090543259557344064386300666065357663220703531534031776815911420

a5 248ucx   0.197 ... 69756051837761778164242337830337018192648802827768629156477871020728799805452914758511130462
a5 shanks   0.197 ... 69756051837761778164242337830337018192648802827768611915098560675901213598556363034321005664
diff1       0.000 ... 00000000000000000000000000000000000000000000000000017241379310344827586206896551724190124797

tn          0.000 ... 37931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034
tn-diff1    0.000 ... 37931034482758620689655172413793103448275862068965499999999999999999999999999999999947806237
diff1       0.000 ... 00000000000000000000000000000000000000000000000000017241379310344827586206896551724190124797


In [56]:
# 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()

Out[56]:
' 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000307596485496974112297044389967563606979040220689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655'

In [57]:
# 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'}

In [58]:
inject_error_terms(1, 5, 709, atan_5_error_dict_248_and_72_skip)

Out[58]:
' 0.1973955598498807583700497651947902934475851037878521015176889402410339699782437857326978280372880441126281180736913601044564798867942393557475654952163032700522107470015645015560061286185526633257318692806643896806189528405825931124251613297313993397113233537821796084176648310525473039665725650488878155309384290579311695934192851806364919697519401708560949527368673738508400812367856158009329822514023246675549211026704574378815474839079978985020075223696837961392278354193255722328413846477441352909705465122438302697560518377617781642423378303370181926488028277686119150985606759012135985563630343731994276670910331615887626717988105544496904233350034648518913474542464252303567548646913690252722465677540'

In [59]:
# 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)

Out[59]:
19739 55598 49880 75837 00497 65194 79029 34475 85103 78785 21015 17688
19739 55598 49880 75837 00497 65194 79029 34475 85103 78785 21015 17688

94024 10339 69978 24378 57326 97828 03728 80441 12628 11807 36913 60104
94024 10339 69978 24378 57326 97828 03728 80441 12628 11807 36913 60104

45647 98867 94239 35574 75654 95216 30327 00522 10747 00156 45015 56006
45647 98867 94239 35574 75654 95216 30327 00522 10747 00156 45015 56006

12861 85526 63325 73186 92806 64389 68061 89528 40582 59311 24251 61329
12861 85526 63325 73186 92806 64389 68061 89528 40582 59311 24251 61329

73139 93397 11323 35378 21796 08417 66483 10525 47303 96657 25650 48887
73139 93397 11323 35378 21796 08417 66483 10525 47303 96657 25650 48887

81553 09384 29057 93116 95934 19285 18063 64919 69751 94017 08560 94952
81553 09384 29057 93116 95934 19285 18063 64919 69751 94017 08560 94952

73686 73738 50840 08123 67856 15800 93298 22514 02324 66755 49211 02670
73686 73738 50840 08123 67856 15800 93298 22514 02324 66755 49211 02670

45743 78815 47483 90799 78985 02007 52236 96837 96139 22783 54193 25572
45743 78815 47483 90799 78985 02007 52236 96837 96139 22783 54193 25572

23284 13846 47744 13529 09705 46512 24383 02697 56051 83776 17781 64242
23284 13846 47744 13529 09705 46512 24383 02697 56051 83776 17781 64242

33783 03370 18192 64880 28277 68611 91509 85606 75901 21359 85563 63034
33783 03370 18192 64880 28277 68611 91509 85606 75901 21359 85563 63034

37319 94276 67091 03316 15887 62671 79881 05544 49690 42333 50034 64851
32100 56649 97826 76360 88711 52327 56610 84900 93773 38023 19504 70687

89134 74542 46425 23035 67548 64691 36902 52722 46567 7540
65729 38513 59243 19759 37947 36057 50636 20935 07853 2833



## Summary

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.

## Update: Further forensics¶

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!

In [60]:
# 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]

In [61]:
# Nope, this one doesn't work.

cyclic_diffs('123456')

111105
222156
332667
437778
488889


In [62]:
# Neither does this.

cyclic_diffs('910672')

803943
843381
237762
181566
619605


In [63]:
# But this one is perfect.

cyclic_diffs('142857')

285714
142857
714285
428571
571428


In [64]:
# And so does the 28-digit sequence plucked out of term 72.

cyclic_diffs('6896551724137931034482758620')

2068965517241379310344827586
2758620689655172413793103448
3448275862068965517241379310
1379310344827586206896551724
1724137931034482758620689655
5172413793103448275862068965
3448275862068965517241379310
4482758620689655172413793103
2758620689655172413793103448
5517241379310344827586206896
3103448275862068965517241379
1034482758620689655172413793
2413793103448275862068965517
3793103448275862068965517241
5862068965517241379310344827
6551724137931034482758620689
3448275862068965517241379310
2413793103448275862068965517
2068965517241379310344827586
1379310344827586206896551724
4137931034482758620689655172
6896551724137931034482758620
1034482758620689655172413793
1724137931034482758620689655
6896551724137931034482758620
4827586206896551724137931034
6206896551724137931034482758



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 [65]:
# 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()

Out[65]:
' 0.0000000000000000000001851279006896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000'

In [66]:
# 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'}

In [67]:
inject_error_terms(1, 5, 709, atan_5_error_dict_248_and_14_zeroes)

Out[67]:
' 0.1973955598498807583700497651947902934475851037878521015176889402410339699782437857326978280372880441126281180736913601044564798867942393557475654952163032700522107470015645015560061286185526633257318692806643896806189528405825931124251613297313993397113233537821796084176648310525473039665725650488878155309384290579311695934192851806364919697519401708560949527368673738508400812367856158009329822514023246675549211026704574378815474839079978985020075223696837961392278354193255722328413846477441352909705465122438302697560518377617781642423378303370181926488028277686119150985606759012135985563630343731994276670910331615887626717988105544496904233350034648518913474542464252303567548646913690252722465677540'

In [68]:
print_diff_table(inject_error_terms(1, 5, 709, atan_5_error_dict_248_and_14_zeroes), atan_5_shanks)

Out[68]:
19739 55598 49880 75837 00497 65194 79029 34475 85103 78785 21015 17688
19739 55598 49880 75837 00497 65194 79029 34475 85103 78785 21015 17688

94024 10339 69978 24378 57326 97828 03728 80441 12628 11807 36913 60104
94024 10339 69978 24378 57326 97828 03728 80441 12628 11807 36913 60104

45647 98867 94239 35574 75654 95216 30327 00522 10747 00156 45015 56006
45647 98867 94239 35574 75654 95216 30327 00522 10747 00156 45015 56006

12861 85526 63325 73186 92806 64389 68061 89528 40582 59311 24251 61329
12861 85526 63325 73186 92806 64389 68061 89528 40582 59311 24251 61329

73139 93397 11323 35378 21796 08417 66483 10525 47303 96657 25650 48887
73139 93397 11323 35378 21796 08417 66483 10525 47303 96657 25650 48887

81553 09384 29057 93116 95934 19285 18063 64919 69751 94017 08560 94952
81553 09384 29057 93116 95934 19285 18063 64919 69751 94017 08560 94952

73686 73738 50840 08123 67856 15800 93298 22514 02324 66755 49211 02670
73686 73738 50840 08123 67856 15800 93298 22514 02324 66755 49211 02670

45743 78815 47483 90799 78985 02007 52236 96837 96139 22783 54193 25572
45743 78815 47483 90799 78985 02007 52236 96837 96139 22783 54193 25572

23284 13846 47744 13529 09705 46512 24383 02697 56051 83776 17781 64242
23284 13846 47744 13529 09705 46512 24383 02697 56051 83776 17781 64242

33783 03370 18192 64880 28277 68611 91509 85606 75901 21359 85563 63034
33783 03370 18192 64880 28277 68611 91509 85606 75901 21359 85563 63034

37319 94276 67091 03316 15887 62671 79881 05544 49690 42333 50034 64851
32100 56649 97826 76360 88711 52327 56610 84900 93773 38023 19504 70687

89134 74542 46425 23035 67548 64691 36902 52722 46567 7540
65729 38513 59243 19759 37947 36057 50636 20935 07853 2833



## Weighing the options

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.