Pi.html

A Sample of some Historical

Calculations

Greg Fee

based on Jonathon Borwein's talk

The Life of

2:00-3:00 pm,March 14, 2003, K9509, Simon Fraser University

>

Archimedes 's Calculation

About 240 BC Archimedes inscribed and circumscribed regular polygons with (6,12,24,48, and 96) sides around a circle and obtained 3+10/71 < < 3+10/70

The following Maple commands generate a picture, illustrating Archimedes polygons.

>    ngon_cir := proc(n::posint)
local cir,t,pn,ngon1,R,ngon2,i;
cir := plot([cos(t),sin(t),t=-Pi..Pi],numpoints=361,colour=red);
pn := evalhf(2*Pi/n);
ngon1 := plots[polygonplot]([seq([cos(pn*i),sin(pn*i)],i=0..n-1)],
  colour=green,style=line);
R := sec(.5*pn);
ngon2 := plots[polygonplot]([seq([R*cos(pn*i),R*sin(pn*i)],i=0..n-1)],
  colour=blue,style=line);
plots[display]({cir,ngon1,ngon2},scaling=constrained);
end:

> ngon_cir(6);

> ngon_cir(12);

> ngon_cir(24);

> ngon_cir(48);

> ngon_cir(96);

> ngon_cir1 := proc(n::posint)
local pn,a1,a2,cir,t,ngon1,R,ngon2;
pn := evalhf(Pi/n);
a1 := evalhf(Pi/2-pn);
a2 := evalhf(Pi/2+pn);
cir := plot([cos(t),sin(t),t=a1..a2],colour=red);
ngon1 := plot([[cos(a1),sin(a1)],[cos(a2),sin(a2)]],colour=green);
R := sec(pn);
ngon2 := plot([[R*cos(a1),1.],[R*cos(a2),1.]],colour=blue,style=line);
plots[display]({cir,ngon1,ngon2},scaling=constrained,axes=none,view=[R*cos(a2)..R*cos(a1),.999*sin(a1)..1.001]);
end:

> ngon_cir1(6);

> ngon_cir1(12);

> ngon_cir1(24);

> ngon_cir1(48);

> ngon_cir1(96);

> ngon_cir_half := proc(n::posint)
local pn,a1,a2,cir,t,ngon1,R,ngon2;
pn := evalhf(Pi/n);
a1 := evalhf(Pi/2-pn);
a2 := evalhf(Pi/2);
cir := plot([cos(t),sin(t),t=a1..a2],colour=red);
ngon1 := plot([[cos(a1),sin(a1)],[cos(a2),sin(a1)]],colour=green);
R := sec(pn);
ngon2 := plot([[R*cos(a1),1.],[R*cos(a2),1.]],colour=blue,style=line);
plots[display]({cir,ngon1,ngon2},scaling=constrained,axes=none,view=[R*cos(a2)..R*cos(a1),.999*sin(a1)..1.001]);
end:

> ngon_cir_half(96);

>

>    Archimedes := proc(n)
local a,b,e,s,eps,k;
a := evalf(2*3^(1/2));
b := 3.;
e := a-b;
s := 6;
userinfo(2,procname,'iteration',1);
userinfo(3,procname,cat("inner polygon has ",s," sides"));
userinfo(4,procname,b<Pi);
userinfo(3,procname,cat("outer polygon has ",s," sides"));
userinfo(4,procname,Pi<a);
eps := Float(1,1-Digits);
for k from 2 to n while eps<e do
  s := 2*s;
  a := 2*a*b/(a+b);
  b := a*b;
  b := b^(1/2);
  e := a-b;
  userinfo(2,procname,'iteration',k);
  userinfo(3,procname,cat("inner polygon has ",s," sides"));
  userinfo(4,procname,b<Pi);
  userinfo(3,procname,cat("outer polygon has ",s," sides"));
  userinfo(4,procname,Pi<a)
od;
userinfo(1,procname,cat(k-1," iterations"));
b+.5*e;
end:

> infolevel[Archimedes] := 5;

> Api := Archimedes(4);

Archimedes:   iteration   1

Archimedes:   "inner polygon has 6 sides"

Archimedes:   3. < Pi

Archimedes:   "outer polygon has 6 sides"

Archimedes:   Pi < 3.464101616

Archimedes:   iteration   2

Archimedes:   "inner polygon has 12 sides"

Archimedes:   3.105828542 < Pi

Archimedes:   "outer polygon has 12 sides"

Archimedes:   Pi < 3.215390310

Archimedes:   iteration   3

Archimedes:   "inner polygon has 24 sides"

Archimedes:   3.132628614 < Pi

Archimedes:   "outer polygon has 24 sides"

Archimedes:   Pi < 3.159659942

Archimedes:   iteration   4

Archimedes:   "inner polygon has 48 sides"

Archimedes:   3.139350204 < Pi

Archimedes:   "outer polygon has 48 sides"

Archimedes:   Pi < 3.146086216

Archimedes:   4 iterations

>

van Ceulen's Calculation

By 1610 Ludolph van Ceulen inscribed and circumscribed regular polygons with 2^62 sides around a circle and obtained Pi to 35 decimal places.

This is similar to Archimedes method except he started with a square instead of a hexagon.

>

>    Ceulen := proc(n)
local a,b,e,s,eps,k;
a := 4.;
b := 2*evalf(sqrt(2));
e := a-b;
s := 4;
userinfo(2,procname,'iteration',1);
userinfo(3,procname,cat("inner polygon has ",s," sides"));
userinfo(4,procname,b<Pi);
userinfo(3,procname,cat("outer polygon has ",s," sides"));
userinfo(4,procname,Pi<a);
eps := Float(1,1-Digits);
for k from 2 to n while eps<e do
  s := 2*s;
  a := 2*a*b/(a+b);
  b := a*b;
  b := b^(1/2);
  e := a-b;
  userinfo(2,procname,'iteration',k);
  userinfo(3,procname,cat("inner polygon has ",s," sides"));
  userinfo(4,procname,b<Pi);
  userinfo(3,procname,cat("outer polygon has ",s," sides"));
  userinfo(4,procname,Pi<a);
od;
userinfo(1,procname,cat(k-1," iterations"));
b+.5*e;
end:

> Digits := 38;

> infolevel[Ceulen] := 5;

> Cpi := Ceulen(61);

Ceulen:   iteration   1

Ceulen:   "inner polygon has 4 sides"

Ceulen:   2.8284271247461900976033774484193961572 < Pi

Ceulen:   "outer polygon has 4 sides"

Ceulen:   Pi < 4.

Ceulen:   iteration   2

Ceulen:   "inner polygon has 8 sides"

Ceulen:   3.0614674589207181738276798722431909341 < Pi

Ceulen:   "outer polygon has 8 sides"

Ceulen:   Pi < 3.3137084989847603904135097936775846286

Ceulen:   iteration   3

Ceulen:   "inner polygon has 16 sides"

Ceulen:   3.1214451522580522855725578956323558548 < Pi

Ceulen:   "outer polygon has 16 sides"

Ceulen:   Pi < 3.1825978780745281105855619623148196574

Ceulen:   iteration   4

Ceulen:   "inner polygon has 32 sides"

Ceulen:   3.1365484905459392638142580444365390675 < Pi

Ceulen:   "outer polygon has 32 sides"

Ceulen:   Pi < 3.1517249074292560984703206813224778332

Ceulen:   iteration   5

Ceulen:   "inner polygon has 64 sides"

Ceulen:   3.1403311569547529123171185243316901321 < Pi

Ceulen:   "outer polygon has 64 sides"

Ceulen:   Pi < 3.1441183852459042627419725613640714930

Ceulen:   iteration   6

Ceulen:   "inner polygon has 128 sides"

Ceulen:   3.1412772509327728680620197707882144083 < Pi

Ceulen:   "outer polygon has 128 sides"

Ceulen:   Pi < 3.1422236299424568453862085069963163126

Ceulen:   iteration   7

Ceulen:   "inner polygon has 256 sides"

Ceulen:   3.1415138011443010763285150594568223079 < Pi

Ceulen:   "outer polygon has 256 sides"

Ceulen:   Pi < 3.1417503691689664591072136279733238900

Ceulen:   iteration   8

Ceulen:   "inner polygon has 512 sides"

Ceulen:   3.1415729403670913841358001102707614295 < Pi

Ceulen:   "outer polygon has 512 sides"

Ceulen:   Pi < 3.1416320807031818057187151878711476690

Ceulen:   iteration   9

Ceulen:   "inner polygon has 1024 sides"

Ceulen:   3.1415877252771597006288542627019187394 < Pi

Ceulen:   "outer polygon has 1024 sides"

Ceulen:   Pi < 3.1416025102568089467636896584926600118

Ceulen:   iteration   10

Ceulen:   "inner polygon has 2048 sides"

Ceulen:   3.1415914215111999739979717637408339557 < Pi

Ceulen:   "outer polygon has 2048 sides"

Ceulen:   Pi < 3.1415951177495890503530922359816959454

Ceulen:   iteration   11

Ceulen:   "inner polygon has 4096 sides"

Ceulen:   3.1415923455701177423403759941573699303 < Pi

Ceulen:   "outer polygon has 4096 sides"

Ceulen:   Pi < 3.1415932696293073107894587868279757644

Ceulen:   iteration   12

Ceulen:   "inner polygon has 8192 sides"

Ceulen:   3.1415925765848726656816060922378753098 < Pi

Ceulen:   "outer polygon has 8192 sides"

Ceulen:   Pi < 3.1415928075996445765282544359605834636

Ceulen:   iteration   13

Ceulen:   "inner polygon has 16384 sides"

Ceulen:   3.1415926343385629890954782636277912940 < Pi

Ceulen:   "outer polygon has 16384 sides"

Ceulen:   Pi < 3.1415926920922543742284195571808838588

Ceulen:   iteration   14

Ceulen:   "inner polygon has 32768 sides"

Ceulen:   3.1415926487769856694851079692771770757 < Pi

Ceulen:   "outer polygon has 32768 sides"

Ceulen:   Pi < 3.1415926632154084162321791900900738404

Ceulen:   iteration   15

Ceulen:   "inner polygon has 65536 sides"

Ceulen:   3.1415926523865913458035255210579638844 < Pi

Ceulen:   "outer polygon has 65536 sides"

Ceulen:   Pi < 3.1415926559961970262692831627712876004

Ceulen:   iteration   16

Ceulen:   "inner polygon has 131072 sides"

Ceulen:   3.1415926532889927652719430421737400035 < Pi

Ceulen:   "outer polygon has 131072 sides"

Ceulen:   Pi < 3.1415926541913941849995693188358437026

Ceulen:   iteration   17

Ceulen:   "inner polygon has 262144 sides"

Ceulen:   3.1415926535145931201633482432810804327 < Pi

Ceulen:   "outer polygon has 262144 sides"

Ceulen:   Pi < 3.1415926537401934750709539916089029610

Ceulen:   iteration   18

Ceulen:   "inner polygon has 524288 sides"

Ceulen:   3.1415926535709932088877183448597721147 < Pi

Ceulen:   "outer polygon has 524288 sides"

Ceulen:   Pi < 3.1415926536273932976131009806397257502

Ceulen:   iteration   19

Ceulen:   "inner polygon has 1048576 sides"

Ceulen:   3.1415926535850932310689057953358123492 < Pi

Ceulen:   "outer polygon has 1048576 sides"

Ceulen:   Pi < 3.1415926535991932532501565291994311718

Ceulen:   iteration   20

Ceulen:   "inner polygon has 2097152 sides"

