//
// C++ transcription of the 
// fortran program piqp (by peter borwein ?)
// with a slight speedup in the function expm()
// added binary printing
//
// done by joerg arndt
//

#include <math.h>
#include <stdlib.h>
#include <iostream.h>


const int ntp=55;
double ic;
int    itp;
double tp[ntp];  //  hold powers of 2


double expm(double &p, double ak);
double series(double m, double &ic);
void   ihex(double x, int bq);



int main(int argc, char **argv)
//
// this program employs the recently discovered digit extraction scheme
// to produce hex digits of pi.
// this code is valid up to ic=2^24 on systems with IEEE arithmetic
//
{
double pid,s1,s2,s3,s4;

    tp[0]=1.0;

    for(int k=1; k<ntp; ++k) // setup table with powers of 2
        tp[k]=2.0*tp[k-1];


    if(argc<=1)
    {
        ic=0;
    }
    else
    {
        ic=atoi(argv[1]); // precision from arg
        ic=fabs(ic);      // force positive
    }
 


    cout<<"\n series 1:"<<endl;
    s1=series(1,ic);
    cout<<" s1="<<s1<<endl;

    cout<<"\n series 2:"<<endl;
    s2=series(4,ic);
    cout<<" s2="<<s2<<endl;

    cout<<"\n series 3:"<<endl;
    s3=series(5,ic);
    cout<<" s3="<<s3<<endl;

    cout<<"\n series 4:"<<endl;
    s4=series(6,ic);
    cout<<" s4="<<s4<<endl;


    pid=fmod(4.0*s1-2.0*s2-s3-s4,1.0)+1.0;
    cout<<"\n sum="<<pid<<endl;


    ihex(pid,0);
    ihex(pid,1);


return 0;
}
// --------------------------



void ihex(double x, int bq)
//
// print the first nhx hex or binary digits of x
//
{
char hx[]="0123456789abcdefghijklmnopqrstuvwxyz";
//    {'0','1','2','3','4','5','6','7','8','9','a','b','c','d','e','f'};
char *bn[]=
    {"0000","0001","0010","0011","0100","0101","0110","0111",
     "1000","1001","1010","1011","1100","1101","1110","1111"};
const int nhx=10; //  number of hex digits to print 
double y;


    if(bq==0)
        cout<<" Pi at hex position "<<ic<<"  = ";
    else
        cout<<" Pi at bin position "<<4*ic<<"  = ";

    y=fabs(x);

    for(int i=1; i<=nhx; ++i)
    {
        y=16.0*(y-floor(y+0.0));

        if(bq==0)
            cout<<hx[(int)floor(y+0.0)];
        else
            cout<<bn[(int)floor(y+0.0)]<<",";
    }

    cout<<endl;

return;
}
// --------------------------


double series(double m, double &ic)
//
// evaluate the series sum_k 16^(ic-k)/(8*k+m)
// using the modular exponentiation technique
//
{
double k;
const double eps=1e-17;
double ak,p,s,t;

    s=0.0;
    itp=ntp-1;

    for(k=0; k<ic; ++k)
    {
        ak=8*k+m;
        p=ic-k;
        t=expm(p,ak);
        s=fmod(s+t/ak,1.0);
    }


    for(k=ic; k<=ic+100; ++k)
    {
        ak=8*k+m;
        t=pow(16.0,(ic-k))/ak;
        if(t<eps)
	    break;
        s=fmod(s+t,1.0);
    }


return s;
}
// --------------------------


double expm(double &p, double ak)
//
// evaluates 16^p mod ak 
// this routine uses left-to-right binary exponentiation scheme
// it is valid for ak <= 2^24
//
{
double p1,r;
double pt;


/*
    for(itp=0; itp<ntp; ++itp)  // find the greatest power less than or equal to p
        if(tp[itp]>p)
            break;
*/

    if(tp[itp]>p)               // slightly faster: update only 
        itp--;


    pt=tp[itp];

    p1=p;
    r=1.0;


x110:                         // perform binary exponentiation algorithm

    if(p1>=pt)
    {
        r=fmod(16.0*r,ak);
        p1-=pt;
    }

    pt*=0.5;

    if(pt>=1.0)
    {
        r=fmod(r*r,ak);
        goto x110;
    }

    r=fmod(r,ak);


return r;
}
// --------------------------