/* tballs -- test file for complex ball arithmetic. Copyright (C) 2018, 2020, 2021, 2022 INRIA This file is part of GNU MPC. GNU MPC is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this program. If not, see http://www.gnu.org/licenses/ . */ #include "mpc-tests.h" #include "mpc-impl.h" /* For the alternative AGM implementation, we need all the power of this include file. */ static int mpc_mpcb_agm (mpc_ptr rop, mpc_srcptr opa, mpc_srcptr opb, mpc_rnd_t rnd) /* Alternative implementation of mpc_agm that uses complex balls. */ { mpfr_prec_t prec; mpc_t b0, diff; mpcb_t a, b, an, bn, anp1, bnp1, res; mpfr_exp_t exp_an, exp_diff; mpcr_t rab; int cmp, equal, re_zero, im_zero, ok, inex; if (!mpc_fin_p (opa) || !mpc_fin_p (opb) || mpc_zero_p (opa) || mpc_zero_p (opb) || mpc_cmp (opa, opb) == 0 || ( mpfr_sgn (mpc_realref (opa)) == -mpfr_sgn (mpc_realref (opb)) && mpfr_sgn (mpc_imagref (opa)) == -mpfr_sgn (mpc_imagref (opb)) && mpfr_cmpabs (mpc_realref (opa), mpc_realref (opb)) == 0 && mpfr_cmpabs (mpc_imagref (opa), mpc_imagref (opb)) == 0) || ( mpfr_zero_p (mpc_imagref (opa)) && mpfr_zero_p (mpc_imagref (opb)) && mpfr_sgn (mpc_realref (opa)) == mpfr_sgn (mpc_realref (opb))) || ( mpfr_zero_p (mpc_realref (opa)) && mpfr_zero_p (mpc_realref (opb)) && mpfr_sgn (mpc_imagref (opa)) == mpfr_sgn (mpc_imagref (opb)))) /* Special cases that are handled separately by mpc_agm; there is no need to rewrite them. */ return mpc_agm (rop, opa, opb, rnd); /* Exclude the case of angle 0, also handled separately by mpc_agm. */ mpc_init2 (b0, 2); mpc_div (b0, opb, opa, MPC_RNDZZ); if (mpfr_zero_p (mpc_imagref (b0)) && mpfr_sgn (mpc_realref (b0)) > 0) { mpc_clear (b0); return mpc_agm (rop, opa, opb, rnd); } mpc_clear (b0); cmp = mpc_cmp_abs (opa, opb); mpcb_init (a); mpcb_init (b); mpcb_init (an); mpcb_init (bn); mpcb_init (anp1); mpcb_init (bnp1); mpcb_init (res); prec = MPC_MAX (MPC_MAX (MPC_MAX_PREC (opa), MPC_MAX_PREC (opb)), MPC_MAX_PREC (rop) + 20); /* So copying opa and opb will be exact, and there is a small safety margin for the result. */ do { mpcb_set_prec (a, prec); mpcb_set_prec (b, prec); mpcb_set_prec (an, prec); mpcb_set_prec (bn, prec); mpcb_set_prec (anp1, prec); mpcb_set_prec (bnp1, prec); mpcb_set_prec (res, prec); /* TODO: Think about the mpcb_set variants; mpcb_set_c, for instance, modifies the precision. It is probably better to add a precision parameter to mpcb_init and potentially round with mpcb_set_xxx. */ mpc_set (a->c, opa, MPC_RNDNN); /* exact */ mpcr_set_zero (a->r); mpc_set (b->c, opb, MPC_RNDNN); mpcr_set_zero (b->r); mpc_set_ui_ui (an->c, 1, 0, MPC_RNDNN); mpcr_set_zero (an->r); if (cmp >= 0) mpcb_div (bn, b, a); else mpcb_div (bn, a, b); /* Iterate until there is a fixed point or (often one iteration earlier) the arithmetic and the geometric mean coincide. */ do { mpcb_add (anp1, an, bn); mpcb_div_2ui (anp1, anp1, 1); mpcb_mul (bnp1, an, bn); mpcb_sqrt (bnp1, bnp1); /* Be aware of the branch cut! The current function does what is needed here. */ equal = mpc_cmp (an->c, bn->c) == 0 || ( mpc_cmp (an->c, anp1->c) == 0 && mpc_cmp (bn->c, bnp1->c) == 0); mpcb_set (an, anp1); mpcb_set (bn, bnp1); } while (!equal); /* Check whether we can conclude, see the error analysis in algorithms.tex. */ if (mpcr_inf_p (anp1->r)) ok = 0; else { mpc_init2 (diff, prec); mpc_sub (diff, an->c, bn->c, MPC_RNDZZ); /* FIXME: We would need to round away, but this is not yet implemented. */ re_zero = mpfr_zero_p (mpc_realref (diff)); if (!re_zero) MPFR_ADD_ONE_ULP (mpc_realref (diff)); im_zero = mpfr_zero_p (mpc_imagref (diff)); if (!im_zero) MPFR_ADD_ONE_ULP (mpc_imagref (diff)); mpcb_set (res, anp1); if (re_zero && im_zero) mpcr_set_zero (rab); else { exp_an = MPC_MIN (mpfr_get_exp (mpc_realref (an->c)), mpfr_get_exp (mpc_imagref (an->c))) - 1; if (re_zero) exp_diff = mpfr_get_exp (mpc_imagref (diff)) + 1; else if (im_zero) exp_diff = mpfr_get_exp (mpc_realref (diff)) + 1; else exp_diff = MPC_MAX (mpfr_get_exp (mpc_realref (diff)), mpfr_get_exp (mpc_imagref (diff)) + 1); mpcr_set_one (rab); (rab->exp) += (exp_diff - exp_an); /* TODO: Should be done by an mpcr function. */ } mpcr_add (rab, rab, an->r); (rab->exp)++; mpcr_add (res->r, rab, bn->r); /* r = 2 * (rab + an->r) + bn->r */ if (cmp >= 0) mpcb_mul (res, res, a); else mpcb_mul (res, res, b); ok = mpcb_can_round (res, MPC_PREC_RE (rop), MPC_PREC_IM (rop), rnd); mpc_clear (diff); } if (!ok) prec += prec + mpcr_get_exp (res->r); } while (!ok); inex = mpcb_round (rop, res, rnd); mpcb_clear (a); mpcb_clear (b); mpcb_clear (an); mpcb_clear (bn); mpcb_clear (anp1); mpcb_clear (bnp1); mpcb_clear (res); return inex; } static int test_agm (void) { mpfr_prec_t prec; mpc_t a, b, agm1, agm2; mpc_rnd_t rnd = MPC_RNDDU; int inex, inexb, ok; prec = 1000; mpc_init2 (a, prec); mpc_init2 (b, prec); mpc_set_si_si (a, 100, 0, MPC_RNDNN); mpc_set_si_si (b, 0, 100, MPC_RNDNN); mpc_init2 (agm1, prec); mpc_init2 (agm2, prec); inex = mpc_agm (agm1, a, b, rnd); inexb = mpc_mpcb_agm (agm2, a, b, rnd); ok = (inex == inexb) && (mpc_cmp (agm1, agm2) == 0); mpc_clear (a); mpc_clear (b); mpc_clear (agm1); mpc_clear (agm2); return !ok; } int main (void) { return test_agm (); }