Ceulen:   3.1415926535886182366142085907724078849 < Pi

Ceulen:   "outer polygon has 2097152 sides"

Ceulen:   Pi < 3.1415926535921432421595153414207270780

Ceulen:   iteration   21

Ceulen:   "inner polygon has 4194304 sides"

Ceulen:   3.1415926535894994880005346604326558615 < Pi

Ceulen:   "outer polygon has 4194304 sides"

Ceulen:   Pi < 3.1415926535903807393868609772936365666

Ceulen:   iteration   22

Ceulen:   "inner polygon has 8388608 sides"

Ceulen:   3.1415926535897198008471162010227865490 < Pi

Ceulen:   "outer polygon has 8388608 sides"

Ceulen:   Pi < 3.1415926535899401136936977570629630320

Ceulen:   iteration   23

Ceulen:   "inner polygon has 16777216 sides"

Ceulen:   3.1415926535897748790587615876187610142 < Pi

Ceulen:   "outer polygon has 16777216 sides"

Ceulen:   Pi < 3.1415926535898299572704069751803633416

Ceulen:   iteration   24

Ceulen:   "inner polygon has 33554432 sides"

Ceulen:   3.1415926535897886486116729343582822426 < Pi

Ceulen:   "outer polygon has 33554432 sides"

Ceulen:   Pi < 3.1415926535898024181645842811581552124

Ceulen:   iteration   25

Ceulen:   "inner polygon has 67108864 sides"

Ceulen:   3.1415926535897920909999007710488205255 < Pi

Ceulen:   "outer polygon has 67108864 sides"

Ceulen:   Pi < 3.1415926535897955333881286077431307922

Ceulen:   iteration   26

Ceulen:   "inner polygon has 134217728 sides"

Ceulen:   3.1415926535897929515969577302218087198 < Pi

Ceulen:   "outer polygon has 134217728 sides"

Ceulen:   Pi < 3.1415926535897938121940146893950326630

Ceulen:   iteration   27

Ceulen:   "inner polygon has 268435456 sides"

Ceulen:   3.1415926535897931667462219700150778698 < Pi

Ceulen:   "outer polygon has 268435456 sides"

Ceulen:   Pi < 3.1415926535897933818954862098083617542

Ceulen:   iteration   28

Ceulen:   "inner polygon has 536870912 sides"

Ceulen:   3.1415926535897932205335380299633965387 < Pi

Ceulen:   "outer polygon has 536870912 sides"

Ceulen:   Pi < 3.1415926535897932743208540899117161284

Ceulen:   iteration   29

Ceulen:   "inner polygon has 1073741824 sides"

Ceulen:   3.1415926535897932339803670449504762923 < Pi

Ceulen:   "outer polygon has 1073741824 sides"

Ceulen:   Pi < 3.1415926535897932474271960599375561034

Ceulen:   iteration   30

Ceulen:   "inner polygon has 2147483648 sides"

Ceulen:   3.1415926535897932373420742986972462361 < Pi

Ceulen:   "outer polygon has 2147483648 sides"

Ceulen:   Pi < 3.1415926535897932407037815524440161834

Ceulen:   iteration   31

Ceulen:   "inner polygon has 4294967296 sides"

Ceulen:   3.1415926535897932381825011121339387223 < Pi

Ceulen:   "outer polygon has 4294967296 sides"

Ceulen:   Pi < 3.1415926535897932390229279255706312088

Ceulen:   iteration   32

Ceulen:   "inner polygon has 8589934592 sides"

Ceulen:   3.1415926535897932383926078154931118438 < Pi

Ceulen:   "outer polygon has 8589934592 sides"

Ceulen:   Pi < 3.1415926535897932386027145188522849654

Ceulen:   iteration   33

Ceulen:   "inner polygon has 17179869184 sides"

Ceulen:   3.1415926535897932384451344913329051242 < Pi

Ceulen:   "outer polygon has 17179869184 sides"

Ceulen:   Pi < 3.1415926535897932384976611671726984046

Ceulen:   iteration   34

Ceulen:   "inner polygon has 34359738368 sides"

Ceulen:   3.1415926535897932384582661602928534443 < Pi

Ceulen:   "outer polygon has 34359738368 sides"

Ceulen:   Pi < 3.1415926535897932384713978292528017644

Ceulen:   iteration   35

Ceulen:   "inner polygon has 68719476736 sides"

Ceulen:   3.1415926535897932384615490775328405243 < Pi

Ceulen:   "outer polygon has 68719476736 sides"

Ceulen:   Pi < 3.1415926535897932384648319947728276044

Ceulen:   iteration   36

Ceulen:   "inner polygon has 137438953472 sides"

Ceulen:   3.1415926535897932384623698068428372944 < Pi

Ceulen:   "outer polygon has 137438953472 sides"

Ceulen:   Pi < 3.1415926535897932384631905361528340644

Ceulen:   iteration   37

Ceulen:   "inner polygon has 274877906944 sides"

Ceulen:   3.1415926535897932384625749891703364869 < Pi

Ceulen:   "outer polygon has 274877906944 sides"

Ceulen:   Pi < 3.1415926535897932384627801714978356794

Ceulen:   iteration   38

Ceulen:   "inner polygon has 549755813888 sides"

Ceulen:   3.1415926535897932384626262847522112851 < Pi

Ceulen:   "outer polygon has 549755813888 sides"

Ceulen:   Pi < 3.1415926535897932384626775803340860832

Ceulen:   iteration   39

Ceulen:   "inner polygon has 1099511627776 sides"

Ceulen:   3.1415926535897932384626391086476799847 < Pi

Ceulen:   "outer polygon has 1099511627776 sides"

Ceulen:   Pi < 3.1415926535897932384626519325431486842

Ceulen:   iteration   40

Ceulen:   "inner polygon has 2199023255552 sides"

Ceulen:   3.1415926535897932384626423146215471595 < Pi

Ceulen:   "outer polygon has 2199023255552 sides"

Ceulen:   Pi < 3.1415926535897932384626455205954143344

Ceulen:   iteration   41

Ceulen:   "inner polygon has 4398046511104 sides"

Ceulen:   3.1415926535897932384626431161150139532 < Pi

Ceulen:   "outer polygon has 4398046511104 sides"

Ceulen:   Pi < 3.1415926535897932384626439176084807470

Ceulen:   iteration   42

Ceulen:   "inner polygon has 8796093022208 sides"

Ceulen:   3.1415926535897932384626433164883806516 < Pi

Ceulen:   "outer polygon has 8796093022208 sides"

Ceulen:   Pi < 3.1415926535897932384626435168617473500

Ceulen:   iteration   43

Ceulen:   "inner polygon has 17592186044416 sides"

Ceulen:   3.1415926535897932384626433665817223262 < Pi

Ceulen:   "outer polygon has 17592186044416 sides"

Ceulen:   Pi < 3.1415926535897932384626434166750640008

Ceulen:   iteration   44

Ceulen:   "inner polygon has 35184372088832 sides"

Ceulen:   3.1415926535897932384626433791050577449 < Pi

Ceulen:   "outer polygon has 35184372088832 sides"

Ceulen:   Pi < 3.1415926535897932384626433916283931636

Ceulen:   iteration   45

Ceulen:   "inner polygon has 70368744177664 sides"

Ceulen:   3.1415926535897932384626433822358915995 < Pi

Ceulen:   "outer polygon has 70368744177664 sides"

Ceulen:   Pi < 3.1415926535897932384626433853667254542

Ceulen:   iteration   46

Ceulen:   "inner polygon has 140737488355328 sides"

Ceulen:   3.1415926535897932384626433830186000631 < Pi

Ceulen:   "outer polygon has 140737488355328 sides"

Ceulen:   Pi < 3.1415926535897932384626433838013085268

Ceulen:   iteration   47

Ceulen:   "inner polygon has 281474976710656 sides"

Ceulen:   3.1415926535897932384626433832142771791 < Pi

Ceulen:   "outer polygon has 281474976710656 sides"

Ceulen:   Pi < 3.1415926535897932384626433834099542950

Ceulen:   iteration   48

Ceulen:   "inner polygon has 562949953421312 sides"

Ceulen:   3.1415926535897932384626433832631964580 < Pi

Ceulen:   "outer polygon has 562949953421312 sides"

Ceulen:   Pi < 3.1415926535897932384626433833121157370

Ceulen:   iteration   49

Ceulen:   "inner polygon has 1125899906842624 sides"

Ceulen:   3.1415926535897932384626433832754262777 < Pi

Ceulen:   "outer polygon has 1125899906842624 sides"

Ceulen:   Pi < 3.1415926535897932384626433832876560974

Ceulen:   iteration   50

Ceulen:   "inner polygon has 2251799813685248 sides"

Ceulen:   3.1415926535897932384626433832784837327 < Pi

Ceulen:   "outer polygon has 2251799813685248 sides"

Ceulen:   Pi < 3.1415926535897932384626433832815411876

Ceulen:   iteration   51

Ceulen:   "inner polygon has 4503599627370496 sides"

Ceulen:   3.1415926535897932384626433832792480965 < Pi

Ceulen:   "outer polygon has 4503599627370496 sides"

Ceulen:   Pi < 3.1415926535897932384626433832800124602

Ceulen:   iteration   52

Ceulen:   "inner polygon has 9007199254740992 sides"

Ceulen:   3.1415926535897932384626433832794391874 < Pi

Ceulen:   "outer polygon has 9007199254740992 sides"

Ceulen:   Pi < 3.1415926535897932384626433832796302784

Ceulen:   iteration   53

Ceulen:   "inner polygon has 18014398509481984 sides"

Ceulen:   3.1415926535897932384626433832794869601 < Pi

Ceulen:   "outer polygon has 18014398509481984 sides"

Ceulen:   Pi < 3.1415926535897932384626433832795347328

Ceulen:   iteration   54

Ceulen:   "inner polygon has 36028797018963968 sides"

Ceulen:   3.1415926535897932384626433832794989033 < Pi

Ceulen:   "outer polygon has 36028797018963968 sides"

Ceulen:   Pi < 3.1415926535897932384626433832795108464

Ceulen:   iteration   55

Ceulen:   "inner polygon has 72057594037927936 sides"

Ceulen:   3.1415926535897932384626433832795018890 < Pi

Ceulen:   "outer polygon has 72057594037927936 sides"

Ceulen:   Pi < 3.1415926535897932384626433832795048748

Ceulen:   iteration   56

Ceulen:   "inner polygon has 144115188075855872 sides"

Ceulen:   3.1415926535897932384626433832795026355 < Pi

Ceulen:   "outer polygon has 144115188075855872 sides"

Ceulen:   Pi < 3.1415926535897932384626433832795033820

Ceulen:   iteration   57

Ceulen:   "inner polygon has 288230376151711744 sides"

Ceulen:   3.1415926535897932384626433832795028222 < Pi

Ceulen:   "outer polygon has 288230376151711744 sides"

Ceulen:   Pi < 3.1415926535897932384626433832795030088

Ceulen:   iteration   58

Ceulen:   "inner polygon has 576460752303423488 sides"

Ceulen:   3.1415926535897932384626433832795028689 < Pi

Ceulen:   "outer polygon has 576460752303423488 sides"

Ceulen:   Pi < 3.1415926535897932384626433832795029156

Ceulen:   iteration   59

Ceulen:   "inner polygon has 1152921504606846976 sides"

Ceulen:   3.1415926535897932384626433832795028806 < Pi

Ceulen:   "outer polygon has 1152921504606846976 sides"

Ceulen:   Pi < 3.1415926535897932384626433832795028922

