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 <float.h>
  16
  17 #include "math.h"
  18 #define INLINE_REM_PIO2
  19 #include "math_private.h"
  20 #include "e_rem_pio2.c"
  21 #include "k_sincos.h"
  22
  23 void
  24 sincos(double x, double *sn, double *cs)
  25 {
  26         double y[2];
  27         int32_t n, ix;
  28
  29         /* High word of x. */
  30         GET_HIGH_WORD(ix, x);
  31
  32         /* |x| ~< pi/4 */
  33         ix &= 0x7fffffff;
  34         if (ix <= 0x3fe921fb) {
  35                 if (ix < 0x3e400000) {          /* |x| < 2**-27 */
  36                         if ((int)x == 0) {      /* Generate inexact. */
  37                                 *sn = x;
  38                                 *cs = 1;
  39                                 return;
  40                         }
  41                 }
  42                 __kernel_sincos(x, 0, 0, sn, cs);
  43                 return;
  44         }
  45
  46         /* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
  47         if (ix >= 0x7ff00000) {
  48                 *sn = x - x;
  49                 *cs = x - x;
  50                 return;
  51         }
  52
  53         /* Argument reduction. */
  54         n = __ieee754_rem_pio2(x, y);
  55
  56         switch(n & 3) {
  57         case 0:
  58                 __kernel_sincos(y[0], y[1], 1, sn, cs);
  59                 break;
  60         case 1:
  61                 __kernel_sincos(y[0], y[1], 1, cs, sn);
  62                 *cs = -*cs;
  63                 break;
  64         case 2:
  65                 __kernel_sincos(y[0], y[1], 1, sn, cs);
  66                 *sn = -*sn;
  67                 *cs = -*cs;
  68                 break;
  69         default:
  70                 __kernel_sincos(y[0], y[1], 1, cs, sn);
  71                 *sn = -*sn;
  72         }
  73 }
  74
  75 #if (LDBL_MANT_DIG == 53)
  76 __weak_reference(sincos, sincosl);
  77 #endif