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