Ceulen:   iteration   60

Ceulen:   "inner polygon has 2305843009213693952 sides"

Ceulen:   3.1415926535897932384626433832795028835 < Pi

Ceulen:   "outer polygon has 2305843009213693952 sides"

Ceulen:   Pi < 3.1415926535897932384626433832795028864

Ceulen:   iteration   61

Ceulen:   "inner polygon has 4611686018427387904 sides"

Ceulen:   3.1415926535897932384626433832795028842 < Pi

Ceulen:   "outer polygon has 4611686018427387904 sides"

Ceulen:   Pi < 3.1415926535897932384626433832795028850

Ceulen:   61 iterations

> Maple_pi := evalf(Pi);

> Cpi_error := Cpi-Maple_pi;

> infolevel[Ceulen] := 0;

> for Digits from 100 by 100 to 500 do
st := time(): CPi := Ceulen(infinity): time_CPi[Digits] := time()-st;
lprint(Digits,time_CPi[Digits]);
od:

100, .196

200, .907

300, 2.003

400, 3.766

500, 6.891

> Digits := 10;

>

Gregory's series for arctan(x)

In 1668 James Gregory used found the following series for arctan(x)

> Gregory_series := arctan(x)=Sum((-1)^n/(2*n+1)*x^(2*n+1),n=0..infinity);

This series may derived by starting with the power series for (1/(1+x^2)) = 1-x^2+x^4-x^6+x^8-x^10+... and integrating term by term on the right side and using the known integral int(1/(1+x^2),x)=arctan(x) on the left side,
If we let x=1 we find a formula for Pi.
Pi/4 = 1 - 1/3 + 1/5 -1/7 + 1/9 - 1/11 + ...

> arctan(1)=Sum((-1)^n/(2*n+1),n=0..infinity);

>    Gregory := proc(n)
local i,s,t,u;
s := -1;
u := -1.;
t := 0.;
for i to n do
  u := u+2;
  s := -s;
  t := t+s/u;
od;
4*t;
end;

> e := 1; st := time(): Gpi||e := evalhf(Gregory(10^e)): timeGpi||e := time()-st; Gpi||e;

> e := 2; st := time(): Gpi||e := evalhf(Gregory(10^e)): timeGpi||e := time()-st; Gpi||e;

> e := 3; st := time(): Gpi||e := evalhf(Gregory(10^e)): timeGpi||e := time()-st; Gpi||e;

> e := 4; st := time(): Gpi||e := evalhf(Gregory(10^e)): timeGpi||e := time()-st; Gpi||e;

> e := 5; st := time(): Gpi||e := evalhf(Gregory(10^e)): timeGpi||e := time()-st; Gpi||e;

> e := 6; st := time(): Gpi||e := evalhf(Gregory(10^e)): timeGpi||e := time()-st; Gpi||e;

> evalf(Pi,18);

>

Newton's Calculation

In 1680 Isaac Newton used the following formula to compute Pi to 17 digits (16 correct)

> Newton_formula := Pi/24=Int(x^(1/2)*(1-x)^(1/2),x=0..1/4)+sqrt(3)/32;

The following Maple commands generate a picture of Newton's formula. The area of the red curved triangle represents the integral Int(x^(1/2)*(1-x)^(1/2),x=0..1/4) and the blue triangle has area sqrt(3)/4*(1/4)*(1/2)=sqrt(3)/32 . The sum of these two is 1/6 of a circle with radius=1/2 which has area 1/6*Pi*(1/2)^2=Pi/24

> cir1 := plot([.5+.5*cos(t),.5*sin(t),t=-Pi..Pi],numpoints=361,colour=green):

> Nt := 40;

> pi3N := evalf(Pi/Nt/3);

> cur_tri1 := plots[polygonplot]([[0.,0.],[.25,0.],seq([.5+.5*cos(pi3N*t),.5*sin(pi3N*t)],t=2*Nt..3*Nt)],colour=red):

> y1 := sqrt(3.)/4;

> tri1 := plots[polygonplot]([[.25,0.],[.5,0.],[.25,y1]],colour=blue):

> Nt := 40;

> plots[display]({cir1,cur_tri1,tri1},scaling=constrained);

> interface(verboseproc=0);

>    Newton := proc()
local d,one,p4,osq3,sq3,ti1,ti2,ti3,c,i,ct,pi;
description "compute Pi using Newton's formula",
"Pi/24=Int(x^(1/2)*(1-x)^(1/2),x=0..1/4)+sqrt(3)/32";
d := Digits;
one := 10^d;
p4 := one;
userinfo(1,Newton,1/24*Pi=
  Int(x^(1/2)*(1-x)^(1/2),x=0..1/4)+1/32*sqrt(3));
userinfo(2,Newton,"use Newton's binomial series");
userinfo(2,Newton,(1-x)^(1/2)=
  Sum(binomial(1/2,n)*(-x)^n,n=0..infinity));
userinfo(2,Newton,
  "use above for the integral and sqrt(1-1/4)=sqrt(3)/2");
userinfo(3,Newton,"use integer arithmetic");
osq3 := 0;
sq3 := one;
userinfo(3,Newton,1/2*sqrt(3),'term',1,sq3);
ti1 := -3;
ti2 := 0;
ti3 := 3;
c := iquo(2*one,3);
userinfo(3,Newton,'integral','term',1,c);
for i while sq3<>osq3 do
  osq3 := sq3;
  p4 := iquo(p4,4);
  userinfo(5,Newton,"divide by 4",p4);
  ti1 := ti1+2;
  ti2 := ti2+2;
  p4 := iquo(ti1*p4,ti2);
  userinfo(4,Newton,1/2*sqrt(3),'term',i+1,p4);
  sq3 := sq3+p4;
  userinfo(4,Newton,1/2*sqrt(3),'sum',p4);
  ti3 := ti3+2;
  ct := iquo(2*p4,ti3);
  userinfo(4,Newton,'integral','term',i+1,ct);
  c := c+ct;
  userinfo(4,Newton,'integral','sum',c);
od;
sq3 := iquo(sq3,2);
userinfo(3,Newton,1/4*sqrt(3)=Float(sq3,-Digits));
userinfo(3,Newton,
  8*Int(x^(1/2)*(1-x)^(1/2),x=0..1/4)=Float(c,-Digits));
pi := c+sq3;
pi := 3*pi;
pi := Float(pi,-d);
userinfo(2,Newton,Pi=pi);
evalf(pi);
end;

> infolevel[Newton] := 5;

> Digits := 17;

> Newton_pi := Newton();

Newton:   1/24*Pi = Int(x^(1/2)*(1-x)^(1/2),x = 0 .. 1/4)+1/32*3^(1/2)

Newton:   "use Newton's binomial series"

Newton:   (1-x)^(1/2) = Sum(binomial(1/2,n)*(-x)^n,n = 0 .. infinity)

Newton:   "use above for the integral and sqrt(1-1/4)=sqrt(3)/2"

Newton:   "use integer arithmetic"

Newton:   1/2*3^(1/2)   term   1   100000000000000000

Newton:   integral   term   1   66666666666666666

Newton:   "divide by 4"   25000000000000000

Newton:   1/2*3^(1/2)   term   2   -12500000000000000

Newton:   1/2*3^(1/2)   sum   -12500000000000000

Newton:   integral   term   2   -5000000000000000

Newton:   integral   sum   61666666666666666

Newton:   "divide by 4"   -3125000000000000

Newton:   1/2*3^(1/2)   term   3   -781250000000000

Newton:   1/2*3^(1/2)   sum   -781250000000000

Newton:   integral   term   3   -223214285714285

Newton:   integral   sum   61443452380952381

Newton:   "divide by 4"   -195312500000000

Newton:   1/2*3^(1/2)   term   4   -97656250000000

Newton:   1/2*3^(1/2)   sum   -97656250000000

Newton:   integral   term   4   -21701388888888

Newton:   integral   sum   61421750992063493

Newton:   "divide by 4"   -24414062500000

Newton:   1/2*3^(1/2)   term   5   -15258789062500

Newton:   1/2*3^(1/2)   sum   -15258789062500

Newton:   integral   term   5   -2774325284090

Newton:   integral   sum   61418976666779403

Newton:   "divide by 4"   -3814697265625

Newton:   1/2*3^(1/2)   term   6   -2670288085937

Newton:   1/2*3^(1/2)   sum   -2670288085937

Newton:   integral   term   6   -410813551682

Newton:   integral   sum   61418565853227721

Newton:   "divide by 4"   -667572021484

Newton:   1/2*3^(1/2)   term   7   -500679016113

Newton:   1/2*3^(1/2)   sum   -500679016113

Newton:   integral   term   7   -66757202148

Newton:   integral   sum   61418499096025573

Newton:   "divide by 4"   -125169754028

Newton:   1/2*3^(1/2)   term   8   -98347663879

Newton:   1/2*3^(1/2)   sum   -98347663879

Newton:   integral   term   8   -11570313397

Newton:   integral   sum   61418487525712176

Newton:   "divide by 4"   -24586915969

Newton:   1/2*3^(1/2)   term   9   -19976869224

Newton:   1/2*3^(1/2)   sum   -19976869224

Newton:   integral   term   9   -2102828339

Newton:   integral   sum   61418485422883837

Newton:   "divide by 4"   -4994217306

Newton:   1/2*3^(1/2)   term   10   -4161847755

Newton:   1/2*3^(1/2)   sum   -4161847755

Newton:   integral   term   10   -396366452

Newton:   integral   sum   61418485026517385

Newton:   "divide by 4"   -1040461938

Newton:   1/2*3^(1/2)   term   11   -884392647

Newton:   1/2*3^(1/2)   sum   -884392647

Newton:   integral   term   11   -76903708

Newton:   integral   sum   61418484949613677

Newton:   "divide by 4"   -221098161

Newton:   1/2*3^(1/2)   term   12   -190948411

Newton:   1/2*3^(1/2)   sum   -190948411

Newton:   integral   term   12   -15275872

Newton:   integral   sum   61418484934337805

Newton:   "divide by 4"   -47737102

Newton:   1/2*3^(1/2)   term   13   -41769964

Newton:   1/2*3^(1/2)   sum   -41769964

Newton:   integral   term   13   -3094071

Newton:   integral   sum   61418484931243734

Newton:   "divide by 4"   -10442491

Newton:   1/2*3^(1/2)   term   14   -9237588

Newton:   1/2*3^(1/2)   sum   -9237588

Newton:   integral   term   14   -637075

Newton:   integral   sum   61418484930606659

Newton:   "divide by 4"   -2309397

Newton:   1/2*3^(1/2)   term   15   -2061961

Newton:   1/2*3^(1/2)   sum   -2061961

Newton:   integral   term   15   -133029

Newton:   integral   sum   61418484930473630

Newton:   "divide by 4"   -515490

Newton:   1/2*3^(1/2)   term   16   -463941

Newton:   1/2*3^(1/2)   sum   -463941

Newton:   integral   term   16   -28117

Newton:   integral   sum   61418484930445513

Newton:   "divide by 4"   -115985

Newton:   1/2*3^(1/2)   term   17   -105111

Newton:   1/2*3^(1/2)   sum   -105111

Newton:   integral   term   17   -6006

Newton:   integral   sum   61418484930439507

Newton:   "divide by 4"   -26277

Newton:   1/2*3^(1/2)   term   18   -23958

Newton:   1/2*3^(1/2)   sum   -23958

Newton:   integral   term   18   -1295

Newton:   integral   sum   61418484930438212

Newton:   "divide by 4"   -5989

Newton:   1/2*3^(1/2)   term   19   -5489

