(root)/
mpfr-4.2.1/
tests/
tsin_cos.c
       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
      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
      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
      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
      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
     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
     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
     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
     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
     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
     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
     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
     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
     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
     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
     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
     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  }