home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cmath / sin.c < prev    next >
Encoding:
C/C++ Source or Header  |  1983-07-02  |  1.4 KB  |  75 lines
  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. static double twoopi    = 0.63661977236758134308;
  9. static double p0    =  .1357884097877375669092680e8;
  10. static double p1    = -.4942908100902844161158627e7;
  11. static double p2    =  .4401030535375266501944918e6;
  12. static double p3    = -.1384727249982452873054457e5;
  13. static double p4    =  .1459688406665768722226959e3;
  14. static double q0    =  .8644558652922534429915149e7;
  15. static double q1    =  .4081792252343299749395779e6;
  16. static double q2    =  .9463096101538208180571257e4;
  17. static double q3    =  .1326534908786136358911494e3;
  18.  
  19. double
  20. cos(arg)
  21. double arg;
  22. {
  23.     double sinus();
  24.     if(arg<0)
  25.         arg = -arg;
  26.     return(sinus(arg, 1));
  27. }
  28.  
  29. double
  30. sin(arg)
  31. double arg;
  32. {
  33.     double sinus();
  34.     return(sinus(arg, 0));
  35. }
  36.  
  37. static double
  38. sinus(arg, quad)
  39. double arg;
  40. int quad;
  41. {
  42.     double modf();
  43.     double e, f;
  44.     double ysq;
  45.     double x,y;
  46.     int k;
  47.     double temp1, temp2;
  48.  
  49.     x = arg;
  50.     if(x<0) {
  51.         x = -x;
  52.         quad = quad + 2;
  53.     }
  54.     x = x*twoopi;    /*underflow?*/
  55.     if(x>32764){
  56.         y = modf(x,&e);
  57.         e = e + quad;
  58.         modf(0.25*e,&f);
  59.         quad = e - 4*f;
  60.     }else{
  61.         k = x;
  62.         y = x - k;
  63.         quad = (quad + k) & 03;
  64.     }
  65.     if (quad & 01)
  66.         y = 1-y;
  67.     if(quad > 1)
  68.         y = -y;
  69.  
  70.     ysq = y*y;
  71.     temp1 = ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y;
  72.     temp2 = ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0);
  73.     return(temp1/temp2);
  74. }
  75.