For n = 4 these two polygons look like:
The outer approximation would have a circumference U
and the inner approximation
Principally that does not look too interesting, but that’s roughly the limits for our Pi computation and if we increase the number of edges it becomes smoother and the circumference of both approach each other.
Now the idea is to increase the number of edges in a loop and to calculate the approximated circumference of the inner and outer polygon in each iteration from the circumferences found in the last iteration. The number of edges is doubled each time. Therefore we first have a look at a polygon inside the circle:
A segment of the polygon inside the circle looks like:
If we double the number n of edges from one step to the next we get a new circumference u2n. This we have to calculate out of the last circumference un. Therefor the length x has to be calculated by the half part of the segment un/n first.
With this the segment for the next step u2n/2n can be calculated. This segment has the divisor 2n instead of n because it is meant for the next step having doubled as much edges already..
Polygon outside the circle
The segment of the polygon outside the circle looks like:
he triangle BDC and therefore
There the approach is a bit different. The triangle MAC is similar to
or
and
and with this
With this we can calculate the inner and the outer polygon and compare them. And if they are close enough to each other, we have Pi
static void Main(string[] args)
{
double n = 4.0;
double r = 0.5;
double un = n * r * Math.Sqrt(2.0);
double u2n = 0;
double Un = r * 8.0;
double U2n = 0;
double x;
bool run = true;
while (run)
{
u2n = 2.0 * n * Math.Sqrt(2.0 * r * r - 2.0 * r * Math.Sqrt(r * r - un * un / 4.0 / n / n));
U2n = 8.0 * n * n * r * (Math.Sqrt(r * r + Un * Un / 4.0 / n / n) - r) / Un;
run = (n < 1E7);
un = u2n;
Un = U2n;
n = n * 2;
Console.WriteLine("n= " + n.ToString() + " u = " + un.ToString() + " U = " + Un.ToString());
}
Console.ReadKey();
}
Gives the output
| n | u | U | U-u |
| 8 | 3.06146745892072 | 3.31370849898476 | 0.252241040064042 |
| 16 | 3.12144515225805 | 3.18259787807453 | 0.0611527258164766 |
| 32 | 3.13654849054594 | 3.15172490742926 | 0.0151764168833175 |
| 64 | 3.14033115695474 | 3.14411838524587 | 0.00378722829112732 |
| 128 | 3.14127725093276 | 3.14222362994234 | 0.00094637900958805 |
| 256 | 3.14151380114415 | 3.1417503691697 | 0.000236568025558714 |
| 512 | 3.14157294036788 | 3.14163208070225 | 5.9140334366159E-05 |
| 1024 | 3.14158772527996 | 3.14160251024197 | 1.47849620111273E-05 |
| 2048 | 3.14159142150464 | 3.14159511771837 | 3.69621373019768E-06 |
| 4096 | 3.14159234561108 | 3.14159326963175 | 9.24020668424674E-07 |
| 8192 | 3.141592576545 | 3.14159280837968 | 2.31834680342047E-07 |
| 16384 | 3.14159263346325 | 3.14159269096297 | 5.74997209845662E-08 |
| 32768 | 3.14159265480759 | 3.14159267557045 | 2.07628629844692E-08 |
| 65536 | 3.14159264532122 | 3.14159269096297 | 4.56417539318466E-08 |
| 131072 | 3.14159260737572 | 3.14159252378847 | -8.35872469018284E-08 |
| 262144 | 3.14159291093967 | 3.14159086957912 | -2.04136055748094E-06 |
| 524288 | 3.14159412519519 | 3.14159252378847 | -1.60140671834341E-06 |
| 1048576 | 3.14159655370482 | 3.14160058362634 | 4.02992151649073E-06 |
| 2097152 | 3.14159655370482 | 3.14159252378847 | -4.02991634684824E-06 |
| 4194304 | 3.14167426502176 | 3.1414451588708 | -0.000229106150956504 |
| 8388608 | 3.1418296818892 | 3.14097079560249 | -0.000858886286714711 |
| 16777216 | 3.14245127249413 | 3.13398329388536 | -0.00846797860877446 |
| 33554432 | 3.14245127249413 | 3.11105678802532 | -0.0313944844688132 |
Even if the delta between U and u is the smallest at n= 32768, it reaches its accuracy at n= 65536 and fits up to 7 decimal places.
Both calculations have an accuracy problem:
For the inner polygon the expression
runs against 0 if n rises and this very small value is multiplied by 2n which becomes quite big. That results in an error by extinction.
To avoid that we can use the fact that
with this
and
Same is for the outer polygon:
becomes very small too if n becomes big and can by written as
and
With this the sequence
static void Main(string[] args)
{
double n = 4.0;
double r = 0.5;
double un = n * r * Math.Sqrt(2.0);
double u2n = 0;
double Un = r * 8.0;
double U2n = 0;
bool run = true;
while (run)
{
u2n = un * r * Math.Sqrt(2.0) /
Math.Sqrt(r *
r + r * Math.Sqrt(r *
r - un / n * un / n / 4));
U2n
= 2.0
* r * Un / (Math.Sqrt(r *
r + Un * Un / 4.0
/ n / n) + r);
run = (n
< 1E8);un = u2n;
Un = U2n;
n = n * 2;
Console.WriteLine("n= " + n.ToString() + " u = "
+
un.ToString() + "
U = " + Un.ToString());
}Console.ReadKey();
}
| n | u | U | U-u |
| 8 | 3.06146745892072 | 3.31370849898476 | 0.252241040064042 |
| 16 | 3.12144515225805 | 3.18259787807453 | 0.0611527258164757 |
| 32 | 3.13654849054594 | 3.15172490742926 | 0.0151764168833162 |
| 64 | 3.14033115695475 | 3.1441183852459 | 0.00378722829115086 |
| 128 | 3.14127725093277 | 3.14222362994246 | 0.000946379009683085 |
| 256 | 3.1415138011443 | 3.14175036916897 | 0.000236568024664319 |
| 512 | 3.14157294036709 | 3.14163208070318 | 5.91403360896692E-05 |
| 1024 | 3.14158772527716 | 3.14160251025681 | 1.47849796485744E-05 |
| 2048 | 3.1415914215112 | 3.14159511774959 | 3.69623838869515E-06 |
| 4096 | 3.14159234557012 | 3.14159326962931 | 9.24059189166826E-07 |
| 8192 | 3.14159257658487 | 3.14159280759965 | 2.31014770868399E-07 |
| 16384 | 3.14159263433856 | 3.14159269209226 | 5.7753690274609E-08 |
| 32768 | 3.14159264877699 | 3.14159266321541 | 1.44384211253623E-08 |
| 65536 | 3.14159265238659 | 3.1415926559962 | 3.60960417111755E-09 |
| 131072 | 3.141592653289 | 3.1415926541914 | 9.02399044377944E-10 |
| 262144 | 3.1415926535146 | 3.14159265374019 | 2.25598206782252E-10 |
| 524288 | 3.141592653571 | 3.14159265362739 | 5.63975532941186E-11 |
| 1048576 | 3.1415926535851 | 3.14159265359919 | 1.40971678774804E-11 |
| 2097152 | 3.14159265358862 | 3.14159265359214 | 3.52207152332085E-12 |
| 4194304 | 3.1415926535895 | 3.14159265359038 | 8.78852546293274E-13 |
| 8388608 | 3.14159265358972 | 3.14159265358994 | 2.17603712826531E-13 |
| 16777216 | 3.14159265358978 | 3.14159265358983 | 5.24025267623074E-14 |
| 33554432 | 3.14159265358979 | 3.1415926535898 | 1.02140518265514E-14 |
| 67108864 | 3.1415926535898 | 3.1415926535898 | -4.44089209850063E-16 |
| 134217728 | 3.1415926535898 | 3.14159265358979 | -2.66453525910038E-15 |
| 268435456 | 3.1415926535898 | 3.14159265358979 | -3.10862446895044E-15 |
Now the accuracy is reached at n=67108864 and there our Pi is
3.1415926535898
and we have an accuracy of 11 decimal places. That’s quite cool
Isaac Newton (1642 - 1727) developed a interesting algorithm to calculate Pi too. He basically calculates the square of 1/6 of a circle with the diameter = 1:
And with d = 1
The square is separated into the 2 areas F and G. For G it’s
For F things are a bit more complicate.
F is calculated as integral:
In which for y is:
Or a bit converted:
Now there is an interesting step:
Newton goes backward. Instead of building the sum function of a series, he builds from the function
The binominal series:
I built this series by a Taylor polynomial. That’s basically the same, but looks a little different.
(For a detailed description of Taylor polynomials see Taylor Polynomials)
To implement an algorithm this version is a bit better
Now the integral of the entire function must be built.
And therefore it must be split into one integral for each part of the series:
And so on.
I intentionally separated one 4 in each denominator as there will be a multiplication by 4 at the end.
So the first integral is a bit different, but for the others it’s:
And with the above formulation the multiplication by 4 can be seen:
The entire formulation becomes:
Or if we separate a part Pn build it from Pn-1:
For n > 1
and
That can be implemented in an iterative loop like:
double n = 1.0;
double F = 1.0/3.0; // initialize with the first element
double G = 3.0 * Math.Sqrt(3.0) / 4.0;
double pi = 0;
double Pn;
double In;
Pn = 1.0 / 2.0 / 4.0; // second element is different
In = Pn / 5.0;
F = F - In;
bool run = true;
while (run)
{
n++;
Pn = Pn * (2.0 * n - 3.0) / 2.0 / 4.0 / n;
In = Pn / (2.0 * n + 3.0);
F = F - In;
run = (In > 1E-15);
pi = G + 6.0 * F;
Console.WriteLine(n.ToString() + " " + pi.ToString() + " " + In.ToString());
}Pn = Pn * (2.0 * n - 3.0) / 2.0 / 4.0 / n;
In = Pn / (2.0 * n + 3.0);
F = F - In;
run = (In > 1E-15);
pi = G + 6.0 * F;
Console.WriteLine(n.ToString() + " " + pi.ToString() + " " + In.ToString());
The value of In can be used as indicator for the accuracy of the calculation and therefore its value is used as breaking condition.With this loop me little console application makes the following output:
| n | PI | In |
| 2 | 3.14234167710523 | 0.00111607142857143 |
| 3 | 3.14169063543856 | 0.000108506944444444 |
| 4 | 3.14160740568004 | 1.38716264204545E-05 |
| 5 | 3.14159508127349 | 2.05406775841346E-06 |
| 6 | 3.14159307855742 | 3.33786010742187E-07 |
| 7 | 3.14159273144802 | 5.78515669878791E-08 |
| 8 | 3.14159266836317 | 1.05141416976326E-08 |
| 9 | 3.14159265647218 | 1.98183226443472E-09 |
| 10 | 3.14159265416507 | 3.84518542610433E-10 |
| 11 | 3.14159265370679 | 7.63793650548905E-11 |
| 12 | 3.14159265361397 | 1.54703575053308E-11 |
| 13 | 3.14159265359486 | 3.18537533514935E-12 |
| 14 | 3.14159265359087 | 6.65148939511314E-13 |
| 15 | 3.14159265359002 | 1.40588298578528E-13 |
| 16 | 3.14159265358984 | 3.00319200311721E-14 |
| 17 | 3.1415926535898 | 6.47548355203134E-15 |
| 18 | 3.1415926535898 | 1.40786421243416E-15 |
| 19 | 3.14159265358979 | 3.08365637030269E-16 |
After 19 iterations the resolution of the double value is reached and Pi is computed with an accuracy of 14 decimal places (o.k. the last one is for sure not reliable
)Newton calculated Pi with the same accuracy 300 years ago. He did all by hand. An admirable effort.
Today a computer usually computes mathematical functions with an accuracy of 15 decimal places. If one want’s to compute with more than these 15 decimal places, a variable type carrying more than 15 decimal places must be introduced and al the needed mathematical operations must be implemented for this type. I implemented these operations in library ExtFloat (see To compute with higher precision ) and with this library, computed with 1100 decimal places the following sequence
ExtFloat ONE = new ExtFloat(1);
ExtFloat TWO = new ExtFloat(2);
ExtFloat THREE = new ExtFloat(3);
ExtFloat FOUR = new ExtFloat(4);
ExtFloat FIVE = new ExtFloat(5);
ExtFloat SIX = new ExtFloat(6);
ExtFloat DELTA = new ExtFloat(ExtFloat.StrToNumber("1E-1021"));
ExtFloat t = new ExtFloat(ONE);
ExtFloat n = new ExtFloat(ONE);
ExtFloat F = new ExtFloat(ONE / THREE);
ExtFloat G = new ExtFloat(THREE);
ExtFloat pi = new ExtFloat();
ExtFloat Pn = new ExtFloat(ExtFloat.StrToNumber("1.25E-1"));
ExtFloat In = new ExtFloat();
G = THREE * ExtFloat.Root(G, 2);
G = G / FOUR;
In = Pn / FIVE;
F = F - In;
bool run = true;
int i=0;
while (run)
{
n = n + ONE;
Pn = Pn * ((TWO * n) - THREE) / TWO / FOUR / n;
In = Pn / (TWO * n + THREE);
F = F - In;
run = (ExtFloat.Cmp(In.Abs(), DELTA) > 0) && (i < 10000);
pi = G + (SIX * F);
}
Returns after 1681 cycles
3.
1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 3421170679 8214808651 3282306647 0938446095 5058223172 5359408128 4811174502 8410270193 8521105559 6446229489 5493038196 4428810975 6659334461 2847564823 3786783165 2712019091 4564856692 3460348610 4543266482 1339360726 0249141273 7245870066 0631558817 4881520920 9628292540 9171536436 7892590360 0113305305 4882046652 1384146951 9415116094 3305727036 5759591953 0921861173 8193261179 3105118548 0744623799 6274956735 1885752724 8912279381 8301194912 9833673362 4406566430 8602139494 6395224737 1907021798 6094370277 0539217176 2931767523 8467481846 7669405132 0005681271 4526356082 7785771342 7577896091 7363717872 1468440901 2249534301 4654958537 1050792279 6892589235 4201995611 2129021960 8640344181 5981362977 4771309960 5187072113 4999999837 2978049951 0597317328 1609631859 5024459455 3469083026 4252230825 3344685035 2619311881 7101000313 7838752886 5875332083 8142061717 7669147303 5982534904 2875546873 1159562863 8823537875 9375195778 1857780532 1712268066 1300192787 6611195909 2164201989
Pi with an accuracy of 1000 decimal places…jappa dappa doooo
The demo project consists of one main window:
A click on the "Compute Pi" button starts the computation. It takes quite some time. Therefore I add a progress bar on the bottom to show the progress.
For evaluation Pi can be divided by 4 and the square of the sine of this value can be computed. The square of sin( Pi / 4 ) is exactly 0.5.
In the text box on top the actual value of In is recorded for each iteration. That shows the information about the change from one iteration to the next and therefore the accuracy of the calculation.
The sine function is basically part of the ExtFlot library. I moved it into the main window in order to have the possibility to set the progress bar value from inside this function to.
C# Demo Project Newtons Pi