1 /* Test file for mpfr_sin.
2 
3 Copyright 2001-2020 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
test_sin(mpfr_ptr a,mpfr_srcptr b,mpfr_rnd_t rnd_mode)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
check53(const char * xs,const char * sin_xs,mpfr_rnd_t rnd_mode)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
check53b(const char * xs,const char * sin_xs,mpfr_rnd_t rnd_mode)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
test_sign(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
check_nans(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
check_regression(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
check_tiny(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
check_binary128(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
check_212(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
main(int argc,char * argv[])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 }
466