lib/msun/src/s_sincos.c

   1 /*-
   2  * ====================================================
   3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   4  *
   5  * Developed at SunPro, a Sun Microsystems, Inc. business.
   6  * Permission to use, copy, modify, and distribute this
   7  * software is freely granted, provided that this notice
   8  * is preserved.
   9  * ====================================================
  10  *
  11  * s_sin.c and s_cos.c merged by Steven G. Kargl.  Descriptions of the
  12  * algorithms are contained in the original files.
  13  */
  14
  15 #include <sys/cdefs.h>
  16 #include <float.h>
  17
  18 #include "math.h"
  19 #define INLINE_REM_PIO2
  20 #include "math_private.h"
  21 #include "e_rem_pio2.c"
  22 #include "k_sincos.h"
  23
  24 void
  25 sincos(double x, double *sn, double *cs)
  26 {
  27         double y[2];
  28         int32_t n, ix;
  29
  30         /* High word of x. */
  31         GET_HIGH_WORD(ix, x);
  32
  33         /* |x| ~< pi/4 */
  34         ix &= 0x7fffffff;
  35         if (ix <= 0x3fe921fb) {
  36                 if (ix < 0x3e400000) {          /* |x| < 2**-27 */
  37                         if ((int)x == 0) {      /* Generate inexact. */
  38                                 *sn = x;
  39                                 *cs = 1;
  40                                 return;
  41                         }
  42                 }
  43                 __kernel_sincos(x, 0, 0, sn, cs);
  44                 return;
  45         }
  46
  47         /* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
  48         if (ix >= 0x7ff00000) {
  49                 *sn = x - x;
  50                 *cs = x - x;
  51                 return;
  52         }
  53
  54         /* Argument reduction. */
  55         n = __ieee754_rem_pio2(x, y);
  56
  57         switch(n & 3) {
  58         case 0:
  59                 __kernel_sincos(y[0], y[1], 1, sn, cs);
  60                 break;
  61         case 1:
  62                 __kernel_sincos(y[0], y[1], 1, cs, sn);
  63                 *cs = -*cs;
  64                 break;
  65         case 2:
  66                 __kernel_sincos(y[0], y[1], 1, sn, cs);
  67                 *sn = -*sn;
  68                 *cs = -*cs;
  69                 break;
  70         default:
  71                 __kernel_sincos(y[0], y[1], 1, cs, sn);
  72                 *sn = -*sn;
  73         }
  74 }
  75
  76 #if (LDBL_MANT_DIG == 53)
  77 __weak_reference(sincos, sincosl);
  78 #endif