Newton:   1/2*3^(1/2)   sum   -5489

Newton:   integral   term   19   -281

Newton:   integral   sum   61418484930437931

Newton:   "divide by 4"   -1372

Newton:   1/2*3^(1/2)   term   20   -1263

Newton:   1/2*3^(1/2)   sum   -1263

Newton:   integral   term   20   -61

Newton:   integral   sum   61418484930437870

Newton:   "divide by 4"   -315

Newton:   1/2*3^(1/2)   term   21   -291

Newton:   1/2*3^(1/2)   sum   -291

Newton:   integral   term   21   -13

Newton:   integral   sum   61418484930437857

Newton:   "divide by 4"   -72

Newton:   1/2*3^(1/2)   term   22   -66

Newton:   1/2*3^(1/2)   sum   -66

Newton:   integral   term   22   -2

Newton:   integral   sum   61418484930437855

Newton:   "divide by 4"   -16

Newton:   1/2*3^(1/2)   term   23   -14

Newton:   1/2*3^(1/2)   sum   -14

Newton:   integral   term   23   0

Newton:   integral   sum   61418484930437855

Newton:   "divide by 4"   -3

Newton:   1/2*3^(1/2)   term   24   -2

Newton:   1/2*3^(1/2)   sum   -2

Newton:   integral   term   24   0

Newton:   integral   sum   61418484930437855

Newton:   "divide by 4"   0

Newton:   1/2*3^(1/2)   term   25   0

Newton:   1/2*3^(1/2)   sum   0

Newton:   integral   term   25   0

Newton:   integral   sum   61418484930437855

Newton:   1/4*3^(1/2) = .43301270189221943

Newton:   8*Int(x^(1/2)*(1-x)^(1/2),x = 0 .. 1/4) = .61418484930437855

Newton:   Pi = 3.14159265358979394

> Maple_pi := evalf(Pi);

> Newton_error := Newton_pi-Maple_pi;

> infolevel[Newton] := 0;

> for Digits from 100 by 100 to 1000 do
st := time(): Npi := Newton(): time_Npi[Digits] := time()-st;
lprint(Digits,time_Npi[Digits]);
od:

100, .45e-1

200, .104

300, .185

400, .648

500, 1.347

600, .993

700, 1.214

800, 1.402

900, 1.892

1000, 1.985

> Digits := 10;

>

Sharp's Calculation

In 1699 Abraham Sharp used Halley's formula

> Halley_formula := arctan(1/sqrt(3))='arctan'(1/sqrt(3));

and Gregory's series for arctan(x)

> Gregory_series := arctan(x)=Sum((-1)^n*x^(2*n+1)/(2*n+1),n=0..infinity);

to compute to 72 (71 correct) decimal places.
He factored out the sqrt(3) from each term, so he only needed to compute sqrt(3) once.

> Sharp_sum := Pi/6=sqrt(3)/3*Sum((-1)^n/(2*n+1)/3^n,n=0..infinity);

>    Sharp := proc()
local t,p3,v,s,j,i;
description "compute Pi=6*arctan(1/sqrt(3)) as Sharp did";
t := 3.;
p3 := sqrt(t);
userinfo(1,Sharp,sqrt(3)=p3);
p3 := p3/t;
userinfo(2,Sharp,1/3*sqrt(3)=p3);
userinfo(3,Sharp,'using','integer','arithmetic');
p3 := op(1,p3)*10^(Digits+op(2,p3));
userinfo(3,Sharp,'term',1,p3);
v := p3;
t := p3;
s := 1;
j := 1;
for i while v<>0 do
  p3 := iquo(p3,3);
  userinfo(4,Sharp,power(3,-i-1/2)=p3);
  j := j+2;
  v := iquo(p3,j);
  s := -s;
  v := s*v;
  userinfo(3,Sharp,'term',i+1,v);
  t := t+v;
  userinfo(4,Sharp,'total',t);
od;
t := 6*t;
userinfo(2,Sharp,'times',6,t);
t := Float(t,-Digits);
evalf(t,Digits);
end;

> infolevel[Sharp] := 4;

> Digits := 72;

> Sharp_pi := Sharp();

Sharp:   3^(1/2) = 1.73205080756887729352744634150587236694280525381038062805580697945193302

Sharp:   1/3*3^(1/2) = .577350269189625764509148780501957455647601751270126876018602326483977673

Sharp:   using   integer   arithmetic

Sharp:   term   1   577350269189625764509148780501957455647601751270126876018602326483977673

Sharp:   power(3,-3/2) = 192450089729875254836382926833985818549200583756708958672867442161325891

Sharp:   term   2   -64150029909958418278794308944661939516400194585569652890955814053775297

Sharp:   total   513200239279667346230354471557295516131201556684557223127646512430202376

Sharp:   power(3,-5/2) = 64150029909958418278794308944661939516400194585569652890955814053775297

Sharp:   term   3   12830005981991683655758861788932387903280038917113930578191162810755059

Sharp:   total   526030245261659029886113333346227904034481595601671153705837675240957435

Sharp:   power(3,-7/2) = 21383343303319472759598102981553979838800064861856550963651938017925099

Sharp:   term   4   -3054763329045638965656871854507711405542866408836650137664562573989299

Sharp:   total   522975481932613390920456461491720192628938729192834503568173112666968136

Sharp:   power(3,-9/2) = 7127781101106490919866034327184659946266688287285516987883979339308366

Sharp:   term   5   791975677900721213318448258576073327362965365253946331987108815478707

Sharp:   total   523767457610514112133774909750296265956301694558088449900160221482446843

Sharp:   power(3,-11/2) = 2375927033702163639955344775728219982088896095761838995961326446436122

Sharp:   term   6   -215993366700196694541394979611656362008081463251076272360120586039647

Sharp:   total   523551464243813915439233514770684609594293613094837373627800100896407196

Sharp:   power(3,-13/2) = 791975677900721213318448258576073327362965365253946331987108815478707

Sharp:   term   7   60921205992363170255265250659697948258689643481072794768239139652208

Sharp:   total   523612385449806278609488780021344307542552302738318446422568340036059404

Sharp:   power(3,-15/2) = 263991892633573737772816086192024442454321788417982110662369605159569

Sharp:   term   8   -17599459508904915851521072412801629496954785894532140710824640343971

Sharp:   total   523594785990297373693637258948931505913055347952423914281857515395715433

Sharp:   power(3,-17/2) = 87997297544524579257605362064008147484773929472660703554123201719856

Sharp:   term   9   5176311620266151721035609533176949852045525263097688444360188336462

Sharp:   total   523599962301917639845358294558464682862907393477687011970301875584051895

Sharp:   power(3,-19/2) = 29332432514841526419201787354669382494924643157553567851374400573285

Sharp:   term   10   -1543812237623238232589567755508914868153928587239661465861810556488

Sharp:   total   523598418489680016607125704990709173948039239549099772308836013773495407

Sharp:   power(3,-21/2) = 9777477504947175473067262451556460831641547719184522617124800191095

Sharp:   term   11   465594166902246451098441069121736230078168939008786791291657151956

Sharp:   total   523598884083846918853576803431778295684269317718038781095627305430647363

Sharp:   power(3,-23/2) = 3259159168315725157689087483852153610547182573061507539041600063698

Sharp:   term   12   -141702572535466311203873368863137113502051416220065545175721741899

Sharp:   total   523598742381274383387265599558409432547155815666622561030082129708905464

Sharp:   power(3,-25/2) = 1086386389438575052563029161284051203515727524353835846347200021232

Sharp:   term   13   43455455577543002102521166451362048140629100974153433853888000849

Sharp:   total   523598785836729960930267702079575883909203956295723535183515983596906313

Sharp:   power(3,-27/2) = 362128796479525017521009720428017067838575841451278615449066673744

Sharp:   term   14   -13412177647389815463741100756593224734762068201899207979595061990

Sharp:   total   523598772424552313540452238338475127315979221533655333284308004001844323

Sharp:   power(3,-29/2) = 120709598826508339173669906809339022612858613817092871816355557914

Sharp:   term   15   4162399959534770316333445062391000779753745304037685235046743376

Sharp:   total   523598776586952273075222554671920189706980001287400637321993239048587699

Sharp:   power(3,-31/2) = 40236532942169446391223302269779674204286204605697623938785185971

Sharp:   term   16   -1297952675553853109394300073218699167880200148570891094799522128

Sharp:   total   523598775288999597521369445277620116488280833407200488751102144249065571

Sharp:   power(3,-33/2) = 13412177647389815463741100756593224734762068201899207979595061990

Sharp:   term   17   406429625678479256477003053230097719235214187936339635745304908

Sharp:   total   523598775695429223199848701754623169718378552642414676687441779994370479

Sharp:   power(3,-35/2) = 4470725882463271821247033585531074911587356067299735993198353996

Sharp:   term   18   -127735025213236337749915245300887854616781601922849599805667257

Sharp:   total   523598775567694197986612364004707924417490698025633074764592180188703222

Sharp:   power(3,-37/2) = 1490241960821090607082344528510358303862452022433245331066117998

Sharp:   term   19   40276809751921367758982284554334008212498703309006630569354540

Sharp:   total   523598775607971007738533731763690208971824706238131778073598810758057762

Sharp:   power(3,-39/2) = 496747320273696869027448176170119434620817340811081777022039332

Sharp:   term   20   -12737110776248637667370466055644088067200444636181584026206136

Sharp:   total   523598775595233896962285094096319742916180618170931333437417226731851626

Sharp:   power(3,-41/2) = 165582440091232289675816058723373144873605780270360592340679777

Sharp:   term   21   4038596099786153406727220944472515728624531226106355910748287

Sharp:   total   523598775599272493062071247503046963860653133899555864663523582642599913

Sharp:   power(3,-43/2) = 55194146697077429891938686241124381624535260090120197446893259

Sharp:   term   22   -1283584806908777439347411307933125154058959536979539475509145

Sharp:   total   523598775597988908255162470063699552552720008745496905126544043167090768

Sharp:   power(3,-45/2) = 18398048899025809963979562080374793874845086696706732482297753

Sharp:   term   23   408845531089462443643990268452773197218779704371260721828838

Sharp:   total   523598775598397753786251932507343542821172781942715684830915303888919606

Sharp:   power(3,-47/2) = 6132682966341936654659854026791597958281695565568910827432584

Sharp:   term   24   -130482616305147588397018170782799956559185012033381081434735

Sharp:   total   523598775598267271169946784918946524650389981986156499818881922807484871

Sharp:   power(3,-49/2) = 2044227655447312218219951342263865986093898521856303609144194

Sharp:   term   25   41718931743822698331019415148242162981508133099108236921310

Sharp:   total   523598775598308990101690607617277544065538224149138007951981031044406181

Sharp:   power(3,-51/2) = 681409218482437406073317114087955328697966173952101203048064

Sharp:   term   26   -13360965068283086393594453217410888797999336744158847118589

Sharp:   total   523598775598295629136622324530883949612320813260340008615236872197287592

Sharp:   power(3,-53/2) = 227136406160812468691105704695985109565988724650700401016021

Sharp:   term   27   4285592569071933371530296315018586972943183483975479264453

Sharp:   total   523598775598299914729191396464255479908635831847312951798720847676552045

Sharp:   power(3,-55/2) = 75712135386937489563701901565328369855329574883566800338673

Sharp:   term   28   -1376584279762499810249125483005970361005992270610305460703

Sharp:   total   523598775598298538144911633964445230783152825876951945806450237371091342

Sharp:   power(3,-57/2) = 25237378462312496521233967188442789951776524961188933446224

