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
  1. /* Computing the value of pi by integrating the area of a (half of a) circle.
  2. // by Joel Matthew Rees, 1 August 2013 through 19 November 2016
  3. // Copyright Joel Matthew Rees
  4. //
  5. // Fair use acknowledged.
  6. //
  7. // All rights to this expression of the method reserved.
  8. //
  9. // (Starting from scratch isn't that hard,
  10. // and you'll like the results better.)
  11. //
  12. // Explanations here:
  13. // http://joels-programming-fun.blogspot.com/2016/11/using-gmp-to-approximate-pi.html
  14. //
  15. // This is based on a program in bc:
  16. // https://osdn.net/users/reiisi/pastebin/4144
  17. // and a program using the standard C floating point math:
  18. // https://osdn.net/users/reiisi/pastebin/4098
  19. //
  20. // Compile with the right libraries:
  21. // cc -Wall -lgmp -o halfpiArea_gmp halfpiArea_gmp.c
  22. */
  23. #include <stdio.h>
  24. #include <stdlib.h>
  25. #include <float.h>
  26. #include <math.h>
  27. #include <gmp.h> /* NOT a standard header! */
  28. /* https://gmplib.org/manual/Initializing-Floats.html
  29. // https://gmplib.org/manual/GMP-Basics.html#GMP-Basics
  30. */
  31. int main( int argc, char * argv[] )
  32. {
  33. /* Declare the gmp/mpf variables and allocate the pointers: */
  34. mpf_t area;
  35. mpf_t deltaX;
  36. mpf_t lastX;
  37. mpf_t x;
  38. mpf_t almostPi;
  39. mpf_t limit;
  40. mpf_t xsquared, ysquared, y, deltaA;
  41. /* Set the default precision:
  42. // log2(10^20) == 66.44. (About 20 digits, like bc -l, I guess.)
  43. // But gmp handles rounding by default,
  44. // so we get better results than with bc.
  45. */
  46. mpf_set_default_prec ( 67 ); /* At least 67 good bits (but actually more). */
  47. /* Actually allocate and initialize the values of the gmp/mpf variables: */
  48. mpf_init_set_str( area, "0.0", 10 );
  49. mpf_init_set_str( deltaX, "0.025", 10 );
  50. mpf_init_set_str( lastX, "-1.0", 10 );
  51. mpf_init_set( x, lastX );
  52. mpf_init_set( almostPi, area );
  53. mpf_init_set_str( limit, "1.0", 10 );
  54. /* Allocations only: */
  55. mpf_init( xsquared );
  56. mpf_init( ysquared );
  57. mpf_init( y );
  58. mpf_init( deltaA );
  59. /* Get the width argument, if there is one.
  60. // A width argument was not the best design choice, really.
  61. */
  62. if ( argc > 1 )
  63. { if ( mpf_set_str( deltaX, argv[ 1 ], 10 ) != 0 ) /* Had to specify it, base 10 only. */
  64. { printf( "unrecognized deltaX: %s\n", argv[ 1 ] );
  65. return EXIT_FAILURE;
  66. }
  67. }
  68. /* Caculate the area of the unit half-circle from 0 to 180 degrees:
  69. // (From 0 radians to pi radians, but saying so might be confusing.)
  70. //
  71. // (Note that starting from the half delta converges
  72. // at least one decimal digit faster for small counts.
  73. //
  74. // That's because the height of the rectangle crosses the circle at midpoint,
  75. // rather than left edge or right. Draw the pictures to see why.
  76. // Except, since we go from -1 to 1 and do the full half-circle,
  77. // we actually get back some of the loss of precision, by accident.
  78. //
  79. // This is right edge, but changing to midpoint or left edge is trivial.)
  80. */
  81. for ( ; mpf_cmp( x, limit ) <= 0; mpf_add( x, x, deltaX ) )
  82. {
  83. mpf_mul( xsquared, x, x );
  84. mpf_ui_sub( ysquared, 1L, xsquared );
  85. mpf_sqrt( y, ysquared );
  86. /* deltaX == x - lastX; */
  87. mpf_set( lastX, x );
  88. mpf_mul( deltaA, deltaX, y );
  89. mpf_add( area, area, deltaA );
  90. if ( mpf_cmp_d( deltaX, 0.0005 ) > 0 )
  91. { gmp_printf( "(%21.20Ff,%21.20Ff): %21.20Ff\n", x, y, area );
  92. }
  93. }
  94. /* pi is the area of the whole circle */
  95. mpf_mul_ui( almostPi, area, 2L );
  96. gmp_printf( "Almost Pi: %25.20Ff\n", almostPi );
  97. /* No need to clear stuff before ending program. */
  98. return EXIT_SUCCESS;
  99. }
Download Printable view

URL of this paste

Embed with JavaScript

Embed with iframe

Raw text