einzelne Stelle von Pi berechnen

[MaR-V-iN]

Hallo leute...

Sitz hier grad an etwas für die schule und hab so nen paar probleme.
Wie die überschrift schon sagt geht es darum einzelne Stellen von Pi zu berechen. Eigentliche Aufgabe ist zwar Pi auf 100 Nachkommastellen zu berechen, aber da Delphi keinen Datentyp hat der 100 nachkommastellen bietet, und ich es gleich auch für theoretisch 1000 und mehr machen will, bräuchte ich mal ne berechnungsmethode für nur eine einzelne stelle.
Die würde ich dann nacheinander für alle gewünschten stellen ausführen...
Bei Wikipedia hab ich diese Formel gefunden:



Man kann damit irgenwie einzelne Stellen von Pi berechnen aber ich hab keine idee wie...

Hoffe ihr könnt mir weiterhelfen...

MaR-V-iN

PS: auch ein verständlicher c++/basic code wäre OK. Könnte den selber umsetzen

[GrafZahl]

http://www.experimentalmath.info/bbp-codes/bbp-alg.pdf

seite 2 ab ca 1/3 der seite ... für die unendliche summe der Sj nehmen die im beispiel code die ersten 100 summanden ... der fehler der dadurch entsteht macht sich erst einige stellen später bemerkbar, so dass man bei "normalen" 64bit floating point operationen im rechner in der regel 9 korrekte stellen bekommt ... wenn du dich darauf verlässt dass 6 stellen korrekt sind, solltest du auf der sicheren seite sein...


beispiel implementationen finden sich hier:
http://www.experimentalmath.info/bbp-codes/

[MaR-V-iN]

Das was ich verstanden hab sah sehr vielversprechend aus...
aber ich hab immernoch nich so recht verstanden wie ich das jetzt anwende...
Ach und wenn es mir jetz jemand erklärt, bitte in "deutsch", d.h. so das ich es verstehe (bin in der 9. Klasse)...

MaR-V-iN

[GrafZahl]

ok, schau dir das PDF nochmal an. dort gibt es 2 für deine anwendung relevante gleichungen auf seite 2 ... die gleichungen (4) und (6)

gleichung 4 sagt dir, wie du dein ergebnis aus den teilergebnissen zusammensetzt: du berechnest die {16^d Sj} und setzt sie in diese gleichung ein... d ist die stelle ab der du berechnen willst, j ist ein index ... du brauchst die werte für die indizies 1, 4, 5 und 6

gleichung 6 sagt dir, wie die {16^d Sj} berechnet werden (das forum hier schein kein latex zu können, daher musst du für die schöne formel ins pdf schauen)

gemeint ist die letzte zeile auf seite 2:

die formel besteht im wesentlichen aus 2 summen (das große griechische Sigma) unterhalb eines Summenzeichens steht der zugehörige index ... in diesem fall k. die summe bezieht sich auf den bruch hinter dem summenzeichen. das summen zeichen darfst du dir analog zur for schleife in C vorstellen: die schleife läuft von k=0 bis k=d, und wertet immer wieder den bruch aus, nur halt mit dem jeweils aktuellen k ... d ist wieder die stelle die du ausrechnen willst, j ist wieder der jeweilige wert aus {1,4,5,6}

danach folgt eine zweite summe, die eigentlich unendlich viele summanden hat, sie funktioniert genauso wie die erste ... für die praktische anwendung denk dir anstelle von "unendlich"(die liegende 8 über dem Sigma) einfach "d+100" ... da du diese summe nicht vollständig berechnest, entsteht ein fehler, der dir bei "d+100" das ergebnis ab ca. der 9 stelle verfälscht (einen derartigen nachteil haben alle numerischen verfahren) ... du kannst aber normalerweise davon ausgehen, dass alles bis zur 6. stelle richtig ist, und dann das ganze für die d+6. stelle wieder neu ansetzen ...

das größte problem beim rechnen wird vermutlich sein dass du große exponenten bekommst ... du möchtest kein 16 hoch 10000 berechnen wenn du danach sowieso wieder mod 8k+j reduzierst ... das würde a) viel zu lange dauern, und b) dein ergebnis verfälschen, da derartig große zahlen nicht mehr vernünftig im rechner dargestellt werden können (die zahl hat 40000 bit) ... als fließkomma zahl kann man die zwar darstellen, aber nicht exakt, was wiederum einen fehler bedeutet, der nach der modulo reduktion nicht mehr unbedeutend ist ...

die lösung dieses problems heißt square-and-multiply algorithmus ... der ist bei wikipedia zu finden und beruht darauf, dass z.b. 16^100 == (((((((16^2)*16)^2)^2)^2)*16)^2)^2

im dummen fall multipliziert man 100 mal 16 mit sich selbst ... im schlauen fall hat man 8 multiplikationen... das ist definitiv schneller und hat den vorteil, dass man zwischendurch an jeder beliebigen stelle immer wieder die modulo reduktion durchführen kann, ohne dass sich am endergebnis etwas verändern würde ... die zahlen bleiben "relativ" klein, und der darstellungsfehler einer fließkommazahl wird minimiert ...

so weit so toll?

[MaR-V-iN]

Zuerst einmal Danke... das war sehr informativ...
Ich hab versucht das so umzusetzen, allerdings hab ich das Gefühl das noch was nicht will...
für die Berechnung von {16^d Sj}:

Code:

function formel6(j,d:Integer):Extended;
var k: Integer;
var res: Extended;
var x: Variant;
begin
res:=0;
    for k := 0 to d do;
    begin
         x:=Power(16, d-k);
         res := res + ((x mod (8*k+j))/(8*k+j));
    end;
    for k := d+1 to d+100 do;
    begin
        res := res + ((Power(16, d-k)/(8*k+j)));
    end;
    Result:=res;
