One of the most famous series to approximate Pi was created by Gottfried Leibniz [see Wikipedia]. It states that
PI/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - ...
That series is nice looking but has a serious drawback: it converges extremely slowly: calculating PI to 10 correct decimal places requires summation of about 5 billion terms.
However, that slow confergence is a great example to show the difference between regular float mathematics and fp64lib. Implementation is straight forward:
/* Example program for using fp64lib.
* Leibniz formular to approximate PI
*/
#include <fp64lib.h>
#define MAX_COUNT 1000 // # of iterations before printing another approximation
float64_t x, n, two, pi, pi4, prev;
float xf, nf, pif, pi4f, prevf;
bool add = false;
void setup() {
prev = x = fp64_sd(1.0);
n = fp64_sd(3.0);
two = fp64_sd(2);
prevf = xf = 1.0;
nf = 3.0;
pi = float64_NUMBER_PI;
pi4 = fp64_ldexp( pi, -2 ); // fast division by 4
pif = M_PI;
pi4f = ldexp( pif, -2 );
Serial.begin(57600);
}
char *fixedLength( char *s, float64_t x ) {
// return a string always with the same length, including trailing zeros
char *t = fp64_to_string( x, 17, 15 );
sprintf( s, "%s%s00000000000000000", t[0] == '-' ? " " : " ", t );
s[18] = '\0';
return s;
}
void printLine( float nf, float64_t x, float xf, const char* s, bool usePi ) {
char buf[40];
float64_t delta = fp64_sub( x, usePi ? pi : pi4 );
float deltaf = xf - (usePi ? pif : pi4f );
char *t = fp64_to_string( x, 17, 15 );
Serial.print( nf, 0 ); Serial.print( s ); Serial.print( ":" );
Serial.print( fixedLength( buf, x ) );
Serial.print( fixedLength( buf, delta ) );
Serial.print( "\t" );
Serial.print( " " ); Serial.print( xf, 7 );
Serial.print( " " );
if( deltaf >= 0.0 )
Serial.print( " ");
Serial.println( deltaf, 7 );
}
void loop() {
int cnt = 0;
while( 1 ) {
prev = x;
prevf = xf;
if( add ) {
float64_t t = fp64_inverse(n);
x = fp64_add( x, t );
xf = xf + 1.0/nf;
} else {
float64_t t = fp64_inverse(n);
x = fp64_sub( x, t );
xf = xf - 1.0/nf;
}
add = !add;
++cnt;
if( cnt == MAX_COUNT ) {
float64_t approx = fp64_ldexp( fp64_add( prev, x ), 1);
float approxf = ldexp( xf + prevf, 1 );
printLine( nf-2, prev, prevf, " low", false );
printLine( nf, x, xf, " high", false );
printLine( nf, approx, approxf," avg", true );
printLine( nf, pi, pif, " pi", true );
Serial.println();
cnt = 0;
}
n = fp64_add( n, two );
nf = nf + 2.0;
}
}
Let the Arduino run for a few minutes and take a look on the output. Here every 50000 iterations were calculated and the output shows the following:
99999 low: 0.78539316339745 -0.00000500000000 0.7853940 -0.0000042
100001 high: 0.78540316329745 0.00000499990000 0.7854040 0.0000058
100001 avg: 3.14159265338980 -0.00000000019999 3.1415958 0.0000031
100001 pi: 3.14159265358979 0000000000000000 3.1415927 0.0000000
199999 low: 0.78539566339744 -0.00000250000001 0.7853965 -0.0000017
200001 high: 0.78540066337244 0.00000249997499 0.7854015 0.0000033
200001 avg: 3.14159265353976 -0.00000000005003 3.1415958 0.0000031
200001 pi: 3.14159265358979 0000000000000000 3.1415927 0.0000000
299999 low: 0.78539649673077 -0.00000166666667 0.7853974 -0.0000008
300001 high: 0.78539983005299 0.00000166665554 0.7854008 0.0000026
300001 avg: 3.14159265356753 -0.00000000002225 3.1415963 0.0000036
300001 pi: 3.14159265358979 0000000000000000 3.1415927 0.0000000
The Leibniz series always creates a value that is either lower than PI/4 (indicated with “low” in the output) or higher than Pi/4 (indicated with “high”). So, ideally taking the average of the low and the high value, multiplying that by 4, should give a good approximation of PI, listed in the line with “avg”. And finally, the value to compare with is in the line tagged “pi”.
Already after 50000 iterations (n=100001), the value of Pi for fp64lib is 3.14159265338980, precise for 9 digits after the decimal point. The low and the high value are nicely distributed on both sides of Pi/4.
On the other side, our float calculation results in PI being 3.1415958, which is accurate to 5 digits only. And you can see, that all the small errors due to the reduced precisions are adding up, as both the low and the high value are above Pi/4, so no chance to come to a good approximation of PI.
This can be verified by taking a look on the result ofter after 500000 iterations:
999999 low: 0.78539766339745 -0.00000050000000 0.7853985 0.0000003
1000001 high: 0.78539866339645 0.00000049999900 0.7853995 0.0000013
1000001 avg: 3.14159265358781 -0.00000000000198 3.1415958 0.0000031
1000001 pi: 3.14159265358979 0000000000000000 3.1415927 0.0000000
Precision for fp64lib has improved from 9 digits to 11 digits, but it stayed at 5 digits for float.
Having said that, I just can repeat my inital comment: the above algorithm is not very well suited to calculate PI. And, fp64lib will also face its challenges as float did, only some million calculations later:
2043999 low: 0.78539791877916 -0.00000024461828 0.7853989 0.0000007
2044001 high: 0.78539840801571 0.00000024461827 0.7853993 0.0000011
2044001 avg: 3.14159265358976 -0.00000000000002 3.1415963 0.0000036
2044001 pi: 3.14159265358979 0000000000000000 3.1415927 0.0000000
This is the point where the best precision is achieved, with PI being accurately approximated to 13 digits. After that, rounding errors in fp64lib will also sum up, slowly, but noticeably:
4999999 low: 0.78539806339805 -0.00000009999939 0.7853991 0.0000009
5000001 high: 0.78539826339801 0.00000010000056 0.7853993 0.0000011
5000001 avg: 3.14159265359213 0.00000000000234 3.1415967 0.0000041
5000001 pi: 3.14159265358979 0000000000000000 3.1415927 0.0000000
5999999 low: 0.78539808006491 -0.00000008333254 0.7853991 0.0000009
6000001 high: 0.78539824673155 0.00000008333410 0.7853993 0.0000011
6000001 avg: 3.14159265359292 0.00000000000313 3.1415967 0.0000041
6000001 pi: 3.14159265358979 0000000000000000 3.1415927 0.0000000