# Pastebin: halfpiArea_gmp.c

A C program using the GMP library to demonstrate simple numerical calculus methods for generating a poor approximation of π, without a lot of command-line control.

Format
C
Post date
2016-11-19 11:19
Publication Period
Unlimited
` /* 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; }`