Sharp:   term   29   442761025654605202126911705060399823715377630898051463968

Sharp:   total   523598775598298980905937288569647357694857886276775661184081135422555310

Sharp:   power(3,-59/2) = 8412459487437498840411322396147596650592174987062977815408

Sharp:   term   30   -142584059109110149837480040612671129671053813340050471447

Sharp:   total   523598775598298838321878179459497520214817273605645990130267795372083863

Sharp:   power(3,-61/2) = 2804153162479166280137107465382532216864058329020992605136

Sharp:   term   31   45969723975068299674378810907910364210886202115098239428

Sharp:   total   523598775598298884291602154527797194593628181516010201016469910470323291

Sharp:   power(3,-63/2) = 934717720826388760045702488460844072288019443006997535045

Sharp:   term   32   -14836789219466488254693690293029270988698721317571389445

Sharp:   total   523598775598298869454812935061308939899937888486739212317748592898933846

Sharp:   power(3,-65/2) = 311572573608796253348567496153614690762673147668999178348

Sharp:   term   33   4793424209366096205362576863901764473271894579523064282

Sharp:   total   523598775598298874248237144427405145262514752388503685589643172421998128

Sharp:   power(3,-67/2) = 103857524536265417782855832051204896920891049222999726116

Sharp:   term   34   -1550112306511424146012773612704550700311806704820891434

Sharp:   total   523598775598298872698124837915980999249741139683952985277836467601106694

Sharp:   power(3,-69/2) = 34619174845421805927618610683734965640297016407666575372

Sharp:   term   35   501727171672779796052443633097608197685464005908211237

Sharp:   total   523598775598298873199852009588760795302184772781561182963300473509317931

Sharp:   power(3,-71/2) = 11539724948473935309206203561244988546765672135888858457

Sharp:   term   36   -162531337302449793087411317764013923193882706139279696

Sharp:   total   523598775598298873037320672286311002214773455017547259769417767370038235

Sharp:   power(3,-73/2) = 3846574982824645103068734520414996182255224045296286152

Sharp:   term   37   52692807983899247987242938635821865510345534867072413

Sharp:   total   523598775598298873090013480270210250202016393653369125279763302237110648

Sharp:   power(3,-75/2) = 1282191660941548367689578173471665394085074681765428717

Sharp:   term   38   -17095888812553978235861042312955538587800995756872382

Sharp:   total   523598775598298873072917591457656271966155351340413586691962306480238266

Sharp:   power(3,-77/2) = 427397220313849455896526057823888464695024893921809572

Sharp:   term   39   5550613250829213712941896854855694346688634985997526

Sharp:   total   523598775598298873078468204708485485679097248195269281038650941466235792

Sharp:   power(3,-79/2) = 142465740104616485298842019274629488231674964640603190

Sharp:   term   40   -1803363798792613737960025560438347952299683096716496

Sharp:   total   523598775598298873076664840909692871941137222634830933086351258369519296

Sharp:   power(3,-81/2) = 47488580034872161766280673091543162743891654880201063

Sharp:   term   41   586278765862619281065193494957322996838168578767914

Sharp:   total   523598775598298873077251119675555491222202416129788256083189426948287210

Sharp:   power(3,-83/2) = 15829526678290720588760224363847720914630551626733687

Sharp:   term   42   -190717188895068922756147281492141215838922308755827

Sharp:   total   523598775598298873077060402486660422299446268848296114867350504639531383

Sharp:   power(3,-85/2) = 5276508892763573529586741454615906971543517208911229

Sharp:   term   43   62076575208983217995138134760187140841688437751896

Sharp:   total   523598775598298873077122479061869405517441406983056302008192193077283279

Sharp:   power(3,-87/2) = 1758836297587857843195580484871968990514505736303743

Sharp:   term   44   -20216509167676526933282534308873206787523054440272

Sharp:   total   523598775598298873077102262552701728990508124448747428801404670022843007

Sharp:   power(3,-89/2) = 586278765862619281065193494957322996838168578767914

Sharp:   term   45   6587401863624935742305544887160932548743467177167

Sharp:   total   523598775598298873077108849954565353926250429993634589733953413490020174

Sharp:   power(3,-91/2) = 195426255287539760355064498319107665612722859589304

Sharp:   term   46   -2147541266896041322583126355155029292447503951530

Sharp:   total   523598775598298873077106702413298457884927846867279434704660965986068644

Sharp:   power(3,-93/2) = 65142085095846586785021499439702555204240953196434

Sharp:   term   47   700452527912328890161521499351640378540225303187

Sharp:   total   523598775598298873077107402865826370213818008388778786345039506211371831

Sharp:   power(3,-95/2) = 21714028365282195595007166479900851734746984398811

Sharp:   term   48   -228568719634549427315864910314745807734178783145

Sharp:   total   523598775598298873077107174297106735664390692523868471599231772032588686

Sharp:   power(3,-97/2) = 7238009455094065198335722159966950578248994799603

Sharp:   term   49   74618654176227476271502290308937634827309224738

Sharp:   total   523598775598298873077107248915760911891866964026158780536866599341813424

Sharp:   power(3,-99/2) = 2412669818364688399445240719988983526082998266534

Sharp:   term   50   -24370402205703923226719603232211954808919174409

Sharp:   total   523598775598298873077107224545358706187943737306555548324911790422639015

Sharp:   power(3,-101/2) = 804223272788229466481746906662994508694332755511

Sharp:   term   51   7962606661269598678037098085772222858359730252

Sharp:   total   523598775598298873077107232507965367457542415343653634097134648782369267

Sharp:   power(3,-103/2) = 268074424262743155493915635554331502898110918503

Sharp:   term   52   -2602664313230516072756462481110014591243795325

Sharp:   total   523598775598298873077107229905301054227026342587191152987120057538573942

Sharp:   power(3,-105/2) = 89358141420914385164638545184777167632703639501

Sharp:   term   53   851029918294422715853700430331211120311463233

Sharp:   total   523598775598298873077107230756330972521449058440891583318331177850037175

Sharp:   power(3,-107/2) = 29786047140304795054879515061592389210901213167

Sharp:   term   54   -278374272339297150045602944500863450569170216

Sharp:   total   523598775598298873077107230477956700182151908395288638817467727280866959

Sharp:   power(3,-109/2) = 9928682380101598351626505020530796403633737722

Sharp:   term   55   91088829175244021574555091931475196363612272

Sharp:   total   523598775598298873077107230569045529357395929969843730748942923644479231

Sharp:   power(3,-111/2) = 3309560793367199450542168340176932134544579240

Sharp:   term   56   -29815863003308103158037552614206595806707921

Sharp:   total   523598775598298873077107230539229666354087826811806178134736327837771310

Sharp:   power(3,-113/2) = 1103186931122399816847389446725644044848193080

Sharp:   term   57   9762716204623007228737959705536672963258345

Sharp:   total   523598775598298873077107230548992382558710834040544137840273000801029655

Sharp:   power(3,-115/2) = 367728977040799938949129815575214681616064360

Sharp:   term   58   -3197643278615651643035911439784475492313603

Sharp:   total   523598775598298873077107230545794739280095182397508226400488525308716052

Sharp:   power(3,-117/2) = 122576325680266646316376605191738227205354786

Sharp:   term   59   1047660903250142276208347052920839548763716

Sharp:   total   523598775598298873077107230546842400183345324673716573453409364857479768

Sharp:   power(3,-119/2) = 40858775226755548772125535063912742401784928

Sharp:   term   60   -343351052325676880438029706419434810099033

Sharp:   total   523598775598298873077107230546499049131019647793278543746989930047380735

Sharp:   power(3,-121/2) = 13619591742251849590708511687970914133928309

Sharp:   term   61   112558609440097930501723237090668711850647

Sharp:   total   523598775598298873077107230546611607740459745723780266984080598759231382

Sharp:   power(3,-123/2) = 4539863914083949863569503895990304711309436

Sharp:   term   62   -36909462716129673687556942243823615539101

Sharp:   total   523598775598298873077107230546574698277743616050092710041836775143692281

Sharp:   power(3,-125/2) = 1513287971361316621189834631996768237103145

Sharp:   term   63   12106303770890532969518677055974145896825

Sharp:   total   523598775598298873077107230546586804581514506583062228718892749289589106

Sharp:   power(3,-127/2) = 504429323787105540396611543998922745701048

Sharp:   term   64   -3971884439268547562178043653534824769299

Sharp:   total   523598775598298873077107230546582832697075238035500050675239214464819807

Sharp:   power(3,-129/2) = 168143107929035180132203847999640915233682

Sharp:   term   65   1303434945186319225831037581392565234369

Sharp:   total   523598775598298873077107230546584136132020424354725881712820607030054176

Sharp:   power(3,-131/2) = 56047702643011726710734615999880305077894

Sharp:   term   66   -427845058343600967257516152670842023495

Sharp:   total   523598775598298873077107230546583708286962080753758624196667936188030681

Sharp:   power(3,-133/2) = 18682567547670575570244871999960101692631

Sharp:   term   67   140470432689252447896577984962106027764

Sharp:   total   523598775598298873077107230546583848757394770006206520774652898294058445

Sharp:   power(3,-135/2) = 6227522515890191856748290666653367230877

Sharp:   term   68   -46129796414001421161098449382617535043

Sharp:   total   523598775598298873077107230546583802627598356004785359676203515676523402

Sharp:   power(3,-137/2) = 2075840838630063952249430222217789076959

Sharp:   term   69   15152122909708496001820658556334226839

Sharp:   total   523598775598298873077107230546583817779721265713281361496862072010750241

Sharp:   power(3,-139/2) = 691946946210021317416476740739263025653

Sharp:   term   70   -4978035584244757679255228350642180040

Sharp:   total   523598775598298873077107230546583812801685681468523682241633721368570201

Sharp:   power(3,-141/2) = 230648982070007105805492246913087675217

Sharp:   term   71   1635808383475227700748171963922607625

Sharp:   total   523598775598298873077107230546583814437494064943751382989805685291177826

Sharp:   power(3,-143/2) = 76882994023335701935164082304362558405

Sharp:   term   72   -537643314848501412134014561568968939

Sharp:   total   523598775598298873077107230546583813899850750095249970855791123722208887

Sharp:   power(3,-145/2) = 25627664674445233978388027434787519468

Sharp:   term   73   176742514996174027437158809895086341

Sharp:   total   523598775598298873077107230546583814076593265091423998292949933617295228

Sharp:   power(3,-147/2) = 8542554891481744659462675811595839822

Sharp:   term   74   -58112618309399623533759699398611155

Sharp:   total   523598775598298873077107230546583814018480646782024374759190234218684073

Sharp:   power(3,-149/2) = 2847518297160581553154225270531946607

Sharp:   term   75   19110861054769003712444464902898970

Sharp:   total   523598775598298873077107230546583814037591507836793378471634699121583043

Sharp:   power(3,-151/2) = 949172765720193851051408423510648869

Sharp:   term   76   -6285912355762873185770916711991052

Sharp:   total   523598775598298873077107230546583814031305595481030505285863782409591991

Sharp:   power(3,-153/2) = 316390921906731283683802807836882956

Sharp:   term   77   2067914522266217540417011815927339

Sharp:   total   523598775598298873077107230546583814033373510003296722826280794225519330

Sharp:   power(3,-155/2) = 105463640635577094561267602612294318

Sharp:   term   78   -680410584745658674588823242659963

Sharp:   total   523598775598298873077107230546583814032693099418551064151691970982859367

Sharp:   power(3,-157/2) = 35154546878525698187089200870764772

