newlib/libm/mathfp/sf_sine.c

   1
   2 /* @(#)z_sinef.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  * sine generator
  11  *
  12  * Input:
  13  *   x - floating point value
  14  *   cosine - indicates cosine value
  15  *
  16  * Output:
  17  *   Sine of x.
  18  *
  19  * Description:
  20  *   This routine calculates sines and cosines.
  21  *
  22  *****************************************************************/
  23
  24 #include "fdlibm.h"
  25 #include "zmath.h"
  26
  27 static const float HALF_PI = 1.570796326;
  28 static const float ONE_OVER_PI = 0.318309886;
  29 static const float r[] = { -0.1666665668,
  30                             0.8333025139e-02,
  31                            -0.1980741872e-03,
  32                             0.2601903036e-5 };
  33
  34 float
  35 sinef (float x,
  36         int cosine)
  37 {
  38   int sgn, N;
  39   float y, XN, g, R, res;
  40   float YMAX = 210828714.0;
  41
  42   switch (numtestf (x))
  43     {
  44       case NAN:
  45         errno = EDOM;
  46         return (x);
  47       case INF:
  48         errno = EDOM;
  49         return (z_notanum_f.f);
  50     }
  51
  52   /* Use sin and cos properties to ease computations. */
  53   if (cosine)
  54     {
  55       sgn = 1;
  56       y = fabsf (x) + HALF_PI;
  57     }
  58   else
  59     {
  60       if (x < 0.0)
  61         {
  62           sgn = -1;
  63           y = -x;
  64         }
  65       else
  66         {
  67           sgn = 1;
  68           y = x;
  69         }
  70     }
  71
  72   /* Check for values of y that will overflow here. */
  73   if (y > YMAX)
  74     {
  75       errno = ERANGE;
  76       return (x);
  77     }
  78
  79   /* Calculate the exponent. */
  80   if (y < 0.0)
  81     N = (int) (y * ONE_OVER_PI - 0.5);
  82   else
  83     N = (int) (y * ONE_OVER_PI + 0.5);
  84   XN = (float) N;
  85
  86   if (N & 1)
  87     sgn = -sgn;
  88
  89   if (cosine)
  90     XN -= 0.5;
  91
  92   y = fabsf (x) - XN * __PI;
  93
  94   if (-z_rooteps_f < y && y < z_rooteps_f)
  95     res = y;
  96
  97   else
  98     {
  99       g = y * y;
 100
 101       /* Calculate the Taylor series. */
 102       R = (((r[3] * g + r[2]) * g + r[1]) * g + r[0]) * g;
 103
 104       /* Finally, compute the result. */
 105       res = y + y * R;
 106     }
 107
 108   res *= sgn;
 109
 110   return (res);
 111 }