BorweinPlouffePi_nth.cpp

/***************************************************************
 * Computing nth hex digit of Pi with BBP formula.
 * Bailey-Borwein-Plouffe_formula for Pi
 * (original https://gist.github.com/komasaru/9e418eb666ab649ef589)
 * removed class declarations, uses only functions
 * Juhani Kaukoranta 31.7.2019
 **************************************************************/
#include <math.h>
#include <iostream>
#include <stdio.h>
#include <stdlib.h>

using namespace std;

// Declaration
int     n;                          // nth Digit to compute
double  pi;                         // Pi
char    pi_hex[14];                 // Pi(Hex)
clock_t t0, t1;                     // Time


/*
 * Convert Pi to Hex-strings
 */
void convHex(double x, char chx[]) {
    double y;
    int i;
    const char hx[] = "0123456789ABCDEF";

    y = fabs(x);
    for (i = 0; i < 16; i++) {
        y = 16.0 * (y - floor(y));
        chx[i] = hx[(int)(y)];
    }
}

/*
 * Compute Modular Exponentiation
 */
long compModExp(int b, int e, int m)
{
    long ans;

    if (e == 0) return 1;

    ans = compModExp(b, e / 2, m);
    ans = (ans * ans) % m;
    if ((e % 2) == 1) ans = (ans * b) % m;

    return ans;
}


/*
 * Compute S
 */
double S(int j)
{
    double s = 0.0;        // Summation of Total, Left
    double t;              // Each term of right summation
    int    r;              // Denominator
    int    k;              // Loop index
    double EPS = 1.0e-17;  // Loop-exit accuration of the right summation

    // Left Sum (0 ... d)
    for (k = 0; k <= n; k++) {
        r = 8 * k + j;
        t = (double)compModExp(16, n - k, r);
        t /= r;
        s += t - (int)t;
        s -= (int)s;
    }

    // Right sum (n + 1 ...)
    while (1) {
        r = 8 * k + j;
        t = pow(16.0, (double)(n - k));
        t /= (double)r;
        if (t < EPS) break;
        s += t;
        s -= (int)s;
        k ++;
    }

    return s;
}

/*
 * Compute PI
 */
void compPi()
{
    // Time (start)
    t0 = clock();

    // Compute Pi
    pi = 4.0 * S(1) - 2.0 * S(4) - S(5) - S(6);
    pi = pi - (int)pi + 1.0;
    convHex(pi, pi_hex);
    cout << "Hex digit: " << (char)pi_hex[0] << endl;
    printf("FRACTION  : %.15f\n",   pi);
    printf("10 HEX DIGITS: %10.10s\n",pi_hex);


    // Time (end of computation)
    t1 = clock();
    cout << "( TIME: " << (double)(t1 - t0) / CLOCKS_PER_SEC
         << " seconds )" << endl;
}

// paaohjelma

int main()
{

	cout << "kuinka mones Piin desimaali?";
	cin >> n;
	cout << "**** PI Computation ( digit: " << n << " )" << endl;
	// Compute PI
	n = n - 1;
	compPi();

}