Using π to Check Your Floating Point Accuracy

In my post on using areas of rectangles to approximate π in C, I blithely asserted that standard C floating point math would lose accuracy compared to bc. But I didn't explain how I came to that conclusion.

I'll show you a little bit about how I got there, tiptoeing around some of the steps where it's easy to get lost.

First, refresh your memory about the bc program:

• The theory
• Working around loss of precision

And here's how we'll call the program from the command line:

echo "q=1;m=32;scale=41;" | cat - halfpiArea_by_m.bc | bc -l

The name for the number of rectangles, "m" was rather awkward, but "q" for quiet should be easy to remember.

I suppose you may be wondering how that works.

First, echo just echoes what it is given into the standard output stream.

Second, the output stream of echo is piped to the standard input stream for cat and the command line tells cat to continue with the contents of the source code of the bc program. This way, cat concatenates things so that the result is a program that sets q, m, and scale appropriately before it begins:

q=1;m=32;scale=41;
/* ... */
if ( m < 1 ) { m = 11; }
s=0;
w = 1/m;

Unless we set them first, q, m, and d start at zero, and scale starts at 20 when we give bc the "-l" option.

So cat assembles the lines to set q, m, and scale with the program source and sends it all to the standard output stream.

Then the command line pipes the output of "cat" to "bc -l", and the program starts with scale, q, and m set to 41, 32, and 1.

Setting m changes the number of rectangles to divide the circle into, thus changing the width of the rectangles. We'll change that and examine the resulting attempt to calculate π.

The C program, halfpiArea, accepts a width as a parameter.

(That was not the best choice. You can change it if you want. I'm using the program as it is.)

The width is the inverse of the number of rectangles because the radius of the circle is 1. So we can calculate the width using bc:

echo "1/32" | bc -l

which tell us that one thirty-second is .03125000000000000000 (0.03125). We can select the decimal fraction that bc gives us on the terminal screen and copy it. In a Linux OS shell, we can usually use right-click after selecting the text. In a Mac OS X shell, the usual copy keyboard shortcut works.

In MSWindows, I think you will have to click the application icon in the left of the title bar and select "Edit" from the menu that pops up, then select "Copy" from the submenu. ;-*

Then we can run the C program by typing the program name and pasting the fraction in after as follows:

./halfpiArea .03125000000000000000

Except we want to run the bc program after the C program, so that we can get the output of both on the terminal screen at the same time:



Now we can select the output estimates for π and copy them to a text editor, where we can compare them at leisure.

We will notice that the two programs produce nearly the same results for counts of rectangles that are powers of two. Two rectangles, for example:

jwr@fun:~/work/mathgames/piArea$ ./halfpiArea .5
(             -0.5,0.866025403784439): 0.433012701892219
(                0,                1): 0.933012701892219
(              0.5,0.866025403784439):  1.36602540378444
(                1,                0):  1.36602540378444
Almost Pi:  2.73205080756888
jwr@fun:~/work/mathgames/piArea$ echo "q=1;m=2;scale=41;" | cat - halfpiArea_by_m.bc | bc -l
3: (1.00000000000000000000000000000000000000000, 0): 1.3660254037844\
3864676372317075293618347140
2.73205080756887729352744634150587236694280

These are the same, except for the last digit of the C program, which is actually correctly rounded, rather than just different:

2.73205080756888
2.73205080756887729352744634150587236694280

Similarly, 2048 rectangles:

jwr@fun:~/work/mathgames/piArea$ echo "1/2048" | bc -l
.00048828125000000000
jwr@fun:~/work/mathgames/piArea$ ./halfpiArea .00048828125000000000
(   -0.99951171875,0.0312461850698753): 1.52569263036501e-05
...

(                1,                0):  1.57078998270572
Almost Pi:  3.14157996541144
jwr@fun:~/work/mathgames/piArea$ echo "q=1;m=2048;scale=41;" | cat - halfpiArea_by_m.bc | bc -l
2049: (1.00000000000000000000000000000000000000000, 0): 1.5707899827\
0572260019288624147880158332158
3.14157996541144520038577248295760316664316

which are the same to the last digit of the C output:

3.14157996541144
3.14157996541144520038577248295760316664316

In fact, the C program is well behaved all the way to 2¹⁶:

jwr@fun:~/work/mathgames/piArea$ echo "2^16" | bc -l
65536
jwr@fun:~/work/mathgames/piArea$ echo "1/(2^16)" | bc -l
.00001525878906250000
jwr@fun:~/work/mathgames/piArea$ ./halfpiArea .00001525878906250000
Almost Pi:  3.14159258349584
jwr@fun:~/work/mathgames/piArea$ echo "q=1;m=2^16;scale=41;" | cat - halfpiArea_by_m.bc | bc -l
65537: (1.00000000000000000000000000000000000000000, 0): 1.570796291\
74791608251698437807118686209362
3.14159258349583216503396875614237372418724
jwr@fun:~/work/mathgames/piArea$ echo "q=1;m=2^16;scale=81;" | cat - halfpiArea_by_m.bc | bc -l
65537: (1.0000000000000000000000000000000000000000000000000000000000\
00000000000000000000000, 0): 1.5707962917479160825169843780711868627\
48822128541338008872254934835925225596242002
3.141592583495832165033968756142373725497644257082676017744509869671\
850451192484004

(Going to scale=81 did take a minute or four. I should time it:

scale=81
a(1)*4
3.141592653589793238462643383279502884197169399375105820974944592307816406286208996

Okay, three minutes, forty seconds.)

Here are the estimates:

3.14159258349584
3.14159258349583216503396875614237372418724
3.141592583495832165033968756142373725497644257082676017744509869671850451192484004 

The C program produces the same as the bc program to the C program's second-to-the-last digit.

Note that the bc program is showing precision loss out in the 37th digit. Again, remember why, from the fly-in-the-pi post.

And, just as a reminder (in bc -l),

scale=81
a(1)*4
3.141592653589793238462643383279502884197169399375105820974944592307816406286208996

So, even with all this effort, we have only calculated π to the sixth digit past the decimal point, which is also expected, see the theory post.

That's not bad results, really, but if we don't choose the number of rectangles carefully, the C program output starts varying from the bc program quite quickly.

Ten, 100, and 1000 are good, but 10,000 is not. 11 is good, but 23 is not. Try a few, and see if you can guess what numbers work better and which work worse, and why.

Now, if our only purpose were to calculate π to six decimal digits, we'd be satisfied with using 2¹⁶ for the count and be satisfied with standard C floating point math.

But, of course, that's not the real purpose.

What we are really talking about here is a way to get good solutions to equations involving integral calculus. (Was I being sneaky for not warning you in advance?)

So we do care about the differences, sometimes. We can't always select the convenient counts of rectangles to estimate with. Which is why I'm going to demonstrate using the Gnu Multi-Precision library real-soon-now.