lib/ack/libp/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
   8 /* $Header$ */
   9
  10 #define __NO_DEFS
  11 #include <math.h>
  12
  13 #if __STDC__
  14 #include <pc_math.h>
  15 #endif
  16
  17 static double
  18 sinus(x, cos_flag)
  19         double x;
  20 {
  21         /*      Algorithm and coefficients from:
  22                         "Software manual for the elementary functions"
  23                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
  24         */
  25
  26         static double r[] = {
  27                 -0.16666666666666665052e+0,
  28                  0.83333333333331650314e-2,
  29                 -0.19841269841201840457e-3,
  30                  0.27557319210152756119e-5,
  31                 -0.25052106798274584544e-7,
  32                  0.16058936490371589114e-9,
  33                 -0.76429178068910467734e-12,
  34                  0.27204790957888846175e-14
  35         };
  36
  37         double  xsqr;
  38         double  y;
  39         int     neg = 0;
  40
  41         if (x < 0) {
  42                 x = -x;
  43                 neg = 1;
  44         }
  45         if (cos_flag) {
  46                 neg = 0;
  47                 y = M_PI_2 + x;
  48         }
  49         else    y = x;
  50
  51         /* ??? avoid loss of significance, if y is too large, error ??? */
  52
  53         y = y * M_1_PI + 0.5;
  54
  55         /*      Use extended precision to calculate reduced argument.
  56                 Here we used 12 bits of the mantissa for a1.
  57                 Also split x in integer part x1 and fraction part x2.
  58         */
  59 #define A1 3.1416015625
  60 #define A2 -8.908910206761537356617e-6
  61         {
  62                 double x1, x2;
  63                 extern double   _fif();
  64
  65                 _fif(y, 1.0,  &y);
  66                 if (_fif(y, 0.5, &x1)) neg = !neg;
  67                 if (cos_flag) y -= 0.5;
  68                 x2 = _fif(x, 1.0, &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 ? -x : x;
  86 }
  87
  88 double
  89 _sin(x)
  90         double x;
  91 {
  92         return sinus(x, 0);
  93 }
  94
  95 double
  96 _cos(x)
  97         double x;
  98 {
  99         if (x < 0) x = -x;
 100         return sinus(x, 1);
 101 }