1 /* mpfr_agm -- arithmetic-geometric mean of two floating-point numbers
2
3 Copyright 1999-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
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 3 of the License, or (at your
11 option) any later version.
12
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.
17
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.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
25
26 /* agm(x,y) is between x and y, so we don't need to save exponent range */
27 int
28 mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mpfr_rnd_t rnd_mode)
29 {
30 int compare, inexact;
31 mp_size_t s;
32 mpfr_prec_t p, q;
33 mp_limb_t *up, *vp, *ufp, *vfp;
34 mpfr_t u, v, uf, vf, sc1, sc2;
35 mpfr_exp_t scaleop = 0, scaleit;
36 unsigned long n; /* number of iterations */
37 MPFR_ZIV_DECL (loop);
38 MPFR_TMP_DECL(marker);
39 MPFR_SAVE_EXPO_DECL (expo);
40
41 MPFR_LOG_FUNC
42 (("op2[%Pd]=%.*Rg op1[%Pd]=%.*Rg rnd=%d",
43 mpfr_get_prec (op2), mpfr_log_prec, op2,
44 mpfr_get_prec (op1), mpfr_log_prec, op1, rnd_mode),
45 ("r[%Pd]=%.*Rg inexact=%d",
46 mpfr_get_prec (r), mpfr_log_prec, r, inexact));
47
48 /* Deal with special values */
49 if (MPFR_ARE_SINGULAR (op1, op2))
50 {
51 /* If a or b is NaN, the result is NaN */
52 if (MPFR_IS_NAN(op1) || MPFR_IS_NAN(op2))
53 {
54 MPFR_SET_NAN(r);
55 MPFR_RET_NAN;
56 }
57 /* now one of a or b is Inf or 0 */
58 /* If a and b is +Inf, the result is +Inf.
59 Otherwise if a or b is -Inf or 0, the result is NaN */
60 else if (MPFR_IS_INF(op1) || MPFR_IS_INF(op2))
61 {
62 if (MPFR_IS_STRICTPOS(op1) && MPFR_IS_STRICTPOS(op2))
63 {
64 MPFR_SET_INF(r);
65 MPFR_SET_SAME_SIGN(r, op1);
66 MPFR_RET(0); /* exact */
67 }
68 else
69 {
70 MPFR_SET_NAN(r);
71 MPFR_RET_NAN;
72 }
73 }
74 else /* a and b are neither NaN nor Inf, and one is zero */
75 { /* If a or b is 0, the result is +0, in particular because the
76 result is always >= 0 with our definition (Maple sometimes
77 chooses a different sign for GaussAGM, but it uses another
78 definition, with possible negative results). */
79 MPFR_ASSERTD (MPFR_IS_ZERO (op1) || MPFR_IS_ZERO (op2));
80 MPFR_SET_POS (r);
81 MPFR_SET_ZERO (r);
82 MPFR_RET (0); /* exact */
83 }
84 }
85
86 /* If a or b is negative (excluding -Infinity), the result is NaN */
87 if (MPFR_UNLIKELY(MPFR_IS_NEG(op1) || MPFR_IS_NEG(op2)))
88 {
89 MPFR_SET_NAN(r);
90 MPFR_RET_NAN;
91 }
92
93 /* Precision of the following calculus */
94 q = MPFR_PREC(r);
95 p = q + MPFR_INT_CEIL_LOG2(q) + 15;
96 MPFR_ASSERTD (p >= 7); /* see algorithms.tex */
97 s = MPFR_PREC2LIMBS (p);
98
99 /* b (op2) and a (op1) are the 2 operands but we want b >= a */
100 compare = mpfr_cmp (op1, op2);
101 if (MPFR_UNLIKELY( compare == 0 ))
102 return mpfr_set (r, op1, rnd_mode);
103 else if (compare > 0)
104 {
105 mpfr_srcptr t = op1;
106 op1 = op2;
107 op2 = t;
108 }
109
110 /* Now b (=op2) > a (=op1) */
111
112 MPFR_SAVE_EXPO_MARK (expo);
113
114 MPFR_TMP_MARK(marker);
115
116 /* Main loop */
117 MPFR_ZIV_INIT (loop, p);
118 for (;;)
119 {
120 mpfr_prec_t eq;
121 unsigned long err = 0; /* must be set to 0 at each Ziv iteration */
122 MPFR_BLOCK_DECL (flags);
123
124 /* Init temporary vars */
125 MPFR_TMP_INIT (up, u, p, s);
126 MPFR_TMP_INIT (vp, v, p, s);
127 MPFR_TMP_INIT (ufp, uf, p, s);
128 MPFR_TMP_INIT (vfp, vf, p, s);
129
130 /* Calculus of un and vn */
131 retry:
132 MPFR_BLOCK (flags,
133 mpfr_mul (u, op1, op2, MPFR_RNDN);
134 /* mpfr_mul(...): faster since PREC(op) < PREC(u) */
135 mpfr_add (v, op1, op2, MPFR_RNDN);
136 /* mpfr_add with !=prec is still good */);
137 if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags)))
138 {
139 mpfr_exp_t e1, e2;
140
141 MPFR_ASSERTN (scaleop == 0);
142 e1 = MPFR_GET_EXP (op1);
143 e2 = MPFR_GET_EXP (op2);
144
145 /* Let's determine scaleop to avoid an overflow/underflow. */
146 if (MPFR_OVERFLOW (flags))
147 {
148 /* Let's recall that emin <= e1 <= e2 <= emax.
149 There has been an overflow. Thus e2 >= emax/2.
150 If the mpfr_mul overflowed, then e1 + e2 > emax.
151 If the mpfr_add overflowed, then e2 = emax.
152 We want: (e1 + scale) + (e2 + scale) <= emax,
153 i.e. scale <= (emax - e1 - e2) / 2. Let's take
154 scale = min(floor((emax - e1 - e2) / 2), -1).
155 This is OK, as:
156 1. emin <= scale <= -1.
157 2. e1 + scale >= emin. Indeed:
158 * If e1 + e2 > emax, then
159 e1 + scale >= e1 + (emax - e1 - e2) / 2 - 1
160 >= (emax + e1 - emax) / 2 - 1
161 >= e1 / 2 - 1 >= emin.
162 * Otherwise, mpfr_mul didn't overflow, therefore
163 mpfr_add overflowed and e2 = emax, so that
164 e1 > emin (see restriction below).
165 e1 + scale > emin - 1, thus e1 + scale >= emin.
166 3. e2 + scale <= emax, since scale < 0. */
167 if (e1 + e2 > MPFR_EMAX_MAX)
168 {
169 scaleop = - (((e1 + e2) - MPFR_EMAX_MAX + 1) / 2);
170 MPFR_ASSERTN (scaleop < 0);
171 }
172 else
173 {
174 /* The addition necessarily overflowed. */
175 MPFR_ASSERTN (e2 == MPFR_EMAX_MAX);
176 /* The case where e1 = emin and e2 = emax is not supported
177 here. This would mean that the precision of e2 would be
178 huge (and possibly not supported in practice anyway). */
179 MPFR_ASSERTN (e1 > MPFR_EMIN_MIN);
180 /* Note: this case is probably impossible to have in practice
181 since we need e2 = emax, and no overflow in the product.
182 Since the product is >= 2^(e1+e2-2), it implies
183 e1 + e2 - 2 <= emax, thus e1 <= 2. Now to get an overflow
184 we need op1 >= 1/2 ulp(op2), which implies that the
185 precision of op2 should be at least emax-2. On a 64-bit
186 computer this is impossible to have, and would require
187 a huge amount of memory on a 32-bit computer. */
188 scaleop = -1;
189 }
190
191 }
192 else /* underflow only (in the multiplication) */
193 {
194 /* We have e1 + e2 <= emin (so, e1 <= e2 <= 0).
195 We want: (e1 + scale) + (e2 + scale) >= emin + 1,
196 i.e. scale >= (emin + 1 - e1 - e2) / 2. let's take
197 scale = ceil((emin + 1 - e1 - e2) / 2). This is OK, as:
198 1. 1 <= scale <= emax.
199 2. e1 + scale >= emin + 1 >= emin.
200 3. e2 + scale <= scale <= emax. */
201 MPFR_ASSERTN (e1 <= e2 && e2 <= 0);
202 scaleop = (MPFR_EMIN_MIN + 2 - e1 - e2) / 2;
203 MPFR_ASSERTN (scaleop > 0);
204 }
205
206 MPFR_ALIAS (sc1, op1, MPFR_SIGN (op1), e1 + scaleop);
207 MPFR_ALIAS (sc2, op2, MPFR_SIGN (op2), e2 + scaleop);
208 op1 = sc1;
209 op2 = sc2;
210 MPFR_LOG_MSG (("Exception in pre-iteration, scale = %"
211 MPFR_EXP_FSPEC "d\n", scaleop));
212 goto retry;
213 }
214
215 MPFR_CLEAR_FLAGS ();
216 mpfr_sqrt (u, u, MPFR_RNDN);
217 mpfr_div_2ui (v, v, 1, MPFR_RNDN);
218
219 scaleit = 0;
220 n = 1;
221 while (mpfr_cmp2 (u, v, &eq) != 0 && eq <= p - 2)
222 {
223 MPFR_BLOCK_DECL (flags2);
224
225 MPFR_LOG_MSG (("Iteration n = %lu\n", n));
226
227 retry2:
228 mpfr_add (vf, u, v, MPFR_RNDN); /* No overflow? */
229 mpfr_div_2ui (vf, vf, 1, MPFR_RNDN);
230 /* See proof in algorithms.tex */
231 if (eq > p / 4)
232 {
233 mpfr_t w;
234 MPFR_BLOCK_DECL (flags3);
235
236 MPFR_LOG_MSG (("4*eq > p\n", 0));
237
238 /* vf = V(k) */
239 mpfr_init2 (w, (p + 1) / 2);
240 MPFR_BLOCK
241 (flags3,
242 mpfr_sub (w, v, u, MPFR_RNDN); /* e = V(k-1)-U(k-1) */
243 mpfr_sqr (w, w, MPFR_RNDN); /* e = e^2 */
244 mpfr_div_2ui (w, w, 4, MPFR_RNDN); /* e*= (1/2)^2*1/4 */
245 mpfr_div (w, w, vf, MPFR_RNDN); /* 1/4*e^2/V(k) */
246 );
247 if (MPFR_LIKELY (! MPFR_UNDERFLOW (flags3)))
248 {
249 mpfr_sub (v, vf, w, MPFR_RNDN);
250 err = MPFR_GET_EXP (vf) - MPFR_GET_EXP (v); /* 0 or 1 */
251 mpfr_clear (w);
252 break;
253 }
254 /* There has been an underflow because of the cancellation
255 between V(k-1) and U(k-1). Let's use the conventional
256 method. */
257 MPFR_LOG_MSG (("4*eq > p -> underflow\n", 0));
258 mpfr_clear (w);
259 MPFR_CLEAR_UNDERFLOW ();
260 }
261 /* U(k) increases, so that U.V can overflow (but not underflow). */
262 MPFR_BLOCK (flags2, mpfr_mul (uf, u, v, MPFR_RNDN););
263 if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags2)))
264 {
265 mpfr_exp_t scale2;
266
267 scale2 = - (((MPFR_GET_EXP (u) + MPFR_GET_EXP (v))
268 - MPFR_EMAX_MAX + 1) / 2);
269 MPFR_EXP (u) += scale2;
270 MPFR_EXP (v) += scale2;
271 scaleit += scale2;
272 MPFR_LOG_MSG (("Overflow in iteration n = %lu, scaleit = %"
273 MPFR_EXP_FSPEC "d (%" MPFR_EXP_FSPEC "d)\n",
274 n, scaleit, scale2));
275 MPFR_CLEAR_OVERFLOW ();
276 goto retry2;
277 }
278 mpfr_sqrt (u, uf, MPFR_RNDN);
279 mpfr_swap (v, vf);
280 n ++;
281 }
282
283 MPFR_LOG_MSG (("End of iterations (n = %lu)\n", n));
284
285 /* the error on v is bounded by (18n+51) ulps, or twice if there
286 was an exponent loss in the final subtraction */
287 err += MPFR_INT_CEIL_LOG2(18 * n + 51); /* 18n+51 should not overflow
288 since n is about log(p) */
289 /* we should have n+2 <= 2^(p/4) [see algorithms.tex] */
290 if (MPFR_LIKELY (MPFR_INT_CEIL_LOG2(n + 2) <= p / 4 &&
291 MPFR_CAN_ROUND (v, p - err, q, rnd_mode)))
292 break; /* Stop the loop */
293
294 /* Next iteration */
295 MPFR_ZIV_NEXT (loop, p);
296 s = MPFR_PREC2LIMBS (p);
297 }
298 MPFR_ZIV_FREE (loop);
299
300 if (MPFR_UNLIKELY ((__gmpfr_flags & (MPFR_FLAGS_ALL ^ MPFR_FLAGS_INEXACT))
301 != 0))
302 {
303 MPFR_ASSERTN (! mpfr_overflow_p ()); /* since mpfr_clear_flags */
304 MPFR_ASSERTN (! mpfr_underflow_p ()); /* since mpfr_clear_flags */
305 MPFR_ASSERTN (! mpfr_divby0_p ()); /* since mpfr_clear_flags */
306 MPFR_ASSERTN (! mpfr_nanflag_p ()); /* since mpfr_clear_flags */
307 }
308
309 /* Setting of the result */
310 inexact = mpfr_set (r, v, rnd_mode);
311 MPFR_EXP (r) -= scaleop + scaleit;
312
313 /* Let's clean */
314 MPFR_TMP_FREE(marker);
315
316 MPFR_SAVE_EXPO_FREE (expo);
317 /* From the definition of the AGM, underflow and overflow
318 are not possible. */
319 return mpfr_check_range (r, inexact, rnd_mode);
320 /* agm(u,v) can be exact for u, v rational only for u=v.
321 Proof (due to Nicolas Brisebarre): it suffices to consider
322 u=1 and v<1. Then 1/AGM(1,v) = 2F1(1/2,1/2,1;1-v^2),
323 and a theorem due to G.V. Chudnovsky states that for x a
324 non-zero algebraic number with |x|<1, then
325 2F1(1/2,1/2,1;x) and 2F1(-1/2,1/2,1;x) are algebraically
326 independent over Q. */
327 }