#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;
}