Sharp:   term   79   223914311328189160427319750769202

Sharp:   total   523598775598298873077107230546583814032917013729879253312119290733628569

Sharp:   power(3,-159/2) = 11718182292841899395696400290254924

Sharp:   term   80   -73699259703408172299977360316068

Sharp:   total   523598775598298873077107230546583814032843314470175845139819313373312501

Sharp:   power(3,-161/2) = 3906060764280633131898800096751641

Sharp:   term   81   24261246983109522558377640352494

Sharp:   total   523598775598298873077107230546583814032867575717158954662377691013664995

Sharp:   power(3,-163/2) = 1302020254760211043966266698917213

Sharp:   term   82   -7987854323682276343351329441209

Sharp:   total   523598775598298873077107230546583814032859587862835272386034339684223786

Sharp:   power(3,-165/2) = 434006751586737014655422232972404

Sharp:   term   83   2630343949010527361548013533166

Sharp:   total   523598775598298873077107230546583814032862218206784282913395887697756952

Sharp:   power(3,-167/2) = 144668917195579004885140744324134

Sharp:   term   84   -866280941290892244821202061821

Sharp:   total   523598775598298873077107230546583814032861351925842992021151066495695131

Sharp:   power(3,-169/2) = 48222972398526334961713581441378

Sharp:   term   85   285343031943942810424340718588

Sharp:   total   523598775598298873077107230546583814032861637268874935963961490836413719

Sharp:   power(3,-171/2) = 16074324132842111653904527147126

Sharp:   term   86   -94001895513696559379558638287

Sharp:   total   523598775598298873077107230546583814032861543266979422267402111277775432

Sharp:   power(3,-173/2) = 5358108044280703884634842382375

Sharp:   term   87   30971722799310427078814117817

Sharp:   total   523598775598298873077107230546583814032861574238702221577829190091893249

Sharp:   power(3,-175/2) = 1786036014760234628211614127458

Sharp:   term   88   -10205920084344197875494937871

Sharp:   total   523598775598298873077107230546583814032861564032782137233631314596955378

Sharp:   power(3,-177/2) = 595345338253411542737204709152

Sharp:   term   89   3363532984482551088910761068

Sharp:   total   523598775598298873077107230546583814032861567396315121716182403507716446

Sharp:   power(3,-179/2) = 198448446084470514245734903050

Sharp:   term   90   -1108650536784751476233155882

Sharp:   total   523598775598298873077107230546583814032861566287664584931430927274560564

Sharp:   power(3,-181/2) = 66149482028156838081911634350

Sharp:   term   91   365466751536778110949788035

Sharp:   total   523598775598298873077107230546583814032861566653131336468209038224348599

Sharp:   power(3,-183/2) = 22049827342718946027303878116

Sharp:   term   92   -120490859796278393591824470

Sharp:   total   523598775598298873077107230546583814032861566532640476671930644632524129

Sharp:   power(3,-185/2) = 7349942447572982009101292705

Sharp:   term   93   39729418635529632481628609

Sharp:   total   523598775598298873077107230546583814032861566572369895307460277114152738

Sharp:   power(3,-187/2) = 2449980815857660669700430901

Sharp:   term   94   -13101501689078399303210860

Sharp:   total   523598775598298873077107230546583814032861566559268393618381877810941878

Sharp:   power(3,-189/2) = 816660271952553556566810300

Sharp:   term   95   4320953819854780722575715

Sharp:   total   523598775598298873077107230546583814032861566563589347438236658533517593

Sharp:   power(3,-191/2) = 272220090650851185522270100

Sharp:   term   96   -1425236076706027149331257

Sharp:   total   523598775598298873077107230546583814032861566562164111361530631384186336

Sharp:   power(3,-193/2) = 90740030216950395174090033

Sharp:   term   97   470155596979017591575596

Sharp:   total   523598775598298873077107230546583814032861566562634266958509648975761932

Sharp:   power(3,-195/2) = 30246676738983465058030011

Sharp:   term   98   -155111162764017769528359

Sharp:   total   523598775598298873077107230546583814032861566562479155795745631206233573

Sharp:   power(3,-197/2) = 10082225579661155019343337

Sharp:   term   99   51178810049041396037275

Sharp:   total   523598775598298873077107230546583814032861566562530334605794672602270848

Sharp:   power(3,-199/2) = 3360741859887051673114445

Sharp:   term   100   -16888150049683676749318

Sharp:   total   523598775598298873077107230546583814032861566562513446455744988925521530

Sharp:   power(3,-201/2) = 1120247286629017224371481

Sharp:   term   101   5573369585218991166027

Sharp:   total   523598775598298873077107230546583814032861566562519019825330207916687557

Sharp:   power(3,-203/2) = 373415762209672408123827

Sharp:   term   102   -1839486513348139941496

Sharp:   total   523598775598298873077107230546583814032861566562517180338816859776746061

Sharp:   power(3,-205/2) = 124471920736557469374609

Sharp:   term   103   607180101153938874998

Sharp:   total   523598775598298873077107230546583814032861566562517787518918013715621059

Sharp:   power(3,-207/2) = 41490640245519156458203

Sharp:   term   104   -200437875582218147141

Sharp:   total   523598775598298873077107230546583814032861566562517587081042431497473918

Sharp:   power(3,-209/2) = 13830213415173052152734

Sharp:   term   105   66173269929057665802

Sharp:   total   523598775598298873077107230546583814032861566562517653254312360555139720

Sharp:   power(3,-211/2) = 4610071138391017384244

Sharp:   term   106   -21848678381000082389

Sharp:   total   523598775598298873077107230546583814032861566562517631405633979555057331

Sharp:   power(3,-213/2) = 1536690379463672461414

Sharp:   term   107   7214508823773110147

Sharp:   total   523598775598298873077107230546583814032861566562517638620142803328167478

Sharp:   power(3,-215/2) = 512230126487890820471

Sharp:   term   108   -2382465704594841025

Sharp:   total   523598775598298873077107230546583814032861566562517636237677098733326453

Sharp:   power(3,-217/2) = 170743375495963606823

Sharp:   term   109   786835831778634132

Sharp:   total   523598775598298873077107230546583814032861566562517637024512930511960585

Sharp:   power(3,-219/2) = 56914458498654535607

Sharp:   term   110   -259883372139975048

Sharp:   total   523598775598298873077107230546583814032861566562517636764629558371985537

Sharp:   power(3,-221/2) = 18971486166218178535

Sharp:   term   111   85843828806417097

Sharp:   total   523598775598298873077107230546583814032861566562517636850473387178402634

Sharp:   power(3,-223/2) = 6323828722072726178

Sharp:   term   112   -28357976332164691

Sharp:   total   523598775598298873077107230546583814032861566562517636822115410846237943

Sharp:   power(3,-225/2) = 2107942907357575392

Sharp:   term   113   9368635143811446

Sharp:   total   523598775598298873077107230546583814032861566562517636831484045990049389

Sharp:   power(3,-227/2) = 702647635785858464

Sharp:   term   114   -3095364034298935

Sharp:   total   523598775598298873077107230546583814032861566562517636828388681955750454

Sharp:   power(3,-229/2) = 234215878595286154

Sharp:   term   115   1022776762424830

Sharp:   total   523598775598298873077107230546583814032861566562517636829411458718175284

Sharp:   power(3,-231/2) = 78071959531762051

Sharp:   term   116   -337973850786848

Sharp:   total   523598775598298873077107230546583814032861566562517636829073484867388436

Sharp:   power(3,-233/2) = 26023986510587350

Sharp:   term   117   111690929229988

Sharp:   total   523598775598298873077107230546583814032861566562517636829185175796618424

Sharp:   power(3,-235/2) = 8674662170195783

Sharp:   term   118   -36913456043386

Sharp:   total   523598775598298873077107230546583814032861566562517636829148262340575038

Sharp:   power(3,-237/2) = 2891554056731927

Sharp:   term   119   12200650028404

Sharp:   total   523598775598298873077107230546583814032861566562517636829160462990603442

Sharp:   power(3,-239/2) = 963851352243975

Sharp:   term   120   -4032850846209

Sharp:   total   523598775598298873077107230546583814032861566562517636829156430139757233

Sharp:   power(3,-241/2) = 321283784081325

Sharp:   term   121   1333127734777

Sharp:   total   523598775598298873077107230546583814032861566562517636829157763267492010

Sharp:   power(3,-243/2) = 107094594693775

Sharp:   term   122   -440718496682

Sharp:   total   523598775598298873077107230546583814032861566562517636829157322548995328

Sharp:   power(3,-245/2) = 35698198231258

Sharp:   term   123   145706931556

Sharp:   total   523598775598298873077107230546583814032861566562517636829157468255926884

Sharp:   power(3,-247/2) = 11899399410419

Sharp:   term   124   -48175706115

Sharp:   total   523598775598298873077107230546583814032861566562517636829157420080220769

Sharp:   power(3,-249/2) = 3966466470139

Sharp:   term   125   15929584217

Sharp:   total   523598775598298873077107230546583814032861566562517636829157436009804986

Sharp:   power(3,-251/2) = 1322155490046

Sharp:   term   126   -5267551753

Sharp:   total   523598775598298873077107230546583814032861566562517636829157430742253233

Sharp:   power(3,-253/2) = 440718496682

Sharp:   term   127   1741970342

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432484223575

Sharp:   power(3,-255/2) = 146906165560

Sharp:   term   128   -576102610

Sharp:   total   523598775598298873077107230546583814032861566562517636829157431908120965

Sharp:   power(3,-257/2) = 48968721853

Sharp:   term   129   190539773

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432098660738

Sharp:   power(3,-259/2) = 16322907284

Sharp:   term   130   -63022808

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432035637930

Sharp:   power(3,-261/2) = 5440969094

Sharp:   term   131   20846624

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432056484554

Sharp:   power(3,-263/2) = 1813656364

Sharp:   term   132   -6896031

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432049588523

Sharp:   power(3,-265/2) = 604552121

Sharp:   term   133   2281328

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432051869851

Sharp:   power(3,-267/2) = 201517373

Sharp:   term   134   -754746

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432051115105

Sharp:   power(3,-269/2) = 67172457

Sharp:   term   135   249711

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432051364816

Sharp:   power(3,-271/2) = 22390819

Sharp:   term   136   -82622

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432051282194

Sharp:   power(3,-273/2) = 7463606

Sharp:   term   137   27339

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432051309533

Sharp:   power(3,-275/2) = 2487868

Sharp:   term   138   -9046

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432051300487

Sharp:   power(3,-277/2) = 829289

Sharp:   term   139   2993

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432051303480

Sharp:   power(3,-279/2) = 276429

Sharp:   term   140   -990

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432051302490

Sharp:   power(3,-281/2) = 92143

Sharp:   term   141   327

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432051302817

Sharp:   power(3,-283/2) = 30714

Sharp:   term   142   -108

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432051302709

Sharp:   power(3,-285/2) = 10238

Sharp:   term   143   35

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432051302744

Sharp:   power(3,-287/2) = 3412

Sharp:   term   144   -11

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432051302733

Sharp:   power(3,-289/2) = 1137

Sharp:   term   145   3

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432051302736

Sharp:   power(3,-291/2) = 379

Sharp:   term   146   -1

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432051302735

Sharp:   power(3,-293/2) = 126

Sharp:   term   147   0

Sharp:   total   523598775598298873077107230546583814032861566562517636829157432051302735

Sharp:   times   6   3141592653589793238462643383279502884197169399375105820974944592307816410

> Maple_pi := evalf(Pi);

