Giu/analistica
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

ex-2: implement correct series estimation in fancier.c

Browse Source
This commit is contained in:
Michele Guerini Rocco 2020-05-24 12:49:47 +00:00
parent 49e9fbd5b7
commit c8f34f1822
1 changed files with 57 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

63
ex-2/fancier.c
View File

@ -146,20 +146,42 @@ void mpq_log2(mpq_t rop, size_t digits) {
} }
fprintf(stderr, "[mpq_log2] series terms: %ld\n", n); fprintf(stderr, "[mpq_log2] series terms: %ld\n", n);
/* calculate series */ /* Compute the nth partial sum */
for (size_t i=1; i<n+1; i++) { for (size_t i=1; i<=n; i++) {
mpq_set_ui(term, 1, 1); mpq_set_ui(term, 1, 1);
mpz_ui_pow_ui(mpq_denref(term), 2, i); mpz_ui_pow_ui(mpq_denref(term), 2, i);
mpz_mul_ui(mpq_denref(term), mpq_denref(term), i); mpz_mul_ui(mpq_denref(term), mpq_denref(term), i);
mpq_add(sum, sum, term); mpq_add(sum, sum, term);
} }
mpq_fprintf(stderr, digits > 50 ? 50 : digits,
"[mpq_log2] log(2)=%q...\n", &sum);
/* Compute upper/lower bounds */
mpq_t upp, low;
mpq_inits(upp, low, NULL);
// upp = sum + term
mpq_add(upp, sum, term);
// low = sum + term * n/(n+2)
mpz_mul_ui(mpq_numref(low), mpq_numref(term), n);
mpz_mul_ui(mpq_denref(low), mpq_denref(term), n+2);
mpq_canonicalize(low);
mpq_add(low, sum, low);
/* Average the two bounds */
mpq_add(sum, upp, low);
mpz_mul_ui(mpq_denref(sum), mpq_denref(sum), 2);
mpq_canonicalize(sum);
size_t d = digits > 50 ? 50 : digits;
mpq_fprintf(stderr, d+2, "[mpq_log2] upper bound: %q...\n", &upp);
mpq_fprintf(stderr, d+2, "[mpq_log2] lower bound: %q...\n", &low);
mpq_fprintf(stderr, d, "[mpq_log2] log(2)=%q...\n", &sum);
// return sum
mpq_set(rop, sum); mpq_set(rop, sum);
// free memory // free memory
mpq_clears(sum, term, NULL); mpq_clears(sum, term, upp, low, NULL);
mpz_clears(temp, prec, width, NULL); mpz_clears(temp, prec, width, NULL);
} }
@ -278,7 +300,7 @@ void mpq_log_ui(mpq_t rop, size_t x, size_t digits) {
fprintf(stderr, "[mpq_log_ui] decimal digits: %ld\n", digits); fprintf(stderr, "[mpq_log_ui] decimal digits: %ld\n", digits);
/* calculate the series of /* Compute the nth partial sum of
* log(A) = log((1+y)/(1-y)) = 2 Σ_{k=0} y^(2k+1)/(2k+1) * log(A) = log((1+y)/(1-y)) = 2 Σ_{k=0} y^(2k+1)/(2k+1)
*/ */
mpq_t sum, term; mpq_inits(sum, term, NULL); mpq_t sum, term; mpq_inits(sum, term, NULL);
@ -299,6 +321,35 @@ void mpq_log_ui(mpq_t rop, size_t x, size_t digits) {
mpq_add(sum, sum, term); mpq_add(sum, sum, term);
} }
/* Compute the upper and lower bounds */
mpq_t upp, low;
mpq_inits(upp, low, NULL);
// upp = sum + term * y2/(1-y2)
mpq_mul(upp, term, y2);
mpq_set_ui(temp, 1, 1);
mpq_sub(temp, temp, y2);
mpq_div(upp, upp, temp);
mpq_add(upp, sum, upp);
// low = sum + term / ((2*n + 3)/(2*n + 1)*y^-2 - 1)
mpq_set_ui(temp, 2*n + 3, 2*n + 1);
mpq_div(temp, temp, y2);
mpz_sub(mpq_numref(temp), mpq_numref(temp), mpq_denref(temp));
mpq_canonicalize(temp);
mpq_div(low, term, temp);
mpq_add(low, sum, low);
/* Average the two bounds */
mpq_add(sum, low, upp);
mpz_div_ui(mpq_numref(sum), mpq_numref(sum), 2);
mpq_canonicalize(sum);
size_t d = digits > 50 ? 50 : digits;
mpq_fprintf(stderr, d+2, "[mpq_log_ui] upper bound: %q...\n", upp);
mpq_fprintf(stderr, d+2, "[mpq_log_ui] lower bound: %q...\n", low);
mpq_fprintf(stderr, d, "[mpq_log_ui] series: %q...\n", sum);
// multiply by 2 // multiply by 2
mpz_mul_ui(mpq_numref(sum), mpq_numref(sum), 2); mpz_mul_ui(mpq_numref(sum), mpq_numref(sum), 2);
mpq_canonicalize(sum); mpq_canonicalize(sum);