xref: /netbsd-src/external/lgpl3/mpfr/dist/tests/tgamma_inc.c (revision f3cfa6f6ce31685c6c4a758bc430e69eb99f50a4)
1 /* Test file for mpfr_gamma_inc
2 
3 Copyright 2016-2018 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 http://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 TEST_FUNCTION mpfr_gamma_inc
26 #define TWO_ARGS
27 #define TEST_RANDOM_POS2 0 /* the 2nd argument is never negative */
28 #define TGENERIC_NOWARNING 1
29 #define TEST_RANDOM_EMAX 8
30 #define TEST_RANDOM_EMIN -32
31 #define REDUCE_EMAX TEST_RANDOM_EMAX
32 #define REDUCE_EMIN TEST_RANDOM_EMIN
33 #include "tgeneric.c"
34 
35 /* do k random tests at precision p */
36 static void
37 test_random (mpfr_prec_t p, int k)
38 {
39   mpfr_t a, x, y, z, t;
40 
41   mpfr_inits2 (p, a, x, y, z, (mpfr_ptr) 0);
42   mpfr_init2 (t, p + 20);
43   while (k--)
44     {
45       do mpfr_urandomb (a, RANDS); while (mpfr_zero_p (a));
46       if (randlimb () & 1)
47         mpfr_neg (a, a, MPFR_RNDN);
48       do mpfr_urandomb (x, RANDS); while (mpfr_zero_p (x));
49       mpfr_gamma_inc (y, a, x, MPFR_RNDN);
50       mpfr_gamma_inc (t, a, x, MPFR_RNDN);
51       if (mpfr_can_round (t, mpfr_get_prec (z), MPFR_RNDN, MPFR_RNDN, p))
52         {
53           mpfr_set (z, t, MPFR_RNDN);
54           if (mpfr_cmp (y, z))
55             {
56               printf ("mpfr_gamma_inc failed for a=");
57               mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
58               printf (" x=");
59               mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
60               printf ("\nexpected ");
61               mpfr_out_str (stdout, 10, 0, z, MPFR_RNDN);
62               printf ("\ngot      ");
63               mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN);
64               printf ("\n");
65               exit (1);
66             }
67         }
68     }
69   mpfr_clears (a, x, y, z, (mpfr_ptr) 0);
70   mpfr_clear (t);
71 }
72 
73 static void
74 specials (void)
75 {
76   mpfr_t a, x;
77   int inex;
78 
79   mpfr_init2 (a, 2);
80   mpfr_init2 (x, 2);
81 
82   /* check gamma_inc(2,0) = 1 is exact */
83   mpfr_set_ui (a, 2, MPFR_RNDN);
84   mpfr_set_ui (x, 0, MPFR_RNDN);
85   mpfr_clear_inexflag ();
86   inex = mpfr_gamma_inc (a, a, x, MPFR_RNDN);
87   if (mpfr_cmp_ui (a, 1))
88     {
89       printf ("Error for gamma_inc(2,0)\n");
90       printf ("expected 1\n");
91       printf ("got      ");
92       mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
93       printf ("\n");
94       exit (1);
95     }
96   if (inex != 0)
97     {
98       printf ("Wrong ternary value for gamma_inc(2,0)\n");
99       printf ("expected 0\n");
100       printf ("got      %d\n", inex);
101       exit (1);
102     }
103   if (mpfr_inexflag_p ())
104     {
105       printf ("Wrong inexact flag for gamma_inc(2,0)\n");
106       printf ("expected 0\n");
107       printf ("got      1\n");
108       exit (1);
109     }
110 
111   /* check gamma_inc(0,1) = 0.219383934395520 */
112   mpfr_set_ui (a, 0, MPFR_RNDN);
113   mpfr_set_ui (x, 1, MPFR_RNDN);
114   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
115   if (mpfr_cmp_ui_2exp (a, 1, -2))
116     {
117       printf ("Error for gamma_inc(0,1)\n");
118       printf ("expected 0.25\n");
119       printf ("got      ");
120       mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
121       printf ("\n");
122       exit (1);
123     }
124 
125   /* check gamma_inc(-1,1) = 0.148495506775922 */
126   mpfr_set_si (a, -1, MPFR_RNDN);
127   mpfr_set_ui (x, 1, MPFR_RNDN);
128   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
129   if (mpfr_cmp_ui_2exp (a, 1, -3))
130     {
131       printf ("Error for gamma_inc(-1,1)\n");
132       printf ("expected 0.125\n");
133       printf ("got      ");
134       mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
135       printf ("\n");
136       exit (1);
137     }
138 
139   /* check gamma_inc(-2,1) = 0.109691967197760 */
140   mpfr_set_si (a, -2, MPFR_RNDN);
141   mpfr_set_ui (x, 1, MPFR_RNDN);
142   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
143   if (mpfr_cmp_ui_2exp (a, 1, -3))
144     {
145       printf ("Error for gamma_inc(-2,1)\n");
146       printf ("expected 0.125\n");
147       printf ("got      ");
148       mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
149       printf ("\n");
150       exit (1);
151     }
152 
153   /* check gamma_inc(-3,1) = 0.109691967197760 */
154   mpfr_set_si (a, -3, MPFR_RNDN);
155   mpfr_set_ui (x, 1, MPFR_RNDN);
156   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
157   if (mpfr_cmp_ui_2exp (a, 3, -5))
158     {
159       printf ("Error for gamma_inc(-3,1)\n");
160       printf ("expected 3/32\n");
161       printf ("got      ");
162       mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
163       printf ("\n");
164       exit (1);
165     }
166 
167   /* check gamma_inc(-100,1) = 0.00364201018241046 */
168   mpfr_set_si (a, -100, MPFR_RNDN);
169   mpfr_set_ui (x, 1, MPFR_RNDN);
170   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
171   if (mpfr_cmp_ui_2exp (a, 1, -8))
172     {
173       printf ("Error for gamma_inc(-100,1)\n");
174       printf ("expected 1/256\n");
175       printf ("got      ");
176       mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
177       printf ("\n");
178       exit (1);
179     }
180 
181   /* TODO: Once internal overflow is supported, add new tests with
182      larger exponents, e.g. 64 (in addition to 25). */
183   mpfr_set_prec (a, 1);
184   mpfr_set_prec (x, 1);
185   mpfr_set_ui_2exp (a, 1, 25, MPFR_RNDN);
186   mpfr_set_ui_2exp (x, 1, -25, MPFR_RNDN);
187   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
188 
189   mpfr_clear (a);
190   mpfr_clear (x);
191 }
192 
193 /* tests for negative integer a: for -n <= a <= -1, perform k tests
194    with random x in 0..|a| and precision 'prec' */
195 static void
196 test_negint (long n, unsigned long k, mpfr_prec_t prec)
197 {
198   long i, j;
199   mpfr_t a, x, y;
200 
201   mpfr_init2 (a, prec);
202   mpfr_init2 (x, prec);
203   mpfr_init2 (y, prec);
204   for (i = 1; i <= n; i++)
205     {
206       mpfr_set_si (a, -i, MPFR_RNDN);
207       for (j = 1; j <= k; j++)
208         {
209           mpfr_urandomb (x, RANDS);
210           mpfr_mul_si (x, x, j, MPFR_RNDN);
211           mpfr_set_prec (y, prec + 20);
212           mpfr_gamma_inc (y, a, x, MPFR_RNDZ);
213           mpfr_gamma_inc (x, a, x, MPFR_RNDZ);
214           mpfr_prec_round (y, prec, MPFR_RNDZ);
215           if (mpfr_cmp (x, y))
216             {
217               printf ("Error in mpfr_gamma_inc(%ld,%ld) with MPFR_RNDZ\n",
218                       -i, j);
219               printf ("expected ");
220               mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN);
221               printf ("\ngot      ");
222               mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
223               printf ("\n");
224               exit (1);
225             }
226         }
227     }
228   mpfr_clear (a);
229   mpfr_clear (x);
230   mpfr_clear (y);
231 }
232 
233 int
234 main (int argc, char *argv[])
235 {
236   mpfr_prec_t p;
237 
238   tests_start_mpfr ();
239 
240   if (argc == 4) /* tgamma_inc a x prec: print gamma_inc(a,x) to prec bits */
241     {
242       mpfr_prec_t p = atoi (argv[3]);
243       mpfr_t a, x;
244       mpfr_init2 (a, p);
245       mpfr_init2 (x, p);
246       mpfr_set_str (a, argv[1], 10, MPFR_RNDN);
247       mpfr_set_str (x, argv[2], 10, MPFR_RNDN);
248       mpfr_gamma_inc (x, a, x, MPFR_RNDN);
249       mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
250       printf ("\n");
251       mpfr_clear (a);
252       mpfr_clear (x);
253       return 0;
254     }
255 
256   specials ();
257 
258   test_negint (30, 10, 53);
259 
260   for (p = MPFR_PREC_MIN; p < 100; p++)
261     test_random (p, 10);
262 
263   test_generic (MPFR_PREC_MIN, 100, 5);
264 
265   tests_end_mpfr ();
266   return 0;
267 }
268