1/* tballs -- test file for complex ball arithmetic.
2
3Copyright (C) 2018, 2020, 2021, 2022 INRIA
4
5This file is part of GNU MPC.
6
7GNU MPC is free software; you can redistribute it and/or modify it under
8the terms of the GNU Lesser General Public License as published by the
9Free Software Foundation; either version 3 of the License, or (at your
10option) any later version.
11
12GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
13WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15more details.
16
17You should have received a copy of the GNU Lesser General Public License
18along 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
26static int
27mpc_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
180static int
181test_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
211int
212main (void)
213{
214   return test_agm ();
215}
216
217