lib/math/sin.c

   1 /*
   2  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
   3  * See the copyright notice in the ACK home directory, in the file "Copyright".
   4  *
   5  * Author: Ceriel J.H. Jacobs
   6  */
   7 /* $Header$ */
   8
   9 #include        <math.h>
  10 #include        <float.h>
  11 #include        <errno.h>
  12 #include        "localmath.h"
  13
  14 static double
  15 sinus(double x, int cos_flag)
  16 {
  17         /*      Algorithm and coefficients from:
  18                         "Software manual for the elementary functions"
  19                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
  20         */
  21
  22         static double r[] = {
  23                 -0.16666666666666665052e+0,
  24                  0.83333333333331650314e-2,
  25                 -0.19841269841201840457e-3,
  26                  0.27557319210152756119e-5,
  27                 -0.25052106798274584544e-7,
  28                  0.16058936490371589114e-9,
  29                 -0.76429178068910467734e-12,
  30                  0.27204790957888846175e-14
  31         };
  32
  33         double  y;
  34         int     neg = 1;
  35
  36         if (__IsNan(x)) {
  37                 errno = EDOM;
  38                 return x;
  39         }
  40         if (x < 0) {
  41                 x = -x;
  42                 neg = -1;
  43         }
  44         if (cos_flag) {
  45                 neg = 1;
  46                 y = M_PI_2 + x;
  47         }
  48         else    y = x;
  49
  50         /* ??? avoid loss of significance, if y is too large, error ??? */
  51
  52         y = y * M_1_PI + 0.5;
  53
  54         if (y >= DBL_MAX/M_PI) return 0.0;
  55
  56         /*      Use extended precision to calculate reduced argument.
  57                 Here we used 12 bits of the mantissa for a1.
  58                 Also split x in integer part x1 and fraction part x2.
  59         */
  60 #define A1 3.1416015625
  61 #define A2 -8.908910206761537356617e-6
  62         {
  63                 double x1, x2;
  64
  65                 modf(y, &y);
  66                 if (modf(0.5*y, &x1)) neg = -neg;
  67                 if (cos_flag) y -= 0.5;
  68                 x2 = modf(x, &x1);
  69                 x = x1 - y * A1;
  70                 x += x2;
  71                 x -= y * A2;
  72 #undef A1
  73 #undef A2
  74         }
  75
  76         if (x < 0) {
  77                 neg = -neg;
  78                 x = -x;
  79         }
  80
  81         /* ??? avoid underflow ??? */
  82
  83         y = x * x;
  84         x += x * y * POLYNOM7(y, r);
  85         return neg==-1 ? -x : x;
  86 }
  87
  88 double
  89 sin(double x)
  90 {
  91         return sinus(x, 0);
  92 }
  93
  94 double
  95 cos(double x)
  96 {
  97         if (x < 0) x = -x;
  98         return sinus(x, 1);
  99 }