Frage der Präzision? (pi, gmp)

Swordfish

Hi Volks!

Ich hab' mir gestern und heute auch zufällig den Kopf über die Berechnung von pi zerbrochen und das Resultat meiner Bemühungen war die Erkenntnis:

pi = 16 arctan(1/5) - 4 arctan(1/239)

Danach hab' ich mir 'ne mp Library für C gesucht und bin auf GMP 4.1.2 gestoßen. Mein Problem ist jedoch, daß ab

3,141592653589793        2384626433832795... // pi
3,141592653589793 [hier] 4107032954041479... // 'mein' pi

nicht's mehr stimmt.
Ich denk mal, es liegt an der Prezision meiner variablen.

Meine Fragen daher:

Welche Prezision für C brauche ich, wenn C das Resultat einer Potenz (64 bit float)^(32 bit integer) ist?
Welche Prezision für T brauche ich, wenn T das Resultat einer Division (Prezision C?) / (32 bit integer) ist?
Welche Prezision für ArcTan brauch ich, wenn ArcTan das Resultat einer Addition (Prezision ArcTan?) / (Prezision T) ist?

Vielleicht kann jemand licht in's Dunkel bringen.

THXIA

Swordfish

Walli

Welches arctan benutzt du denn?

Swordfish

mein eigenes!?

arctan(x) definiert sich doch gerade durch

x^3   x^5   x^7   x^9
  arctan(x) = x - --- + --- - --- + --- ...
                   3     5     7     9

Walli

Dann liegts evtl. an deiner Implementierung. Wie siehts aus wenn du ein paar mehr Iterationsschritte machst?

Btw: Hat gmp keine eigene arctan Implementierung?

Swordfish

Ja, gmp müsst schon 'ne Implementierung dafür haben, aber das nutzt mir nichts wenn ich 16 arctan(1/5) - 4 arctan(1/239) beim konvergieren gegen pi zusehen möchte.

Ich hab's mal so implementiert:

dStartTime = GetTime();

// make sure that strPi is different from strOldPi
// ...preliminary code and topic to change ;)
strPi[0] = 1;	strOldPi[0] = 2;
strPi[1] = 0;	strOldPi[1] = 0;

for(__int32 i = 0, z = 3; strcmp(strPi, strOldPi); i++, z += 2)
{
  free(strPi);
  free(strOldPi);

  mpf_pow_ui(mpfC1, mpfX1, z);
  mpf_pow_ui(mpfC2, mpfX2, z);
  mpf_div_ui(mpfT1, mpfC1, z);
  mpf_div_ui(mpfT2, mpfC2, z);

  if(blnSign)
  {
    mpf_add(mpfArcTan1, mpfArcTan1, mpfT1);
    mpf_add(mpfArcTan2, mpfArcTan2, mpfT2);
  } else {
    mpf_sub(mpfArcTan1, mpfArcTan1, mpfT1);
    mpf_sub(mpfArcTan2, mpfArcTan2, mpfT2);			
  }

  blnSign = !blnSign;

  // Save old pi for comparison:
  strOldPi = mpf_get_str(NULL, &eExp, 10, LONG_MAX, mpfPi);

  mpf_mul_ui(mpfPP1, mpfArcTan1, 16);
  mpf_mul_ui(mpfPP2, mpfArcTan2, 4);
  mpf_sub(mpfPi, mpfPP1, mpfPP2);

  strPi = mpf_get_str(NULL, &eExp, 10, LONG_MAX, mpfPi);

  szPiLength     = strlen(strPi);
  szOldPiLength  = strlen(strOldPi);
  szOldLastPiPos = szLastPiPos;

  for(szLastPiPos = szOldLastPiPos; szLastPiPos < szPiLength && szLastPiPos < szOldPiLength; szLastPiPos++)
    if(strPi[szLastPiPos] != strOldPi[szLastPiPos]) break;

  if(szLastPiPos)
  {
    printf("New digits found: %i\n", szLastPiPos - szOldLastPiPos);
    for(p = szOldLastPiPos; p < szLastPiPos; p++)
    {
      printf("Digit %i: %c  Elapsed time: %.3f\n", p, strPi[p], GetTime()-dStartTime);
    }
  }
}