Startseite → Pi berechnen → Kettenbruch in C programmiert

Kettenbruch in C programmiert

Vielen Dank an Andreas Noe für die folgenden Ausführungen und die beiden C-Programme.

Das Programm verwendet für die Näherung von Pi den Kettenbruch, den Johann Heinrich Lambert 1770 veröffentlichte.

Ausblenden

\frac{4}{\pi}=1+\frac{1^2}{3+\frac{2^2}{5+\frac{3^2}{\ddots}}}

Ein Teilbruch besteht also aus der Summe der n-ten ungeraden Zahl und dem Quotienten von n2 und dem nächsten Teilbruch. Für n = 1:

\frac{4}{\pi}=2\,n-1+\frac{n^2}{2\,(n+1)-1+\frac{(n+1)^2}{2\,(n+2)-1+\frac{(n+2)^2}{\ddots}}}

Diese rekursive Vorschrift kann mit einer rekursiven Funktion berechnet werden. Mit wachsender Rekursionstiefe werden die Änderungen am Gesamtergebnis stetig kleiner. Laut Wikipedia wächst die Genauigkeit des Ergebnisses mit jedem Schritt um ca. 0,76555 Stellen. Damit kann die Berechnung abgebrochen werden, wenn die gewünschte Genauigkeit erreicht ist.

In dieser Version des Programms wird die Berechnung mit dem Datentyp double durchgeführt. Ein double-Wert besteht aus 64Bit. Davon entfallen 52Bit auf die Mantisse, ein Bit auf das Vorzeichen und elf Bit auf den Exponenten. Die für double erreichbare Genauigkeit liegt also bei etwa:

\frac{\ln(2^{52})}{\ln(10)}=15,654

Dezimalstellen.

Für knapp 16 Dezimalstellen reicht es also aus, die Rekursionstiefe auf etwa:

\frac{15,654}{0,76555}=20,448

also 21 Teilbrüche zu begrenzen.

C-Code
  1. #include <stdio.h>
  2. #include <math.h>


  5. /**
  6.  @ brief *
  7.  Berechnet den n-ten Teilbruch des Kettenbruchs.
  8.  
  9.  @param n
  10.  Gibt den aktuell zu berechnenden Teilbruch an.
  11.  @param max_n
  12.  Gibt die maximale Rekursionstiefe an.
  13.  @return
  14.  Fuer n=1 ist der Rueckgabewert rund 4/Pi.
  15.  */
  16. double kettenbruch( int n, int max_n );


  19. //main_______________________________________________________________
  20. int main( int argc, char** argv )
  21.   // Pi naehern //////////////////////////////////////////////////////
  22.   double Pi = 4.0 / kettenbruch( 1, 21 );
  23.  
  24.   // Ergebnis mit 15 Dezimalstellen (14 Nachkommastellen)
  25.   // ausgeben.
  26.   printf( "%.14lf\n", Pi );
  27.  
  28.   return 0;
  29. }//main


  32. //kettenbruch________________________________________________________
  33. double kettenbruch( int n, int max_n ) {
  34.   double Result;
  35.  
  36.   if ( n > max_n ) {
  37.   // Die maximale Rekursiostiefe ist erreicht
  38.   // -> Abbruch
  39.   Result = 1.0;
  40.   }
  41.   else {
  42.   Result = 2*n-1 + n*n / kettenbruch( n+1, max_n );
  43.   }
  44.  
  45.   return Result;
  46. }//kettenbruch

Das zweite Programm funktioniert nach dem gleichen Prinzip, nur, dass der Datentyp double durch den Gleitkommatyp aus der GMP Bibliothek ersetzt wird.

In C kann man keine Strukturen als Rückgabetyp verwenden. Deshalb erfolgt die Rückgabe per Referenzparameter (Zeiger). Beim Aufruf des Programms gibt man an, wie viele Dezimalstellen berechnet werden sollen. Die letzte Stelle wird gerundet (leider). Unter Linux kann man etwa 40000 Dezimalstellen berechnen. Auf meinem guten alten Core2 Duo (Das Programm kann nur eine CPU verwenden.) braucht es dafür etwa 200s. In dem Bild habe ich ein paar Versuche aufgezeichnet.

Versuche

Unter Windows kann man etwa 19000 Stellen berechnen, bevor der Stack zu Ende geht. Unter Linux ist die Größe des Stacks ein Systemparameter und kann zur Laufzeit durch einen Systemaufruf geändert werden. Unter Windows gibt man die gewünschte Größe des Stacks dem Compiler als Parameter mit.

