diff --git a/ex-2/fancier.c b/ex-2/fancier.c index a08e315..c2bc256 100644 --- a/ex-2/fancier.c +++ b/ex-2/fancier.c @@ -3,6 +3,7 @@ #include #include #include +#include /* The Euler-Mascheroni constant is here computed to * arbitrary precision using GMP rational numbers @@ -12,8 +13,8 @@ * * with: * - * A(N) = Σ_(k = 1)^(k_max) (N^k/k!) * H(k) - * B(N) = Σ_(k = 0)^(k_max) (N^k/k!) + * A(N) = Σ_(k = 1)^(k_max) (N^k/k!)² * H(k) + * B(N) = Σ_(k = 0)^(k_max) (N^k/k!)² * C(N) = 1/4N Σ_(k = 0)^(2N) ((2k)!)^3/((k!)^4*(16N))^2k * H(k) = Σ_(j = 1)^(k) (1/k), the k-th harmonic number * @@ -332,9 +333,13 @@ void a_series(mpq_t rop, size_t N) { mpz_inits(fact, power, stop, NULL); mpz_set_ui(stop, N); + mpz_set_ui(fact, 1); + mpz_set_ui(power, 1); for (size_t k = 0; k < 5*N; k++) { - mpz_fac_ui(fact, k); - mpz_pow_ui(power, stop, k); + if (k > 0) { + mpz_mul_ui(fact, fact, k); + mpz_mul(power, power, stop); + } mpq_set_z(term, power); mpq_set_z(factq, fact); @@ -365,9 +370,13 @@ void b_series(mpq_t rop, size_t N){ mpz_inits(fact, power, stop, NULL); mpz_set_ui(stop, N); + mpz_set_ui(fact, 1); + mpz_set_ui(power, 1); for (size_t k = 0; k < 5*N; k++) { - mpz_fac_ui(fact, k); - mpz_pow_ui(power, stop, k); + if (k > 0) { + mpz_mul_ui(fact, fact, k); + mpz_mul(power, power, stop); + } mpq_set_z(term, power); mpq_set_z(factq, fact); @@ -388,25 +397,29 @@ void b_series(mpq_t rop, size_t N){ /* C series */ void c_series(mpq_t rop, size_t N) { mpq_t sum, term, stop2, fact2q; - mpz_t fact1, fact2, power, stop1; + mpz_t fact1, pfact1, fact2, pfact2, power, stop1; mpq_inits(sum, term, stop2, fact2q, NULL); - mpz_inits(fact1, fact2, power, stop1, NULL); + mpz_inits(fact1, pfact1, fact2, pfact2, power, stop1, NULL); mpq_set_ui(stop2, 4*N, 1); + mpz_set_ui(fact1, 1); + mpz_set_ui(fact2, 1); + mpz_set_ui(stop1, 1); + for (size_t k = 0; k <= 2*N; k++) { - mpz_fac_ui(fact1, 2*k); - mpz_pow_ui(fact1, fact1, 3); + if (k > 0) { + mpz_mul_ui(fact1, fact1, 2*k * (2*k - 1)); + mpz_mul_ui(fact2, fact2, k); + mpz_mul_ui(stop1, stop1, 16*N); + mpz_mul_ui(stop1, stop1, 16*N); + } + mpz_pow_ui(pfact1, fact1, 3); + mpz_pow_ui(pfact2, fact2, 4); + mpz_mul(pfact2, pfact2, stop1); - mpz_fac_ui(fact2, k); - mpz_pow_ui(fact2, fact2, 4); - - mpz_set_ui(stop1, 16*N); - mpz_pow_ui(stop1, stop1, 2*k); - mpz_mul(fact2, fact2, stop1); - - mpq_set_z(term, fact1); - mpq_set_z(fact2q, fact2); + mpq_set_z(term, pfact1); + mpq_set_z(fact2q, pfact2); mpq_div(term, term, fact2q); mpq_add(sum, sum, term); @@ -417,7 +430,7 @@ void c_series(mpq_t rop, size_t N) { // free memory mpq_clears(sum, term, stop2, fact2q, NULL); - mpz_clears(fact1, fact2, power, stop1, NULL); + mpz_clears(fact1, pfact1, fact2, pfact2, power, stop1, NULL); } @@ -467,20 +480,48 @@ int main(int argc, char** argv) { // number of decimal places size_t D = atol(argv[1]); // number of terms - size_t N = floor(1.0/4 * log(10) * (double)D); + size_t N = 10 + floor(1.0/8 * log(10) * (double)D); + + fprintf(stderr, "[main] N: %ld\n", N); /* calculate series */ mpq_t A, B, C; mpq_inits(A, B, C, NULL); + + clock_t begin, end; + double time_a, time_b, time_c, time_log; + + begin = clock(); a_series(A, N); + end = clock(); + time_a = (double)(end - begin)/CLOCKS_PER_SEC; + fprintf(stderr, "[main] a series time: %gs\n", time_a); + + begin = clock(); b_series(B, N); + end = clock(); + time_b = (double)(end - begin)/CLOCKS_PER_SEC; + fprintf(stderr, "[main] b series time: %gs\n", time_b); + + begin = clock(); c_series(C, N); + end = clock(); + time_c = (double)(end - begin)/CLOCKS_PER_SEC; + fprintf(stderr, "[main] c series time: %gs\n", time_c); /* calculate γ... */ mpq_t gamma, corr, logn; mpq_inits(gamma, corr, logn, NULL); + begin = clock(); mpq_log_ui(logn, N, D); + end = clock(); + time_log = (double)(end - begin)/CLOCKS_PER_SEC; + fprintf(stderr, "[main] log(n) time: %gs\n", time_log); + + fprintf(stderr, "[main] total time: %gs\n", + time_a + time_b + time_c + time_log); + mpq_div(gamma, A, B); mpq_div(corr, C, B); mpq_div(corr, corr, B); @@ -493,15 +534,14 @@ int main(int argc, char** argv) { mpq_sub(diff, gamma, exact); /* print results */ - size_t d = D+1 > 50 ? 50 : D+1; + size_t d = D+1 > 100 ? 100 : D+1; - fprintf(stderr, "[main] N: %ld\n", N); mpq_fprintf(stderr, d, "[main] approx:\t%q...\n", gamma); mpq_fprintf(stderr, d, "[main] exact:\t%q...\n", exact); mpq_fprintf(stderr, d, "[main] diff:\t%q...\n", diff); mpq_fprintf(stderr, d, "[main] digits:\t%q\n", indicator); - mpq_fprintf(stdout, D+1, "%q\n", gamma); + mpq_fprintf(stdout, D, "%q\n", gamma); // free memory mpq_clears(A, B, C, gamma, corr, logn,