newlib/libm/mathfp/s_sine.c

   1
   2 /* @(#)z_sine.c 1.0 98/08/13 */
   3 /******************************************************************
   4  * The following routines are coded directly from the algorithms
   5  * and coefficients given in "Software Manual for the Elementary
   6  * Functions" by William J. Cody, Jr. and William Waite, Prentice
   7  * Hall, 1980.
   8  ******************************************************************/
   9
  10 /*
  11 FUNCTION
  12         <<sin>>, <<cos>>, <<sine>>, <<sinf>>, <<cosf>>, <<sinef>>---sine or cosine
  13 INDEX
  14 sin
  15 INDEX
  16 sinf
  17 INDEX
  18 cos
  19 INDEX
  20 cosf
  21 SYNOPSIS
  22         #include <math.h>
  23         double sin(double <[x]>);
  24         float  sinf(float <[x]>);
  25         double cos(double <[x]>);
  26         float cosf(float <[x]>);
  27
  28 DESCRIPTION
  29         <<sin>> and <<cos>> compute (respectively) the sine and cosine
  30         of the argument <[x]>.  Angles are specified in radians.
  31 RETURNS
  32         The sine or cosine of <[x]> is returned.
  33
  34 PORTABILITY
  35         <<sin>> and <<cos>> are ANSI C.
  36         <<sinf>> and <<cosf>> are extensions.
  37
  38 QUICKREF
  39         sin ansi pure
  40         sinf - pure
  41 */
  42
  43 /******************************************************************
  44  * sine
  45  *
  46  * Input:
  47  *   x - floating point value
  48  *   cosine - indicates cosine value
  49  *
  50  * Output:
  51  *   Sine of x.
  52  *
  53  * Description:
  54  *   This routine calculates sines and cosines.
  55  *
  56  *****************************************************************/
  57
  58 #include "fdlibm.h"
  59 #include "zmath.h"
  60
  61 #ifndef _DOUBLE_IS_32BITS
  62
  63 static const double HALF_PI = 1.57079632679489661923;
  64 static const double ONE_OVER_PI = 0.31830988618379067154;
  65 static const double r[] = { -0.16666666666666665052,
  66                              0.83333333333331650314e-02,
  67                             -0.19841269841201840457e-03,
  68                              0.27557319210152756119e-05,
  69                             -0.25052106798274584544e-07,
  70                              0.16058936490371589114e-09,
  71                             -0.76429178068910467734e-12,
  72                              0.27204790957888846175e-14 };
  73
  74 double
  75 sine (double x,
  76         int cosine)
  77 {
  78   int sgn, N;
  79   double y, XN, g, R, res;
  80   double YMAX = 210828714.0;
  81
  82   switch (numtest (x))
  83     {
  84       case NAN:
  85         errno = EDOM;
  86         return (x);
  87       case INF:
  88         errno = EDOM;
  89         return (z_notanum.d);
  90     }
  91
  92   /* Use sin and cos properties to ease computations. */
  93   if (cosine)
  94     {
  95       sgn = 1;
  96       y = fabs (x) + HALF_PI;
  97     }
  98   else
  99     {
 100       if (x < 0.0)
 101         {
 102           sgn = -1;
 103           y = -x;
 104         }
 105       else
 106         {
 107           sgn = 1;
 108           y = x;
 109         }
 110     }
 111
 112   /* Check for values of y that will overflow here. */
 113   if (y > YMAX)
 114     {
 115       errno = ERANGE;
 116       return (x);
 117     }
 118
 119   /* Calculate the exponent. */
 120   if (y < 0.0)
 121     N = (int) (y * ONE_OVER_PI - 0.5);
 122   else
 123     N = (int) (y * ONE_OVER_PI + 0.5);
 124   XN = (double) N;
 125
 126   if (N & 1)
 127     sgn = -sgn;
 128
 129   if (cosine)
 130     XN -= 0.5;
 131
 132   y = fabs (x) - XN * __PI;
 133
 134   if (-z_rooteps < y && y < z_rooteps)
 135     res = y;
 136
 137   else
 138     {
 139       g = y * y;
 140
 141       /* Calculate the Taylor series. */
 142       R = (((((((r[6] * g + r[5]) * g + r[4]) * g + r[3]) * g + r[2]) * g + r[1]) * g + r[0]) * g);
 143
 144       /* Finally, compute the result. */
 145       res = y + y * R;
 146     }
 147
 148   res *= sgn;
 149
 150   return (res);
 151 }
 152
 153 #endif /* _DOUBLE_IS_32BITS */