/******************************************************************************\
*                                                                              *
*   Advanced Problem 2 - The Most Amazing Mathematical Formula                 *
*                                                                              *
*   flamingspinach                                                             *
*                                                                              *
*   Started:  2006-10-06                                                       *
*   Finished: 2006-10-06                                                       *
*                                                                              *
*   NOTE: Made use of Euler's Formula to eliminate the need for complex<>.     *
*         coprime() makes use of the Euclidean algorithm. Reformulated the     *
*         recursive function according to my understanding of it rather than   *
*         by the given function.                                               *
*                                                                              *
\******************************************************************************/

#include <iostream>
#include <iomanip>
#include <ctime>
#include <cmath>
#define _PI 3.14159265358979
using namespace std;

bool coprime (int a, int b) {
    if (!b) return !(a-1);
    return coprime(b, a % b);
}

int p_recursive (int n, int k) {
    if (n < k) return 0;
    if (n == k) return 1;
    return p_recursive(n, k+1) + p_recursive(n-k, k);
}    

double p_Rademacher (int n) {
    double a, b, c, d, h, k, r, x;
    for (k = 1; k*k <= n; ++k) {
        a = sqrt(k)*pow(n-1/24.,-3/2.);
        b = sqrt(2/3.*_PI*_PI/k/k*(n-1/24.));
        b = b * cosh(b) - sinh(b);
        for (c = 0, h = 1; h <= k; ++h)
            if (coprime(h,k)) {
                for (d = 0, r = 1; r <= k-1; ++r)
                    d += r/k*(h*r/k - (int)(h*r/k) - 1/2.);
                c += cos(_PI*(d - 2.*h*n/k));
            }
        x += a * b * c;
    }
    return floor(x/_PI/sqrt(8)+0.5);
}

double p(int n, bool useRademacher) { return useRademacher ? p_Rademacher(n) : p_recursive(n,1); }

int main () {
    double timer1, timer2, x;
    int i;
    cout << setprecision(10);
    
    cout << "p(n), by recursion: " << endl;
    timer1 = -clock();
    for (i = 0; timer1 + clock() < 5 * CLK_TCK; ++i) {                          // replace constraint with "i < somenumber" if required
        timer2 = -clock(); x = p(i, false); timer2 += clock();                  /* remove comment marker for better timing --> */
        cout << "p(" << setw(3) << i << ") = " << setw(25) << x
             << "     Time: " << setw(10) << timer2/CLK_TCK << "s" << endl;     /* <-- remove comment marker for better timing */
    }
    timer1 += clock();
    cout << "p(n), by recursion." << endl
         << "Time elapsed total: " << timer1/CLK_TCK << "s" << endl
         << "Total n reached: " << i - 1 << endl
         << endl;
    
    system("pause");
    cout << endl;
    
    cout << "p(n), by Rademacher series: " << endl;
    timer1 = -clock();
    for (i = 0; timer1 + clock() < 5 * CLK_TCK; ++i) {                          // replace constraint with "i < somenumber" if required
        timer2 = -clock(); x = p(i, true); timer2 += clock();                   /* remove comment marker for better timing --> */
        cout << "p(" << setw(3) << i << ") = " << setw(25) << x
             << "     Time: " << setw(10) << timer2/CLK_TCK << "s" << endl;     /* <-- remove comment marker for better timing */
    }
    timer1 += clock();
    cout << "p(n), by Rademacher series." << endl
         << "Time elapsed total: " << timer1/CLK_TCK << "s" << endl
         << "Total n reached: " << i - 1 << endl
         << endl;
    
    system("pause");
    return 0;
}