sysdeps/ieee754/dbl-64/s_sincos.c

   1 /* Compute sine and cosine of argument.
   2    Copyright (C) 1997-2024 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