This article illustrates the convergence - or lack of it - of a numeric sequence in the world of digital computers, and the need to consider such not unfrequent cases when designing algorithms in scientific computing.
The subject is a very old and known one: to find the limit of the ratio between the perimeter of a regular polygon with \(n\) sides inscript in a circle, and the diameter of the circle, when \(n\) tends to infinity, but the concern of this article is, on the road, to get acquainted with some particularities of numerical computing in digital computers.
It is not dificult to find a formula for that perimeter using elementary trigonometry: beeing \(r\) the circle’s radius, the side of the polygon is \(2 \times r \times sin(\frac{\theta}{2})\), where \(\theta\) is the central angle subtended by the polygon’s side, equal to \(\frac{2\pi}{n}\). So the polygon’s perimeter is \(2n \times r \times sin(\frac{\theta}{2})\).
The ratio searched, perimeter / diameter is, then, \(n \times sin(\frac{\pi}{n})\)
Moving from integers to real numbers, let \(f(x) = x.sin(\frac{\pi}{x})\). We are interested in \(L = \lim_{x\to\infty} f(x)\), and we’ll use some elementary calculus to find it:
\[
L = \lim_{x\to\infty} x.sin(\frac{\pi}{x}) = \lim_{x\to\infty}\frac{sin(\frac{\pi}{x})}{\frac{1}{x}} = \lim_{x\to\infty}\frac{cos(\frac{\pi}{x})(-\frac{\pi}{x^2})}{-\frac{1}{x^2}} = \lim_{x\to\infty}\pi. cos(\frac{\pi}{x}) = \pi\times1 = \pi
\]
Let’s do some calculations, using the R open source programming language, to watch the results the formula \(n.sin(\frac{\pi}{n})\) throws when we take several values of \(n\):
for (n in seq(3,1003,50)){
cat(sprintf('%10d %23.14f\n',n,n*sin(pi/n)))
}Trigonometric formula:
==========================================
Number of sides Perimeter/diameter ratio
========================================== 3 2.59807621135332
53 3.13975327836686
103 3.14110556975456
153 3.14137190072373
203 3.14146725256340
253 3.14151192005327
303 3.14153636624320
353 3.14155118232917
403 3.14156083455835
453 3.14156747096938
503 3.14157222861400
553 3.14157575511904
603 3.14157844132946
653 3.14158053445816
703 3.14158219706711
753 3.14158353961050
803 3.14158463926489
853 3.14158555127347
903 3.14158631602220
953 3.14158696358895
1003 3.14158751674684From this answer table is strikingly evident that the ratio tends indeed to the irrational number \(\pi\), as we expected. Even from the geometrical point of view, this is not surprising, since a regular polygon with many sides tends to confound itself with the circumscribing circle, so its perimeter approximates the circumference (the circle’s perimeter), and, hence, its relation to the diameter gets close to \(\pi\).
But the rate of convergency to \(\pi\) is rather low, and that’s why we have computed only one of each fifty consecutive numbers, as shown in the former table.
Many centuries ago, Archimedes (287 - 212 BCE), one of the three greatest matematicians ever born, developed an algorithm for faster search of more digits of \(\pi\), duplicating the number of sides on each step, progressing, thus, much faster. We can now easily use our trigonometry to do something similar:
for (n in seq(2,19)){
m = 2**n
cat(sprintf('%10d %23.14e\n',m,m*sin(pi/m)))
}Archimedes's fashion:
==========================================
Number of sides Perimeter/diameter ratio
========================================== 4 2.82842712474619
8 3.06146745892072
16 3.12144515225805
32 3.13654849054594
64 3.14033115695475
128 3.14127725093277
256 3.14151380114430
512 3.14157294036709
1024 3.14158772527716
2048 3.14159142151120
4096 3.14159234557012
8192 3.14159257658487
16384 3.14159263433856
32768 3.14159264877699
65536 3.14159265238659
131072 3.14159265328899
262144 3.14159265351459
524288 3.14159265357099But there was not yet trigonometry on those old times, so Archimedes went on in a purely pythagorical fashion, beginning from the inscript hexagon:
Using the Pythagorical theorem, he reached a recursive formula relating the side length of a polygon of \(2n\) sides to the side length of a polygon of \(n\) sides, and, consequently, the respective perimeters:
\[
l_{2n} = \sqrt{2r^2 - r \sqrt{4r^2 - {l_n}^2}}
\] If we take a unitary radius, without loosing generality we can simplify this formula. Furthermore, in that case the diameter’s length is 2, and the searched ratio is just the semiperimeter, \(SP_{2n}\): \[
l_{2n} = \sqrt{2 - \sqrt{4 - {l_n}^2}}
\] \[
SP_{2n} = n \times \sqrt{2 - \sqrt{4 - {l_n}^2}}
\] Let’s apply this to our search:
m = 3; lp = 1 # Initialization: 3 sides, side length of the hexagon
for (n in seq(6,34)){
m = 2*m # Number of sides
p = lp * m # Perimeter
cat(sprintf('%10d %14.8e %24.16e\n',m,lp,p/2))
lp = sqrt(2 - sqrt(4 - lp**2))
} Archimedes without trigonometry (only Pythagoras):
=======================================================
Number of sides Side length Semiperimeter
======================================================= 6 1.00000000e+00 3.0000000000000000
12 5.17638090e-01 3.1058285412302498
24 2.61052384e-01 3.1326286132812369
48 1.30806258e-01 3.1393502030468721
96 6.54381656e-02 3.1410319508905298
192 3.27234633e-02 3.1414524722853443
384 1.63622792e-02 3.1415576079116221
768 8.18120805e-03 3.1415838921489359
1536 4.09061258e-03 3.1415904632367617
3072 2.04530736e-03 3.1415921060430483
6144 1.02265381e-03 3.1415925165881546
12288 5.11326924e-04 3.1415926186407894
24576 2.55663464e-04 3.1415926453212157
49152 1.27831732e-04 3.1415926453212157
98304 6.39158660e-05 3.1415926453212157
196608 3.19579330e-05 3.1415926453212157
393216 1.59789648e-05 3.1415923038117377
786432 7.98948238e-06 3.1415923038117377
1572864 3.99473424e-06 3.1415868396550413
3145728 1.99736712e-06 3.1415868396550413
6291456 9.98711352e-07 3.1416742650217575
12582912 4.99355676e-07 3.1416742650217575
25165824 2.49344118e-07 3.1374750995027831
50331648 1.24672059e-07 3.1374750995027831
100663296 6.32202728e-08 3.1819805153394638
201326592 2.98023224e-08 3.0000000000000000
402653184 1.49011612e-08 3.0000000000000000
805306368 0.00000000e+00 0.0000000000000000
1610612736 0.00000000e+00 0.0000000000000000He arrived - manually and using attic numerals (a predecessor notation of roman numerals)! - as far as 96 sides, reaching 4 correct digits, a better approximation to the real \(\pi\) that the known fraction \(\frac{22}{7}\), used in western mathematics during many centuries after Archimedes. As a matter of fact, he concluded correctly that \(\frac{223}{71} \lt \pi \lt \frac{22}{7}\).
If we presently compute further approximations, we can see that the ratio perimeter/diameter (or, simply, the semiperimeter of a circle with unit radius) converges towards \(\pi\) until 24576 sides.
There apparently stops until 196608 sides (within the 16 decimal digits used in the table). Afterwards, a clear divergence begins, arriving to the value 3 at the 201326592 sides case. Two more steps, and the side length, and thus the semiperimeter, drop suddenly to zero, where it finally remains!
The reason of this weird behaviour is numerical instability of the algorithm in the floating point number system used by digital computers. Concretely, due to what is known as “catastrophic cancellation”. But, beeing aware of this condition, is not difficult to avoid it:
To understand what’s happening, look at the formula \(l_{2n} = \sqrt{2 - \sqrt{4 - {l_n}^2}}\). As the number of sides \(n\) grows, the length of each side shrinks. When it has approached zero enough, the result of the difference \(4 - {l_n}^2\) is almost 4. To notice that it’s a little less than 4, it’d be required a number of decimal digits greater than those kept within the computer floating point numeric representation (presently, around 15 decimal digits according to IEEE 754 standard, in what is called double precision, that keeps 8 bytes for each number). Even if a difference is noticeable, the number of correct significant digits in that difference gets shorter when \(n\) grows further than some value (in our example, 24576). It’s said that the less significant digits get cancelled, and this is the harmful cancelative substraction.
At certain point in the algorithm, (when \(n = 805306368\)), all the difference digits have been cancelled, and we get \(l_{2n} = \sqrt{2 - \sqrt{4}} = 0\).
So, the ultimate issue’s cause is that the computing number system is not the good and trustable mathematical real number field, which is infinite and continuous, but a finite and discrete subset of reals, called floating point system, which all digital computers, not having infinite speed nor storage capacity, use to represent real numbers, and fails to strictly obey most properties of the real number arithmetic. The specially dangerous arithmetic operation is substraction, and you should always be aware of not substracting too similar numbers.
To avoid in our formula that harmful operation it’s enough to multiply and divide it by its “conjugate”, employing a word reminiscent of the complex number’s conjugate:
\[ l_{2n}^2 = 2 - \sqrt{4 - {l_n}^2}\times\frac{2 + \sqrt{4 - {l_n}^2}}{2 + \sqrt{4 - {l_n}^2}} = \frac{4-(4-{l_n}^2)}{2 + \sqrt{4 - {l_n}^2}} = \frac{{l_n}^2}{2 + \sqrt{4 - {l_n}^2}} \] And extracting the square root, we reach the stabilized formula for \(l_{2n}\): \[ l_{2n} = \frac{l_n}{\sqrt{2 + \sqrt{4 - {l_n}^2}}} \] Where the harmful substraction has been changed into an innocuous sum.
Now see the effect of this transformation on the results:
m = 3; lp = 1 # initialize to 3 sides and hexagon's side length
for (n in seq(6,32)){
m = 2*m # Number of sides
p = lp * m # Perimeter
cat(sprintf('%10d %14.8e %24.16f\n',m,lp,p/2))
lp = lp/sqrt(2 + sqrt(4 - lp**2))
}Pythagorical Archimedes stabilized:
=======================================================
Number of sides Side length Semiperimeter
======================================================= 6 1.00000000e+00 3.0000000000000000
12 5.17638090e-01 3.1058285412302489
24 2.61052384e-01 3.1326286132812378
48 1.30806258e-01 3.1393502030468667
96 6.54381656e-02 3.1410319508905093
192 3.27234633e-02 3.1414524722854615
384 1.63622792e-02 3.1415576079118575
768 8.18120805e-03 3.1415838921483177
1536 4.09061258e-03 3.1415904632280496
3072 2.04530736e-03 3.1415921059992709
6144 1.02265381e-03 3.1415925166921563
12288 5.11326924e-04 3.1415926193653831
24576 2.55663464e-04 3.1415926450336897
49152 1.27831732e-04 3.1415926514507664
98304 6.39158662e-05 3.1415926530550360
196608 3.19579331e-05 3.1415926534561032
393216 1.59789665e-05 3.1415926535563701
786432 7.98948327e-06 3.1415926535814371
1572864 3.99474164e-06 3.1415926535877032
3145728 1.99737082e-06 3.1415926535892700
6291456 9.98685409e-07 3.1415926535896617
12582912 4.99342704e-07 3.1415926535897594
25165824 2.49671352e-07 3.1415926535897842
50331648 1.24835676e-07 3.1415926535897905
100663296 6.24178380e-08 3.1415926535897913
201326592 3.12089190e-08 3.1415926535897922
402653184 1.56044595e-08 3.1415926535897922We can appreciate now how the theorical convergence to \(\pi\) has been installed into the squence, at least as far as the 15 decimal digit’s double precision allows. (Decimal digit in the sense of beeing one of {0,1,2,3,4,5,6,7,8,9}, not meaning fractional digits).
Working with quadruple precision (16 bytes per number), the same behaviour is observed, but with far more sides in the polygon and significant digits in the final answer.
Just for the sake of it, let’s run the stabilized algorithm with quadruple precision. To avoid an undesirable flood of numbers, we’ll only print the final results.
# initialize to 3 sides and hexagon's side length
(m <- mpfr(3, 128))
(lp <- mpfr(1, 128))
for (n in seq(6,68)){
m = 2*m # Number of sides
p = lp * m # Perimeter
lp = lp/sqrt(2 + sqrt(4 - lp**2))
}Pythagorical Archimedes stabilized in quadruple precision:
==========================================================Number of sides: 27670116110564327424.Side length: 1.135373860028851529778e-19Semiperimeter: 3.141592653589793238462643383279502884254Embedded pi: 3.141592653589793238462643383279502884195The comparison between the found semiperimeter and the embedded \(\pi\) the R language brings incorporated, allows us to ascertain to have gotten 37 correct significative digits. Thus, the sequence has continued converging numerically, through the stabilized algorithm.
To end with a historical reference, we’ll say that Archimede’s method for finding \(\pi\) was widely used by everyone until the development of infinite series expansions (in India during the 15th century and in Europe during the 17th century). To know the present methods and state of the hunting-for-\(\pi\) subject, see http://www.numberworld.org/misc_runs/pi-5t/details.html.
And if you want to know more about these kind of themes, and scientific computing in general, learn Numerical Analysis!
DORN S. y Mc.CRACKEN D. Numerical Methods with FORTRAN IV case studies. Wiley International, USA, 1972.
FORSYTHE George, MOLER Cleve and MALCOLM Michael. Computer Methods for mathematical computations. Prentice-Hall, Englewood Cliffs, N.J., 1977.
OVERTON, Michael: Numerical Computing with IEEE Floating Point Arithmetic, SIAM, Philadelphia, 2001.
ANDREE, Richard, et al.: Computer programming: Techniques, Analysis and Mathematics. Prentice-Hall, Englewood Cliffs, N.J., 1973.
MÄCHLER, Martin: Arbitrarily Accurate Computation with R: The Rmpfr Package.
https://cran.r-project.org/web/packages/Rmpfr/vignettes/Rmpfr-pkg.pdf