Annotation of OpenXM_contrib/gmp/mpfr/agm.c, Revision 1.1.1.2
1.1 maekawa 1: /* mpfr_agm -- arithmetic-geometric mean of two floating-point numbers
2:
1.1.1.2 ! ohara 3: Copyright 1999, 2000, 2001, 2002 Free Software Foundation.
1.1 maekawa 4:
5: This file is part of the MPFR Library.
6:
7: The MPFR Library is free software; you can redistribute it and/or modify
1.1.1.2 ! ohara 8: it under the terms of the GNU Lesser General Public License as published by
! 9: the Free Software Foundation; either version 2.1 of the License, or (at your
1.1 maekawa 10: option) any later version.
11:
12: The MPFR Library is distributed in the hope that it will be useful, but
13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.1.1.2 ! ohara 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 15: License for more details.
16:
1.1.1.2 ! ohara 17: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 18: along with the MPFR Library; see the file COPYING.LIB. If not, write to
19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20: MA 02111-1307, USA. */
21:
22: #include "gmp.h"
23: #include "gmp-impl.h"
24: #include "mpfr.h"
1.1.1.2 ! ohara 25: #include "mpfr-impl.h"
1.1 maekawa 26:
27: void
1.1.1.2 ! ohara 28: mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode)
1.1 maekawa 29: {
1.1.1.2 ! ohara 30: int s, go_on, compare;
! 31: mp_prec_t p, q;
1.1 maekawa 32: double uo, vo;
33: mp_limb_t *up, *vp, *tmpp, *tmpup, *tmpvp, *ap, *bp;
34: mpfr_t u, v, tmp, tmpu, tmpv, a, b;
1.1.1.2 ! ohara 35: TMP_DECL(marker);
1.1 maekawa 36:
37: /* If a or b is NaN, the result is NaN */
1.1.1.2 ! ohara 38: if (MPFR_IS_NAN(op1) || MPFR_IS_NAN(op2))
! 39: {
! 40: MPFR_SET_NAN(r);
! 41: __mpfr_flags |= MPFR_FLAGS_NAN;
! 42: return;
! 43: }
1.1 maekawa 44:
1.1.1.2 ! ohara 45: /* If a or b is negative (including -Infinity), the result is NaN */
! 46: if ((MPFR_SIGN(op1) < 0) || (MPFR_SIGN(op2) < 0))
! 47: {
! 48: MPFR_SET_NAN(r);
! 49: __mpfr_flags |= MPFR_FLAGS_NAN;
! 50: return;
! 51: }
1.1 maekawa 52:
1.1.1.2 ! ohara 53: MPFR_CLEAR_NAN(r);
1.1 maekawa 54:
1.1.1.2 ! ohara 55: /* If a or b is +Infinity, the result is +Infinity */
! 56: if (MPFR_IS_INF(op1) || MPFR_IS_INF(op2))
! 57: {
! 58: MPFR_SET_INF(r);
! 59: MPFR_SET_SAME_SIGN(r, op1);
! 60: return;
! 61: }
1.1 maekawa 62:
1.1.1.2 ! ohara 63: MPFR_CLEAR_INF(r);
1.1 maekawa 64:
65: /* If a or b is 0, the result is 0 */
1.1.1.2 ! ohara 66: if ((MPFR_NOTZERO(op1) && MPFR_NOTZERO(op2)) == 0)
! 67: {
! 68: MPFR_SET_ZERO(r);
! 69: return;
1.1 maekawa 70: }
71:
72: /* precision of the following calculus */
1.1.1.2 ! ohara 73: q = MPFR_PREC(r);
1.1 maekawa 74: p = q + 15;
75:
76: /* Initialisations */
77: go_on=1;
78:
1.1.1.2 ! ohara 79: TMP_MARK(marker);
1.1 maekawa 80: s=(p-1)/BITS_PER_MP_LIMB+1;
1.1.1.2 ! ohara 81: MPFR_INIT(ap, a, p, s);
! 82: MPFR_INIT(bp, b, p, s);
! 83: MPFR_INIT(up, u, p, s);
! 84: MPFR_INIT(vp, v, p, s);
! 85: MPFR_INIT(tmpup, tmpu, p, s);
! 86: MPFR_INIT(tmpvp, tmpv, p, s);
! 87: MPFR_INIT(tmpp, tmp, p, s);
1.1 maekawa 88:
89:
90:
1.1.1.2 ! ohara 91: /* b and a are the 2 operands but we want b >= a */
! 92: if ((compare = mpfr_cmp (op1,op2)) > 0)
! 93: {
! 94: mpfr_set (b,op1,GMP_RNDN);
! 95: mpfr_set (a,op2,GMP_RNDN);
! 96: }
! 97: else if (compare == 0)
! 98: {
! 99: mpfr_set (r, op1, rnd_mode);
! 100: TMP_FREE(marker);
! 101: return;
! 102: }
! 103: else
! 104: {
! 105: mpfr_set (b,op2,GMP_RNDN);
! 106: mpfr_set (a,op1,GMP_RNDN);
! 107: }
1.1 maekawa 108:
1.1.1.2 ! ohara 109: vo = mpfr_get_d1 (b);
! 110: uo = mpfr_get_d1 (a);
1.1 maekawa 111:
112: mpfr_set(u,a,GMP_RNDN);
113: mpfr_set(v,b,GMP_RNDN);
114:
115: /* Main loop */
116:
117: while (go_on) {
1.1.1.2 ! ohara 118: int err, can_round;
! 119: mp_prec_t eq;
1.1 maekawa 120:
1.1.1.2 ! ohara 121: err=1 + (int) ((3.0/2.0*(double)_mpfr_ceil_log2((double)p)+1.0)*_mpfr_ceil_exp2(-(double)p)
! 122: +3.0*_mpfr_ceil_exp2(-2.0*(double)p*uo/(vo-uo)));
1.1 maekawa 123: if(p-err-3<=q) {
124: p=q+err+4;
1.1.1.2 ! ohara 125: err= 1 +
! 126: (int) ((3.0/2.0*_mpfr_ceil_log2((double)p)+1.0)*_mpfr_ceil_exp2(-(double)p)
! 127: +3.0*_mpfr_ceil_exp2(-2.0*(double)p*uo/(vo-uo)));
1.1 maekawa 128: }
129:
130: /* Calculus of un and vn */
1.1.1.2 ! ohara 131: do
! 132: {
! 133: mpfr_mul(tmp, u, v, GMP_RNDN);
! 134: mpfr_sqrt (tmpu, tmp, GMP_RNDN);
! 135: mpfr_add(tmp, u, v, GMP_RNDN);
! 136: mpfr_div_2ui(tmpv, tmp, 1, GMP_RNDN);
! 137: mpfr_set(u, tmpu, GMP_RNDN);
! 138: mpfr_set(v, tmpv, GMP_RNDN);
! 139: }
! 140: while (mpfr_cmp2(u, v, &eq) != 0 && eq <= p - 2);
1.1 maekawa 141:
142: /* Roundability of the result */
143: can_round=mpfr_can_round(v,p-err-3,GMP_RNDN,rnd_mode,q);
144:
145: if (can_round)
146: go_on=0;
147:
148: else {
149: go_on=1;
150: p+=5;
151: s=(p-1)/BITS_PER_MP_LIMB+1;
1.1.1.2 ! ohara 152: MPFR_INIT(up, u, p, s);
! 153: MPFR_INIT(vp, v, p, s);
! 154: MPFR_INIT(tmpup, tmpu, p, s);
! 155: MPFR_INIT(tmpvp, tmpv, p, s);
! 156: MPFR_INIT(tmpp, tmp, p, s);
1.1 maekawa 157: mpfr_set(u,a,GMP_RNDN);
158: mpfr_set(v,b,GMP_RNDN);
159: }
160: }
161: /* End of while */
162:
163: /* Setting of the result */
164:
1.1.1.2 ! ohara 165: mpfr_set(r,v,rnd_mode);
1.1 maekawa 166:
167: /* Let's clean */
1.1.1.2 ! ohara 168: TMP_FREE(marker);
! 169:
1.1 maekawa 170: return ;
171: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to agm.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr