=================================================================== RCS file: /home/cvs/OpenXM_contrib/gmp/mpfr/Attic/agm.c,v retrieving revision 1.1.1.1 retrieving revision 1.1.1.2 diff -u -p -r1.1.1.1 -r1.1.1.2 --- OpenXM_contrib/gmp/mpfr/Attic/agm.c 2000/09/09 14:12:19 1.1.1.1 +++ OpenXM_contrib/gmp/mpfr/Attic/agm.c 2003/08/25 16:06:06 1.1.1.2 @@ -1,141 +1,144 @@ /* mpfr_agm -- arithmetic-geometric mean of two floating-point numbers -Copyright (C) 1999 PolKA project, Inria Lorraine and Loria +Copyright 1999, 2000, 2001, 2002 Free Software Foundation. This file is part of the MPFR Library. The MPFR Library is free software; you can redistribute it and/or modify -it under the terms of the GNU Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. The MPFR Library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. -You should have received a copy of the GNU Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include -#include #include "gmp.h" #include "gmp-impl.h" #include "mpfr.h" +#include "mpfr-impl.h" - -/*Memory gestion */ -#define MON_INIT(xp, x, p, s) xp = (mp_ptr) TMP_ALLOC(s*BYTES_PER_MP_LIMB); x -> _mp_prec = p; x -> _mp_d = xp; x -> _mp_size = s; x -> _mp_exp = 0; - - - - void -#ifdef __STDC__ -mpfr_agm(mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, unsigned char rnd_mode) -#else -mpfr_agm(r, a, b, rnd_mode) - mpfr_ptr r; - mpfr_srcptr a; - mpfr_srcptr b; - unsigned char rnd_mode; -#endif +mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode) { - int s, p, q, 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; mpfr_t u, v, tmp, tmpu, tmpv, a, b; - TMP_DECL(marker1); - TMP_DECL(marker2); + TMP_DECL(marker); - /* If a or b is NaN, the result is NaN */ - if (FLAG_NAN(op1) || FLAG_NAN(op2)) - { SET_NAN(r); return; } + if (MPFR_IS_NAN(op1) || MPFR_IS_NAN(op2)) + { + MPFR_SET_NAN(r); + __mpfr_flags |= MPFR_FLAGS_NAN; + 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); + __mpfr_flags |= MPFR_FLAGS_NAN; + return; + } - /* If a or b is negative, the result is NaN */ - if ((SIGN(op1)<0)||(SIGN(op2)<0)) - { 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_CLEAR_INF(r); /* If a or b is 0, the result is 0 */ - if ((SIGN(op1)==0)||(SIGN(op2)==0)) - { SET_ZERO(r); - return; + if ((MPFR_NOTZERO(op1) && MPFR_NOTZERO(op2)) == 0) + { + MPFR_SET_ZERO(r); + return; } /* precision of the following calculus */ - q = PREC(r); + q = MPFR_PREC(r); p = q + 15; - /* Initialisations */ go_on=1; - TMP_MARK(marker1); + TMP_MARK(marker); s=(p-1)/BITS_PER_MP_LIMB+1; - MON_INIT(ap, a, p, s); - MON_INIT(bp, b, p, s); - TMP_MARK(marker2); - MON_INIT(up, u, p, s); - MON_INIT(vp, v, p, s); - MON_INIT(tmpup, tmpu, p, s); - MON_INIT(tmpvp, tmpv, p, s); - MON_INIT(tmpp, tmp, p, s); + MPFR_INIT(ap, a, p, s); + MPFR_INIT(bp, b, p, s); + MPFR_INIT(up, u, p, s); + MPFR_INIT(vp, v, p, s); + MPFR_INIT(tmpup, tmpu, p, s); + MPFR_INIT(tmpvp, tmpv, p, s); + MPFR_INIT(tmpp, tmp, p, s); - /* 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); - } + /* b and a are the 2 operands but we want b >= a */ + 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); + TMP_FREE(marker); + return; + } + else + { + mpfr_set (b,op2,GMP_RNDN); + mpfr_set (a,op1,GMP_RNDN); + } - vo=mpfr_get_d(b); - uo=mpfr_get_d(a); + vo = mpfr_get_d1 (b); + uo = mpfr_get_d1 (a); mpfr_set(u,a,GMP_RNDN); mpfr_set(v,b,GMP_RNDN); - /* Main loop */ while (go_on) { - int err, eq, can_round; + int err, can_round; + mp_prec_t eq; - eq=0; - - err=ceil((3.0/2.0*log((double)p)/log(2.0)+1.0)*exp(-(double)p*log(2.0))+3.0*exp(-2.0*(double)p*uo*log(2.0)/(vo-uo))); + err=1 + (int) ((3.0/2.0*(double)_mpfr_ceil_log2((double)p)+1.0)*_mpfr_ceil_exp2(-(double)p) + +3.0*_mpfr_ceil_exp2(-2.0*(double)p*uo/(vo-uo))); if(p-err-3<=q) { p=q+err+4; - err=ceil((3.0/2.0*log((double)p)/log(2.0)+1.0)*exp(-(double)p*log(2.0))+3.0*exp(-2.0*(double)p*uo*log(2.0)/(vo-uo))); + err= 1 + + (int) ((3.0/2.0*_mpfr_ceil_log2((double)p)+1.0)*_mpfr_ceil_exp2(-(double)p) + +3.0*_mpfr_ceil_exp2(-2.0*(double)p*uo/(vo-uo))); } /* Calculus of un and vn */ - while (eq<=p-2) { - mpfr_mul(tmp,u,v,GMP_RNDN); - mpfr_sqrt(tmpu,tmp,GMP_RNDN); - mpfr_add(tmp,u,v,GMP_RNDN); - 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); - else - eq=mpfr_cmp2(u,v); - } + do + { + mpfr_mul(tmp, u, v, GMP_RNDN); + mpfr_sqrt (tmpu, tmp, GMP_RNDN); + mpfr_add(tmp, u, v, GMP_RNDN); + mpfr_div_2ui(tmpv, tmp, 1, GMP_RNDN); + mpfr_set(u, tmpu, GMP_RNDN); + mpfr_set(v, tmpv, GMP_RNDN); + } + while (mpfr_cmp2(u, v, &eq) != 0 && eq <= p - 2); - /* printf("avant can_round %i bits faux\n v : ",err+3); - mpfr_print_raw(v); printf("\n u : "); - mpfr_print_raw(u);printf("\n");*/ - - /* Roundability of the result */ can_round=mpfr_can_round(v,p-err-3,GMP_RNDN,rnd_mode,q); @@ -145,14 +148,12 @@ mpfr_agm(r, a, b, rnd_mode) else { go_on=1; p+=5; - TMP_FREE(marker2); - TMP_MARK(marker2); s=(p-1)/BITS_PER_MP_LIMB+1; - MON_INIT(up, u, p, s); - MON_INIT(vp, v, p, s); - MON_INIT(tmpup, tmpu, p, s); - MON_INIT(tmpvp, tmpv, p, s); - MON_INIT(tmpp, tmp, p, s); + MPFR_INIT(up, u, p, s); + MPFR_INIT(vp, v, p, s); + MPFR_INIT(tmpup, tmpu, p, s); + MPFR_INIT(tmpvp, tmpv, p, s); + MPFR_INIT(tmpp, tmp, p, s); mpfr_set(u,a,GMP_RNDN); mpfr_set(v,b,GMP_RNDN); } @@ -161,12 +162,10 @@ mpfr_agm(r, a, b, rnd_mode) /* Setting of the result */ - mpfr_set(r,v,rnd_mode); + mpfr_set(r,v,rnd_mode); - /* Let's clean */ - TMP_FREE(marker1); - + TMP_FREE(marker); + return ; } -