xref: /netbsd-src/external/lgpl3/mpfr/dist/tests/tsin_cos.c (revision ba125506a622fe649968631a56eba5d42ff57863)
1 /* Test file for mpfr_sin_cos.
2 
3 Copyright 2000-2004, 2006-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
large_test(char * X,int prec,int N)26 large_test (char *X, int prec, int N)
27 {
28   int i;
29   mpfr_t x, s, c;
30 
31   mpfr_init2 (x, prec);
32   mpfr_init2 (s, prec);
33   mpfr_init2 (c, prec);
34   mpfr_set_str (x, X, 10, MPFR_RNDN);
35 
36   for (i = 0; i < N; i++)
37     mpfr_sin_cos (s, c, x, MPFR_RNDN);
38 
39   mpfr_clear (x);
40   mpfr_clear (s);
41   mpfr_clear (c);
42 }
43 
44 static void
check53(const char * xs,const char * sin_xs,const char * cos_xs,mpfr_rnd_t rnd_mode)45 check53 (const char *xs, const char *sin_xs, const char *cos_xs,
46          mpfr_rnd_t rnd_mode)
47 {
48   mpfr_t xx, s, c;
49 
50   mpfr_inits2 (53, xx, s, c, (mpfr_ptr) 0);
51   mpfr_set_str1 (xx, xs); /* should be exact */
52   mpfr_sin_cos (s, c, xx, rnd_mode);
53   if (mpfr_cmp_str1 (s, sin_xs))
54     {
55       printf ("mpfr_sin_cos failed for x=%s, rnd=%s\n",
56               xs, mpfr_print_rnd_mode (rnd_mode));
57       printf ("mpfr_sin_cos gives sin(x)=");
58       mpfr_out_str(stdout, 10, 0, s, MPFR_RNDN);
59       printf(", expected %s\n", sin_xs);
60       exit (1);
61     }
62   if (mpfr_cmp_str1 (c, cos_xs))
63     {
64       printf ("mpfr_sin_cos failed for x=%s, rnd=%s\n",
65               xs, mpfr_print_rnd_mode (rnd_mode));
66       printf ("mpfr_sin_cos gives cos(x)=");
67       mpfr_out_str(stdout, 10, 0, c, MPFR_RNDN);
68       printf(", expected %s\n", cos_xs);
69       exit (1);
70     }
71   mpfr_clears (xx, s, c, (mpfr_ptr) 0);
72 }
73 
74 static void
check53sin(const char * xs,const char * sin_xs,mpfr_rnd_t rnd_mode)75 check53sin (const char *xs, const char *sin_xs, mpfr_rnd_t rnd_mode)
76 {
77   mpfr_t xx, s, c;
78 
79   mpfr_inits2 (53, xx, s, c, (mpfr_ptr) 0);
80   mpfr_set_str1 (xx, xs); /* should be exact */
81   mpfr_sin_cos (s, c, xx, rnd_mode);
82   if (mpfr_cmp_str1 (s, sin_xs))
83     {
84       printf ("mpfr_sin_cos failed for x=%s, rnd=%s\n",
85               xs, mpfr_print_rnd_mode (rnd_mode));
86       printf ("mpfr_sin_cos gives sin(x)=");
87       mpfr_out_str(stdout, 10, 0, s, MPFR_RNDN);
88       printf(", expected %s\n", sin_xs);
89       exit (1);
90     }
91   mpfr_clears (xx, s, c, (mpfr_ptr) 0);
92 }
93 
94 static void
check53cos(const char * xs,const char * cos_xs,mpfr_rnd_t rnd_mode)95 check53cos (const char *xs, const char *cos_xs, mpfr_rnd_t rnd_mode)
96 {
97   mpfr_t xx, c, s;
98 
99   mpfr_inits2 (53, xx, s, c, (mpfr_ptr) 0);
100   mpfr_set_str1 (xx, xs); /* should be exact */
101   mpfr_sin_cos (s, c, xx, rnd_mode);
102   if (mpfr_cmp_str1 (c, cos_xs))
103     {
104       printf ("mpfr_sin_cos failed for x=%s, rnd=%s\n",
105               xs, mpfr_print_rnd_mode (rnd_mode));
106       printf ("mpfr_sin_cos gives cos(x)=");
107       mpfr_out_str(stdout, 10, 0, c, MPFR_RNDN);
108       printf(", expected %s\n", cos_xs);
109       exit (1);
110     }
111   mpfr_clears (xx, s, c, (mpfr_ptr) 0);
112 }
113 
114 static void
check_nans(void)115 check_nans (void)
116 {
117   mpfr_t  x, s, c;
118 
119   mpfr_init2 (x, 123L);
120   mpfr_init2 (s, 123L);
121   mpfr_init2 (c, 123L);
122 
123   /* sin(NaN)==NaN, cos(NaN)==NaN */
124   mpfr_set_nan (x);
125   mpfr_sin_cos (s, c, x, MPFR_RNDN);
126   MPFR_ASSERTN (mpfr_nan_p (s));
127   MPFR_ASSERTN (mpfr_nan_p (c));
128 
129   /* sin(+Inf)==NaN, cos(+Inf)==NaN */
130   mpfr_set_inf (x, 1);
131   mpfr_sin_cos (s, c, x, MPFR_RNDN);
132   MPFR_ASSERTN (mpfr_nan_p (s));
133   MPFR_ASSERTN (mpfr_nan_p (c));
134 
135   /* sin(-Inf)==NaN, cos(-Inf)==NaN */
136   mpfr_set_inf (x, -1);
137   mpfr_sin_cos (s, c, x, MPFR_RNDN);
138   MPFR_ASSERTN (mpfr_nan_p (s));
139   MPFR_ASSERTN (mpfr_nan_p (c));
140 
141   /* check zero */
142   mpfr_set_ui  (x, 0, MPFR_RNDN);
143   mpfr_sin_cos (s, c, x, MPFR_RNDN);
144   MPFR_ASSERTN (mpfr_cmp_ui (s, 0) == 0 && MPFR_IS_POS (s));
145   MPFR_ASSERTN (mpfr_cmp_ui (c, 1) == 0);
146   mpfr_neg (x, x, MPFR_RNDN);
147   mpfr_sin_cos (s, c, x, MPFR_RNDN);
148   MPFR_ASSERTN (mpfr_cmp_ui (s, 0) == 0 && MPFR_IS_NEG (s));
149   MPFR_ASSERTN (mpfr_cmp_ui (c, 1) == 0);
150 
151   /* coverage test */
152   mpfr_set_prec (x, 2);
153   mpfr_set_ui (x, 4, MPFR_RNDN);
154   mpfr_set_prec (s, 2);
155   mpfr_set_prec (c, 2);
156   mpfr_sin_cos (s, c, x, MPFR_RNDN);
157   MPFR_ASSERTN (mpfr_cmp_si_2exp (s, -3, -2) == 0);
158   MPFR_ASSERTN (mpfr_cmp_si_2exp (c, -3, -2) == 0);
159 
160   mpfr_clear (x);
161   mpfr_clear (s);
162   mpfr_clear (c);
163 }
164 
165 static void
overflowed_sin_cos0(void)166 overflowed_sin_cos0 (void)
167 {
168   mpfr_t x, y, z;
169   int emax, inex, rnd, err = 0;
170   mpfr_exp_t old_emax;
171 
172   old_emax = mpfr_get_emax ();
173 
174   mpfr_init2 (x, 8);
175   mpfr_init2 (y, 8);
176   mpfr_init2 (z, 8);
177 
178   for (emax = -1; emax <= 0; emax++)
179     {
180       mpfr_set_ui_2exp (z, 1, emax, MPFR_RNDN);
181       mpfr_nextbelow (z);
182       set_emax (emax);  /* 1 is not representable. */
183       /* and if emax < 0, 1 - eps is not representable either. */
184       RND_LOOP (rnd)
185         {
186           mpfr_set_si (x, 0, MPFR_RNDN);
187           mpfr_neg (x, x, MPFR_RNDN);
188           mpfr_clear_flags ();
189           inex = mpfr_sin_cos (x, y, x, (mpfr_rnd_t) rnd);
190           if (! mpfr_overflow_p ())
191             {
192               printf ("Error in overflowed_sin_cos0 (rnd = %s):\n"
193                       "  The overflow flag is not set.\n",
194                       mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
195               err = 1;
196             }
197           if (! (mpfr_zero_p (x) && MPFR_IS_NEG (x)))
198             {
199               printf ("Error in overflowed_sin_cos0 (rnd = %s):\n"
200                       "  Got sin =  ",
201                       mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
202               mpfr_dump (x);
203               printf ("  instead of -0.\n");
204               err = 1;
205             }
206           if (rnd == MPFR_RNDZ || rnd == MPFR_RNDD)
207             {
208               if (inex == 0)
209                 {
210                   printf ("Error in overflowed_sin_cos0 (rnd = %s):\n"
211                           "  The inexact value must be non-zero.\n",
212                           mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
213                   err = 1;
214                 }
215               if (! mpfr_equal_p (y, z))
216                 {
217                   printf ("Error in overflowed_sin_cos0 (rnd = %s):\n"
218                           "  Got cos =  ",
219                           mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
220                   mpfr_dump (y);
221                   printf ("  instead of 0.11111111E%d.\n", emax);
222                   err = 1;
223                 }
224             }
225           else
226             {
227               if (inex == 0)
228                 {
229                   printf ("Error in overflowed_sin_cos0 (rnd = %s):\n"
230                           "  The inexact value must be non-zero.\n",
231                           mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
232                   err = 1;
233                 }
234               if (! (mpfr_inf_p (y) && MPFR_IS_POS (y)))
235                 {
236                   printf ("Error in overflowed_sin_cos0 (rnd = %s):\n"
237                           "  Got cos =  ",
238                           mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
239                   mpfr_dump (y);
240                   printf ("  instead of +Inf.\n");
241                   err = 1;
242                 }
243             }
244         }
245       set_emax (old_emax);
246     }
247 
248   if (err)
249     exit (1);
250   mpfr_clear (x);
251   mpfr_clear (y);
252   mpfr_clear (z);
253 }
254 
255 static void
tiny(void)256 tiny (void)
257 {
258   mpfr_t x, s, c;
259   int i, inex;
260 
261   mpfr_inits2 (64, x, s, c, (mpfr_ptr) 0);
262 
263   for (i = -1; i <= 1; i += 2)
264     {
265       mpfr_set_si (x, i, MPFR_RNDN);
266       mpfr_set_exp (x, mpfr_get_emin ());
267       inex = mpfr_sin_cos (s, c, x, MPFR_RNDN);
268       MPFR_ASSERTN (inex != 0);
269       MPFR_ASSERTN (mpfr_equal_p (s, x));
270       MPFR_ASSERTN (!mpfr_nan_p (c) && mpfr_cmp_ui (c, 1) == 0);
271     }
272 
273   mpfr_clears (x, s, c, (mpfr_ptr) 0);
274 }
275 
276 /* bug found in nightly tests */
277 static void
test20071214(void)278 test20071214 (void)
279 {
280   mpfr_t a, b;
281   int inex;
282 
283   mpfr_init2 (a, 4);
284   mpfr_init2 (b, 4);
285 
286   mpfr_set_ui_2exp (a, 3, -4, MPFR_RNDN);
287   inex = mpfr_sin_cos (a, b, a, MPFR_RNDD);
288   MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 11, -6) == 0);
289   MPFR_ASSERTN(mpfr_cmp_ui_2exp (b, 15, -4) == 0);
290   MPFR_ASSERTN(inex == 10);
291 
292   mpfr_set_ui_2exp (a, 3, -4, MPFR_RNDN);
293   inex = mpfr_sin_cos (a, b, a, MPFR_RNDU);
294   MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 3, -4) == 0);
295   MPFR_ASSERTN(mpfr_cmp_ui (b, 1) == 0);
296   MPFR_ASSERTN(inex == 5);
297 
298   mpfr_set_ui_2exp (a, 3, -4, MPFR_RNDN);
299   inex = mpfr_sin_cos (a, b, a, MPFR_RNDN);
300   MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 3, -4) == 0);
301   MPFR_ASSERTN(mpfr_cmp_ui (b, 1) == 0);
302   MPFR_ASSERTN(inex == 5);
303 
304   mpfr_clear (a);
305   mpfr_clear (b);
306 }
307 
308 /* check that mpfr_sin_cos and test_mpfr_sincos_fast agree */
309 static void
test_mpfr_sincos_fast(void)310 test_mpfr_sincos_fast (void)
311 {
312   mpfr_t x, y, z, yref, zref, h;
313   mpfr_prec_t p = 1000;
314   int i, inex, inexref;
315   mpfr_rnd_t r;
316 
317   mpfr_init2 (x, p);
318   mpfr_init2 (y, p);
319   mpfr_init2 (z, p);
320   mpfr_init2 (yref, p);
321   mpfr_init2 (zref, p);
322   mpfr_init2 (h, p);
323   mpfr_set_ui (x, 0, MPFR_RNDN);
324   /* we generate a random value x, compute sin(x) and cos(x) with both
325      mpfr_sin_cos and mpfr_sincos_fast, and check the values and the flags
326      agree */
327   for (i = 0; i < 100; i++)
328     {
329       mpfr_urandomb (h, RANDS);
330       mpfr_add (x, x, h, MPFR_RNDN);
331       r = RND_RAND ();
332       inexref = mpfr_sin_cos (yref, zref, x, r);
333       inex = mpfr_sincos_fast (y, z, x, r);
334       if (mpfr_cmp (y, yref))
335         {
336           printf ("mpfr_sin_cos and mpfr_sincos_fast disagree\n");
337           printf ("x="); mpfr_dump (x);
338           printf ("rnd=%s\n", mpfr_print_rnd_mode (r));
339           printf ("yref="); mpfr_dump (yref);
340           printf ("y="); mpfr_dump (y);
341           exit (1);
342         }
343       if (mpfr_cmp (z, zref))
344         {
345           printf ("mpfr_sin_cos and mpfr_sincos_fast disagree\n");
346           printf ("x="); mpfr_dump (x);
347           printf ("rnd=%s\n", mpfr_print_rnd_mode (r));
348           printf ("zref="); mpfr_dump (zref);
349           printf ("z="); mpfr_dump (z);
350           exit (1);
351         }
352       if (inex != inexref)
353         {
354           printf ("mpfr_sin_cos and mpfr_sincos_fast disagree\n");
355           printf ("x="); mpfr_dump (x);
356           printf ("rnd=%s\n", mpfr_print_rnd_mode (r));
357           printf ("inexref=%d inex=%d\n", inexref, inex);
358           exit (1);
359         }
360     }
361   mpfr_clear (x);
362   mpfr_clear (y);
363   mpfr_clear (z);
364   mpfr_clear (yref);
365   mpfr_clear (zref);
366   mpfr_clear (h);
367 }
368 
369 static void
bug20091007(void)370 bug20091007 (void)
371 {
372   mpfr_t x, y, z, yref, zref;
373   mpfr_prec_t p = 1000;
374   int inex, inexref;
375   mpfr_rnd_t r = MPFR_RNDZ;
376 
377   mpfr_init2 (x, p);
378   mpfr_init2 (y, p);
379   mpfr_init2 (z, p);
380   mpfr_init2 (yref, p);
381   mpfr_init2 (zref, p);
382 
383   mpfr_set_str (x, "1.9ecdc22ba77a5ab2560f7e84289e2a328906f47377ea3fd4c82d1bb2f13ee05c032cffc1933eadab7b0a5498e03e3bd0508968e59c25829d97a0b54f20cd4662c8dfffa54e714de41fc8ee3e0e0b244d110a194db05b70022b7d77f88955d415b09f17dd404576098dc51a583a3e49c35839551646e880c7eb790a01a4@1", 16, MPFR_RNDN);
384   inexref = mpfr_sin_cos (yref, zref, x, r);
385   inex = mpfr_sincos_fast (y, z, x, r);
386 
387   if (mpfr_cmp (y, yref))
388     {
389       printf ("mpfr_sin_cos and mpfr_sincos_fast disagree (bug20091007)\n");
390       printf ("yref="); mpfr_dump (yref);
391       printf ("y="); mpfr_dump (y);
392       exit (1);
393     }
394   if (mpfr_cmp (z, zref))
395     {
396       printf ("mpfr_sin_cos and mpfr_sincos_fast disagree (bug20091007)\n");
397       printf ("zref="); mpfr_dump (zref);
398       printf ("z="); mpfr_dump (z);
399       exit (1);
400     }
401   if (inex != inexref)
402     {
403       printf ("mpfr_sin_cos and mpfr_sincos_fast disagree (bug20091007)\n");
404       printf ("inexref=%d inex=%d\n", inexref, inex);
405       exit (1);
406     }
407 
408   mpfr_clear (x);
409   mpfr_clear (y);
410   mpfr_clear (z);
411   mpfr_clear (yref);
412   mpfr_clear (zref);
413 }
414 
415 /* Note: with the sin_cos.c code before r6507, the disagreement occurs
416    only on the return ("inexact") value, which is new in r6444. */
417 static void
bug20091008(void)418 bug20091008 (void)
419 {
420   mpfr_t x, y, z, yref, zref;
421   mpfr_prec_t p = 1000;
422   int inex, inexref;
423   mpfr_rnd_t r = MPFR_RNDN;
424 
425   mpfr_init2 (x, p);
426   mpfr_init2 (y, p);
427   mpfr_init2 (z, p);
428   mpfr_init2 (yref, p);
429   mpfr_init2 (zref, p);
430 
431   mpfr_set_str (x, "c.91813724e28ef6a711d33e6505984699daef7fe93636c1ed5d0168bc96989cc6802f7f9e405c902ec62fb90cd39c9d21084c8ad8b5af4c4aa87bf402e2e4a78e6fe1ffeb6dbbbdbbc2983c196c518966ccc1e094ed39ee77984ef2428069d65de37928e75247edbe7007245e682616b5ebbf05f2fdefc74ad192024f10", 16, MPFR_RNDN);
432   inexref = mpfr_sin_cos (yref, zref, x, r);
433   inex = mpfr_sincos_fast (y, z, x, r);
434 
435   if (mpfr_cmp (y, yref))
436     {
437       printf ("mpfr_sin_cos and mpfr_sincos_fast disagree (bug20091008)\n");
438       printf ("yref="); mpfr_dump (yref);
439       printf ("y="); mpfr_dump (y);
440       exit (1);
441     }
442   if (mpfr_cmp (z, zref))
443     {
444       printf ("mpfr_sin_cos and mpfr_sincos_fast disagree (bug20091008)\n");
445       printf ("zref="); mpfr_dump (zref);
446       printf ("z="); mpfr_dump (z);
447       exit (1);
448     }
449   /* sin(x) is rounded up, cos(x) is rounded up too, thus we should get 5
450      for the return value */
451   if (inex != inexref)
452     {
453       printf ("mpfr_sin_cos and mpfr_sincos_fast disagree (bug20091008)\n");
454       printf ("inexref=%d inex=%d (5 expected)\n", inexref, inex);
455       exit (1);
456     }
457 
458   mpfr_clear (x);
459   mpfr_clear (y);
460   mpfr_clear (z);
461   mpfr_clear (yref);
462   mpfr_clear (zref);
463 }
464 
465 static void
bug20091013(void)466 bug20091013 (void)
467 {
468   mpfr_t x, y, z, yref, zref;
469   mpfr_prec_t p = 1000;
470   int inex, inexref;
471   mpfr_rnd_t r = MPFR_RNDN;
472 
473   mpfr_init2 (x, p);
474   mpfr_init2 (y, p);
475   mpfr_init2 (z, p);
476   mpfr_init2 (yref, p);
477   mpfr_init2 (zref, p);
478 
479   mpfr_set_str (x, "3.240ff3fdcb1ee7cd667b96287593ae24e20fb63ed7c2d5bf4bd0f2cc5509283b04e7628e66382605f14ed5967cef15296041539a1bdaa626c777c7fbb6f2068414759b78cee14f37848689b3a170f583656be4e0837f464d8210556a3a822d4ecfdd59f4e0d5fdb76bf7e15b8a57234e2160b98e14c17bbdf27c4643b8@1", 16, MPFR_RNDN);
480   inexref = mpfr_sin_cos (yref, zref, x, r);
481   inex = mpfr_sincos_fast (y, z, x, r);
482 
483   if (mpfr_cmp (y, yref))
484     {
485       printf ("mpfr_sin_cos and mpfr_sincos_fast disagree (bug20091013)\n");
486       printf ("yref="); mpfr_dump (yref);
487       printf ("y="); mpfr_dump (y);
488       exit (1);
489     }
490   if (mpfr_cmp (z, zref))
491     {
492       printf ("mpfr_sin_cos and mpfr_sincos_fast disagree (bug20091013)\n");
493       printf ("zref="); mpfr_dump (zref);
494       printf ("z="); mpfr_dump (z);
495       exit (1);
496     }
497   /* sin(x) is rounded down and cos(x) is rounded down, thus we should get
498      2+4*2 = 10 as return value */
499   if (inex != inexref)
500     {
501       printf ("mpfr_sin_cos and mpfr_sincos_fast disagree (bug20091013)\n");
502       printf ("inexref=%d inex=%d (10 expected)\n", inexref, inex);
503       exit (1);
504     }
505 
506   mpfr_clear (x);
507   mpfr_clear (y);
508   mpfr_clear (z);
509   mpfr_clear (yref);
510   mpfr_clear (zref);
511 }
512 
513 /* Bug reported by Laurent Fousse for the 2.4 branch.
514    No problem in the trunk.
515    https://sympa.inria.fr/sympa/arc/mpfr/2009-11/msg00044.html */
516 static void
bug20091122(void)517 bug20091122 (void)
518 {
519   mpfr_t x, y, z, yref, zref;
520   mpfr_prec_t p = 3;
521   mpfr_rnd_t r = MPFR_RNDN;
522 
523   mpfr_init2 (x, 5);
524   mpfr_init2 (y, p);
525   mpfr_init2 (z, p);
526   mpfr_init2 (yref, p);
527   mpfr_init2 (zref, p);
528 
529   mpfr_set_str (x, "0.11111E49", 2, MPFR_RNDN);
530   mpfr_sin_cos (yref, zref, x, r);
531 
532   mpfr_sin (y, x, r);
533   mpfr_cos (z, x, r);
534 
535   if (! mpfr_equal_p (y, yref))
536     {
537       printf ("mpfr_sin_cos and mpfr_sin disagree (bug20091122)\n");
538       printf ("yref = "); mpfr_dump (yref);
539       printf ("y    = "); mpfr_dump (y);
540       exit (1);
541     }
542   if (! mpfr_equal_p (z, zref))
543     {
544       printf ("mpfr_sin_cos and mpfr_cos disagree (bug20091122)\n");
545       printf ("zref = "); mpfr_dump (zref);
546       printf ("z    = "); mpfr_dump (z);
547       exit (1);
548     }
549 
550   mpfr_clear (x);
551   mpfr_clear (y);
552   mpfr_clear (z);
553   mpfr_clear (yref);
554   mpfr_clear (zref);
555 }
556 
557 static void
consistency(void)558 consistency (void)
559 {
560   mpfr_t x, s1, s2, c1, c2;
561   mpfr_exp_t emin, emax;
562   mpfr_rnd_t rnd;
563   unsigned int flags_sin, flags_cos, flags, flags_before, flags_ref;
564   int inex_sin, is, inex_cos, ic, inex, inex_ref;
565   int i;
566 
567   emin = mpfr_get_emin ();
568   emax = mpfr_get_emax ();
569 
570   for (i = 0; i <= 10000; i++)
571     {
572       mpfr_init2 (x, MPFR_PREC_MIN + (randlimb () % 8));
573       mpfr_inits2 (MPFR_PREC_MIN + (randlimb () % 8), s1, s2, c1, c2,
574                    (mpfr_ptr) 0);
575       if (i < 8 * MPFR_RND_MAX)
576         {
577           int j = i / MPFR_RND_MAX;
578           if (j & 1)
579             set_emin (MPFR_EMIN_MIN);
580           mpfr_set_si (x, (j & 2) ? 1 : -1, MPFR_RNDN);
581           mpfr_set_exp (x, mpfr_get_emin ());
582           rnd = (mpfr_rnd_t) (i % MPFR_RND_MAX);
583           if (rnd == MPFR_RNDF)
584             goto end;
585           flags_before = 0;
586           if (j & 4)
587             set_emax (-17);
588         }
589       else
590         {
591           tests_default_random (x, 256, -5, 50, 0);
592           rnd = RND_RAND_NO_RNDF ();
593           flags_before = RAND_BOOL () ?
594             (unsigned int) (MPFR_FLAGS_ALL ^ MPFR_FLAGS_ERANGE) :
595             (unsigned int) 0;
596         }
597       __gmpfr_flags = flags_before;
598       inex_sin = mpfr_sin (s1, x, rnd);
599       is = inex_sin < 0 ? 2 : inex_sin > 0 ? 1 : 0;
600       flags_sin = __gmpfr_flags;
601       __gmpfr_flags = flags_before;
602       inex_cos = mpfr_cos (c1, x, rnd);
603       ic = inex_cos < 0 ? 2 : inex_cos > 0 ? 1 : 0;
604       flags_cos = __gmpfr_flags;
605       __gmpfr_flags = flags_before;
606       inex = mpfr_sin_cos (s2, c2, x, rnd);
607       flags = __gmpfr_flags;
608       inex_ref = is + 4 * ic;
609       flags_ref = flags_sin | flags_cos;
610       if (!(mpfr_equal_p (s1, s2) && mpfr_equal_p (c1, c2)) ||
611           inex != inex_ref || flags != flags_ref)
612         {
613           printf ("mpfr_sin_cos and mpfr_sin/mpfr_cos disagree on %s,"
614                   " i = %d\nx = ", mpfr_print_rnd_mode (rnd), i);
615           mpfr_dump (x);
616           printf ("s1 = ");
617           mpfr_dump (s1);
618           printf ("s2 = ");
619           mpfr_dump (s2);
620           printf ("c1 = ");
621           mpfr_dump (c1);
622           printf ("c2 = ");
623           mpfr_dump (c2);
624           printf ("inex_sin = %d (s = %d), inex_cos = %d (c = %d), "
625                   "inex = %d (expected %d)\n",
626                   inex_sin, is, inex_cos, ic, inex, inex_ref);
627           printf ("flags_sin = 0x%x, flags_cos = 0x%x, "
628                   "flags = 0x%x (expected 0x%x)\n",
629                   flags_sin, flags_cos, flags, flags_ref);
630           exit (1);
631         }
632     end:
633       mpfr_clears (x, s1, s2, c1, c2, (mpfr_ptr) 0);
634       set_emin (emin);
635       set_emax (emax);
636     }
637 }
638 
639 static void
coverage_01032011(void)640 coverage_01032011 (void)
641 {
642   mpfr_t val, cval, sval, svalf;
643   int status_f, status;
644 
645   mpfr_init2 (val, MPFR_PREC_MIN);
646   mpfr_init2 (cval, MPFR_PREC_MIN);
647   mpfr_init2 (sval, MPFR_PREC_MIN);
648   mpfr_init2 (svalf, MPFR_PREC_MIN);
649 
650   mpfr_set_str1 (val, "-0.7");
651 
652   status_f = mpfr_sincos_fast (svalf, NULL, val, MPFR_RNDN);
653   status = mpfr_sin_cos (sval, cval, val, MPFR_RNDN);
654   if (! mpfr_equal_p (svalf, sval) || VSIGN (status_f) != VSIGN (status))
655     {
656       printf ("mpfr_sincos_fast differ from mpfr_sin_cos result:\n"
657               " sin fast is ");
658       mpfr_dump (svalf);
659       printf (" sin is ");
660       mpfr_dump (sval);
661       printf ("status_f = %d, status = %d\n", status_f, status);
662       exit (1);
663     }
664 
665   mpfr_clears(val, cval, sval, svalf, (mpfr_ptr) 0);
666 }
667 
668 /* tsin_cos prec [N] performs N tests with prec bits */
669 int
main(int argc,char * argv[])670 main (int argc, char *argv[])
671 {
672   tests_start_mpfr ();
673 
674   if (argc > 1)
675     {
676       if (argc != 3 && argc != 4)
677         {
678           fprintf (stderr, "Usage: tsin_cos x prec [n]\n");
679           exit (1);
680         }
681       large_test (argv[1], atoi (argv[2]), (argc > 3) ? atoi (argv[3]) : 1);
682       goto end;
683     }
684 
685   bug20091013 ();
686   bug20091008 ();
687   bug20091007 ();
688   bug20091122 ();
689   consistency ();
690 
691   test_mpfr_sincos_fast ();
692 
693   check_nans ();
694 
695   /* worst case from PhD thesis of Vincent Lefe`vre: x=8980155785351021/2^54 */
696   check53 ("4.984987858808754279e-1", "4.781075595393330379e-1",
697            "8.783012931285841817e-1", MPFR_RNDN);
698   check53 ("4.984987858808754279e-1", "4.781075595393329824e-1",
699            "8.783012931285840707e-1", MPFR_RNDD);
700   check53 ("4.984987858808754279e-1", "4.781075595393329824e-1",
701            "8.783012931285840707e-1", MPFR_RNDZ);
702   check53 ("4.984987858808754279e-1", "4.781075595393330379e-1",
703            "8.783012931285841817e-1", MPFR_RNDU);
704   check53 ("1.00031274099908640274",  "8.416399183372403892e-1",
705            "0.540039116973283217504", MPFR_RNDN);
706   check53 ("1.00229256850978698523",  "8.427074524447979442e-1",
707            "0.538371757797526551137", MPFR_RNDZ);
708   check53 ("1.00288304857059840103",  "8.430252033025980029e-1",
709            "0.537874062022526966409", MPFR_RNDZ);
710   check53 ("1.00591265847407274059",  "8.446508805292128885e-1",
711            "0.53531755997839769456",  MPFR_RNDN);
712 
713   /* check one argument only */
714   check53sin ("1.00591265847407274059", "8.446508805292128885e-1", MPFR_RNDN);
715   check53cos ("1.00591265847407274059", "0.53531755997839769456",  MPFR_RNDN);
716 
717   overflowed_sin_cos0 ();
718   tiny ();
719   test20071214 ();
720 
721   coverage_01032011 ();
722 
723  end:
724   tests_end_mpfr ();
725   return 0;
726 }
727