1  /* Compute sine and cosine of argument.
       2     Copyright (C) 1997-2023 Free Software Foundation, Inc.
       3     This file is part of the GNU C Library.
       4  
       5     The GNU C Library is free software; you can redistribute it and/or
       6     modify it under the terms of the GNU Lesser General Public
       7     License as published by the Free Software Foundation; either
       8     version 2.1 of the License, or (at your option) any later version.
       9  
      10     The GNU C Library is distributed in the hope that it will be useful,
      11     but WITHOUT ANY WARRANTY; without even the implied warranty of
      12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      13     Lesser General Public License for more details.
      14  
      15     You should have received a copy of the GNU Lesser General Public
      16     License along with the GNU C Library; if not, see
      17     <https://www.gnu.org/licenses/>.  */
      18  
      19  #include <errno.h>
      20  #include <math.h>
      21  
      22  #include <math_private.h>
      23  #include <fenv_private.h>
      24  #include <math-underflow.h>
      25  #include <libm-alias-double.h>
      26  
      27  #ifndef SECTION
      28  # define SECTION
      29  #endif
      30  
      31  #define IN_SINCOS
      32  #include "s_sin.c"
      33  
      34  void
      35  SECTION
      36  __sincos (double x, double *sinx, double *cosx)
      37  {
      38    mynumber u;
      39    int k;
      40  
      41    SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
      42  
      43    u.x = x;
      44    k = u.i[HIGH_HALF] & 0x7fffffff;
      45  
      46    if (k < 0x400368fd)
      47      {
      48        double a, da, y;
      49        /* |x| < 2^-27 => cos (x) = 1, sin (x) = x.  */
      50        if (k < 0x3e400000)
      51  	{
      52  	  if (k < 0x3e500000)
      53  	    math_check_force_underflow (x);
      54  	  *sinx = x;
      55  	  *cosx = 1.0;
      56  	  return;
      57  	}
      58        /* |x| < 0.855469.  */
      59        else if (k < 0x3feb6000)
      60  	{
      61  	  *sinx = do_sin (x, 0);
      62  	  *cosx = do_cos (x, 0);
      63  	  return;
      64  	}
      65  
      66        /* |x| < 2.426265.  */
      67        y = hp0 - fabs (x);
      68        a = y + hp1;
      69        da = (y - a) + hp1;
      70        *sinx = copysign (do_cos (a, da), x);
      71        *cosx = do_sin (a, da);
      72        return;
      73      }
      74    /* |x| < 2^1024.  */
      75    if (k < 0x7ff00000)
      76      {
      77        double a, da, xx;
      78        unsigned int n;
      79  
      80        /* If |x| < 105414350 use simple range reduction.  */
      81        n = k < 0x419921FB ? reduce_sincos (x, &a, &da) : __branred (x, &a, &da);
      82        n = n & 3;
      83  
      84        if (n == 1 || n == 2)
      85  	{
      86  	  a = -a;
      87  	  da = -da;
      88  	}
      89  
      90        if (n & 1)
      91  	{
      92  	  double *temp = cosx;
      93  	  cosx = sinx;
      94  	  sinx = temp;
      95  	}
      96  
      97        *sinx = do_sin (a, da);
      98        xx = do_cos (a, da);
      99        *cosx = (n & 2) ? -xx : xx;
     100        return;
     101      }
     102  
     103    if (isinf (x))
     104      __set_errno (EDOM);
     105  
     106    *sinx = *cosx = x / x;
     107  }
     108  #ifndef __sincos
     109  libm_alias_double (__sincos, sincos)
     110  #endif