1 /* Test file for mpfr_sin.
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 #ifdef CHECK_EXTERNAL
26 static int
27 test_sin (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode)
28 {
29 int res;
30 int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_get_prec (a)>=53;
31 if (ok)
32 {
33 mpfr_print_raw (b);
34 }
35 res = mpfr_sin (a, b, rnd_mode);
36 if (ok)
37 {
38 printf (" ");
39 mpfr_print_raw (a);
40 printf ("\n");
41 }
42 return res;
43 }
44 #else
45 #define test_sin mpfr_sin
46 #endif
47
48 static void
49 check53 (const char *xs, const char *sin_xs, mpfr_rnd_t rnd_mode)
50 {
51 mpfr_t xx, s;
52
53 mpfr_init2 (xx, 53);
54 mpfr_init2 (s, 53);
55 mpfr_set_str1 (xx, xs); /* should be exact */
56 test_sin (s, xx, rnd_mode);
57 if (mpfr_cmp_str1 (s, sin_xs))
58 {
59 printf ("mpfr_sin failed for x=%s, rnd=%s\n",
60 xs, mpfr_print_rnd_mode (rnd_mode));
61 printf ("mpfr_sin gives sin(x)=");
62 mpfr_out_str (stdout, 10, 0, s, MPFR_RNDN);
63 printf (", expected %s\n", sin_xs);
64 exit (1);
65 }
66 mpfr_clear (xx);
67 mpfr_clear (s);
68 }
69
70 static void
71 check53b (const char *xs, const char *sin_xs, mpfr_rnd_t rnd_mode)
72 {
73 mpfr_t xx, s;
74
75 mpfr_init2 (xx, 53);
76 mpfr_init2 (s, 53);
77 mpfr_set_str (xx, xs, 2, MPFR_RNDN); /* should be exact */
78 test_sin (s, xx, rnd_mode);
79 if (mpfr_cmp_str (s, sin_xs, 2, MPFR_RNDN))
80 {
81 printf ("mpfr_sin failed in rounding mode %s for\n x = %s\n",
82 mpfr_print_rnd_mode (rnd_mode), xs);
83 printf (" got ");
84 mpfr_out_str (stdout, 2, 0, s, MPFR_RNDN);
85 printf ("\nexpected %s\n", sin_xs);
86 exit (1);
87 }
88 mpfr_clear (xx);
89 mpfr_clear (s);
90 }
91
92 static void
93 test_sign (void)
94 {
95 mpfr_t pid, piu, x, y;
96 int p, k;
97
98 mpfr_init2 (pid, 4096);
99 mpfr_const_pi (pid, MPFR_RNDD);
100 mpfr_init2 (piu, 4096);
101 mpfr_const_pi (piu, MPFR_RNDU);
102 mpfr_init (x);
103 mpfr_init2 (y, 2);
104 for (p = 8; p <= 128; p++)
105 for (k = 2; k <= 6; k += 2)
106 {
107 mpfr_set_prec (x, p);
108 mpfr_mul_ui (x, pid, k, MPFR_RNDD);
109 test_sin (y, x, MPFR_RNDN);
110 if (MPFR_IS_POS (y))
111 {
112 printf ("Error in test_sign for sin(%dpi-epsilon), prec = %d"
113 " for argument.\nResult should have been negative.\n",
114 k, p);
115 exit (1);
116 }
117 mpfr_mul_ui (x, piu, k, MPFR_RNDU);
118 test_sin (y, x, MPFR_RNDN);
119 if (MPFR_IS_NEG (y))
120 {
121 printf ("Error in test_sign for sin(%dpi+epsilon), prec = %d"
122 " for argument.\nResult should have been positive.\n",
123 k, p);
124 exit (1);
125 }
126 }
127
128 /* worst case on 53 bits */
129 mpfr_set_prec (x, 53);
130 mpfr_set_prec (y, 53);
131 mpfr_set_str (x, "6134899525417045", 10, MPFR_RNDN);
132 test_sin (y, x, MPFR_RNDN);
133 mpfr_set_str_binary (x, "11011010111101011110111100010101010101110000000001011E-106");
134 MPFR_ASSERTN(mpfr_cmp (x, y) == 0);
135
136 /* Bug on Special cases */
137 mpfr_set_str_binary (x, "0.100011011010111101E-32");
138 test_sin (y, x, MPFR_RNDN);
139 if (mpfr_cmp_str (y, "0.10001101101011110100000000000000000000000000000000000E-32", 2, MPFR_RNDN))
140 {
141 printf("sin special 97 error:\nx=");
142 mpfr_dump (x); printf("y=");
143 mpfr_dump (y);
144 exit (1);
145 }
146
147 mpfr_set_prec (x, 53);
148 mpfr_set_prec (y, 53);
149 mpfr_set_str_binary (x, "1.1001001000011111101101010100010001000010110100010011");
150 mpfr_set_str_binary (y, "1.1111111111111111111111111111111111111111111111111111e-1");
151 test_sin (x, x, MPFR_RNDZ);
152 MPFR_ASSERTN(mpfr_cmp (x, y) == 0);
153
154 mpfr_clear (pid);
155 mpfr_clear (piu);
156 mpfr_clear (x);
157 mpfr_clear (y);
158 }
159
160 static void
161 check_nans (void)
162 {
163 mpfr_t x, y;
164
165 mpfr_init2 (x, 123L);
166 mpfr_init2 (y, 123L);
167
168 mpfr_set_nan (x);
169 test_sin (y, x, MPFR_RNDN);
170 if (! mpfr_nan_p (y))
171 {
172 printf ("Error: sin(NaN) != NaN\n");
173 exit (1);
174 }
175
176 mpfr_set_inf (x, 1);
177 test_sin (y, x, MPFR_RNDN);
178 if (! mpfr_nan_p (y))
179 {
180 printf ("Error: sin(Inf) != NaN\n");
181 exit (1);
182 }
183
184 mpfr_set_inf (x, -1);
185 test_sin (y, x, MPFR_RNDN);
186 if (! mpfr_nan_p (y))
187 {
188 printf ("Error: sin(-Inf) != NaN\n");
189 exit (1);
190 }
191
192 mpfr_clear (x);
193 mpfr_clear (y);
194 }
195
196 #define TEST_FUNCTION test_sin
197 #ifndef MPFR_USE_MINI_GMP
198 #define REDUCE_EMAX 262143 /* otherwise arg. reduction is too expensive */
199 #else
200 #define REDUCE_EMAX 16383 /* reduce further since mini-gmp works in O(n^2) */
201 #endif
202 #include "tgeneric.c"
203
204 const char xs[] = "0.111011111110110000111000001100000111110E-1";
205
206 static void
207 check_regression (void)
208 {
209 mpfr_t x, y;
210 mpfr_prec_t p;
211 int i;
212
213 p = strlen (xs) - 2 - 3;
214 mpfr_inits2 (p, x, y, (mpfr_ptr) 0);
215
216 mpfr_set_str (x, xs, 2, MPFR_RNDN);
217 i = mpfr_sin (y, x, MPFR_RNDN);
218 if (i >= 0
219 || mpfr_cmp_str (y, "0.111001110011110011110001010110011101110E-1",
220 2, MPFR_RNDN))
221 {
222 printf ("Regression test failed (1) i=%d\ny=", i);
223 mpfr_dump (y);
224 exit (1);
225 }
226 mpfr_clears (x, y, (mpfr_ptr) 0);
227 }
228
229 /* Test provided by Christopher Creutzig, 2007-05-21. */
230 static void
231 check_tiny (void)
232 {
233 mpfr_t x, y;
234
235 mpfr_init2 (x, 53);
236 mpfr_init2 (y, 53);
237
238 mpfr_set_ui (x, 1, MPFR_RNDN);
239 mpfr_set_exp (x, mpfr_get_emin ());
240 mpfr_sin (y, x, MPFR_RNDD);
241 if (mpfr_cmp (x, y) < 0)
242 {
243 printf ("Error in check_tiny: got sin(x) > x for x = 2^(emin-1)\n");
244 exit (1);
245 }
246
247 mpfr_sin (y, x, MPFR_RNDU);
248 mpfr_mul_2ui (y, y, 1, MPFR_RNDU);
249 if (mpfr_cmp (x, y) > 0)
250 {
251 printf ("Error in check_tiny: got sin(x) < x/2 for x = 2^(emin-1)\n");
252 exit (1);
253 }
254
255 mpfr_neg (x, x, MPFR_RNDN);
256 mpfr_sin (y, x, MPFR_RNDU);
257 if (mpfr_cmp (x, y) > 0)
258 {
259 printf ("Error in check_tiny: got sin(x) < x for x = -2^(emin-1)\n");
260 exit (1);
261 }
262
263 mpfr_sin (y, x, MPFR_RNDD);
264 mpfr_mul_2ui (y, y, 1, MPFR_RNDD);
265 if (mpfr_cmp (x, y) < 0)
266 {
267 printf ("Error in check_tiny: got sin(x) > x/2 for x = -2^(emin-1)\n");
268 exit (1);
269 }
270
271 mpfr_clear (y);
272 mpfr_clear (x);
273 }
274
275 static void
276 check_binary128 (void)
277 {
278 mpfr_t x, y, z;
279
280 mpfr_init2 (x, 113);
281 mpfr_init2 (y, 113);
282 mpfr_init2 (z, 113);
283
284 /* number closest to a odd multiple of pi/2 in the binary128 format:
285 8794873135033829349702184924722639 * 2^1852 */
286 mpfr_set_str (x, "1b19ee7c329d7d951906d1e11b5cfp1852", 16, MPFR_RNDN);
287 mpfr_cos (y, x, MPFR_RNDN);
288 mpfr_set_str (z, "1.ad1a2037cd7820f748483f5d39c3p-124", 16, MPFR_RNDN);
289 if (! mpfr_equal_p (y, z))
290 {
291 printf ("Error in check_binary128 (cos x)\n");
292 printf ("expected "); mpfr_dump (z);
293 printf ("got "); mpfr_dump (y);
294 exit (1);
295 }
296 mpfr_sin (y, x, MPFR_RNDN);
297 mpfr_set_ui (z, 1, MPFR_RNDN);
298 if (! mpfr_equal_p (y, z))
299 {
300 printf ("Error in check_binary128 (sin x)\n");
301 printf ("expected "); mpfr_dump (z);
302 printf ("got "); mpfr_dump (y);
303 exit (1);
304 }
305
306 /* now multiply x by 2, so that it is near an even multiple of pi/2 */
307 mpfr_mul_2ui (x, x, 1, MPFR_RNDN);
308 mpfr_cos (y, x, MPFR_RNDN);
309 mpfr_set_si (z, -1, MPFR_RNDN);
310 if (! mpfr_equal_p (y, z))
311 {
312 printf ("Error in check_binary128 (cos 2x)\n");
313 printf ("expected "); mpfr_dump (z);
314 printf ("got "); mpfr_dump (y);
315 exit (1);
316 }
317 mpfr_sin (y, x, MPFR_RNDN);
318 mpfr_set_str (z, "3.5a34406f9af041ee90907eba7386p-124", 16, MPFR_RNDN);
319 if (! mpfr_equal_p (y, z))
320 {
321 printf ("Error in check_binary128 (sin 2x)\n");
322 printf ("expected "); mpfr_dump (z);
323 printf ("got "); mpfr_dump (y);
324 exit (1);
325 }
326
327 mpfr_clear (x);
328 mpfr_clear (y);
329 mpfr_clear (z);
330 }
331
332 /* check Ziv's loop with precision 212 bits */
333 static void
334 check_212 (void)
335 {
336 mpfr_t x, y, z;
337
338 mpfr_init2 (x, 212);
339 mpfr_init2 (y, 212);
340 mpfr_init2 (z, 212);
341
342 mpfr_set_str (x, "f.c34b10aa02f796d435a3db0146b4e8a0b2850422f778af06be66p+0", 16, MPFR_RNDN);
343 mpfr_sin (y, x, MPFR_RNDN);
344 mpfr_set_str (z, "-e.0c2d5c189f8a0d185d7036b87b90f3040f4f2aa0f46f901bad44p-8", 16, MPFR_RNDN);
345 if (! mpfr_equal_p (y, z))
346 {
347 printf ("Error in check_212\n");
348 printf ("expected "); mpfr_dump (z);
349 printf ("got "); mpfr_dump (y);
350 exit (1);
351 }
352
353 mpfr_clear (x);
354 mpfr_clear (y);
355 mpfr_clear (z);
356 }
357
358 int
359 main (int argc, char *argv[])
360 {
361 mpfr_t x, c, s, c2, s2;
362
363 tests_start_mpfr ();
364
365 check_regression ();
366 check_nans ();
367
368 /* worst case from PhD thesis of Vincent Lefe`vre: x=8980155785351021/2^54 */
369 check53 ("4.984987858808754279e-1", "4.781075595393330379e-1", MPFR_RNDN);
370 check53 ("4.984987858808754279e-1", "4.781075595393329824e-1", MPFR_RNDD);
371 check53 ("4.984987858808754279e-1", "4.781075595393329824e-1", MPFR_RNDZ);
372 check53 ("4.984987858808754279e-1", "4.781075595393330379e-1", MPFR_RNDU);
373 check53 ("1.00031274099908640274", "8.416399183372403892e-1", MPFR_RNDN);
374 check53 ("1.00229256850978698523", "8.427074524447979442e-1", MPFR_RNDZ);
375 check53 ("1.00288304857059840103", "8.430252033025980029e-1", MPFR_RNDZ);
376 check53 ("1.00591265847407274059", "8.446508805292128885e-1", MPFR_RNDN);
377
378 /* Other worst cases showing a bug introduced on 2005-01-29 in rev 3248 */
379 check53b ("1.0111001111010111010111111000010011010001110001111011e-21",
380 "1.0111001111010111010111111000010011010001101001110001e-21",
381 MPFR_RNDU);
382 check53b ("1.1011101111111010000001010111000010000111100100101101",
383 "1.1111100100101100001111100000110011110011010001010101e-1",
384 MPFR_RNDU);
385
386 mpfr_init2 (x, 2);
387
388 mpfr_set_str (x, "0.5", 10, MPFR_RNDN);
389 test_sin (x, x, MPFR_RNDD);
390 if (mpfr_cmp_ui_2exp (x, 3, -3)) /* x != 0.375 = 3/8 */
391 {
392 printf ("mpfr_sin(0.5, MPFR_RNDD) failed with precision=2\n");
393 exit (1);
394 }
395
396 /* bug found by Kevin Ryde */
397 mpfr_const_pi (x, MPFR_RNDN);
398 mpfr_mul_ui (x, x, 3L, MPFR_RNDN);
399 mpfr_div_ui (x, x, 2L, MPFR_RNDN);
400 test_sin (x, x, MPFR_RNDN);
401 if (mpfr_cmp_ui (x, 0) >= 0)
402 {
403 printf ("Error: wrong sign for sin(3*Pi/2)\n");
404 exit (1);
405 }
406
407 /* Can fail on an assert */
408 mpfr_set_prec (x, 53);
409 mpfr_set_str (x, "77291789194529019661184401408", 10, MPFR_RNDN);
410 mpfr_init2 (c, 4); mpfr_init2 (s, 42);
411 mpfr_init2 (c2, 4); mpfr_init2 (s2, 42);
412
413 test_sin (s, x, MPFR_RNDN);
414 mpfr_cos (c, x, MPFR_RNDN);
415 mpfr_sin_cos (s2, c2, x, MPFR_RNDN);
416 if (mpfr_cmp (c2, c))
417 {
418 printf("cos differs for x=77291789194529019661184401408");
419 exit (1);
420 }
421 if (mpfr_cmp (s2, s))
422 {
423 printf("sin differs for x=77291789194529019661184401408");
424 exit (1);
425 }
426
427 mpfr_set_str_binary (x, "1.1001001000011111101101010100010001000010110100010011");
428 test_sin (x, x, MPFR_RNDZ);
429 if (mpfr_cmp_str (x, "1.1111111111111111111111111111111111111111111111111111e-1", 2, MPFR_RNDN))
430 {
431 printf ("Error for x= 1.1001001000011111101101010100010001000010110100010011\nGot ");
432 mpfr_dump (x);
433 exit (1);
434 }
435
436 mpfr_set_prec (s, 9);
437 mpfr_set_prec (x, 190);
438 mpfr_const_pi (x, MPFR_RNDN);
439 mpfr_sin (s, x, MPFR_RNDZ);
440 if (mpfr_cmp_str (s, "0.100000101e-196", 2, MPFR_RNDN))
441 {
442 printf ("Error for x ~= pi\n");
443 mpfr_dump (s);
444 exit (1);
445 }
446
447 mpfr_clear (s2);
448 mpfr_clear (c2);
449 mpfr_clear (s);
450 mpfr_clear (c);
451 mpfr_clear (x);
452
453 test_generic (MPFR_PREC_MIN, 100, 15);
454 test_generic (MPFR_SINCOS_THRESHOLD-1, MPFR_SINCOS_THRESHOLD+1, 2);
455 test_sign ();
456 check_tiny ();
457 check_binary128 ();
458 check_212 ();
459
460 data_check ("data/sin", mpfr_sin, "mpfr_sin");
461 bad_cases (mpfr_sin, mpfr_asin, "mpfr_sin", 256, -40, 0, 4, 128, 800, 50);
462
463 tests_end_mpfr ();
464 return 0;
465 }