diff options
| -rw-r--r-- | agm.c | 52 | ||||
| -rw-r--r-- | set_q.c | 22 |
2 files changed, 48 insertions, 26 deletions
@@ -81,7 +81,7 @@ mpfr_agm (r, op2, op1, rnd_mode) mp_rnd_t rnd_mode; #endif { - int s, go_on; + int s, go_on, compare; mp_prec_t p, q; double uo, vo; mp_limb_t *up, *vp, *tmpp, *tmpup, *tmpvp, *ap, *bp; @@ -91,25 +91,36 @@ mpfr_agm (r, op2, op1, rnd_mode) TMP_DECL(marker2); /* If a or b is NaN, the result is NaN */ - if (MPFR_IS_NAN(op1) || MPFR_IS_NAN(op2)) - { MPFR_SET_NAN(r); return; } + if (MPFR_IS_NAN(op1) || MPFR_IS_NAN(op2)) + { + MPFR_SET_NAN(r); + return; + } /* If a or b is negative (including -Infinity), the result is NaN */ if ((MPFR_SIGN(op1) < 0) || (MPFR_SIGN(op2) < 0)) - { MPFR_SET_NAN(r); return; } + { + MPFR_SET_NAN(r); + return; + } MPFR_CLEAR_NAN(r); /* If a or b is +Infinity, the result is +Infinity */ if (MPFR_IS_INF(op1) || MPFR_IS_INF(op2)) - { MPFR_SET_INF(r); MPFR_SET_SAME_SIGN(r, op1); return; } + { + MPFR_SET_INF(r); + MPFR_SET_SAME_SIGN(r, op1); + return; + } MPFR_CLEAR_INF(r); /* If a or b is 0, the result is 0 */ if ((MPFR_NOTZERO(op1) && MPFR_NOTZERO(op2)) == 0) - { MPFR_SET_ZERO(r); - return; + { + MPFR_SET_ZERO(r); + return; } /* precision of the following calculus */ @@ -133,12 +144,21 @@ mpfr_agm (r, op2, op1, rnd_mode) /* b and a will be the 2 operands but I want b>= a */ - if (mpfr_cmp(op1,op2) > 0) { - mpfr_set(b,op1,GMP_RNDN); mpfr_set(a,op2,GMP_RNDN); - } - else { - mpfr_set(b,op2,GMP_RNDN); mpfr_set(a,op1,GMP_RNDN); - } + if ((compare = mpfr_cmp (op1,op2)) > 0) + { + mpfr_set (b,op1,GMP_RNDN); + mpfr_set (a,op2,GMP_RNDN); + } + else if (compare == 0) + { + mpfr_set (r, op1, rnd_mode); + return; + } + else + { + mpfr_set (b,op2,GMP_RNDN); + mpfr_set (a,op1,GMP_RNDN); + } vo=mpfr_get_d(b); uo=mpfr_get_d(a); @@ -171,10 +191,10 @@ mpfr_agm (r, op2, op1, rnd_mode) mpfr_div_2exp(tmpv,tmp,1,GMP_RNDN); mpfr_set(u,tmpu,GMP_RNDN); mpfr_set(v,tmpv,GMP_RNDN); - if (mpfr_cmp(v,u)>=0) - eq=mpfr_cmp2(v,u); + if (mpfr_cmp(v,u) >= 0) + eq = mpfr_cmp2(v,u); else - eq=mpfr_cmp2(u,v); + eq = mpfr_cmp2(u,v); } /* printf("avant can_round %i bits faux\n v : ",err+3); @@ -42,22 +42,24 @@ mpfr_set_q (f, q, rnd) mpfr_t n,d; MPFR_CLEAR_FLAGS(f); - num = mpq_numref(q); - sign = mpz_cmp_ui(num, 0); - if (sign==0) { - MPFR_SET_ZERO(f); - return; - } + num = mpq_numref (q); + sign = mpz_cmp_ui (num, 0); + if (sign == 0) + { + MPFR_SET_ZERO(f); + return; + } den = mpq_denref(q); prec = MPFR_PREC(f); - mpfr_init2(n, mpz_sizeinbase(num, 2)); - mpfr_set_z(n, num, GMP_RNDZ); /* result is exact */ + mpfr_init2 (n, mpz_sizeinbase(num, 2)); + mpfr_set_z (n, num, GMP_RNDZ); /* result is exact */ mpfr_init2(d, mpz_sizeinbase(den, 2)); mpfr_set_z(d, den, GMP_RNDZ); /* result is exact */ - MPFR_PREC(f) = prec; + MPFR_PREC(f) = prec; mpfr_div(f, n, d, rnd); - mpfr_clear(n); mpfr_clear(d); + mpfr_clear(n); + mpfr_clear(d); } |