Sections

Archimedes’ Method
Doing it in Python

This post is part of a series. Other posts are:

One of the many ways in which USA differs from the rest of the world is the way people writes dates. In Turkey and in Chile —and I think in all Europe— we write day/month/year. On USA they write month/day/year, and it is one of the typical causes of misunderstanding. The only positive side effect is that days like today are written as 3/14 and people calls it the Pi day¹.

Based on this nimble premise, every year in this day there are people around the world that does something about the number $π.$ Let’s be some of them, now that we calculated $e$, the square root of two, and even the golden ration with ten thousand digits.

There are several formulas to find the numeric value of $π.$ In general this is a hard problem. The number of digits known increased slowly until the invention of electronic computers. Since then, it has been a tradition for each new high-performance computer to prove its power by breaking the record of the numbers of digits of $π.$

Archimedes’ method

In today’s words, we would say that Archimedes was an engineer and a mathematician that lived between 287 BCE and 211 BCE in Syracuse, Sicily. At that time the city was part of the Greek cultural sphere, as many places in the north of the Mediterranean sea. Archimedes has been famous since then for many inventions and ideas. One of these ideas was to realize that a regular polygon with many sides —let’s say $n$— will be more and more like a circle when $n$ gets bigger. And it is easy to calculate the perimeter of a circle.

If we draw the polygon inside a circle so its vertices lay in the circumference, then each edge will cover 360/$n$ degrees. For example an hexagon has 6 edges, each one covering 360/6 or 60 degrees. Let’s put the name $α(n)$ to this angle 360°/$n.$ In the following figure, $α(n)$ corresponds to the edge AB.

Since we know that CA and CB are radius of the circle, we know that their length is 1. The line CD bisects the angle $α(n),$ and CD forms a right angle with AB. We can therefore use a little trigonometry and find that the edge length is $2\sin(α(n)/2),$ and thus the polygon perimeter is \[\text{Perimeter}(n)=2n\sin\left(\frac{α(n)}{2}\right).\] Archimedes idea was to duplicate the number of edges on every step. In other words, he calculated $\text{Perimeter}(2n),$ then $\text{Perimeter}(4n),$ and so on. Apparently he went as far as 96 edges, that is $6⋅2^4$ so he did four steps.

So we need to calculate $\sin(α(2n)/2)$ on every step. Archimedes probably used Pythagoras theorem. In today’s terms, he knew that \[\sin^2(α(n))+\cos^2(α(n))=1\] and from Euclides he could have knew that \[\cos(α(n))=\cos^2(α(n)/2)-\sin^2(α(n)/2)=\cos^2(α(2n))-\sin^2(α(2n)).\] For the last equation we need to remember that $α(n)=360°/n,$ thus $α(n)/2=α(2n).$ Combining the two equations we have \[\sin^2(α(2n))=\frac{1-\sqrt{1-\sin^2(α(n))}}{2}.\]

Let’s write this in more simple terms. Let’s call $S_n=\sin^2(α(n)).$ Then on each step we calculate $$S_{2n} = .

We just need an initial condition. For an hexagon, $\alpha(6)$ is 60 degrees, so we need to know $\sin(30°).$ But we know that AB las length 1, so $\sin(30°)$ must be 0.5 and $S_6=(0.5)^2=0.25.$

num_decimals = 70
one = 10**num_decimals

def babylon_sqrt(a):
  y = (a + one) // 2
  x = a * one // y
  while abs(x-y)>10:
    y = (x + y) // 2
    x = a * one // y
  return x

n = 6
S = one//4
for i in range(4):
  S = (one - babylon_sqrt(one-S)) // 2
  n *= 2
  print(n * babylon_sqrt(S), n)

31058285412302491481867860514885799401888268158391661657680384877806828 12
31326286132812381971617494694917362446497769154815740453748067600431896 24
31393502030468672071351468212084211891503505893625928559354072358246272 48
31410319508905096381113529264596601070364122161628349002972048407024896 96

This is what Archimedes would have found, had he used decimal representation of numbers. Let’s continue and see if we can get better.

for i in range(50):
  S = (one - babylon_sqrt(one-S)) // 2
  n *= 2
  print(n * babylon_sqrt(S))

31414524722854620754506093089612256452476623045496751758215827748232704
31415576079118576455164633451298595415043764795884965479512616709230336
31415838921483184086689696037211533505200449157810925621400265067056384
31415904632280500957384585059309517235542823086757970500652177683217920
31415921059992715505447766406101173531274972549662846725845323446904832
31415925166921574475928740847688319059677188923681876144715475913502720
31415926193653839551895493120653190422221927255948419354173694893821952
31415926450336908966721415089192384127226230112421013582657082245505024
31415926514507676517042536404922190204484274723854209040478846245879808
31415926530550368416911231804154742022572768058868531003145303754211328
31415926534561041392646431596150783313543414874294138275298342424608768
31415926535563709636628233165541133642749135041841220969851082707173376
31415926535814376697626683659225751788703495073716490005061207210328064
31415926535877043462876483788980471857502408672476606214721562547847168
31415926535892710154188945540564999738004063788315788586955240738127872
31415926535896626827017061710907747199859793791593346574119397561991168
31415926535897605995224090799271347533561151459628822795407744200343552
31415926535897850787275848074223367208751396830888638314880949228142592
31415926535897911985288787393140192102034265840864587927336901813469184
31415926535897927284792022222880574573760314838884579931551179435147264
31415926535897931109667830930316368707217160322740085130963712521273344
31415926535897932065886783107175360897801705020854854565141473413562368
31415926535897932304941521151390111674024112028330539973014455312187392
31415926535897932364705205662443799538615730707258650243223721145794560
31415926535897932379646126790207221515422136434931877180247352014798848
31415926535897932383381357072148077010289894182971508912137531354513408
31415926535897932384315164642633290884048468389738980001448613665308672
31415926535897932384548616535254594352490714114540669280266144008110080
31415926535897932384606979508409920219601438181560456063415522080849920
31415926535897932384621570251698751686379129363053427278369181740826624
31415926535897932384625217937520959553073552793721127010324988565651456
31415926535897932384626129858976511519747158691086061625736709428740096
31415926535897932384626357839340399511415560167850924405121309988618240
31415926535897932384626414834431371509332660537182730638328272527556608
31415926535897932384626429083204114508811935629878326115855774846550016
31415926535897932384626432645397300258681754404025273216786559192268800
31415926535897932384626433535945596696149209099726520231277910305013760
31415926535897932384626433758582670805516072711312115552981291714478080
31415926535897932384626433814241939332857788309436586420331329378844672
31415926535897932384626433828156756464693216405478076505306277859033088
31415926535897932384626433831635460747652073983619226692110789306744832
31415926535897932384626433832505136818391769316055762511195638440394752
31415926535897932384626433832722555836076627540080820753162610713559040
31415926535897932384626433832776910590497962674944305542446383829614592
31415926535897932384626433832790499279102275084810781851510210469822464
31415926535897932384626433832793896451248473289996958689830322640519168
31415926535897932384626433832794745744285930729159631688503107178725376
31415926535897932384626433832794958067474479835392254424402143727321088
31415926535897932384626433832795011148090039538724652374268096771260416
31415926535897932384626433832795024418302034287989994302872405776793600

It is not very fast, we get less than one new decimal on each step. And we may have rounding errors hiding somewhere. So let’s stop here, and let’s be happy for now with only 30 decimals.

π = n * babylon_sqrt(S)
print(str(π)[0:31])

3141592653589793238462643383279