> Sharp_error := Sharp_pi-Maple_pi;

>

> infolevel[Sharp] := 0;

> for Digits from 100 by 100 to 1000 do
st := time(): Spi := Sharp(): time_Spi[Digits] := time()-st;
lprint(Digits,time_Spi[Digits]);
od:

100, .37e-1

200, .79e-1

300, .136

400, .189

500, .608

600, .592

700, .783

800, .907

900, 1.103

1000, 1.296

> Digits := 10;

>

Machin's Calculation

In 1706 John Machin discovered the following formula

> Machin_formula := Pi/4=4*arctan(1/5)-arctan(1/239);

Multiply the above by 4 to find

> Machin_Pi := 4*Machin_formula;

Verification of Machin's formula using the addition theorem for arctan

> arctan_add := proc(a,b) (a+b)/(1-a*b) end;

> a1 := 1/5;

> a2 := arctan_add(a1,a1);

> eq1 := arctan(a1)+arctan(a1)=arctan(a2);

Numerical check

> evalf(eq1,12);

> a4 := arctan_add(a2,a2);

> eq2 := arctan(a2)+arctan(a2)=arctan(a4);

> evalf(eq2,12);

> eq3 := 4*arctan(a1)=arctan(a4);

> evalf(eq3,12);

> b1 := 1/239;

> a5 := arctan_add(a4,-b1);

> eq4 := 4*arctan(a1)-arctan(b1)=arctan(a5);

> evalf(eq4,12);

Derivation of Machin's formula

> Max_N := 5;

>    printlevel := 2;
for N from 2 to Max_N do
Mf1[N,1] := 1/N;
Mf2[N,1] := arctan_add(1,-Mf1[N,1]);
for j while Mf1[N,j]<1 do
  Mf1[N,j+1] := arctan_add(Mf1[N,j],Mf1[N,1]);
  Mf2[N,j+1] := arctan_add(1,-Mf1[N,j+1]);
od;
od;
printlevel:=1;

Machin computed to 100 digits to the right of the decimal point, or 101 significant digits using his formula

> interface(verboseproc=0);

> Machin := proc()
local g,t1,t2,t;
description "computes Pi using Machin's formula",
"Pi/4=4*arctan(1/5)-arctan(1/239)";
g := length(Digits+9)+2;
userinfo(2,Machin,cat("using ",g," extra guard digits"));
Digits := Digits+g;
t1 := arctan1over(5);
userinfo(1,Machin,arctan(1/5)=t1);
t2 := arctan1over(239);
userinfo(1,Machin,arctan(1/239)=t2);
t := 16*t1-4*t2;
userinfo(1,Machin,Pi=16*arctan(1/5)-4*arctan(1/239));
Digits := Digits-g;
evalf(t);
end;

>

>    arctan1over := proc(n)
local one,n2,pn,t,a,s,k,i;
description "computes arctan(1/n) using arctan's Maclaurin series";
one := 10^Digits;
userinfo(1,procname,'computing',arctan(1/n));
userinfo(2,procname,'using','integer','arithmetic');
n2 := n^2;
pn := iquo(one,n);
userinfo(3,procname,'term',1,pn);
t := pn;
a := t;
s := 1;
k := 1;
for i while t<>0 do
  pn := iquo(pn,n2);
  k := k+2;
  userinfo(4,procname,power(n,-k)=pn);
  t := iquo(pn,k);
  s := -s;
  t := s*t;
  userinfo(3,procname,'term',i+1,t);
  a := a+t;
  userinfo(4,procname,total=a);
od;
userinfo(2,procname,'used',i-1,'terms');
Float(a,-Digits);
end;

>

> Digits := 101;

> infolevel[Machin] := 2;

> infolevel[arctan1over] := 4;

> Machin_pi := Machin();

> Maple_pi := evalf(Pi);

> Machin_error := Machin_pi-Maple_pi;

> infolevel[Machin] := 0;

> infolevel[arctan1over] := 0;

In 1873 William Shanks computed Pi to 707 places (527 correct) using Machin's formula

> Digits := 708;

> Spi := Machin();

> for Digits from 1000 by 1000 to 10000 do
st := time(); Mpi := Machin(); time_Mpi[Digits] := time()-st;
lprint(Digits,time_Mpi[Digits]);
od:

1000, .332

2000, 1.536

3000, 3.005

4000, 5.015

5000, 7.420

6000, 10.429

7000, 13.888

8000, 17.644

9000, 22.558

10000, 27.661

> Digits := 10;

>

Euler's Formulas

In 1779 Leonard Euler published the following formulas for Pi

> Euler_formula1 := Pi=4*arctan(1/2)+4*arctan(1/3);

> Euler_formula2 := Pi=20*arctan(1/7)+8*arctan(3/79);

They were also discovered by Hutton in 1776.

Verification of Euler's formulas using the addition theorem for arctan

> combine(Euler_formula1,arctan);

> combine(Euler_formula2,arctan);

Euler also discovered a new rational function series for arctan(x), which may be derived by repeated integration by parts applied to the integral int(1/(1+x^2),x) = x/(1+x^2) - int(-2*x^2/(1+x^2)^2,x)

The following two procedures use Euler's formulas to compute using his arctan series.

> Euler1 := proc()
local p;
Digits := Digits+2;
p := 2.4*`evalf/hypergeom/kernel`([1,1],[3/2],1/10)
+.56*`evalf/hypergeom/kernel`([1,1],[3/2],1/50);
Digits := Digits-2;
evalf(p);
end;

>

> Euler2 := proc()
local p;
Digits := Digits+2;
p := 2.8*`evalf/hypergeom/kernel`([1,1],[3/2],1/50)
+.30336*`evalf/hypergeom/kernel`([1,1],[3/2],9/6250);
Digits := Digits-2;
evalf(p);
end;

> Digits := 1000;

> st := time(): Epi1 := Euler1(): time_Epi1 := time()-st;

> st := time(): Epi2 := Euler2(): time_Epi2 := time()-st;

> Maple_pi := evalf(Pi):

> Echeck1 := Epi1-Maple_pi;

> Echeck2 := Epi2-Maple_pi;

> Digits := 10;

>

Rutherford's Calculation

In 1841 William Rutherford using the following formula to computed Pi to 208 places (152 correct)

> Rutherford_formula := Pi/4=4*arctan(1/5)-arctan(1/70)+arctan(1/99);

Verification of Rutherford's formula using the addition theorem for arctan

> combine(Rutherford_formula,arctan);

>

> Rutherford := proc()
local g,t1,t2,t3,t;
g := length(Digits+9)+2;
userinfo(2,procname,cat("using ",g," extra guard digits"));
Digits := Digits+g;
t1 := arctan1over(5);
userinfo(1,procname,arctan(1/5)=t1);
t2 := arctan1over(70);
userinfo(1,procname,arctan(1/70)=t2);
t3 := arctan1over(99);
userinfo(1,procname,arctan(1/99)=t3);
t := 16*t1-4*t2+4*t3;
userinfo(1,procname,Pi=16*arctan(1/5)-4*arctan(1/70)+4*arctan(1/99));
Digits := Digits-g;
evalf(t);
end:

> Digits := 208;

> Rpi := Rutherford();

> Maple_pi := evalf(Pi);

> R_error := Rpi-Maple_pi;

>

Dase's Calculation

In 1844 Zacharias Dase used the following formula to compute Pi to 205 places (200 correct) in his head.

> Dase_formula := Pi/4=arctan(1/2)+arctan(1/5)+arctan(1/8);

Verification of Dase's formula using the addition theorem for arctan

> combine(Dase_formula,arctan);

>

> Dase := proc()
local g,t1,t2,t3,t;
g := length(Digits+9)+2;
userinfo(2,procname,cat("using ",g," extra guard digits"));
Digits := Digits+g;
t1 := arctan1over(2);
userinfo(1,procname,arctan(1/2)=t1);
t2 := arctan1over(5);
userinfo(1,procname,arctan(1/5)=t2);
t3 := arctan1over(8);
userinfo(1,procname,arctan(1/8)=t3);
t := t1+t2+t3;
t := 4*t;
userinfo(1,procname,Pi=4*arctan(1/2)+4*arctan(1/5)+4*arctan(1/8));
Digits := Digits-g;
evalf(t);
end:

> Digits := 205;

> Dpi := Dase();

> Maple_pi := evalf(Pi);

> D_error := Dpi-Maple_pi;

> Digits := 10;

>

Clausen's Calculation

In 1847 Thomas Clausen used the following formula to compute Pi to 250 places (248 correct).

> Clausen_formula := Pi/4=2*arctan(1/3)+arctan(1/7);

Verification of Clausen's formula using the addition theorem for arctan

> combine(Clausen_formula,arctan);

>

> Clausen := proc()
local g,t1,t2,t;
g := length(Digits+9)+2;
userinfo(2,procname,cat("using ",g," extra guard digits"));
Digits := Digits+g;
t1 := arctan1over(3);
userinfo(1,procname,arctan(1/3)=t1);
t2 := arctan1over(7);
userinfo(1,procname,arctan(1/7)=t2);
t := 8*t1+4*t2;
userinfo(1,procname,Pi=8*arctan(1/3)+4*arctan(1/7));
Digits := Digits-g;
evalf(t);
end:

> Digits := 250;

> Cpi := Clausen();

> Maple_pi := evalf(Pi);

> C_error := Cpi-Maple_pi;

> Digits := 10;

>

Ramanujan's Formulas

Before 1927 Srinivasa Ramaujan discovered several infinite series for 1/Pi .

Here are a few of them.

>

> SR1_series := 1/Pi=
Sum(binomial(2*n,n)^3*(42*n+5)/2^(12*n+4),
n=0..infinity);

> SR2_series := 1/Pi=
1/4*
Sum((-1)^n*(1123+21460*n)*(4*n)!
/(882^(2*n+1)*n!^4*2^(8*n)),
n=0..infinity);

> SR3_series := 1/Pi=
8^(1/2)/9801*
Sum((4*n)!*(1103+26390*n)
/(396^(4*n)*n!^4),
n=0..infinity);

>

>    SR1 := proc()
local t;
Digits := Digits+4;
t := 5*`evalf/hypergeom/kernel`(
  [1/2, 1/2, 1/2, 47/42],
  [1, 1, 5/42],1/64);
t := 16/t;
Digits := Digits-4;
evalf(t);
end;

>    SR2 := proc()
local t;
Digits := Digits+4;
t := 1123*`evalf/hypergeom/kernel`(
  [3/4, 1/2, 1/4, 22583/21460],
  [1123/21460, 1, 1],-1/777924);
t := 3528/t;
Digits := Digits-4;
evalf(t);
end;

>    SR3 := proc()
local t;
Digits := Digits+4;
t := 1103*sqrt(8.)*
`evalf/hypergeom/kernel`(
  [3/4, 1/2, 1/4, 27493/26390],
  [1, 1, 1103/26390],1/96059601);
t := 9801/t;
Digits := Digits-4;
evalf(t);
end;

> Digits := 10000;

> st := time(): Rp1 := SR1(): time_Rp1 := time()-st;

> st := time(): Rp2 := SR2(): time_Rp2 := time()-st;

> st := time(): Rp3 := SR3(): time_Rp3 := time()-st;

> Rp1-evalf(Pi);

> Rp2-evalf(Pi);

> Rp3-evalf(Pi);

> Digits := 10;

>

Shanks's and Wrench's Calculation

In 1961 Daniel Shanks and John W. Wrench, Jr. computed Pi to 100,000 decimals
The calculation took about 8 hours.

They used an IMB7090 computer and the following formulas for Pi.

