sys/src/libc/port/sin.c

   1 /*
   2         C program for floating point sin/cos.
   3         Calls modf.
   4         There are no error exits.
   5         Coefficients are #3370 from Hart & Cheney (18.80D).
   6 */
   7
   8 #include <u.h>
   9 #include <libc.h>
  10
  11 #define p0      .1357884097877375669092680e8
  12 #define p1     -.4942908100902844161158627e7
  13 #define p2      .4401030535375266501944918e6
  14 #define p3     -.1384727249982452873054457e5
  15 #define p4      .1459688406665768722226959e3
  16 #define q0      .8644558652922534429915149e7
  17 #define q1      .4081792252343299749395779e6
  18 #define q2      .9463096101538208180571257e4
  19 #define q3      .1326534908786136358911494e3
  20
  21 static
  22 double
  23 sinus(double arg, int quad)
  24 {
  25         double e, f, ysq, x, y, temp1, temp2;
  26         int k;
  27
  28         x = arg;
  29         if(x < 0) {
  30                 x = -x;
  31                 quad += 2;
  32         }
  33         x *= 1/PIO2;    /* underflow? */
  34         if(x > 32764) {
  35                 y = modf(x, &e);
  36                 e += quad;
  37                 modf(0.25*e, &f);
  38                 quad = e - 4*f;
  39         } else {
  40                 k = x;
  41                 y = x - k;
  42                 quad += k;
  43                 quad &= 3;
  44         }
  45         if(quad & 1)
  46                 y = 1-y;
  47         if(quad > 1)
  48                 y = -y;
  49
  50         ysq = y*y;
  51         temp1 = ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y;
  52         temp2 = ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0);
  53         return temp1/temp2;
  54 }
  55
  56 double
  57 cos(double arg)
  58 {
  59         if(arg < 0)
  60                 arg = -arg;
  61         return sinus(arg, 1);
  62 }
  63
  64 double
  65 sin(double arg)
  66 {
  67         return sinus(arg, 0);
  68 }