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

ex-2/fancier: fix crash in mpq_log_ui on power of 2

Browse Source
This commit is contained in:
Michele Guerini Rocco 2020-05-25 23:30:03 +02:00
parent 5c7f46ce2f
commit ce91072040
1 changed files with 24 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
ex-2/fancier.c
View File

@ -85,6 +85,10 @@ char* mpq_dec_str(mpq_t z, size_t mantissa_size) {
}
/* printf() like function with special conversion %q, which
* prints the deciamal expansion of mpq_t rational numbers.
* The second argument is the number of digits shown.
*/
void mpq_fprintf(FILE *fp, size_t digits, char* fmt, ...) {
va_list args;
va_start(args, fmt);
@ -243,6 +247,26 @@ void mpq_log_ui(mpq_t rop, size_t x, size_t digits) {
fprintf(stderr, "[mpq_log_ui] binary digits: %ld\n", m);
fprintf(stderr, "[mpq_log_ui] a=%g\n", mpq_get_d(a));
// calculate log(2)*(m-1)
mpq_t log2; mpq_init(log2);
mpq_log2(log2, digits);
mpz_mul_ui(mpq_numref(log2), mpq_numref(log2), m-1);
mpq_canonicalize(log2);
/* When x is a power of two a=1 and y=0.
* The series for log(1+y/1-y) is not defined for y=0
* so this case must be handled separately.
*/
if (mpq_cmp_ui(a, 1, 1) == 0) {
mpq_set(rop, log2);
// free memory
mpq_clears(a, log2, NULL);
mpz_clears(xz, power, NULL);
return;
}
/* calculate y = (a-1)/(a+1) */
mpq_t y, temp; mpq_inits(y, temp, NULL);
@ -366,13 +390,6 @@ void mpq_log_ui(mpq_t rop, size_t x, size_t digits) {
mpq_fprintf(stderr, d+2, "[mpq_log_ui] lower bound: %q...\n", low);
mpq_fprintf(stderr, d, "[mpq_log_ui] series: %q...\n", sum);
// calculate log(2)*(m-1)
mpq_t log2; mpq_init(log2);
mpq_log2(log2, digits);
mpz_mul_ui(mpq_numref(log2), mpq_numref(log2), m-1);
mpq_canonicalize(log2);
// series += log(2)*(m-1)
mpq_add(sum, sum, log2);
mpq_set(rop, sum);