1 /* tj1 -- test file for the Bessel function of first kind (order 1)
2
3 Copyright 2007-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 TEST_FUNCTION mpfr_j1
26 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 8, RANDS)
27 #define REDUCE_EMAX 262143 /* otherwise arg. reduction is too expensive */
28 #include "tgeneric.c"
29
30 /* test for small input, where j1(x) = x/2 - x^3/16 + ... */
31 static void
test_small(void)32 test_small (void)
33 {
34 mpfr_t x, y;
35 int inex, sign;
36 mpfr_exp_t e, emin;
37
38 mpfr_init2 (x, 10);
39 mpfr_init2 (y, 9);
40 emin = mpfr_get_emin ();
41 for (e = -4; e >= -30; e--)
42 {
43 if (e == -30)
44 {
45 e = mpfr_get_emin_min () - 1;
46 set_emin (e + 1);
47 }
48 for (sign = -1; sign <= 1; sign += 2)
49 {
50 mpfr_set_si_2exp (x, sign, e, MPFR_RNDN);
51 mpfr_nexttoinf (x);
52 inex = mpfr_j1 (y, x, MPFR_RNDN);
53 if (e >= -29)
54 {
55 /* since |x| is just above 2^e, |j1(x)| is just above 2^(e-1),
56 thus y should be 2^(e-1) and the inexact flag should be
57 of opposite sign of x */
58 MPFR_ASSERTN(mpfr_cmp_si_2exp0 (y, sign, e - 1) == 0);
59 MPFR_ASSERTN(VSIGN (inex) * sign < 0);
60 }
61 else
62 {
63 /* here |y| should be 0.5*2^emin and the inexact flag should
64 have the sign of x */
65 MPFR_ASSERTN(mpfr_cmp_si_2exp0 (y, sign, e) == 0);
66 MPFR_ASSERTN(VSIGN (inex) * sign > 0);
67 }
68 }
69 }
70 set_emin (emin);
71 mpfr_clear (x);
72 mpfr_clear (y);
73 }
74
75 /* a test that fails with GMP_CHECK_RANDOMIZE=1613146232984428
76 on revision 14429 */
77 static void
bug20210215(void)78 bug20210215 (void)
79 {
80 mpfr_t x, y;
81 int inex;
82
83 mpfr_init2 (x, 221);
84 mpfr_init2 (y, 1);
85 mpfr_set_str (x, "1.6484611511696130037307738844228498447763863563070374544054791168614e+01", 10, MPFR_RNDN);
86 mpfr_clear_flags ();
87 inex = mpfr_j1 (y, x, MPFR_RNDZ);
88 MPFR_ASSERTN (mpfr_cmp_si_2exp0 (y, -1, -9) == 0);
89 MPFR_ASSERTN (inex > 0);
90 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_INEXACT);
91 mpfr_clear (x);
92 mpfr_clear (y);
93 }
94
95 int
main(int argc,char * argv[])96 main (int argc, char *argv[])
97 {
98 mpfr_t x, y;
99
100 tests_start_mpfr ();
101
102 bug20210215 ();
103
104 test_small ();
105
106 mpfr_init (x);
107 mpfr_init (y);
108
109 /* special values */
110 mpfr_set_nan (x);
111 mpfr_j1 (y, x, MPFR_RNDN);
112 MPFR_ASSERTN(mpfr_nan_p (y));
113
114 mpfr_set_inf (x, 1); /* +Inf */
115 mpfr_j1 (y, x, MPFR_RNDN);
116 MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_POS (y));
117
118 mpfr_set_inf (x, -1); /* -Inf */
119 mpfr_j1 (y, x, MPFR_RNDN);
120 MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_POS (y));
121
122 mpfr_set_ui (x, 0, MPFR_RNDN); /* +0 */
123 mpfr_j1 (y, x, MPFR_RNDN);
124 MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_POS (y)); /* j1(+0)=+0 */
125
126 mpfr_set_ui (x, 0, MPFR_RNDN);
127 mpfr_neg (x, x, MPFR_RNDN); /* -0 */
128 mpfr_j1 (y, x, MPFR_RNDN);
129 MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_NEG (y)); /* j1(-0)=-0 */
130
131 mpfr_set_prec (x, 53);
132 mpfr_set_prec (y, 53);
133
134 mpfr_set_ui (x, 1, MPFR_RNDN);
135 mpfr_j1 (y, x, MPFR_RNDN);
136 mpfr_set_str_binary (x, "0.0111000010100111001001111011101001011100001100011011");
137 if (mpfr_cmp (x, y))
138 {
139 printf ("Error in mpfr_j1 for x=1, rnd=MPFR_RNDN\n");
140 printf ("Expected "); mpfr_dump (x);
141 printf ("Got "); mpfr_dump (y);
142 exit (1);
143 }
144 mpfr_clear (x);
145 mpfr_clear (y);
146
147 test_generic (MPFR_PREC_MIN, 100, 10);
148
149 data_check ("data/j1", mpfr_j1, "mpfr_j1");
150
151 tests_end_mpfr ();
152
153 return 0;
154 }
155