xref: /netbsd-src/external/lgpl3/mpfr/dist/tests/taway.c (revision ba125506a622fe649968631a56eba5d42ff57863)
1 /* Test file for round away.
2 
3 Copyright 2000-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 DISP(s,t)                                       \
26   do                                                    \
27     {                                                   \
28       printf (s);                                       \
29       mpfr_out_str (stdout, 2, 0, t, MPFR_RNDN);        \
30     }                                                   \
31   while (0)
32 
33 #define DISP2(s,t) do { DISP(s,t); putchar ('\n'); } while (0)
34 
35 #define SPECIAL_MAX 12
36 
37 static void
set_special(mpfr_ptr x,unsigned int select)38 set_special (mpfr_ptr x, unsigned int select)
39 {
40   MPFR_ASSERTN (select < SPECIAL_MAX);
41   switch (select)
42     {
43     case 0:
44       MPFR_SET_NAN (x);
45       break;
46     case 1:
47       MPFR_SET_INF (x);
48       MPFR_SET_POS (x);
49       break;
50     case 2:
51       MPFR_SET_INF (x);
52       MPFR_SET_NEG (x);
53       break;
54     case 3:
55       MPFR_SET_ZERO (x);
56       MPFR_SET_POS  (x);
57       break;
58     case 4:
59       MPFR_SET_ZERO (x);
60       MPFR_SET_NEG  (x);
61       break;
62     case 5:
63       mpfr_set_str_binary (x, "1");
64       break;
65     case 6:
66       mpfr_set_str_binary (x, "-1");
67       break;
68     case 7:
69       mpfr_set_str_binary (x, "1e-1");
70       break;
71     case 8:
72       mpfr_set_str_binary (x, "1e+1");
73       break;
74     case 9:
75       mpfr_const_pi (x, MPFR_RNDN);
76       break;
77     case 10:
78       mpfr_const_pi (x, MPFR_RNDN);
79       MPFR_SET_EXP (x, MPFR_GET_EXP (x)-1);
80       break;
81     default:
82       mpfr_urandomb (x, RANDS);
83       break;
84     }
85 }
86 
87 /* same as mpfr_cmp, but returns 0 for both NaN's */
88 static int
mpfr_compare(mpfr_srcptr a,mpfr_srcptr b)89 mpfr_compare (mpfr_srcptr a, mpfr_srcptr b)
90 {
91   return (MPFR_IS_NAN(a)) ? !MPFR_IS_NAN(b) :
92     (MPFR_IS_NAN(b) || mpfr_cmp(a, b));
93 }
94 
95 static void
test3(int (* testfunc)(mpfr_ptr,mpfr_srcptr,mpfr_srcptr,mpfr_rnd_t),const char * foo)96 test3 (int (*testfunc)(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t),
97        const char *foo)
98 {
99   mpfr_t ref1, ref2, ref3;
100   mpfr_t res1;
101   mpfr_prec_t p1, p2, p3;
102   int i, inexa, inexd;
103   mpfr_rnd_t r;
104 
105   p1 = (randlimb () % 200) + MPFR_PREC_MIN;
106   p2 = (randlimb () % 200) + MPFR_PREC_MIN;
107   p3 = (randlimb () % 200) + MPFR_PREC_MIN;
108 
109   mpfr_init2 (ref1, p1);
110   mpfr_init2 (ref2, p2);
111   mpfr_init2 (ref3, p3);
112   mpfr_init2 (res1, p1);
113 
114   /* for each variable, consider each of the following 6 possibilities:
115      NaN, +Infinity, -Infinity, +0, -0 or a random number */
116   for (i = 0; i < SPECIAL_MAX * SPECIAL_MAX; i++)
117     {
118       set_special (ref2, i%SPECIAL_MAX);
119       set_special (ref3, i/SPECIAL_MAX);
120 
121       inexa = testfunc (res1, ref2, ref3, MPFR_RNDA);
122       r = MPFR_IS_POS (res1) ? MPFR_RNDU : MPFR_RNDD;
123       inexd = testfunc (ref1, ref2, ref3, r);
124 
125       if (mpfr_compare (res1, ref1) || inexa != inexd)
126         {
127           printf ("Error with RNDA for %s with ", foo);
128           DISP("x=",ref2); DISP2(", y=",ref3);
129           printf ("inexa=%d inexd=%d\n", inexa, inexd);
130           printf ("expected "); mpfr_dump (ref1);
131           printf ("got      "); mpfr_dump (res1);
132           exit (1);
133         }
134     }
135 
136   mpfr_clear (ref1);
137   mpfr_clear (ref2);
138   mpfr_clear (ref3);
139   mpfr_clear (res1);
140 }
141 
142 static void
test4(int (* testfunc)(mpfr_ptr,mpfr_srcptr,mpfr_srcptr,mpfr_srcptr,mpfr_rnd_t),const char * foo)143 test4 (int (*testfunc)(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_srcptr,
144                        mpfr_rnd_t), const char *foo)
145 {
146   mpfr_t ref, op1, op2, op3;
147   mpfr_prec_t pout, p1, p2, p3;
148   mpfr_t res;
149   int i, j, k, inexa, inexd;
150   mpfr_rnd_t r;
151 
152   pout = (randlimb () % 200) + MPFR_PREC_MIN;
153   p1 = (randlimb () % 200) + MPFR_PREC_MIN;
154   p2 = (randlimb () % 200) + MPFR_PREC_MIN;
155   p3 = (randlimb () % 200) + MPFR_PREC_MIN;
156 
157   mpfr_init2 (ref, pout);
158   mpfr_init2 (res, pout);
159   mpfr_init2 (op1, p1);
160   mpfr_init2 (op2, p2);
161   mpfr_init2 (op3, p3);
162 
163   /* for each variable, consider each of the following 6 possibilities:
164      NaN, +Infinity, -Infinity, +0, -0 or a random number */
165 
166   for (i = 0; i < SPECIAL_MAX; i++)
167     {
168       set_special (op1, i);
169       for (j = 0; j < SPECIAL_MAX; j++)
170         {
171           set_special (op2, j);
172           for (k = 0; k < SPECIAL_MAX; k++)
173             {
174               set_special (op3, k);
175 
176               inexa = testfunc (res, op1, op2, op3, MPFR_RNDA);
177               r = MPFR_IS_POS (res) ? MPFR_RNDU : MPFR_RNDD;
178               inexd = testfunc (ref, op1, op2, op3, r);
179 
180               if (mpfr_compare (res, ref) || inexa != inexd)
181                 {
182                   printf ("Error with RNDA for %s with ", foo);
183                   DISP("a=", op1); DISP(", b=", op2); DISP2(", c=", op3);
184                   printf ("inexa=%d inexd=%d\n", inexa, inexd);
185                   DISP("expected ", ref); DISP2(", got ", res);
186                   exit (1);
187                 }
188             }
189         }
190     }
191 
192   mpfr_clear (ref);
193   mpfr_clear (op1);
194   mpfr_clear (op2);
195   mpfr_clear (op3);
196   mpfr_clear (res);
197 }
198 
199 static void
test2ui(int (* testfunc)(mpfr_ptr,mpfr_srcptr,unsigned long int,mpfr_rnd_t),const char * foo)200 test2ui (int (*testfunc)(mpfr_ptr, mpfr_srcptr, unsigned long int, mpfr_rnd_t),
201          const char *foo)
202 {
203   mpfr_t ref1, ref2;
204   unsigned int ref3;
205   mpfr_t res1;
206   mpfr_prec_t p1, p2;
207   int i, inexa, inexd;
208   mpfr_rnd_t r;
209 
210   p1 = (randlimb () % 200) + MPFR_PREC_MIN;
211   p2 = (randlimb () % 200) + MPFR_PREC_MIN;
212 
213   mpfr_init2 (ref1, p1);
214   mpfr_init2 (ref2, p2);
215   mpfr_init2 (res1, p1);
216 
217   /* ref2 can be NaN, +Inf, -Inf, +0, -0 or any number
218      ref3 can be 0 or any number */
219   for (i = 0; i < SPECIAL_MAX * 2; i++)
220     {
221       set_special (ref2, i % SPECIAL_MAX);
222       ref3 = i / SPECIAL_MAX == 0 ? 0 : randlimb ();
223 
224       inexa = testfunc (res1, ref2, ref3, MPFR_RNDA);
225       r = MPFR_IS_POS (res1) ? MPFR_RNDU : MPFR_RNDD;
226       inexd = testfunc (ref1, ref2, ref3, r);
227 
228       if (mpfr_compare (res1, ref1) || inexa != inexd)
229         {
230           printf ("Error with RNDA for %s for c=%u\n", foo, ref3);
231           DISP2("a=",ref2);
232           printf ("inexa=%d inexd=%d\n", inexa, inexd);
233           printf ("expected "); mpfr_dump (ref1);
234           printf ("got      "); mpfr_dump (res1);
235           exit (1);
236         }
237     }
238 
239   mpfr_clear (ref1);
240   mpfr_clear (ref2);
241   mpfr_clear (res1);
242 }
243 
244 static void
testui2(int (* testfunc)(mpfr_ptr,unsigned long int,mpfr_srcptr,mpfr_rnd_t),const char * foo)245 testui2 (int (*testfunc)(mpfr_ptr, unsigned long int, mpfr_srcptr, mpfr_rnd_t),
246          const char *foo)
247 {
248   mpfr_t ref1, ref3;
249   unsigned int ref2;
250   mpfr_t res1;
251   mpfr_prec_t p1, p3;
252   int i, inexa, inexd;
253   mpfr_rnd_t r;
254 
255   p1 = (randlimb () % 200) + MPFR_PREC_MIN;
256   p3 = (randlimb () % 200) + MPFR_PREC_MIN;
257 
258   mpfr_init2 (ref1, p1);
259   mpfr_init2 (ref3, p3);
260   mpfr_init2 (res1, p1);
261 
262   for (i = 0; i < SPECIAL_MAX * 2; i++)
263     {
264       set_special (ref3, i % SPECIAL_MAX);
265       ref2 = i / SPECIAL_MAX == 0 ? 0 : randlimb ();
266 
267       inexa = testfunc (res1, ref2, ref3, MPFR_RNDA);
268       r = MPFR_IS_POS (res1) ? MPFR_RNDU : MPFR_RNDD;
269       inexd = testfunc (ref1, ref2, ref3, r);
270 
271       if (mpfr_compare (res1, ref1) || inexa != inexd)
272         {
273           printf ("Error with RNDA for %s for b=%u\n", foo, ref2);
274           DISP2("a=", ref3);
275           printf ("inexa=%d inexd=%d\n", inexa, inexd);
276           DISP("expected ", ref1); DISP2(", got ", res1);
277           exit (1);
278         }
279     }
280 
281   mpfr_clear (ref1);
282   mpfr_clear (ref3);
283   mpfr_clear (res1);
284 }
285 
286 /* foo(mpfr_ptr, mpfr_srcptr, mp_rndt) */
287 static void
test2(int (* testfunc)(mpfr_ptr,mpfr_srcptr,mpfr_rnd_t),const char * foo)288 test2 (int (*testfunc)(mpfr_ptr, mpfr_srcptr, mpfr_rnd_t), const char *foo)
289 {
290   mpfr_t ref1, ref2;
291   mpfr_t res1;
292   mpfr_prec_t p1, p2;
293   int i, inexa, inexd;
294   mpfr_rnd_t r;
295 
296   p1 = (randlimb () % 200) + MPFR_PREC_MIN;
297   p2 = (randlimb () % 200) + MPFR_PREC_MIN;
298 
299   mpfr_init2 (ref1, p1);
300   mpfr_init2 (ref2, p2);
301   mpfr_init2 (res1, p1);
302 
303   for (i = 0; i < SPECIAL_MAX; i++)
304     {
305       set_special (ref2, i);
306 
307       /* first round to away */
308       inexa = testfunc (res1, ref2, MPFR_RNDA);
309 
310       r = MPFR_IS_POS (res1) ? MPFR_RNDU : MPFR_RNDD;
311       inexd = testfunc (ref1, ref2, r);
312       if (mpfr_compare (res1, ref1) || inexa != inexd)
313         {
314           printf ("Error with RNDA for %s with ", foo);
315           DISP2("x=", ref2);
316           printf ("inexa=%d inexd=%d\n", inexa, inexd);
317           DISP("expected ", ref1); DISP2(", got ", res1);
318           exit (1);
319         }
320     }
321 
322   mpfr_clear (ref1);
323   mpfr_clear (ref2);
324   mpfr_clear (res1);
325 }
326 
327 /* one operand, two results, like mpfr_sin_cos */
328 static void
test3a(int (* testfunc)(mpfr_ptr,mpfr_ptr,mpfr_srcptr,mpfr_rnd_t),const char * foo)329 test3a (int (*testfunc)(mpfr_ptr, mpfr_ptr, mpfr_srcptr, mpfr_rnd_t),
330         const char *foo)
331 {
332   mpfr_t ref1, ref2, ref3;
333   mpfr_t res1, res2;
334   mpfr_prec_t p1, p2, p3;
335   int i, inexa, inexd;
336   mpfr_rnd_t r;
337 
338   p1 = (randlimb () % 200) + MPFR_PREC_MIN;
339   p2 = (randlimb () % 200) + MPFR_PREC_MIN;
340   p3 = (randlimb () % 200) + MPFR_PREC_MIN;
341 
342   mpfr_init2 (ref1, p1);
343   mpfr_init2 (ref2, p2);
344   mpfr_init2 (ref3, p3);
345   mpfr_init2 (res1, p1);
346   mpfr_init2 (res2, p2);
347 
348   for (i = 0; i < SPECIAL_MAX; i++)
349     {
350       set_special (ref3, i);
351 
352       inexa = testfunc (res1, res2, ref3, MPFR_RNDA);
353 
354       /* first check wrt the first operand */
355       r = MPFR_IS_POS (res1) ? MPFR_RNDU : MPFR_RNDD;
356       inexd = testfunc (ref1, ref2, ref3, r);
357       /* the low 2 bits of the inexact flag concern the 1st operand */
358       if (mpfr_compare (res1, ref1) || (inexa & 3) != (inexd & 3))
359         {
360           printf ("Error with RNDA for %s (1st operand)\n", foo);
361           DISP2("a=",ref3);
362           DISP("expected ", ref1); printf ("\n");
363           DISP("got      ", res1); printf ("\n");
364           printf ("inexa=%d inexd=%d\n", inexa & 3, inexd & 3);
365           exit (1);
366         }
367 
368       /* now check wrt the second operand */
369       r = MPFR_IS_POS (res2) ? MPFR_RNDU : MPFR_RNDD;
370       inexd = testfunc (ref1, ref2, ref3, r);
371       /* bits 2..3 of the inexact flag concern the 2nd operand */
372       if (mpfr_compare (res2, ref2) || (inexa >> 2) != (inexd >> 2))
373         {
374           printf ("Error with RNDA for %s (2nd operand)\n", foo);
375           DISP2("a=",ref3);
376           DISP("expected ", ref2); printf ("\n");
377           DISP("got      ", res2); printf ("\n");
378           printf ("inexa=%d inexd=%d\n", inexa >> 2, inexd >> 2);
379           exit (1);
380         }
381 
382     }
383 
384   mpfr_clear (ref1);
385   mpfr_clear (ref2);
386   mpfr_clear (ref3);
387   mpfr_clear (res1);
388   mpfr_clear (res2);
389 }
390 
391 int
main(void)392 main (void)
393 {
394   int N = 20;
395 
396   tests_start_mpfr ();
397 
398   while (N--)
399     {
400       /* no need to test mpfr_round, mpfr_ceil, mpfr_floor, mpfr_trunc since
401          they take no rounding mode */
402 
403       test2ui (mpfr_add_ui, "mpfr_add_ui");
404       test2ui (mpfr_div_2exp, "mpfr_div_2exp");
405       test2ui (mpfr_div_2ui, "mpfr_div_2ui");
406       test2ui (mpfr_div_ui, "mpfr_div_ui");
407       test2ui (mpfr_mul_2exp, "mpfr_mul_2exp");
408       test2ui (mpfr_mul_2ui, "mpfr_mul_2ui");
409       test2ui (mpfr_mul_ui, "mpfr_mul_ui");
410       test2ui (mpfr_pow_ui, "mpfr_pow_ui");
411       test2ui (mpfr_sub_ui, "mpfr_sub_ui");
412 
413       testui2 (mpfr_ui_div, "mpfr_ui_div");
414       testui2 (mpfr_ui_sub, "mpfr_ui_sub");
415       testui2 (mpfr_ui_pow, "mpfr_ui_pow");
416 
417       test2 (mpfr_sqr, "mpfr_sqr");
418       test2 (mpfr_sqrt, "mpfr_sqrt");
419       test2 (mpfr_abs, "mpfr_abs");
420       test2 (mpfr_neg, "mpfr_neg");
421 
422       test2 (mpfr_log, "mpfr_log");
423       test2 (mpfr_log2, "mpfr_log2");
424       test2 (mpfr_log10, "mpfr_log10");
425       test2 (mpfr_log1p, "mpfr_log1p");
426 
427       test2 (mpfr_exp, "mpfr_exp");
428       test2 (mpfr_exp2, "mpfr_exp2");
429       test2 (mpfr_exp10, "mpfr_exp10");
430       test2 (mpfr_expm1, "mpfr_expm1");
431       test2 (mpfr_eint, "mpfr_eint");
432 
433       test2 (mpfr_sinh, "mpfr_sinh");
434       test2 (mpfr_cosh, "mpfr_cosh");
435       test2 (mpfr_tanh, "mpfr_tanh");
436       test2 (mpfr_asinh, "mpfr_asinh");
437       test2 (mpfr_acosh, "mpfr_acosh");
438       test2 (mpfr_atanh, "mpfr_atanh");
439       test2 (mpfr_sech, "mpfr_sech");
440       test2 (mpfr_csch, "mpfr_csch");
441       test2 (mpfr_coth, "mpfr_coth");
442 
443       test2 (mpfr_asin, "mpfr_asin");
444       test2 (mpfr_acos, "mpfr_acos");
445       test2 (mpfr_atan, "mpfr_atan");
446       test2 (mpfr_cos, "mpfr_cos");
447       test2 (mpfr_sin, "mpfr_sin");
448       test2 (mpfr_tan, "mpfr_tan");
449       test2 (mpfr_sec, "mpfr_sec");
450       test2 (mpfr_csc, "mpfr_csc");
451       test2 (mpfr_cot, "mpfr_cot");
452 
453       test2 (mpfr_erf,  "mpfr_erf");
454       test2 (mpfr_erfc, "mpfr_erfc");
455       test2 (mpfr_j0,   "mpfr_j0");
456       test2 (mpfr_j1,   "mpfr_j1");
457       test2 (mpfr_y0,   "mpfr_y0");
458       test2 (mpfr_y1,   "mpfr_y1");
459       test2 (mpfr_zeta, "mpfr_zeta");
460       test2 (mpfr_gamma, "mpfr_gamma");
461       test2 (mpfr_lngamma, "mpfr_lngamma");
462 
463       test2 (mpfr_rint, "mpfr_rint");
464       test2 (mpfr_rint_ceil, "mpfr_rint_ceil");
465       test2 (mpfr_rint_floor, "mpfr_rint_floor");
466       test2 (mpfr_rint_round, "mpfr_rint_round");
467       test2 (mpfr_rint_trunc, "mpfr_rint_trunc");
468       test2 (mpfr_frac, "mpfr_frac");
469 
470       test3 (mpfr_add, "mpfr_add");
471       test3 (mpfr_sub, "mpfr_sub");
472       test3 (mpfr_mul, "mpfr_mul");
473       test3 (mpfr_div, "mpfr_div");
474 
475       test3 (mpfr_agm, "mpfr_agm");
476       test3 (mpfr_min, "mpfr_min");
477       test3 (mpfr_max, "mpfr_max");
478 
479       /* we don't test reldiff since it does not guarantee correct rounding,
480          thus we can get different results with RNDA and RNDU or RNDD. */
481       test3 (mpfr_dim, "mpfr_dim");
482 
483       test3 (mpfr_remainder, "mpfr_remainder");
484       test3 (mpfr_pow, "mpfr_pow");
485       test3 (mpfr_atan2, "mpfr_atan2");
486       test3 (mpfr_hypot, "mpfr_hypot");
487 
488       test3a (mpfr_sin_cos, "mpfr_sin_cos");
489 
490       test4 (mpfr_fma, "mpfr_fma");
491       test4 (mpfr_fms, "mpfr_fms");
492 
493       test2 (mpfr_li2, "mpfr_li2");
494       test2 (mpfr_rec_sqrt, "mpfr_rec_sqrt");
495       test3 (mpfr_fmod, "mpfr_fmod");
496       test3a (mpfr_modf, "mpfr_modf");
497       test3a (mpfr_sinh_cosh, "mpfr_sinh_cosh");
498 
499 #if MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)
500       test2 (mpfr_ai, "mpfr_ai");
501       test2 (mpfr_digamma, "mpfr_digamma");
502 #endif
503     }
504 
505   tests_end_mpfr ();
506   return 0;
507 }
508