end;

und dann für stelle 1:

Code:

y := 4 * formel6(1, 1) - 2 * formel6(4, 1) - formel6(5, 1) - formel6(6, 1);

Stimmt das soweit??
das ergebniss ist dann aber:

Code:

5,41306863451018e-127

und das kommt mir zu klein vor...

Also was ist falsch??

MaR-V-iN

[GrafZahl]

ahhh ... mein fehler

ich habe vergessen zu erwähnen was die großen geschweiften klammern bedeuten: es zählt nur der bruchtanteil

soll heißen: {5,8435346543} = 0,8435346543

oder auch: Ergebnis = Ausdruck - (int) Ausdruck ...

edit:

ich war mal so frei das in c zu implementieren und zu kommentieren...

Code:

#include <stdio.h>
#include <math.h>

#define step 1
#define maxcol 50
#define sp 10
  

// funktion pow16
// berechnet 16 hoch x mod n
long double pow16(int x, int n)
{
       int msb;           //position des most significant bit von x
       int mask;          //eine bitmaske
       double foo;        //ohne diesen umweg macht mein compiler leider sche*** ...
       long double result;

       //speziealfälle
       if(x==0&&n==1)return 0;
       if(x==0)return 1;
       if(x==1)
       {
           return 16-n*(int)(16/n);
       }
      
       foo=(log(x)/log(2));        //logarithmus von x zur basis 2, abgerundet:
                                   //höchstwertige bitstelle von x
                                   //(bei 0 anfangen zu zählen da 2 hoch 0 = 1)
       msb=foo;
      
       mask=1<<msb-1;              //wir setzen bit msb-1 von x in unserer maske
       result=16;                  //initialisierung square and multiply algorithmus
      
                                   // es folgt square and multiply...
       for(;mask;mask=mask>>1)
       {
           result*=result;//square ... bei jedem schleifendurchlauf
           if(x&mask)
               result*=16; //multiply ... nur wenn dass entsprechende bit im
                           //exponenten gesetzt ist (siehe mask)
           result-=n*(int)(result/n);//mod
       }
       result-=n*(int)(result/n);//mod
       return result;
}

// funktion reihe
// berechnet die reihenformel für die ziffern die auf die stelle d folgen
// j bezeichnet den index der BBP formel
long double reihe(int d, int j)
{
       int k;
       long double result,jk,term;
      
       result=(long double)0;
      
       for(k=0;k<d;k++)//erste summe
       /*
       anmerkung:
       laut der summenformel sollte hier k<=d als bedinung stehen...
       k<d liefert allerdings das korrekte ergebnis, von daher gehe
       ich davon aus, dass sich in die formel im pdf ein tippfehler
       eingeschlichen hat, und der laufindex nur bis d-1 laufen sollte
       anstatt bis d ...
       */
       {
           jk=8*k+j;
           result += pow16(d-k,jk)/jk;
           result -= (int)result; // auf gebrochenen anteil reduzieren
       }
       term=1.;
       for(;term>1e-20;k++)//zweite summe
       {
           jk=8*k+j;
           term=pow(16.,(double)d-k)/jk;
           result += term;
           result -= (int)result; // auf gebrochenen anteil reduzieren
       }
       return result;
}

int main(int argc, char *argv[])
{
  int d,idx,col;
  long double frac,s1,s4,s5,s6,pi;
  char hex[]="0123456789ABCDEF";
    
  d=0; //wir hätten gerne die berechnung der stellen nach stelle 0
       //also ab der ersten dezimalstelle...
    
  printf("berechne die ersten 1000 hexadezimalen Nachkommastellen von Pi:\n\n");

  putchar('3');  
  putchar(',');  
  col=0;
  while(1)
  {
  
      //die Sj nach gleichung 6
      s1=reihe(d,1);
      s4=reihe(d,4);
      s5=reihe(d,5);
      s6=reihe(d,6);
      
      frac=4*s1-2*s4-s5-s6;//ergebnis zusammenbauen nach gleichung 2
      
      //ergebnis in den richtigen bereich verschieben...
      frac-=(int)frac;
      if(frac<0)
            frac+=1;
        /*
          anmerkung:
          unser ergebnis ist aus der modulorechnung hervorgegangen, repräsentiert also eine restklasse.
          diese restklasse hat unendlich viele mögliche repräsentanten, also werte die bezüglich dieser
          restklasse kongruent (umgangssprachlich bedeutet "kongruent" soviel wie "gleichbedeutend")
          sind ... da wir das suchen was bei pi hinter dem komma steht, wissen wir: das gesuchte ergebnis
          liegt zwischen 0 und 1 ... somit wählen wir den vertreter der restklasse, der zwischen 0 und 1
          liegt ...
        */

      for(idx=0;idx<step;idx++)
      {
                            
          if(col%sp==0&&col!=0&&col!=maxcol)
          {//alle 10 zeichen ein space
              putchar(' ');
              col++;
              fflush(stdin);
          }
          if(col>maxcol)
          {//zeilenumbruch
              putchar('\n');
              putchar(' ');
              putchar(' ');
              col=0;
              fflush(stdin);
          }
          frac*=16;
          putchar(hex[(int)frac]);
          col++;
          frac-=(int)frac;
      }
      d+=step;
      if(d>=1000) break;
  }
  printf("\n\n");
  
  system("PAUSE");    
  return 0;
}

für d > 1400 gibts irgendwo nen überlauf, der zum crash führt ...