1 /* mpc_agm -- AGM of a complex number.
2
3 Copyright (C) 2022 INRIA
4
5 This file is part of GNU MPC.
6
7 GNU MPC is free software; you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
12 GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15 more details.
16
17 You should have received a copy of the GNU Lesser General Public License
18 along with this program. If not, see http://www.gnu.org/licenses/ .
19 */
20
21 #include "mpc-impl.h"
22
23 static int
mpc_agm_angle_zero(mpc_ptr rop,mpc_srcptr a,mpc_srcptr b,mpc_rnd_t rnd,int cmp)24 mpc_agm_angle_zero (mpc_ptr rop, mpc_srcptr a, mpc_srcptr b, mpc_rnd_t rnd,
25 int cmp)
26 /* AGM for angle 0 between a and b, but they are neither real nor
27 purely imaginary. cmp is mpc_cmp_abs (a, b). */
28 {
29 mpfr_prec_t prec;
30 int inex;
31 mpc_t agm;
32 mpfr_t a0, b0;
33
34 prec = MPC_MAX_PREC (rop);
35
36 mpc_init2 (agm, 2);
37 mpfr_init2 (a0, 2);
38 mpfr_set_ui (a0, 1, MPFR_RNDN);
39 mpfr_init2 (b0, 2);
40
41 do {
42 prec += 20;
43 mpc_set_prec (agm, prec);
44 mpfr_set_prec (b0, prec);
45
46 if (cmp >= 0)
47 mpfr_div (b0, mpc_realref (b), mpc_realref (a), MPFR_RNDZ);
48 else
49 mpfr_div (b0, mpc_realref (a), mpc_realref (b), MPFR_RNDZ);
50
51 mpfr_agm (b0, a0, b0, MPFR_RNDN);
52
53 if (cmp >= 0)
54 mpc_mul_fr (agm, a, b0, MPC_RNDNN);
55 else
56 mpc_mul_fr (agm, b, b0, MPC_RNDNN);
57
58 } while (!mpfr_can_round (mpc_realref (agm), prec - 3,
59 MPFR_RNDN, MPFR_RNDU, mpfr_get_prec (mpc_realref (rop)) + 1)
60 || !mpfr_can_round (mpc_imagref (agm), prec - 3,
61 MPFR_RNDN, MPFR_RNDU, mpfr_get_prec (mpc_imagref (rop)) + 1));
62
63 inex = mpc_set (rop, agm, rnd);
64
65 mpc_clear (agm);
66 mpfr_clear (a0);
67 mpfr_clear (b0);
68
69 return inex;
70 }
71
72
73 static int
mpc_agm_general(mpc_ptr rop,mpc_srcptr a,mpc_srcptr b,mpc_rnd_t rnd)74 mpc_agm_general (mpc_ptr rop, mpc_srcptr a, mpc_srcptr b, mpc_rnd_t rnd)
75 /* AGM for not extremely special numbers:
76 Finite, non-zero, and a != -b; if the angle is 0, then we are neither
77 in the real nor in the purely imaginary case.
78 We follow the strategy outlined in algorithms.tex. */
79 {
80 mpc_t b0, a1, an, bn, anp1, bnp1;
81 int cmp, equal, n, i;
82 mpfr_prec_t prec, N, k1, L, exp_diff, kR, kI;
83 mpfr_exp_t exp_re_a1, exp_re_b0, exp_im_b0;
84 int okR, okI, inex;
85
86 /* Determine whether to compute AGM (1, b0) with b0 = a/b or b0 = b/a. */
87 cmp = mpc_cmp_abs (a, b);
88
89 /* Compute an approximation k1 of the precision loss in the first
90 iteration. */
91 mpc_init2 (b0, 2);
92 mpc_init2 (a1, 2);
93 prec = 1;
94 do {
95 prec *= 2;
96 mpc_set_prec (b0, prec);
97 mpc_set_prec (a1, prec);
98 if (cmp >= 0)
99 mpc_div (b0, b, a, MPC_RNDZZ);
100 else
101 mpc_div (b0, a, b, MPC_RNDZZ);
102 if (mpfr_zero_p (mpc_imagref (b0))
103 && mpfr_sgn (mpc_realref (b0)) > 0) {
104 mpc_clear (b0);
105 mpc_clear (a1);
106 return mpc_agm_angle_zero (rop, a, b, rnd, cmp);
107 }
108 mpc_add_ui (a1, b0, 1, MPC_RNDNN);
109 mpc_div_2ui (a1, a1, 1, MPC_RNDNN);
110 exp_re_a1 = mpfr_get_exp (mpc_realref (a1));
111 } while (exp_re_a1 == -prec);
112 exp_re_b0 = mpfr_get_exp (mpc_realref (b0));
113 exp_im_b0 = mpfr_get_exp (mpc_imagref (b0));
114 mpc_clear (a1);
115 mpc_clear (b0);
116 k1 = MPC_MAX (3, - 2 * exp_re_a1 - 2);
117
118 /* Compute the number n of iterations and the target precision. */
119 N = MPC_MAX_PREC (rop) + 20;
120 /* 20 is an arbitrary safety margin. */
121 do {
122 /* With the notation of algorithms.tex, compute 2*L, which is
123 an integer; then correct this when taking the logarithm. */
124 if (exp_im_b0 <= -1)
125 if (exp_re_b0 <= -1)
126 L = MPC_MAX (6, - MPC_MAX (exp_re_b0, exp_im_b0) + 1);
127 else if (exp_re_a1 <= -2)
128 L = - 2 * MPC_MAX (exp_re_a1, exp_im_b0 - 1) + 3;
129 else
130 L = 6;
131 else
132 L = 6;
133 L = MPC_MAX (1, mpc_ceil_log2 (L) - 1);
134 n = L + mpc_ceil_log2 (N + 4) + 3;
135 prec = N + (n + k1 + 7 + 1) / 2;
136
137 mpc_init2 (an, prec);
138 mpc_init2 (bn, prec);
139 mpc_init2 (anp1, prec);
140 mpc_init2 (bnp1, prec);
141
142 /* Compute the argument for AGM (1, b0) at the working precision. */
143 if (cmp >= 0)
144 mpc_div (bn, b, a, MPC_RNDZZ);
145 else
146 mpc_div (bn, a, b, MPC_RNDZZ);
147 mpc_set_ui (an, 1, MPC_RNDNN);
148
149 equal = 0;
150 /* In practice, a fixed point of the AGM operation is reached before
151 the last iteration, so we may stop when an==anp1 and bn==bnp1.
152 Also in practice one observes that often one iteration earlier one
153 has an==bn, which is also tested for an early abort strategy. */
154 /* Execute the AGM iterations. */
155 for (i = 1; !equal && i <= n; i++) {
156 mpc_add (anp1, an, bn, MPC_RNDNN);
157 mpc_div_2ui (anp1, anp1, 1, MPC_RNDNN);
158 mpc_mul (bnp1, an, bn, MPC_RNDNN);
159 mpc_sqrt (bnp1, bnp1, MPC_RNDNN);
160 equal = (mpc_cmp (an, anp1) == 0 && mpc_cmp (bn, bnp1) == 0)
161 || mpc_cmp (anp1, bnp1) == 0;
162 mpc_swap (an, anp1);
163 mpc_swap (bn, bnp1);
164 }
165
166 /* Remultiply. */
167 if (cmp >= 0)
168 mpc_mul (an, an, a, MPC_RNDNN);
169 else
170 mpc_mul (an, an, b, MPC_RNDNN);
171
172 exp_diff = mpfr_get_exp (mpc_imagref (an))
173 - mpfr_get_exp (mpc_realref (an));
174 kR = MPC_MAX (exp_diff + 1, 0);
175 kI = MPC_MAX (-exp_diff + 1, 0);
176 /* Use the trick of asking mpfr_can_round whether it can do a directed
177 rounding at precision + 1; then the whole uncertainty interval is
178 contained in the upper or the lower half of the interval between two
179 representable numbers, and mpc_set reveals the inexact value
180 regardless of the rounding direction. */
181 okR = mpfr_can_round (mpc_realref (an), N - kR,
182 MPFR_RNDN, MPFR_RNDU, mpfr_get_prec (mpc_realref (rop)) + 1);
183 okI = mpfr_can_round (mpc_imagref (an), N - kI,
184 MPFR_RNDN, MPFR_RNDU, mpfr_get_prec (mpc_imagref (rop)) + 1);
185
186 if (!okR || !okI)
187 /* Until a counterexample is found, we assume that this happens
188 only once and increase the precision only moderately. */
189 N += MPC_MAX (kR, kI);
190 } while (!okR || !okI);
191
192 inex = mpc_set (rop, an, rnd);
193
194 mpc_clear (an);
195 mpc_clear (bn);
196 mpc_clear (anp1);
197 mpc_clear (bnp1);
198
199 return inex;
200 }
201
202
203 int
mpc_agm(mpc_ptr rop,mpc_srcptr a,mpc_srcptr b,mpc_rnd_t rnd)204 mpc_agm (mpc_ptr rop, mpc_srcptr a, mpc_srcptr b, mpc_rnd_t rnd)
205 {
206 int inex_re, inex_im;
207
208 if (!mpc_fin_p (a) || !mpc_fin_p (b)) {
209 mpfr_set_nan (mpc_realref (rop));
210 mpfr_set_nan (mpc_imagref (rop));
211 return MPC_INEX (0, 0);
212 }
213 else if (mpc_zero_p (a) || mpc_zero_p (b))
214 return mpc_set_ui_ui (rop, 0, 0, MPC_RNDNN);
215 else if (mpc_cmp (a, b) == 0)
216 return mpc_set (rop, a, rnd);
217 else if ( mpfr_sgn (mpc_realref (a)) == -mpfr_sgn (mpc_realref (b))
218 && mpfr_sgn (mpc_imagref (a)) == -mpfr_sgn (mpc_imagref (b))
219 && mpfr_cmpabs (mpc_realref (a), mpc_realref (b)) == 0
220 && mpfr_cmpabs (mpc_imagref (a), mpc_imagref (b)) == 0)
221 /* a == -b */
222 return mpc_set_ui_ui (rop, 0, 0, MPC_RNDNN);
223 else if (mpfr_zero_p (mpc_imagref (a))
224 && mpfr_zero_p (mpc_imagref (b))
225 && mpfr_sgn (mpc_realref (a)) == mpfr_sgn (mpc_realref (b))) {
226 /* angle 0, real values */
227 inex_re = mpfr_agm (mpc_realref (rop), mpc_realref (a),
228 mpc_realref (b), MPC_RND_RE (rnd));
229 mpfr_set_ui (mpc_imagref (rop), 0, MPFR_RNDN);
230 return MPC_INEX (inex_re, 0);
231 }
232 else if (mpfr_zero_p (mpc_realref (a))
233 && mpfr_zero_p (mpc_realref (b))
234 && mpfr_sgn (mpc_imagref (a)) == mpfr_sgn (mpc_imagref (b))) {
235 /* angle 0, purely imaginary values */
236 inex_im = mpfr_agm (mpc_imagref (rop), mpc_imagref (a),
237 mpc_imagref (b), MPC_RND_IM (rnd));
238 mpfr_set_ui (mpc_realref (rop), 0, MPFR_RNDN);
239 return MPC_INEX (0, inex_im);
240 }
241 else
242 return mpc_agm_general (rop, a, b, rnd);
243 }
244
245