1 /* Test file for mpfr_agm.
2
3 Copyright 1999, 2001-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 #include "mpfr-test.h"
24
25 #define check(a,b,r) check4(a,b,r,0.0)
26
27 static void
28 check4 (const char *as, const char *bs, mpfr_rnd_t rnd_mode,
29 const char *res, int inex)
30 {
31 mpfr_t ta, tb, tc, tres;
32 mpfr_exp_t emin, emax;
33 int i;
34
35 emin = mpfr_get_emin ();
36 emax = mpfr_get_emax ();
37
38 mpfr_inits2 (53, ta, tb, tc, tres, (mpfr_ptr) 0);
39
40 for (i = 0; i <= 2; i++)
41 {
42 unsigned int expflags, newflags;
43 int inex2;
44
45 mpfr_set_str1 (ta, as);
46 mpfr_set_str1 (tb, bs);
47 mpfr_set_str1 (tc, res);
48
49 if (i > 0)
50 {
51 mpfr_exp_t ea, eb, ec, e0;
52
53 set_emin (MPFR_EMIN_MIN);
54 set_emax (MPFR_EMAX_MAX);
55
56 ea = mpfr_get_exp (ta);
57 eb = mpfr_get_exp (tb);
58 ec = mpfr_get_exp (tc);
59
60 e0 = i == 1 ? __gmpfr_emin : __gmpfr_emax;
61 if ((i == 1 && ea < eb) || (i == 2 && ea > eb))
62 {
63 mpfr_set_exp (ta, e0);
64 mpfr_set_exp (tb, e0 + (eb - ea));
65 mpfr_set_exp (tc, e0 + (ec - ea));
66 }
67 else
68 {
69 mpfr_set_exp (ta, e0 + (ea - eb));
70 mpfr_set_exp (tb, e0);
71 mpfr_set_exp (tc, e0 + (ec - eb));
72 }
73 }
74
75 __gmpfr_flags = expflags =
76 RAND_BOOL () ? MPFR_FLAGS_ALL ^ MPFR_FLAGS_ERANGE : 0;
77 inex2 = mpfr_agm (tres, ta, tb, rnd_mode);
78 newflags = __gmpfr_flags;
79 expflags |= MPFR_FLAGS_INEXACT;
80
81 if (VSIGN (inex2) != inex || newflags != expflags ||
82 ! mpfr_equal_p (tres, tc))
83 {
84 printf ("mpfr_agm failed in rnd_mode=%s for\n",
85 mpfr_print_rnd_mode (rnd_mode));
86 printf (" a = ");
87 mpfr_out_str (stdout, 10, 0, ta, MPFR_RNDN);
88 printf ("\n");
89 printf (" b = ");
90 mpfr_out_str (stdout, 10, 0, tb, MPFR_RNDN);
91 printf ("\n");
92 printf ("expected inex = %d, flags = %u,\n"
93 " ", inex, expflags);
94 mpfr_dump (tc);
95 printf ("got inex = %d, flags = %u,\n"
96 " ", inex2, newflags);
97 mpfr_dump (tres);
98 exit (1);
99 }
100
101 set_emin (emin);
102 set_emax (emax);
103 }
104
105 mpfr_clears (ta, tb, tc, tres, (mpfr_ptr) 0);
106 }
107
108 static void
109 check_large (void)
110 {
111 mpfr_t a, b, agm;
112 int inex;
113
114 mpfr_init2 (a, 82);
115 mpfr_init2 (b, 82);
116 mpfr_init2 (agm, 82);
117
118 mpfr_set_ui (a, 1, MPFR_RNDN);
119 mpfr_set_str_binary (b, "0.1111101100001000000001011000110111101000001011111000100001000101010100011111110010E-39");
120 mpfr_agm (agm, a, b, MPFR_RNDN);
121 mpfr_set_str_binary (a, "0.1110001000111101101010101010101101001010001001001011100101111011110101111001111100E-4");
122 if (mpfr_cmp (agm, a))
123 {
124 printf ("mpfr_agm failed for precision 82\n");
125 exit (1);
126 }
127
128 /* problem found by Damien Fischer <damien@maths.usyd.edu.au> 4 Aug 2003:
129 produced a divide-by-zero exception */
130 mpfr_set_prec (a, 268);
131 mpfr_set_prec (b, 268);
132 mpfr_set_prec (agm, 268);
133 mpfr_set_str (a, "703.93543315330225238487276503953366664991725089988315253092140138947103394917006", 10, MPFR_RNDN);
134 mpfr_set_str (b, "703.93543315330225238487279020523738740563816490895994499256063816906728642622316", 10, MPFR_RNDN);
135 mpfr_agm (agm, a, b, MPFR_RNDN);
136
137 mpfr_set_prec (a, 18);
138 mpfr_set_prec (b, 70);
139 mpfr_set_prec (agm, 67);
140 mpfr_set_str_binary (a, "0.111001111100101000e8");
141 mpfr_set_str_binary (b, "0.1101110111100100010100110000010111011011011100110100111001010100100001e10");
142 inex = mpfr_agm (agm, a, b, MPFR_RNDN);
143 mpfr_set_str_binary (b, "0.1111110010011101101100010101011011010010010000001010100011000110011e9");
144 if (mpfr_cmp (agm, b))
145 {
146 printf ("Error in mpfr_agm (1)\n");
147 exit (1);
148 }
149 if (inex >= 0)
150 {
151 printf ("Wrong flag for mpfr_agm (1)\n");
152 exit (1);
153 }
154
155 /* test worst case: 9 consecutive ones after the last bit */
156 mpfr_set_prec (a, 2);
157 mpfr_set_prec (b, 2);
158 mpfr_set_ui (a, 1, MPFR_RNDN);
159 mpfr_set_ui (b, 2, MPFR_RNDN);
160 mpfr_set_prec (agm, 904);
161 mpfr_agm (agm, a, b, MPFR_RNDZ);
162
163 mpfr_clear (a);
164 mpfr_clear (b);
165 mpfr_clear (agm);
166 }
167
168 static void
169 check_eq (void)
170 {
171 mpfr_t a, b, agm;
172 int p;
173
174 mpfr_init2 (a, 17);
175 mpfr_init2 (b, 9);
176
177 mpfr_set_str_binary (b, "0.101000000E-3");
178 mpfr_set (a, b, MPFR_RNDN);
179
180 for (p = MPFR_PREC_MIN; p <= 2; p++)
181 {
182 int inex;
183
184 mpfr_init2 (agm, p);
185 inex = mpfr_agm (agm, a, b, MPFR_RNDU);
186 if (mpfr_cmp_ui_2exp (agm, 5 - p, -5) != 0)
187 {
188 printf ("Error in check_eq for p = %d: expected %d*2^(-5), got ",
189 p, 5 - p);
190 mpfr_dump (agm);
191 exit (1);
192 }
193 if (inex <= 0)
194 {
195 printf ("Wrong ternary value in check_eq for p = %d\n", p);
196 printf ("expected 1\n");
197 printf ("got %d\n", inex);
198 exit (1);
199 }
200 mpfr_clear (agm);
201 }
202
203 mpfr_clear (a);
204 mpfr_clear (b);
205 }
206
207 static void
208 check_special (void)
209 {
210 mpfr_t x, y, m;
211
212 mpfr_init2 (x, 123L);
213 mpfr_init2 (y, 123L);
214 mpfr_init2 (m, 123L);
215
216 /* agm(1,nan) is NaN */
217 mpfr_set_ui (x, 1L, MPFR_RNDN);
218 mpfr_set_nan (y);
219 mpfr_agm (m, x, y, MPFR_RNDN);
220 MPFR_ASSERTN (mpfr_nan_p (m));
221 mpfr_agm (m, y, x, MPFR_RNDN);
222 MPFR_ASSERTN (mpfr_nan_p (m));
223
224 /* agm(1,+inf) == +inf */
225 mpfr_set_ui (x, 1L, MPFR_RNDN);
226 mpfr_set_inf (y, 1);
227 mpfr_agm (m, x, y, MPFR_RNDN);
228 MPFR_ASSERTN (mpfr_inf_p (m));
229 MPFR_ASSERTN (mpfr_sgn (m) > 0);
230 mpfr_agm (m, y, x, MPFR_RNDN);
231 MPFR_ASSERTN (mpfr_inf_p (m));
232 MPFR_ASSERTN (mpfr_sgn (m) > 0);
233
234 /* agm(+inf,+inf) == +inf */
235 mpfr_set_inf (x, 1);
236 mpfr_set_inf (y, 1);
237 mpfr_agm (m, x, y, MPFR_RNDN);
238 MPFR_ASSERTN (mpfr_inf_p (m));
239 MPFR_ASSERTN (mpfr_sgn (m) > 0);
240
241 /* agm(-inf,+inf) is NaN */
242 mpfr_set_inf (x, -1);
243 mpfr_set_inf (y, 1);
244 mpfr_agm (m, x, y, MPFR_RNDN);
245 MPFR_ASSERTN (mpfr_nan_p (m));
246 mpfr_agm (m, y, x, MPFR_RNDN);
247 MPFR_ASSERTN (mpfr_nan_p (m));
248
249 /* agm(+0,+inf) is NaN */
250 mpfr_set_ui (x, 0, MPFR_RNDN);
251 mpfr_set_inf (y, 1);
252 mpfr_agm (m, x, y, MPFR_RNDN);
253 MPFR_ASSERTN (mpfr_nan_p (m));
254 mpfr_agm (m, y, x, MPFR_RNDN);
255 MPFR_ASSERTN (mpfr_nan_p (m));
256
257 /* agm(-0,+inf) is NaN */
258 mpfr_set_ui (x, 0, MPFR_RNDN);
259 mpfr_neg (x, x, MPFR_RNDN);
260 mpfr_set_inf (y, 1);
261 mpfr_agm (m, x, y, MPFR_RNDN);
262 MPFR_ASSERTN (mpfr_nan_p (m));
263 mpfr_agm (m, y, x, MPFR_RNDN);
264 MPFR_ASSERTN (mpfr_nan_p (m));
265
266 /* agm(+0,-inf) is NaN */
267 mpfr_set_ui (x, 0, MPFR_RNDN);
268 mpfr_set_inf (y, -1);
269 mpfr_agm (m, x, y, MPFR_RNDN);
270 MPFR_ASSERTN (mpfr_nan_p (m));
271 mpfr_agm (m, y, x, MPFR_RNDN);
272 MPFR_ASSERTN (mpfr_nan_p (m));
273
274 /* agm(-0,-inf) is NaN */
275 mpfr_set_ui (x, 0, MPFR_RNDN);
276 mpfr_neg (x, x, MPFR_RNDN);
277 mpfr_set_inf (y, -1);
278 mpfr_agm (m, x, y, MPFR_RNDN);
279 MPFR_ASSERTN (mpfr_nan_p (m));
280 mpfr_agm (m, y, x, MPFR_RNDN);
281 MPFR_ASSERTN (mpfr_nan_p (m));
282
283 /* agm(+0,1) == +0 */
284 mpfr_set_ui (x, 0, MPFR_RNDN);
285 mpfr_set_ui (y, 1, MPFR_RNDN);
286 mpfr_agm (m, x, y, MPFR_RNDN);
287 MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));
288 mpfr_agm (m, y, x, MPFR_RNDN);
289 MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));
290
291 /* agm(-0,1) == +0 */
292 mpfr_set_ui (x, 0, MPFR_RNDN);
293 mpfr_neg (x, x, MPFR_RNDN);
294 mpfr_set_ui (y, 1, MPFR_RNDN);
295 mpfr_agm (m, x, y, MPFR_RNDN);
296 MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));
297 mpfr_agm (m, y, x, MPFR_RNDN);
298 MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));
299
300 /* agm(+0,-1) == +0 */
301 mpfr_set_ui (x, 0, MPFR_RNDN);
302 mpfr_set_si (y, -1, MPFR_RNDN);
303 mpfr_agm (m, x, y, MPFR_RNDN);
304 MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));
305 mpfr_agm (m, y, x, MPFR_RNDN);
306 MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));
307
308 /* agm(-0,-1) == +0 */
309 mpfr_set_ui (x, 0, MPFR_RNDN);
310 mpfr_neg (x, x, MPFR_RNDN);
311 mpfr_set_si (y, -1, MPFR_RNDN);
312 mpfr_agm (m, x, y, MPFR_RNDN);
313 MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));
314 mpfr_agm (m, y, x, MPFR_RNDN);
315 MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));
316
317 /* agm(-0,+0) == +0 */
318 mpfr_set_ui (x, 0, MPFR_RNDN);
319 mpfr_neg (x, x, MPFR_RNDN);
320 mpfr_set_ui (y, 0, MPFR_RNDN);
321 mpfr_agm (m, x, y, MPFR_RNDN);
322 MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));
323 mpfr_agm (m, y, x, MPFR_RNDN);
324 MPFR_ASSERTN (MPFR_IS_ZERO (m) && MPFR_IS_POS(m));
325
326 /* agm(1,1) == 1 */
327 mpfr_set_ui (x, 1, MPFR_RNDN);
328 mpfr_set_ui (y, 1, MPFR_RNDN);
329 mpfr_agm (m, x, y, MPFR_RNDN);
330 MPFR_ASSERTN (mpfr_cmp_ui (m, 1) == 0);
331
332 /* agm(-1,-2) is NaN */
333 mpfr_set_si (x, -1, MPFR_RNDN);
334 mpfr_set_si (y, -2, MPFR_RNDN);
335 mpfr_agm (m, x, y, MPFR_RNDN);
336 MPFR_ASSERTN (mpfr_nan_p (m));
337 mpfr_agm (m, y, x, MPFR_RNDN);
338 MPFR_ASSERTN (mpfr_nan_p (m));
339
340 mpfr_clear (x);
341 mpfr_clear (y);
342 mpfr_clear (m);
343 }
344
345 #define TEST_FUNCTION mpfr_agm
346 #define TWO_ARGS
347 #define TEST_RANDOM_POS 4
348 #define TEST_RANDOM_POS2 4
349 #include "tgeneric.c"
350
351 int
352 main (int argc, char* argv[])
353 {
354 tests_start_mpfr ();
355
356 check_special ();
357 check_large ();
358 check_eq ();
359 check4 ("2.0", "1.0", MPFR_RNDN, "1.456791031046906869", -1);
360 check4 ("6.0", "4.0", MPFR_RNDN, "4.949360872472608925", 1);
361 check4 ("62.0", "61.0", MPFR_RNDN, "61.498983718845075902", -1);
362 check4 ("0.5", "1.0", MPFR_RNDN, "0.72839551552345343459", -1);
363 check4 ("1.0", "2.0", MPFR_RNDN, "1.456791031046906869", -1);
364 check4 ("234375765.0", "234375000.0", MPFR_RNDN, "234375382.49984394025", 1);
365 check4 ("8.0", "1.0", MPFR_RNDU, "3.61575617759736274873", 1);
366 check4 ("1.0", "44.0", MPFR_RNDU, "13.3658354512981243907", 1);
367 check4 ("1.0", "3.7252902984619140625e-9", MPFR_RNDU,
368 "0.07553933569711989657765", 1);
369 test_generic (MPFR_PREC_MIN, 300, 17);
370
371 tests_end_mpfr ();
372 return 0;
373 }