1 /* Test file for mpfr_sinu.
2
3 Copyright 2020-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 static void
26 test_singular (void)
27 {
28 mpfr_t x, y;
29 int inexact;
30
31 mpfr_init (x);
32 mpfr_init (y);
33
34 /* check u = 0 */
35 mpfr_set_ui (x, 17, MPFR_RNDN);
36 inexact = mpfr_sinu (y, x, 0, MPFR_RNDN);
37 MPFR_ASSERTN(mpfr_nan_p (y));
38
39 /* check x = NaN */
40 mpfr_set_nan (x);
41 inexact = mpfr_sinu (y, x, 1, MPFR_RNDN);
42 MPFR_ASSERTN(mpfr_nan_p (y));
43
44 /* check x = +Inf */
45 mpfr_set_inf (x, 1);
46 inexact = mpfr_sinu (y, x, 1, MPFR_RNDN);
47 MPFR_ASSERTN(mpfr_nan_p (y));
48
49 /* check x = -Inf */
50 mpfr_set_inf (x, -1);
51 inexact = mpfr_sinu (y, x, 1, MPFR_RNDN);
52 MPFR_ASSERTN(mpfr_nan_p (y));
53
54 /* check x = +0 */
55 mpfr_set_zero (x, 1);
56 inexact = mpfr_sinu (y, x, 1, MPFR_RNDN);
57 MPFR_ASSERTN(mpfr_zero_p (y) && mpfr_signbit (y) == 0);
58 MPFR_ASSERTN(inexact == 0);
59
60 /* check x = -0 */
61 mpfr_set_zero (x, -1);
62 inexact = mpfr_sinu (y, x, 1, MPFR_RNDN);
63 MPFR_ASSERTN(mpfr_zero_p (y) && mpfr_signbit (y) != 0);
64 MPFR_ASSERTN(inexact == 0);
65
66 mpfr_clear (x);
67 mpfr_clear (y);
68 }
69
70 static void
71 test_exact (void)
72 {
73 mpfr_t x, y;
74 int inexact;
75
76 mpfr_init (x);
77 mpfr_init (y);
78
79 /* check 2*pi*x/u = pi/2 thus x/u = 1/4, for example x=1 and u=4 */
80 mpfr_set_ui (x, 1, MPFR_RNDN);
81 inexact = mpfr_sinu (y, x, 4, MPFR_RNDN);
82 MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0 && inexact == 0);
83
84 /* check 2*pi*x/u = pi thus x/u = 1/2, for example x=2 and u=4 */
85 mpfr_set_ui (x, 2, MPFR_RNDN);
86 inexact = mpfr_sinu (y, x, 4, MPFR_RNDN);
87 MPFR_ASSERTN(mpfr_zero_p (y) && mpfr_signbit (y) == 0);
88 MPFR_ASSERTN(inexact == 0);
89
90 /* check 2*pi*x/u = 3*pi/2 thus x/u = 3/4, for example x=3 and u=4 */
91 mpfr_set_ui (x, 3, MPFR_RNDN);
92 inexact = mpfr_sinu (y, x, 4, MPFR_RNDN);
93 MPFR_ASSERTN(mpfr_cmp_si (y, -1) == 0 && inexact == 0);
94
95 /* check 2*pi*x/u = 2*pi thus x/u = 1, for example x=4 and u=4 */
96 mpfr_set_ui (x, 4, MPFR_RNDN);
97 inexact = mpfr_sinu (y, x, 4, MPFR_RNDN);
98 MPFR_ASSERTN(mpfr_zero_p (y) && mpfr_signbit (y) == 0);
99 MPFR_ASSERTN(inexact == 0);
100
101 /* check x/u = 2^16, for example x=3*2^16 and u=3 */
102 mpfr_set_ui_2exp (x, 3, 16, MPFR_RNDN);
103 inexact = mpfr_sinu (y, x, 3, MPFR_RNDN);
104 MPFR_ASSERTN(mpfr_zero_p (y) && mpfr_signbit (y) == 0);
105 MPFR_ASSERTN(inexact == 0);
106
107 /* check 2*pi*x/u = -pi/2 thus x/u = -1/4, for example x=-1 and u=4 */
108 mpfr_set_si (x, -1, MPFR_RNDN);
109 inexact = mpfr_sinu (y, x, 4, MPFR_RNDN);
110 MPFR_ASSERTN(mpfr_cmp_si (y, -1) == 0 && inexact == 0);
111
112 /* check 2*pi*x/u = -pi thus x/u = -1/2, for example x=-2 and u=4 */
113 mpfr_set_si (x, -2, MPFR_RNDN);
114 inexact = mpfr_sinu (y, x, 4, MPFR_RNDN);
115 MPFR_ASSERTN(mpfr_zero_p (y) && mpfr_signbit (y) != 0);
116 MPFR_ASSERTN(inexact == 0);
117
118 /* check 2*pi*x/u = -3*pi/2 thus x/u = -3/4, for example x=-3 and u=4 */
119 mpfr_set_si (x, -3, MPFR_RNDN);
120 inexact = mpfr_sinu (y, x, 4, MPFR_RNDN);
121 MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0 && inexact == 0);
122
123 /* check 2*pi*x/u = -2*pi thus x/u = -1, for example x=-4 and u=4 */
124 mpfr_set_si (x, -4, MPFR_RNDN);
125 inexact = mpfr_sinu (y, x, 4, MPFR_RNDN);
126 MPFR_ASSERTN(mpfr_zero_p (y) && mpfr_signbit (y) != 0);
127 MPFR_ASSERTN(inexact == 0);
128
129 /* check x/u = -2^16, for example x=-3*2^16 and u=3 */
130 mpfr_set_si_2exp (x, -3, 16, MPFR_RNDN);
131 inexact = mpfr_sinu (y, x, 3, MPFR_RNDN);
132 MPFR_ASSERTN(mpfr_zero_p (y) && mpfr_signbit (y) != 0);
133 MPFR_ASSERTN(inexact == 0);
134
135 /* check 2*pi*x/u = pi/6, for example x=1 and u=12 */
136 mpfr_set_ui (x, 1, MPFR_RNDN);
137 inexact = mpfr_sinu (y, x, 12, MPFR_RNDN);
138 MPFR_ASSERTN(mpfr_cmp_ui_2exp (y, 1, -1) == 0 && inexact == 0);
139
140 /* check 2*pi*x/u = 5*pi/6, for example x=5 and u=12 */
141 mpfr_set_ui (x, 5, MPFR_RNDN);
142 inexact = mpfr_sinu (y, x, 12, MPFR_RNDN);
143 MPFR_ASSERTN(mpfr_cmp_ui_2exp (y, 1, -1) == 0 && inexact == 0);
144
145 /* check 2*pi*x/u = 7*pi/6, for example x=5 and u=12 */
146 mpfr_set_ui (x, 7, MPFR_RNDN);
147 inexact = mpfr_sinu (y, x, 12, MPFR_RNDN);
148 MPFR_ASSERTN(mpfr_cmp_si_2exp (y, -1, -1) == 0 && inexact == 0);
149
150 /* check 2*pi*x/u = 11*pi/6, for example x=5 and u=12 */
151 mpfr_set_ui (x, 11, MPFR_RNDN);
152 inexact = mpfr_sinu (y, x, 12, MPFR_RNDN);
153 MPFR_ASSERTN(mpfr_cmp_si_2exp (y, -1, -1) == 0 && inexact == 0);
154
155 mpfr_clear (x);
156 mpfr_clear (y);
157 }
158
159 static void
160 test_regular (void)
161 {
162 mpfr_t x, y, z;
163 int inexact;
164
165 mpfr_init2 (x, 53);
166 mpfr_init2 (y, 53);
167 mpfr_init2 (z, 53);
168
169 mpfr_set_ui (x, 17, MPFR_RNDN);
170 inexact = mpfr_sinu (y, x, 42, MPFR_RNDN);
171 /* y should be sin(2*17*pi/42) rounded to nearest */
172 mpfr_set_str (z, "0x9.035be4a906768p-4", 16, MPFR_RNDN);
173 MPFR_ASSERTN(mpfr_equal_p (y, z));
174 MPFR_ASSERTN(inexact > 0);
175
176 mpfr_clear (x);
177 mpfr_clear (y);
178 mpfr_clear (z);
179 }
180
181 /* Check argument reduction with large hard-coded inputs. The following
182 values were generated with gen_random(sin,10,53,100,20), where the
183 Sage code for gen_random is given in the tcosu.c file. */
184 static void
185 test_large (void)
186 {
187 static struct {
188 const char *x;
189 unsigned long u;
190 const char *y;
191 } t[] = {
192 { "-0x4.338bcaf01cb4p+16", 37, "0xa.a85c8758c1228p-4" },
193 { "-0x6.01aa844d11acp+16", 12, "0x4.46e2cea96ee2cp-4" },
194 { "-0x7.997759de3b3dp+16", 35, "0xf.144b13bc0b008p-4" },
195 { "-0xc.74f72253bc7ep+16", 74, "-0x9.1c24e8ac6fb48p-4" },
196 { "-0x4.ea97564d3ecp+12", 93, "0x3.2d49033cc4556p-4" },
197 { "0xc.1cd85e0766ffp+16", 48, "-0xd.74b642b1c7368p-4" },
198 { "0xe.add2f0441bap+16", 82, "-0xf.fc35c5a5e1448p-4" },
199 { "0xc.8feb943a0fe5p+16", 43, "-0x5.8af3e0e11b794p-4" },
200 { "-0xd.d61b06ace6d9p+16", 55, "0xa.175c9f7e4ede8p-4" },
201 { "0xb.cf036b4e5fap+12", 46, "0x1.b7128e1975d9ap-4" },
202 /* The following case does not come from the random values.
203 Here, x mod u may need more precision than x. */
204 { "0x3.p+16", 99, "-0xb.e4a8a986c9fd03ap-5" }
205 };
206 int i;
207 mpfr_t x, y, z;
208
209 mpfr_inits2 (53, x, y, z, (mpfr_ptr) 0);
210 for (i = 0; i < numberof (t); i++)
211 {
212 mpfr_set_str (x, t[i].x, 0, MPFR_RNDN);
213 mpfr_set_str (y, t[i].y, 0, MPFR_RNDN);
214 mpfr_sinu (z, x, t[i].u, MPFR_RNDN);
215 MPFR_ASSERTN (mpfr_equal_p (y, z));
216 /* Same test after reducing the precision of x to its minimum one. */
217 mpfr_prec_round (x, mpfr_min_prec (x), MPFR_RNDN);
218 mpfr_sinu (z, x, t[i].u, MPFR_RNDN);
219 MPFR_ASSERTN (mpfr_equal_p (y, z));
220 mpfr_prec_round (x, 53, MPFR_RNDN);
221 }
222 mpfr_clears (x, y, z, (mpfr_ptr) 0);
223 }
224
225 #define TEST_FUNCTION mpfr_sinu
226 #define ULONG_ARG2
227 #include "tgeneric.c"
228
229 static int
230 mpfr_sin2pi (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r)
231 {
232 return mpfr_sinu (y, x, 1, r);
233 }
234
235 int
236 main (void)
237 {
238 tests_start_mpfr ();
239
240 test_singular ();
241 test_exact ();
242 test_regular ();
243 test_large ();
244
245 test_generic (MPFR_PREC_MIN, 100, 1000);
246
247 data_check ("data/sin2pi", mpfr_sin2pi, "mpfr_sin2pi");
248
249 tests_end_mpfr ();
250 return 0;
251 }