1 /* Test file for mpfr_hypot.
2
3 Copyright 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 /* Non-zero when extended exponent range */
26 static int ext = 0;
27
28 static void
29 special (void)
30 {
31 mpfr_t x, y, z;
32
33 mpfr_init (x);
34 mpfr_init (y);
35 mpfr_init (z);
36
37 mpfr_set_nan (x);
38 mpfr_set_ui (y, 0, MPFR_RNDN);
39 mpfr_hypot (z, x, x, MPFR_RNDN);
40 MPFR_ASSERTN(mpfr_nan_p (z));
41 mpfr_hypot (z, x, y, MPFR_RNDN);
42 MPFR_ASSERTN(mpfr_nan_p (z));
43 mpfr_hypot (z, y, x, MPFR_RNDN);
44 MPFR_ASSERTN(mpfr_nan_p (z));
45
46 mpfr_set_inf (x, 1);
47 mpfr_set_inf (y, -1);
48 mpfr_hypot (z, x, y, MPFR_RNDN);
49 MPFR_ASSERTN(mpfr_inf_p (z) && mpfr_sgn (z) > 0);
50
51 mpfr_set_inf (x, -1);
52 mpfr_set_nan (y);
53 mpfr_hypot (z, x, y, MPFR_RNDN);
54 MPFR_ASSERTN(mpfr_inf_p (z) && mpfr_sgn (z) > 0);
55
56 mpfr_set_nan (x);
57 mpfr_set_inf (y, -1);
58 mpfr_hypot (z, x, y, MPFR_RNDN);
59 MPFR_ASSERTN(mpfr_inf_p (z) && mpfr_sgn (z) > 0);
60
61 mpfr_clear (x);
62 mpfr_clear (y);
63 mpfr_clear (z);
64 }
65
66 static void
67 test_large (void)
68 {
69 mpfr_t x, y, z, t;
70
71 mpfr_init (x);
72 mpfr_init (y);
73 mpfr_init (z);
74 mpfr_init (t);
75
76 mpfr_set_ui (x, 21, MPFR_RNDN);
77 mpfr_set_ui (y, 28, MPFR_RNDN);
78 mpfr_set_ui (z, 35, MPFR_RNDN);
79
80 mpfr_mul_2ui (x, x, MPFR_EMAX_DEFAULT-6, MPFR_RNDN);
81 mpfr_mul_2ui (y, y, MPFR_EMAX_DEFAULT-6, MPFR_RNDN);
82 mpfr_mul_2ui (z, z, MPFR_EMAX_DEFAULT-6, MPFR_RNDN);
83
84 mpfr_hypot (t, x, y, MPFR_RNDN);
85 if (mpfr_cmp (z, t))
86 {
87 printf ("Error in test_large: got\n");
88 mpfr_out_str (stdout, 2, 0, t, MPFR_RNDN);
89 printf ("\ninstead of\n");
90 mpfr_out_str (stdout, 2, 0, z, MPFR_RNDN);
91 printf ("\n");
92 exit (1);
93 }
94
95 mpfr_set_prec (x, 53);
96 mpfr_set_prec (t, 53);
97 mpfr_set_prec (y, 53);
98 mpfr_set_str_binary (x, "0.11101100011110000011101000010101010011001101000001100E-1021");
99 mpfr_set_str_binary (y, "0.11111001010011000001110110001101011100001000010010100E-1021");
100 mpfr_hypot (t, x, y, MPFR_RNDN);
101 mpfr_set_str_binary (z, "0.101010111100110111101110111110100110010011001010111E-1020");
102 if (mpfr_cmp (z, t))
103 {
104 printf ("Error in test_large: got\n");
105 mpfr_out_str (stdout, 2, 0, t, MPFR_RNDN);
106 printf ("\ninstead of\n");
107 mpfr_out_str (stdout, 2, 0, z, MPFR_RNDN);
108 printf ("\n");
109 exit (1);
110 }
111
112 mpfr_set_prec (x, 240);
113 mpfr_set_prec (y, 22);
114 mpfr_set_prec (z, 2);
115 mpfr_set_prec (t, 2);
116 mpfr_set_str_binary (x, "0.100111011010010010110100000100000001100010011100110101101111111101011110111011011101010110100101111000111100010100110000100101011110111011100110100110100101110101101100011000001100000001111101110100100100011011011010110111100110010101000111e-7");
117 mpfr_set_str_binary (y, "0.1111000010000011000111E-10");
118 mpfr_hypot (t, x, y, MPFR_RNDN);
119 mpfr_set_str_binary (z, "0.11E-7");
120 if (mpfr_cmp (z, t))
121 {
122 printf ("Error in test_large: got\n");
123 mpfr_out_str (stdout, 2, 0, t, MPFR_RNDN);
124 printf ("\ninstead of\n");
125 mpfr_out_str (stdout, 2, 0, z, MPFR_RNDN);
126 printf ("\n");
127 exit (1);
128 }
129
130 mpfr_clear (x);
131 mpfr_clear (y);
132 mpfr_clear (z);
133 mpfr_clear (t);
134 }
135
136 static void
137 test_small (void)
138 {
139 mpfr_t x, y, z1, z2;
140 int inex1, inex2;
141 unsigned int flags;
142
143 /* Test hypot(x,x) with x = 2^(emin-1). Result is x * sqrt(2). */
144 mpfr_inits2 (8, x, y, z1, z2, (mpfr_ptr) 0);
145 mpfr_set_si_2exp (x, 1, mpfr_get_emin () - 1, MPFR_RNDN);
146 mpfr_set_si_2exp (y, 1, mpfr_get_emin () - 1, MPFR_RNDN);
147 mpfr_set_ui (z1, 2, MPFR_RNDN);
148 inex1 = mpfr_sqrt (z1, z1, MPFR_RNDN);
149 inex2 = mpfr_mul (z1, z1, x, MPFR_RNDN);
150 MPFR_ASSERTN (inex2 == 0);
151 mpfr_clear_flags ();
152 inex2 = mpfr_hypot (z2, x, y, MPFR_RNDN);
153 flags = __gmpfr_flags;
154 if (mpfr_cmp (z1, z2) != 0)
155 {
156 printf ("Error in test_small%s\nExpected ",
157 ext ? ", extended exponent range" : "");
158 mpfr_out_str (stdout, 2, 0, z1, MPFR_RNDN);
159 printf ("\nGot ");
160 mpfr_out_str (stdout, 2, 0, z2, MPFR_RNDN);
161 printf ("\n");
162 exit (1);
163 }
164 if (! SAME_SIGN (inex1, inex2))
165 {
166 printf ("Bad ternary value in test_small%s\nExpected %d, got %d\n",
167 ext ? ", extended exponent range" : "", inex1, inex2);
168 exit (1);
169 }
170 if (flags != MPFR_FLAGS_INEXACT)
171 {
172 printf ("Bad flags in test_small%s\nExpected %u, got %u\n",
173 ext ? ", extended exponent range" : "",
174 (unsigned int) MPFR_FLAGS_INEXACT, flags);
175 exit (1);
176 }
177 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0);
178 }
179
180 static void
181 test_large_small (void)
182 {
183 mpfr_t x, y, z;
184 int inexact, inex2, r;
185
186 mpfr_init2 (x, 3);
187 mpfr_init2 (y, 2);
188 mpfr_init2 (z, 2);
189
190 mpfr_set_ui_2exp (x, 1, mpfr_get_emax () / 2, MPFR_RNDN);
191 mpfr_set_ui_2exp (y, 1, -1, MPFR_RNDN);
192 inexact = mpfr_hypot (z, x, y, MPFR_RNDN);
193 if (inexact >= 0 || mpfr_cmp (x, z))
194 {
195 printf ("Error 1 in test_large_small%s\n",
196 ext ? ", extended exponent range" : "");
197 exit (1);
198 }
199
200 mpfr_mul_ui (x, x, 5, MPFR_RNDN);
201 inexact = mpfr_hypot (z, x, y, MPFR_RNDN);
202 if (mpfr_cmp (x, z) >= 0)
203 {
204 printf ("Error 2 in test_large_small%s\n",
205 ext ? ", extended exponent range" : "");
206 printf ("x = ");
207 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
208 printf ("\n");
209 printf ("y = ");
210 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
211 printf ("\n");
212 printf ("z = ");
213 mpfr_out_str (stdout, 2, 0, z, MPFR_RNDN);
214 printf (" (in precision 2) instead of\n ");
215 mpfr_out_str (stdout, 2, 2, x, MPFR_RNDU);
216 printf ("\n");
217 exit (1);
218 }
219
220 RND_LOOP(r)
221 {
222 mpfr_set_ui_2exp (x, 1, mpfr_get_emax () - 1, MPFR_RNDN);
223 mpfr_set_ui_2exp (y, 1, mpfr_get_emin (), MPFR_RNDN);
224 inexact = mpfr_hypot (z, x, y, (mpfr_rnd_t) r);
225 inex2 = mpfr_add_ui (y, x, 1, (mpfr_rnd_t) r);
226 if (! mpfr_equal_p (y, z) || ! SAME_SIGN (inexact, inex2))
227 {
228 printf ("Error 3 in test_large_small, %s%s\n",
229 mpfr_print_rnd_mode ((mpfr_rnd_t) r),
230 ext ? ", extended exponent range" : "");
231 printf ("Expected ");
232 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
233 printf (", inex = %d\n", inex2);
234 printf ("Got ");
235 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
236 printf (", inex = %d\n", inexact);
237 exit (1);
238 }
239 }
240
241 mpfr_clear (x);
242 mpfr_clear (y);
243 mpfr_clear (z);
244 }
245
246 static void
247 check_overflow (void)
248 {
249 mpfr_t x, y;
250 int inex, r;
251
252 mpfr_inits2 (8, x, y, (mpfr_ptr) 0);
253 mpfr_set_ui (x, 1, MPFR_RNDN);
254 mpfr_setmax (x, mpfr_get_emax ());
255
256 RND_LOOP(r)
257 {
258 mpfr_clear_overflow ();
259 inex = mpfr_hypot (y, x, x, (mpfr_rnd_t) r);
260 if (!mpfr_overflow_p ())
261 {
262 printf ("No overflow in check_overflow for %s%s\n",
263 mpfr_print_rnd_mode ((mpfr_rnd_t) r),
264 ext ? ", extended exponent range" : "");
265 exit (1);
266 }
267 MPFR_ASSERTN (MPFR_IS_POS (y));
268 if (r == MPFR_RNDZ || r == MPFR_RNDD)
269 {
270 MPFR_ASSERTN (inex < 0);
271 MPFR_ASSERTN (!mpfr_inf_p (y));
272 mpfr_nexttoinf (y);
273 }
274 else
275 {
276 MPFR_ASSERTN (inex > 0);
277 }
278 MPFR_ASSERTN (mpfr_inf_p (y));
279 }
280
281 mpfr_clears (x, y, (mpfr_ptr) 0);
282 }
283
284 #define TWO_ARGS
285 #define TEST_FUNCTION mpfr_hypot
286 #include "tgeneric.c"
287
288 static void
289 alltst (void)
290 {
291 mpfr_exp_t emin, emax;
292
293 ext = 0;
294 test_small ();
295 test_large_small ();
296 check_overflow ();
297
298 emin = mpfr_get_emin ();
299 emax = mpfr_get_emax ();
300 set_emin (MPFR_EMIN_MIN);
301 set_emax (MPFR_EMAX_MAX);
302 if (mpfr_get_emin () != emin || mpfr_get_emax () != emax)
303 {
304 ext = 1;
305 test_small ();
306 test_large_small ();
307 check_overflow ();
308 set_emin (emin);
309 set_emax (emax);
310 }
311 }
312
313 /* Test failing with GMP_CHECK_RANDOMIZE=1513841234 (overflow flag not set).
314 The bug was in fact in mpfr_nexttoinf which didn't set the overflow flag. */
315 static void
316 bug20171221 (void)
317 {
318 mpfr_t x, u, y;
319 int inex;
320 mpfr_exp_t emax;
321
322 mpfr_init2 (x, 12);
323 mpfr_init2 (u, 12);
324 mpfr_init2 (y, 11);
325 mpfr_set_str_binary (x, "0.111111111110E0");
326 mpfr_set_str_binary (u, "0.111011110100E-177");
327 emax = mpfr_get_emax ();
328 set_emax (0);
329 mpfr_clear_flags ();
330 inex = mpfr_hypot (y, x, u, MPFR_RNDU);
331 MPFR_ASSERTN(mpfr_inf_p (y) && mpfr_sgn (y) > 0);
332 MPFR_ASSERTN(inex > 0);
333 MPFR_ASSERTN(mpfr_inexflag_p ());
334 MPFR_ASSERTN(mpfr_overflow_p ());
335 set_emax (emax);
336 mpfr_clear (x);
337 mpfr_clear (u);
338 mpfr_clear (y);
339 }
340
341 /* check overflow for x=0xf.ffffffffffff8p+1020 and u=0xf.fffffp+124
342 with rounding upwards and emax=1024 (as in binary64) */
343 static void
344 test_overflow (void)
345 {
346 mpfr_t x, u, y;
347 int inex;
348 mpfr_exp_t emax;
349
350 mpfr_init2 (x, 53);
351 mpfr_init2 (u, 53);
352 mpfr_init2 (y, 53);
353 mpfr_set_str (x, "0xf.ffffffffffff8p+1020", 16, MPFR_RNDN);
354 mpfr_set_str (u, "0xf.fffffp+124", 16, MPFR_RNDN);
355 emax = mpfr_get_emax ();
356 set_emax (1024);
357 mpfr_clear_flags ();
358 inex = mpfr_hypot (y, x, u, MPFR_RNDU);
359 MPFR_ASSERTN(mpfr_inf_p (y) && mpfr_sgn (y) > 0);
360 MPFR_ASSERTN(inex > 0);
361 MPFR_ASSERTN(mpfr_inexflag_p ());
362 MPFR_ASSERTN(mpfr_overflow_p ());
363 set_emax (emax);
364 mpfr_clear (x);
365 mpfr_clear (u);
366 mpfr_clear (y);
367 }
368
369 int
370 main (int argc, char *argv[])
371 {
372 tests_start_mpfr ();
373
374 bug20171221 ();
375 special ();
376
377 test_large ();
378 alltst ();
379 test_overflow ();
380
381 test_generic (MPFR_PREC_MIN, 100, 10);
382
383 tests_end_mpfr ();
384 return 0;
385 }