C-Code
  1. #include<stdlib.h>
  2. #include<stdio.h>
  3. #include<string.h>
  4. #include<gmp.h>
  5. #include<math.h>

  7. typedef unsigned long int ulong;

  9. ulong ulBitWidth; //!< Genauigkeit der Gleitkommazahlen.

  11. char *strUsageTemplate = "usage: %s DECIMAL_PLACES\n\n\DECIMAL_PLACES A decimal number greater than 0.\n";


  14. /**
  15.  @ brief *
  16.  Berechnet den n-ten Teilbruch des Kettenbruchs.
  17.  
  18.  @param n
  19.  Gibt den aktuell zu berechnenden Teilbruch an.
  20.  @param max_n
  21.  Gibt die maximale Rekursionstiefe an.
  22.  @param ptrResult
  23.  Referenzparameter fuer den Rueckgabewert. Fuer n=1 ist der
  24.  Rueckgabewert rund 4/Pi.
  25.  */
  26. void continuedFraction( ulong n,
  27.   ulong max_n,
  28.   mpf_t* ptrResult );


  31. //main_______________________________________________________________
  32. int
  33. main( int argc, char **argv ) {
  34.   ulong ulDecimalPlaces; //!< Anzahl der gewuenschten
  35.   //Dezimalstellen.
  36.   ulong ulRecursionDepth; //!< Benoetigte Rekursionstiefe.
  37.   mpf_t mpfPi; //!< Ergebnis.
  38.  
  39.   // Anzahl der gewuenschten Dezimalstellen ermitteln ////////////////
  40.   if ( argc < 2 ) {
  41.   fprintf( stderr,
  42.   strUsageTemplate,
  43.   argv[0]
  44.   );
  45.  
  46.   exit( 1 );
  47.   }
  48.  
  49.   if ( 1 > sscanf( argv[1], "%lu", &ulDecimalPlaces ) ) {
  50.   fprintf( stderr,
  51.   strUsageTemplate,
  52.   argv[0]
  53.   );
  54.  
  55.   exit( 2 );
  56.   }
  57.  
  58.   // Benoetigte Rekursionstiefe berechnen ////////////////////////////
  59.   ulRecursionDepth = 3+(ulong) ceil((double)ulDecimalPlaces/0.76555);
  60.  
  61.   // Benoetigte Genauigkeit der Gleitkommazahlen berechnen ///////////
  62.   ulBitWidth = (ulong) ceil(
  63.   (double)ulDecimalPlaces * log( 10.0 ) / log( 2.0 ) );
  64.  
  65.   /* Die Genauigkeit (Bitbreite der Mantisse) wird auf ein Vielfaches
  66.   * von 32 aufgerundet. **/
  67.   if ( ulBitWidth & 0x0000001F ) { // wenn 0 != ulBitWidth modulo 32
  68.   ulBitWidth >>= 5; // ulBitWidth = ulBitWidth / 32
  69.   ulBitWidth += 1; // ulBitWidth = ulBitWidth + 1
  70.   ulBitWidth <<= 5; // ulBitWidth = ulBitWidth * 32
  71.   }
  72.  
  73.   // Ergebnisvariable initialisieren ////////////////////////////////
  74.   mpf_init2( mpfPi, ulBitWidth );
  75.  
  76.   // Kettenbruch berechnen ( Zwischenergebnis = Pi/4 ) //////////////
  77.   continuedFraction( 1, ulRecursionDepth, &mpfPi );
  78.  
  79.   // Ergebnis = 4 / Zwischenergebnis = Pi ///////////////////////////
  80.   mpf_ui_div( mpfPi, 4, mpfPi );
  81.  
  82.   // Ergebnis ausgeben //////////////////////////////////////////////
  83.   // ( Die letzte Stelle wird von gmp_printf() gerundet. )
  84.   gmp_printf( "%.*Ff\n", ulDecimalPlaces-1, mpfPi );
  85.  
  86.   return 0;
  87. }//main


  90. //continuedFraction__________________________________________________
  91. void
  92. continuedFraction( ulong n,
  93.   ulong max_n,
  94.   mpf_t* ptrResult ) {
  95.  
  96.   if ( n == max_n ) {
  97.   // Die maximale Rekursiostiefe ist erreicht -> Abbruch //////////
  98.   mpf_set_ui( *ptrResult, 1 );
  99.   }
  100.   else {
  101.   // Lokale Variablen anlegen und initialisieren //////////////////
  102.   mpf_t mpfSummand;
  103.   mpf_t mpfFraction;
  104.   mpf_t mpfDenominator;
  105.  
  106.   mpf_init2( mpfSummand, ulBitWidth );
  107.   mpf_init2( mpfFraction, ulBitWidth );
  108.   mpf_init2( mpfDenominator, ulBitWidth );
  109.  
  110.   //Rueckgabewert berechnen/////////////////////////////////////////
  111.  
  112.   /*
  113.   * ptrResult = 2*n - 1 + n^2 / continuedFraction( n+1, max_n )
  114.   * \______________ ______________/
  115.   * V
  116.   * |
  117.   * +--> mpfDenominator
  118.   */
  119.   continuedFraction( n+1, max_n, &mpfDenominator );
  120.  
  121.   /*
  122.   * ptrResult = 2*n - 1 + n^2 / mpfDenominator
  123.   * \ /
  124.   * |
  125.   * +------------------------> mpfFraction
  126.   */
  127.   mpf_set_ui( mpfFraction, n*n );
  128.  
  129.   /*
  130.   * ptrResult = 2*n - 1 + mpfFraction / mpfDenominator
  131.   * \____________ _____________/
  132.   * V
  133.   * |
  134.   * +------------> mpfFraction
  135.   */
  136.   mpf_div( mpfFraction, mpfFraction, mpfDenominator );
  137.  
  138.   /*
  139.   * ptrResult = 2*n - 1 + mpfFraction
  140.   * |
  141.   * +---------------------------------> mpfSummand
  142.   */
  143.   mpf_set_ui( mpfSummand, n );
  144.  
  145.   /*
  146.   * ptrResult = 2*mpfSummand - 1 + mpfFraction
  147.   * \____ _____/
  148.   * V
  149.   * |
  150.   * +------------------------------> mpfSummand
  151.   */
  152.   mpf_mul_ui( mpfSummand, mpfSummand, 2 );
  153.  
  154.   /*
  155.   * ptrResult = mpfSummand - 1 + mpfFraction
  156.   * \_____ ______/
  157.   * V
  158.   * |
  159.   * +-----------------------------> mpfSummand
  160.   */
  161.   mpf_sub_ui( mpfSummand, mpfSummand, 1 );
  162.  
  163.   /*
  164.   * ptrResult = mpfSummand + mpfFraction
  165.   */
  166.   mpf_add( *ptrResult, mpfSummand, mpfFraction );
  167.   }
  168.  
  169. }//continuedFraction