Or maybe math really is this fun.
Okay, where this came from -- I gave my son a challenge, and he wrote a short program in perl to calculate π by the area of a circle.
So I wrote one in my favorite calculator, bc.
Now, its a nice program, and I can check the accuracy with bc itself, using the arctangent function. But bc is not as fast as C, and this old intel processor is not that fast, either.
And the program had a bug. Well, in my original explanation, I forgot to mention the problem of rounding. Bc always rounds down. Or, rather, it just drops fractional bits past the current scale This produced some notable effects, which are partially resolved in bc by doubling the scale.
Perhaps one reason bc is not as well known as it might be is that bc doesn't help you with rounding.
And it gives you scale instead of precision or accuracy.
So the limits of computers bite you pretty quick with bc, where other similar solutions help you ignore the limits -- sometimes until it's too late.
Which is part of the reason I'm writing this rant.
Anyway, I also wrote a nice little program for calculating π in standard C, using the language's ordinary math. It uses the same basic approach as I used in the bc version, except that the bc version I showed has a few more bells and whistles.
(I actually started with a very minimalistic version in bc, then wrote the C program and uploaded it to the pastebin linked above, then added the comments and extra features to the bc version before uploading it to a pastebin, if you are wondering why.)
So, if you want a full explanation of the method C source, look at my rants on the bc version (linked above) and the comments in the bc source. I'll limit my comments here to details of C programming and the standard C run-time environment, and issues of rounding and accuracy.
(Modifying the C source to do midpoint and so on is something I leave for you to play with.)
Here's the source of the standard C program, with extra comments:
/* Computing the value of pi by integrating the area of a (half of a) circle.
// by Joel Matthew Rees, 1 August 2013
// Copyright Joel Matthew Rees
//
// Fair use acknowledged.
//
// All rights to this expression of the method reserved.
//
// (Starting from scratch isn't that hard,
// and you'll like the results better.)
// -- Really, writing small programs like this from scratch
// is an excellent way to learn the basic grammar and runtime
// of a programming language. DIY!
//
//
// Compile with the right libraries:
// cc -Wall -lm -o halfpiArea halfpiArea.c
//
// C language programs have to be compiled. You can't just run them.
// You have to save the source as halfpiArea.c
// and then run the compiler command above in a shell:
// http://joels-programming-fun.blogspot.com/2013/08/get-shell.html
//
// The compiler will catch simple errors in the source,
// and try to give you hints on how to fix them.
// Learning how to read the hints is a topic for another rant or more.
//
// If the compiler succeeds, the command line above will save
// the executable program in a file called halfpiArea,
// which you then run by typing
// ./halfpiArea
// in the shell.
//
//
*/
/* The following lines help the compiler to figure out
// what functions you are using from the standard library
// and how you should be using them:
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* They literally include information about those functions
// into this program.
*/
/* The main function of a C program in a standard runtime environment
// is pretty much always called "main()".
// And it has two or three standard arguments.
// I only tell the program about two:
// argc is the integer count of the command line arguments,
// and argv[] is the array of command line arguments, as character strings.
// I'll tell you more about the command line arguments below.
*/
int main( int argc, char * argv[] )
{
/* All the variables here are declared as
// double length floating point.
// Accuracy is usually more than 11 useful digits,
// Exponent range is usually more than about 10^(+-235).
// Simplify way too much,
// you can write numbers from 1e-235 to 1e235,
// and, if we get the math write,
// we can guarantee at least 11 good digits for pi.
//
// (Actually, the math normally supported in modern C libraries
// is IEEE 754, which would be 10^(+-308) and at least 15 good digits.
// Check C Language Data Types and IEEE 754 on Wikipedia, stackexchange,
// etc., for more information.
//
// If you want to know what your compiler gives,
// #include <limits.h> and play with the math.)
*/
double area = 0.0; /* Sum the area here. */
double deltaX = 0.025; /* Width of the rectangles (by default 1/40). */
double lastX = -1.0; /* Track the left side. */
double x; /* The right side. */
double almostPi; /* The approximation of pi. */
/* You will notice that we declare the names and types
// of the variables before we use them.
// C also allows us to set the initial value as we declare them.
*/
if ( argc > 1 ) /* Command line arguments? */
{ deltaX = strtod( argv[ 1 ], NULL ); /* Assume it's a different width. */
}
/* This allows us to set a different width on the command line.
// 1/11 would be 0.09090909090909090909 to twenty decimal places:
// ./halfpiArea 0.09090909090909090909
// 1/8 == 0.125, 1/10 == 0.1, 1/100 == 0.01, 1/200000 == 0.000005
*/
/* Note that starting from the half delta converges at least one decimal digit faster.
// (That's because the height of the rectangle crosses the circle at midpoint,
// rather than left edge or right. Draw the pictures to see why.
// This is right edge, but changing to midpoint or left edge is trivial.)
*/
for ( x = -1.0 + deltaX; x <= 1.0; x += deltaX )
{
double ysquared = 1.0 - ( x * x );
double y = sqrt( ysquared );
/* deltaX = x - lastX; Reminder for humans, only! */
lastX = x;
area += deltaX * y;
/* Suppress long listings of the areas.*/
if ( deltaX > 0.00005 ) /* 1/20000 */
{ printf( "(%17.15g,%17.15g): %17.15g\n", x, y, area );
/* %17.15g is a floating point field
// of 17 digits with 15 digits to the right of the decimal.
*/
}
}
almostPi = area * 2;
printf( "Almost Pi: %17.15g\n", almostPi );
return EXIT_SUCCESS;
}
Calculating π to eleven digits works okay with C doubles. But if you look at the listing of the areas, you see some odd differences, because C always rounds but bc always does as many exact digits as you tell it to with the scale variable.
Play with different widths, compare to the bc version, and see what you can learn.
(Modern C, from about C99, allows you to declare floating point types that use alternative rounding.)
And you should notice that C is very fast when compared to bc.
When you want to check the real value of π, remember:
echo "scale=41; a(1)*4;" | bc -lfrom the shell command line, or
scale=41in bc, adjust scale to one more than the number of good digits you want.
a(1)*4
And don't try
./halfpiArea 0.0000000000000000000000000000000000000001unless you really want to watch your computer contemplate its navel until you give up and shut it down. Start with
./halfpiArea 0.00001and work your way up to 0.0000001 or so, where you get really noticeable waits.
And see if you can figure out whether midpoint is really better than right side, as I said it would be, when you start using widths smaller than 1/100 or so.
Oh. And play with the source code a bit, too.
[JMR201611182140:
I decided to post a few hints about playing with the programs:
http://joels-programming-fun.blogspot.com/2016/11/using-to-check-your-floating-point.html.
]
[JMR201611292224:
Here it is, demonstrating the principle with GMP:
http://joels-programming-fun.blogspot.com/2016/11/using-gmp-to-approximate-pi.html.
]
