All pastes #244866 Raw Edit

Something

public cpp v1 · immutable
#244866 ·published 2006-11-10 08:22 UTC
rendered paste body
/******************************************************************************\*                                                                              **   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.14159265358979using 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;}