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 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 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 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
Indexes created Sun Sep 13 00:25:22 UTC 2026