home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_10_12 / 1012084c < prev    next >
Text File  |  1992-10-13  |  3KB  |  74 lines
  1. #include <errno.h>
  2. #include "math.h"
  3. #include <float.h>
  4. #include "xmath.h"
  5. #define twobypi .63661977236758134308
  6.  
  7. double _Sin(x, qoff)
  8. double x;
  9. unsigned int qoff;
  10. {                               /* sin(x) or cos(x) */
  11.     double g, g4, g2;
  12.     float gr = 1 / DBL_EPSILON;
  13.     int quad;
  14.  
  15.     if (FLT_ROUNDS > 0) {
  16.         /* Compiler should eliminate dead code if
  17.          * FLT_ROUNDS is compile time constant */
  18.         if ((g = x) < 0)
  19.             gr = -gr;
  20.         if (FLT_ROUNDS == 1)
  21.             g = (g * twobypi + gr) - gr;
  22.         /* ANINT(x*2/pi) */
  23.         else
  24.             g = ((g * twobypi + (FLT_ROUNDS == 2 ?
  25.                         -.5 : .5)) + gr) - gr;
  26.     } else {
  27.         /* assume FLT_ROUNDS ==0 for positive numbers */
  28.         g = x * twobypi;
  29.         g = x < 0 ?
  30.             gr - ((.5 - g) + gr) : ((g + .5) + gr) - gr;
  31.     }
  32. #if 0
  33.     /* Reduce g to [-3,3] by ignoring reductions by 2pi
  34.      * when casting to int, so there won't ever be
  35.      * overflow; if the argument is excessively large, it
  36.      * will be treated as reducing to first or last octant
  37.      * as there is not enough precision remaining to
  38.      * distinguish octants */
  39.     gr *= 4;
  40.     /* qoff = qoff + (g-ANINT(g*4)) */
  41.     if (fabs(g) >= 1 / DBL_EPSILON)
  42.         errno = ERANGE;
  43.     quad = g - (FLT_ROUNDS > 0 || g >= 0 ?
  44.         (g + gr) - gr : gr - (gr - g));
  45. #else
  46. /* much faster on SPARC MIPS and HP but could overflow */
  47.     quad = ((int) g) & 3;
  48. #endif
  49.     /* Get remainder after subtracting nearest multiple of
  50.      * pi/2 */
  51.         g = (x - g * (3294198. / 2097152.)) - g *
  52.             3.139164786504813217e-7;
  53. /* Now calculate +- sin() or cos(), -Pi/4 <= x <= Pi/4 */
  54.     g2 = g * g;
  55.     g4 = g2 * g2;
  56.     if ((qoff += quad) & 1)     /* cosine series */
  57.         g = 1 + g2 * (-.499999999999999994 + g2 *
  58.             (.041666666666666452 + g2 *
  59.                 (-.001388888888886110 + g2 *
  60.                     .000024801587283884))
  61.             + g4 * g4 * (-.000000275573130985 + g2 *
  62.                 (.000000002087558246 - g2 *
  63.                     .000000000011353383)));
  64.     else                        /* sine series */
  65.         g += g * g2 * (-.16666666666666616 + g2 *
  66.             (.00833333333332036 + g2 *
  67.                 (-.00019841269828653 + g2 *
  68.                     .0000027557313377252))
  69.             + g4 * g4 * (-.000000025050717097 + g2 *
  70.                 .00000000015894743));
  71.     return qoff & 2? -g: g;
  72. }
  73. WRAP_EOF
  74.