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
mpc_mpcb_agm(mpc_ptr rop,mpc_srcptr opa,mpc_srcptr opb,mpc_rnd_t rnd)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
test_agm(void)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
main(void)212 main (void)
213 {
214 return test_agm ();
215 }
216
217