/* 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;
}