ex-2/fancier: optimise a little and print execution time
This commit is contained in:
parent
ab05d23c7e
commit
cf0a3216b3
@ -3,6 +3,7 @@
|
||||
#include <stdarg.h>
|
||||
#include <math.h>
|
||||
#include <gmp.h>
|
||||
#include <time.h>
|
||||
|
||||
/* 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,
|
||||
|
||||
Loading…
Reference in New Issue
Block a user