1 /* tballs -- test file for complex ball arithmetic. 2 3 Copyright (C) 2018, 2020, 2021, 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-tests.h" 22 #include "mpc-impl.h" 23 /* For the alternative AGM implementation, we need all the power of 24 this include file. */ 25 26 static int 27 mpc_mpcb_agm (mpc_ptr rop, mpc_srcptr opa, mpc_srcptr opb, mpc_rnd_t rnd) 28 /* Alternative implementation of mpc_agm that uses complex balls. */ 29 { 30 mpfr_prec_t prec; 31 mpc_t b0, diff; 32 mpcb_t a, b, an, bn, anp1, bnp1, res; 33 mpfr_exp_t exp_an, exp_diff; 34 mpcr_t rab; 35 int cmp, equal, re_zero, im_zero, ok, inex; 36 37 if (!mpc_fin_p (opa) || !mpc_fin_p (opb) 38 || mpc_zero_p (opa) || mpc_zero_p (opb) 39 || mpc_cmp (opa, opb) == 0 40 || ( mpfr_sgn (mpc_realref (opa)) == -mpfr_sgn (mpc_realref (opb)) 41 && mpfr_sgn (mpc_imagref (opa)) == -mpfr_sgn (mpc_imagref (opb)) 42 && mpfr_cmpabs (mpc_realref (opa), mpc_realref (opb)) == 0 43 && mpfr_cmpabs (mpc_imagref (opa), mpc_imagref (opb)) == 0) 44 || ( mpfr_zero_p (mpc_imagref (opa)) 45 && mpfr_zero_p (mpc_imagref (opb)) 46 && mpfr_sgn (mpc_realref (opa)) == mpfr_sgn (mpc_realref (opb))) 47 || ( mpfr_zero_p (mpc_realref (opa)) 48 && mpfr_zero_p (mpc_realref (opb)) 49 && mpfr_sgn (mpc_imagref (opa)) == mpfr_sgn (mpc_imagref (opb)))) 50 /* Special cases that are handled separately by mpc_agm; there is 51 no need to rewrite them. */ 52 return mpc_agm (rop, opa, opb, rnd); 53 54 /* Exclude the case of angle 0, also handled separately by mpc_agm. */ 55 mpc_init2 (b0, 2); 56 mpc_div (b0, opb, opa, MPC_RNDZZ); 57 if (mpfr_zero_p (mpc_imagref (b0)) && mpfr_sgn (mpc_realref (b0)) > 0) { 58 mpc_clear (b0); 59 return mpc_agm (rop, opa, opb, rnd); 60 } 61 mpc_clear (b0); 62 63 cmp = mpc_cmp_abs (opa, opb); 64 65 mpcb_init (a); 66 mpcb_init (b); 67 mpcb_init (an); 68 mpcb_init (bn); 69 mpcb_init (anp1); 70 mpcb_init (bnp1); 71 mpcb_init (res); 72 prec = MPC_MAX (MPC_MAX (MPC_MAX_PREC (opa), MPC_MAX_PREC (opb)), 73 MPC_MAX_PREC (rop) + 20); 74 /* So copying opa and opb will be exact, and there is a small safety 75 margin for the result. */ 76 do { 77 mpcb_set_prec (a, prec); 78 mpcb_set_prec (b, prec); 79 mpcb_set_prec (an, prec); 80 mpcb_set_prec (bn, prec); 81 mpcb_set_prec (anp1, prec); 82 mpcb_set_prec (bnp1, prec); 83 mpcb_set_prec (res, prec); 84 /* TODO: Think about the mpcb_set variants; mpcb_set_c, for instance, 85 modifies the precision. It is probably better to add a precision 86 parameter to mpcb_init and potentially round with mpcb_set_xxx. */ 87 mpc_set (a->c, opa, MPC_RNDNN); /* exact */ 88 mpcr_set_zero (a->r); 89 mpc_set (b->c, opb, MPC_RNDNN); 90 mpcr_set_zero (b->r); 91 mpc_set_ui_ui (an->c, 1, 0, MPC_RNDNN); 92 mpcr_set_zero (an->r); 93 if (cmp >= 0) 94 mpcb_div (bn, b, a); 95 else 96 mpcb_div (bn, a, b); 97 98 /* Iterate until there is a fixed point or (often one iteration 99 earlier) the arithmetic and the geometric mean coincide. */ 100 do { 101 mpcb_add (anp1, an, bn); 102 mpcb_div_2ui (anp1, anp1, 1); 103 mpcb_mul (bnp1, an, bn); 104 mpcb_sqrt (bnp1, bnp1); 105 /* Be aware of the branch cut! The current function does 106 what is needed here. */ 107 equal = mpc_cmp (an->c, bn->c) == 0 108 || ( mpc_cmp (an->c, anp1->c) == 0 109 && mpc_cmp (bn->c, bnp1->c) == 0); 110 mpcb_set (an, anp1); 111 mpcb_set (bn, bnp1); 112 } while (!equal); 113 114 /* Check whether we can conclude, see the error analysis in 115 algorithms.tex. */ 116 if (mpcr_inf_p (anp1->r)) 117 ok = 0; 118 else { 119 mpc_init2 (diff, prec); 120 mpc_sub (diff, an->c, bn->c, MPC_RNDZZ); 121 /* FIXME: We would need to round away, but this is not yet 122 implemented. */ 123 re_zero = mpfr_zero_p (mpc_realref (diff)); 124 if (!re_zero) 125 MPFR_ADD_ONE_ULP (mpc_realref (diff)); 126 im_zero = mpfr_zero_p (mpc_imagref (diff)); 127 if (!im_zero) 128 MPFR_ADD_ONE_ULP (mpc_imagref (diff)); 129 130 mpcb_set (res, anp1); 131 132 if (re_zero && im_zero) 133 mpcr_set_zero (rab); 134 else { 135 exp_an = MPC_MIN (mpfr_get_exp (mpc_realref (an->c)), 136 mpfr_get_exp (mpc_imagref (an->c))) - 1; 137 if (re_zero) 138 exp_diff = mpfr_get_exp (mpc_imagref (diff)) + 1; 139 else if (im_zero) 140 exp_diff = mpfr_get_exp (mpc_realref (diff)) + 1; 141 else 142 exp_diff = MPC_MAX (mpfr_get_exp (mpc_realref (diff)), 143 mpfr_get_exp (mpc_imagref (diff)) + 1); 144 mpcr_set_one (rab); 145 (rab->exp) += (exp_diff - exp_an); 146 /* TODO: Should be done by an mpcr function. */ 147 } 148 mpcr_add (rab, rab, an->r); 149 (rab->exp)++; 150 mpcr_add (res->r, rab, bn->r); 151 /* r = 2 * (rab + an->r) + bn->r */ 152 if (cmp >= 0) 153 mpcb_mul (res, res, a); 154 else 155 mpcb_mul (res, res, b); 156 ok = mpcb_can_round (res, MPC_PREC_RE (rop), MPC_PREC_IM (rop), 157 rnd); 158 159 mpc_clear (diff); 160 } 161 162 if (!ok) 163 prec += prec + mpcr_get_exp (res->r); 164 } while (!ok); 165 166 inex = mpcb_round (rop, res, rnd); 167 168 mpcb_clear (a); 169 mpcb_clear (b); 170 mpcb_clear (an); 171 mpcb_clear (bn); 172 mpcb_clear (anp1); 173 mpcb_clear (bnp1); 174 mpcb_clear (res); 175 176 return inex; 177 } 178 179 180 static int 181 test_agm (void) 182 { 183 mpfr_prec_t prec; 184 mpc_t a, b, agm1, agm2; 185 mpc_rnd_t rnd = MPC_RNDDU; 186 int inex, inexb, ok; 187 188 prec = 1000; 189 190 mpc_init2 (a, prec); 191 mpc_init2 (b, prec); 192 mpc_set_si_si (a, 100, 0, MPC_RNDNN); 193 mpc_set_si_si (b, 0, 100, MPC_RNDNN); 194 mpc_init2 (agm1, prec); 195 mpc_init2 (agm2, prec); 196 197 inex = mpc_agm (agm1, a, b, rnd); 198 inexb = mpc_mpcb_agm (agm2, a, b, rnd); 199 200 ok = (inex == inexb) && (mpc_cmp (agm1, agm2) == 0); 201 202 mpc_clear (a); 203 mpc_clear (b); 204 mpc_clear (agm1); 205 mpc_clear (agm2); 206 207 return !ok; 208 } 209 210 211 int 212 main (void) 213 { 214 return test_agm (); 215 } 216 217
Indexes created Sat Sep 12 00:25:55 UTC 2026