1 /* mpfr_agm -- arithmetic-geometric mean of two floating-point numbers
3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
4 Contributed by the Arenaire and Cacao projects, INRIA.
6 This file is part of the GNU MPFR Library.
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 2.1 of the License, or (at your
11 option) any later version.
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LIB. If not, write to
20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
21 MA 02110-1301, USA. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 /* agm(x,y) is between x and y, so we don't need to save exponent range */
28 mpfr_agm (mpfr_ptr r
, mpfr_srcptr op2
, mpfr_srcptr op1
, mp_rnd_t rnd_mode
)
33 mp_limb_t
*up
, *vp
, *tmpp
;
35 unsigned long n
; /* number of iterations */
36 unsigned long err
= 0;
38 MPFR_TMP_DECL(marker
);
40 MPFR_LOG_FUNC (("op2[%#R]=%R op1[%#R]=%R rnd=%d", op2
,op2
,op1
,op1
,rnd_mode
),
41 ("r[%#R]=%R inexact=%d", r
, r
, inexact
));
43 /* Deal with special values */
44 if (MPFR_ARE_SINGULAR (op1
, op2
))
46 /* If a or b is NaN, the result is NaN */
47 if (MPFR_IS_NAN(op1
) || MPFR_IS_NAN(op2
))
52 /* now one of a or b is Inf or 0 */
53 /* If a and b is +Inf, the result is +Inf.
54 Otherwise if a or b is -Inf or 0, the result is NaN */
55 else if (MPFR_IS_INF(op1
) || MPFR_IS_INF(op2
))
57 if (MPFR_IS_STRICTPOS(op1
) && MPFR_IS_STRICTPOS(op2
))
60 MPFR_SET_SAME_SIGN(r
, op1
);
61 MPFR_RET(0); /* exact */
69 else /* a and b are neither NaN nor Inf, and one is zero */
70 { /* If a or b is 0, the result is +0 since a sqrt is positive */
71 MPFR_ASSERTD (MPFR_IS_ZERO (op1
) || MPFR_IS_ZERO (op2
));
74 MPFR_RET (0); /* exact */
79 /* If a or b is negative (excluding -Infinity), the result is NaN */
80 if (MPFR_UNLIKELY(MPFR_IS_NEG(op1
) || MPFR_IS_NEG(op2
)))
86 /* Precision of the following calculus */
88 p
= q
+ MPFR_INT_CEIL_LOG2(q
) + 15;
89 MPFR_ASSERTD (p
>= 7); /* see algorithms.tex */
90 s
= (p
- 1) / BITS_PER_MP_LIMB
+ 1;
92 /* b (op2) and a (op1) are the 2 operands but we want b >= a */
93 compare
= mpfr_cmp (op1
, op2
);
94 if (MPFR_UNLIKELY( compare
== 0 ))
96 mpfr_set (r
, op1
, rnd_mode
);
97 MPFR_RET (0); /* exact */
105 /* Now b(=op2) >= a (=op1) */
107 MPFR_TMP_MARK(marker
);
110 MPFR_ZIV_INIT (loop
, p
);
115 /* Init temporary vars */
116 MPFR_TMP_INIT (up
, u
, p
, s
);
117 MPFR_TMP_INIT (vp
, v
, p
, s
);
118 MPFR_TMP_INIT (tmpp
, tmp
, p
, s
);
120 /* Calculus of un and vn */
121 mpfr_mul (u
, op1
, op2
, GMP_RNDN
); /* Faster since PREC(op) < PREC(u) */
122 mpfr_sqrt (u
, u
, GMP_RNDN
);
123 mpfr_add (v
, op1
, op2
, GMP_RNDN
); /* add with !=prec is still good*/
124 mpfr_div_2ui (v
, v
, 1, GMP_RNDN
);
126 while (mpfr_cmp2 (u
, v
, &eq
) != 0 && eq
<= p
- 2)
128 mpfr_add (tmp
, u
, v
, GMP_RNDN
);
129 mpfr_div_2ui (tmp
, tmp
, 1, GMP_RNDN
);
130 /* See proof in algorithms.tex */
135 mpfr_init2 (w
, (p
+ 1) / 2);
136 mpfr_sub (w
, v
, u
, GMP_RNDN
); /* e = V(k-1)-U(k-1) */
137 mpfr_sqr (w
, w
, GMP_RNDN
); /* e = e^2 */
138 mpfr_div_2ui (w
, w
, 4, GMP_RNDN
); /* e*= (1/2)^2*1/4 */
139 mpfr_div (w
, w
, tmp
, GMP_RNDN
); /* 1/4*e^2/U(k) */
140 mpfr_sub (v
, tmp
, w
, GMP_RNDN
);
141 err
= MPFR_GET_EXP (tmp
) - MPFR_GET_EXP (v
); /* 0 or 1 */
145 mpfr_mul (u
, u
, v
, GMP_RNDN
);
146 mpfr_sqrt (u
, u
, GMP_RNDN
);
150 /* the error on v is bounded by (18n+51) ulps, or twice if there
151 was an exponent loss in the final subtraction */
152 err
+= MPFR_INT_CEIL_LOG2(18 * n
+ 51); /* 18n+51 should not overflow
153 since n is about log(p) */
154 /* we should have n+2 <= 2^(p/4) [see algorithms.tex] */
155 if (MPFR_LIKELY (MPFR_INT_CEIL_LOG2(n
+ 2) <= p
/ 4 &&
156 MPFR_CAN_ROUND (v
, p
- err
, q
, rnd_mode
)))
157 break; /* Stop the loop */
160 MPFR_ZIV_NEXT (loop
, p
);
161 s
= (p
- 1) / BITS_PER_MP_LIMB
+ 1;
163 MPFR_ZIV_FREE (loop
);
165 /* Setting of the result */
166 inexact
= mpfr_set (r
, v
, rnd_mode
);
169 MPFR_TMP_FREE(marker
);
171 return inexact
; /* agm(u,v) can be exact for u, v rational only for u=v.
172 Proof (due to Nicolas Brisebarre): it suffices to consider
173 u=1 and v<1. Then 1/AGM(1,v) = 2F1(1/2,1/2,1;1-v^2),
174 and a theorem due to G.V. Chudnovsky states that for x a
175 non-zero algebraic number with |x|<1, then
176 2F1(1/2,1/2,1;x) and 2F1(-1/2,1/2,1;x) are algebraically
177 independent over Q. */