In 1896 Stormer found the following arctan formula for Pi.

> Stormer_formula := Pi=24*arctan(1/8)+8*arctan(1/57)+4*arctan(1/239);

Before 1863 Gauss found the following arctan formula for Pi

> Gauss_formula := Pi=48*arctan(1/18)+32*arctan(1/57)-20*arctan(1/239);

Verification of the formulas using the addition theorem for arctan

> combine(Stormer_formula,arctan);

> combine(Gauss_formula,arctan);

> Wrench1 := proc()
local g,t1,t2,t3,t;
g := length(Digits+9)+2;
userinfo(2,procname,cat("using ",g," extra guard digits"));
Digits := Digits+g;
t1 := arctan1over(8);
userinfo(1,procname,arctan(1/8)=t1);
t2 := arctan1over(57);
userinfo(1,procname,arctan(1/57)=t2);
t3 := arctan1over(239);
userinfo(1,procname,arctan(1/239)=t2);
t := 24*t1+8*t2+4*t3;
userinfo(1,procname,Pi=24*arctan(1/8)+8*arctan(1/57)+4*arctan(1/239));
Digits := Digits-g;
evalf(t);
end:

> Wrench2 := proc()
local g,t1,t2,t3,t;
g := length(Digits+9)+2;
userinfo(2,procname,cat("using ",g," extra guard digits"));
Digits := Digits+g;
t1 := arctan1over(18);
userinfo(1,procname,arctan(1/18)=t1);
t2 := arctan1over(57);
userinfo(1,procname,arctan(1/57)=t2);
t3 := arctan1over(239);
userinfo(1,procname,arctan(1/239)=t2);
t := 48*t1+32*t2-20*t3;
userinfo(1,procname,Pi=24*arctan(1/8)+8*arctan(1/57)+4*arctan(1/239));
Digits := Digits-g;
evalf(t);
end:

> W1 := Wrench1();

> W2 := Wrench2();

> Digits := 1000;

> st := time(): W1 := Wrench1(): time_W1 := time()-st;

> st := time(): W2 := Wrench2(): time_W2 := time()-st;

> W1-W2;

> Digits := 10;

>

Bailey's Calculation using Borwein's quartic formula

Before 1987, Jonathan Borwein and Peter Borwein had discovered a quartically convergent algorithm for computing .

In 1986 David Bailey used the quartic algorithm to compute to 29,360,000 decimal digits. This took 28 hours CPU time on a Cray 2 computer.

The algorithm is as follows

> y[0] = sqrt(2)-1;

> a[0]=6-4*sqrt(2);

> y[k+1]=(1-(1-y[k]^4)^(1/4))/(1+(1-y[k]^4)^(1/4));

> a[k+1]=a[k]*(1+y[k+1])^4 - 2^(2*k+3)*y[k+1]*(1+y[k+1]+y[k+1]^2);

Then a[k] converges quartically to 1/ .

Here is a procedure called BB4 which performs the Borwein quartic algorithm

>    BB4 := proc()
local g,t,a0,y0,a,y,err,eps,p2,k,oa0;
g := length(Digits+9)+2;
userinfo(2,BB4,'using',g,'guard','digits');
Digits := Digits+g;
t := evalf(2^(1/2));
a0 := 6-4*t;
y0 := t-1;
userinfo(3,BB4,y[0]=y0);
userinfo(3,BB4,a[0]=a0);
err := 1.;
eps := Float(1,-iquo(Digits+4,4));
p2 := 2;
for k while eps<err do
  oa0 := a0;
  t := 1-y0^4;
  t := t^(1/2);
  t := t^(1/2);
  y0 := (1-t)/(1+t);
  userinfo(3,BB4,y[k]=y0);
  p2 := 4*p2;
  t := 1+y0;
  a0 := a0*t^4-p2*y0*(t+y0^2);
  userinfo(3,BB4,a[k]=a0);
  err := abs(a0-oa0);
  userinfo(4,BB4,Error=err);
od;
userinfo(1,BB4,
  cat(`used `,k-1,` iterations of the Borwein's quartic algorithm`));
t := evalf(1/a0);
Digits := Digits-g;
evalf(t);
end:

> Digits := 50;

> infolevel[BB4] := 4;

> BBpi := BB4();

BB4:   using   4   guard   digits

BB4:   y[0] = .41421356237309504880168872420969807856967187537694807

BB4:   a[0] = .34314575050761980479324510316120768572131249849220772

BB4:   y[1] = .373488546332513358613022412855034494969481731845898753e-2

BB4:   a[1] = .318309886931161151029133158227185981842532714542415883

BB4:   Error = .24835863576458653764111944934021703878779783949791837e-1

BB4:   y[2] = .243231134188186167948111009505177035432623395686094427e-10

BB4:   a[2] = .318309886183790671537767526745028724068924835886669518

BB4:   Error = .747370479491365631482157257773607878655746365e-9

BB4:   y[3] = .437508679045000000000000000000000000000000019141384424e-43

BB4:   a[3] = .318309886183790671537767526745028724068919291480912869

BB4:   Error = .5544405756649e-41

BB4:   used 3 iterations of the Borwein's quartic algorithm

> Maple_pi := evalf(Pi);

> BB_error := BBpi-Maple_pi;

> infolevel[BB4] := 0;

> for Digits from 1000 by 1000 to 10000 do
st := time(): BBpi := BB4(): time_BBpi[Digits] := time()-st;
lprint(Digits,time_BBpi[Digits]);
od:

1000, .152

2000, .402

3000, .939

4000, 1.368

5000, 1.912

6000, 2.480

7000, 3.120

8000, 3.977

9000, 4.839

10000, 5.460

> Digits := 10;

>

Bailey, Borwein and Plouffe binary digit extraction formula

In 1996, David Bailey, Peter Borwein and Simon Plouffe had discovered an algorithm that could compute a base 16 digit of without computing the previous digits.

The algorithm is based on the following formula for Pi.

> BBP_formula := Pi=
Sum(1/16^n*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)),
n=0..infinity);

Numerically check the formula

> Digits := 40;

> check1 := evalf(BBP_formula);

> evalb(check1);

> Digits := 10;

BBP(n) will find a few base 16 digits of Pi starting with the n'th to the right of the decimal place

>    BBP := proc(N::posint)
local N1,s,pk,k,nu,os,b,bp;
Digits := 16;
N1 := N-1;
s := 0.;
pk := -7;
for k from 0 to N1 do
  pk := pk+8;
  nu := modp(4*power(16,N1-k),pk);
  s := s+nu/Float(pk,0);
end do;
s := s-trunc(s);
os := s+1;
b := 1.;
while os<>s do
  os := s;
  pk := pk+8.;
  b := 16.*b;
  bp := b*pk;
  s := s+4./bp;
end do;
pk := -4;
for k from 0 to N1 do
  pk := pk+8;
  nu := modp(2*power(16,N1-k),pk);
  s := s-nu/Float(pk,0);
end do;
os := s+1;
b := 1.;
while os<>s do
  os := s;
  pk := pk+8.;
  b := 16.*b;
  bp := b*pk;
  s := s-2./bp;
end do;
pk := -3;
for k from 0 to N1 do
  pk := pk+8;
  nu := modp(power(16,N1-k),pk);
  s := s-nu/Float(pk,0);
end do;
os := s+1;
b := 1.;
while os<>s do
  os := s;
  pk := pk+8.;
  b := 16.*b;
  bp := b*pk;
  s := s-1/bp;
end do;
pk := -2;
for k from 0 to N1 do
  pk := pk+8;
  nu := modp(power(16,N1-k),pk);
  s := s-nu/Float(pk,0);
end do;
os := s+1;
b := 1.;
while os<>s do
  os := s;
  pk := pk+8.;
  b := 16.*b;
  bp := b*pk;
  s := s-1/bp;
end do;
s := s-trunc(s);
if s<0 then s := s+1; fi;
cfb(s,16);
end:

>    cfb := proc(aa::numeric,b::posint)
local a,sa,La,d,n,ba,i,q,s;
if b=1 then ERROR(`the base must be > 1`) end if;
a := aa;
sa := sign(a);
if type(a,integer) then a := Float(sa*a,0);
elif type(a,fraction) then a := evalf(sa*a);
else a := Float(sa*op(1,a),op(2,a));
end if;
if a=0 then RETURN([1, [0], b, 0]);
end if;
La := length(op(1,a));
d := ilog10(a);
n := trunc(evalhf(d*ln(10)/ln(b)+1));
Digits := La+9;
a := evalf(a/(Float(b,0)^n));
while 1.<=a do
  a := a/Float(b,0);
  n := 1+n;
end do;
ba := b*a;
while ba<1. do
  a := ba;
  ba := b*a;
  n := n-1;
end do;
d := ceil(evalhf(La*ln(10)/ln(b)+1));
for i to d do
  a := b*a;
  q := trunc(a);
  a := a-q;
  s[i] := q;
end do;
[sa,[seq(s[i],i=1..d)],b,n];
end:

Try it out

> b19 := BBP(19);

> Digits := 50;

> pi := evalf(Pi);

> pi16 := cfb(pi,16);

> pi2 := pi16[2];

> pi2[19+1..34];

>

Kanada's Calculation

In 2002 Kanada used the following formulas to compute to a trillion digits (10^12 digits).

> Kpi1 := 48*arctan(1/49)+128*arctan(1/57)
-20*arctan(1/239)+48*arctan(1/110443);

> Kpi2 := 176*arctan(1/57)+28*arctan(1/239)
-48*arctan(1/682)+96*arctan(1/12943);

Verification of Rutherford's formula using the addition theorem for arctan

> combine(Kpi1,arctan);

> combine(Kpi2,arctan);

>

> Kanada1 := proc()
local g,t1,t2,t3,t4,t;
g := length(Digits+9)+2;
userinfo(2,procname,cat("using ",g," extra guard digits"));
Digits := Digits+g;
t1 := arctan1over(49);
userinfo(1,procname,arctan(1/49)=t1);
t2 := arctan1over(57);
userinfo(1,procname,arctan(1/57)=t2);
t3 := arctan1over(239);
userinfo(1,procname,arctan(1/239)=t3);
t4 := arctan1over(110443);
userinfo(1,procname,arctan(1/110443)=t4);
t := 48*t1+128*t2-20*t3+48*t4;
userinfo(1,procname,Pi=48*arctan(1/49)+128*arctan(1/57)
-20*arctan(1/239)+48*arctan(1/110443));
Digits := Digits-g;
evalf(t);
end:

> Kanada2 := proc()
local g,t1,t2,t3,t4,t;
g := length(Digits+9)+2;
userinfo(2,procname,cat("using ",g," extra guard digits"));
Digits := Digits+g;
t1 := arctan1over(57);
userinfo(1,procname,arctan(1/57)=t1);
t2 := arctan1over(239);
userinfo(1,procname,arctan(1/239)=t2);
t3 := arctan1over(682);
userinfo(1,procname,arctan(1/682)=t3);
t4 := arctan1over(12943);
userinfo(1,procname,arctan(1/12943)=t4);
t := 176*t1+28*t2-48*t3+96*t4;
userinfo(1,procname,Pi=176*arctan(1/57)+28*arctan(1/239)
-48*arctan(1/682)+96*arctan(1/12943));
Digits := Digits-g;
evalf(t);
end:

> Digits := 100;

> K_pi1 := Kanada1();

> Maple_pi := evalf(Pi);

> K_p11_error := K_pi1-Maple_pi;

> K_pi2 := Kanada2();

> K_pi2_error := K_pi2-Maple_pi;

> Digits := 10;

>