/* Computing the value of pi by integrating the area of a (half of a) circle. // by Joel Matthew Rees, 1 August 2013 through 19 November 2016 // 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.) // // Explanations here: // http://joels-programming-fun.blogspot.com/2016/11/using-gmp-to-approximate-pi.html // // This is based on a program in bc: // https://osdn.net/users/reiisi/pastebin/4144 // and a program using the standard C floating point math: // https://osdn.net/users/reiisi/pastebin/4098 // // Compile with the right libraries: // cc -Wall -lgmp -o halfpiArea_gmp halfpiArea_gmp.c */ #include <stdio.h> #include <stdlib.h> #include <float.h> #include <math.h> #include <gmp.h> /* NOT a standard header! */ /* https://gmplib.org/manual/Initializing-Floats.html // https://gmplib.org/manual/GMP-Basics.html#GMP-Basics */ int main( int argc, char * argv[] ) { /* Declare the gmp/mpf variables and allocate the pointers: */ mpf_t area; mpf_t deltaX; mpf_t lastX; mpf_t x; mpf_t almostPi; mpf_t limit; mpf_t xsquared, ysquared, y, deltaA; /* Set the default precision: // log2(10^20) == 66.44. (About 20 digits, like bc -l, I guess.) // But gmp handles rounding by default, // so we get better results than with bc. */ mpf_set_default_prec ( 67 ); /* At least 67 good bits (but actually more). */ /* Actually allocate and initialize the values of the gmp/mpf variables: */ mpf_init_set_str( area, "0.0", 10 ); mpf_init_set_str( deltaX, "0.025", 10 ); mpf_init_set_str( lastX, "-1.0", 10 ); mpf_init_set( x, lastX ); mpf_init_set( almostPi, area ); mpf_init_set_str( limit, "1.0", 10 ); /* Allocations only: */ mpf_init( xsquared ); mpf_init( ysquared ); mpf_init( y ); mpf_init( deltaA ); /* Get the width argument, if there is one. // A width argument was not the best design choice, really. */ if ( argc > 1 ) { if ( mpf_set_str( deltaX, argv[ 1 ], 10 ) != 0 ) /* Had to specify it, base 10 only. */ { printf( "unrecognized deltaX: %s\n", argv[ 1 ] ); return EXIT_FAILURE; } } /* Caculate the area of the unit half-circle from 0 to 180 degrees: // (From 0 radians to pi radians, but saying so might be confusing.) // // (Note that starting from the half delta converges // at least one decimal digit faster for small counts. // // 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. // Except, since we go from -1 to 1 and do the full half-circle, // we actually get back some of the loss of precision, by accident. // // This is right edge, but changing to midpoint or left edge is trivial.) */ for ( ; mpf_cmp( x, limit ) <= 0; mpf_add( x, x, deltaX ) ) { mpf_mul( xsquared, x, x ); mpf_ui_sub( ysquared, 1L, xsquared ); mpf_sqrt( y, ysquared ); /* deltaX == x - lastX; */ mpf_set( lastX, x ); mpf_mul( deltaA, deltaX, y ); mpf_add( area, area, deltaA ); if ( mpf_cmp_d( deltaX, 0.0005 ) > 0 ) { gmp_printf( "(%21.20Ff,%21.20Ff): %21.20Ff\n", x, y, area ); } } /* pi is the area of the whole circle */ mpf_mul_ui( almostPi, area, 2L ); gmp_printf( "Almost Pi: %25.20Ff\n", almostPi ); /* No need to clear stuff before ending program. */ return EXIT_SUCCESS; }