1 /* origin: FreeBSD /usr/src/lib/msun/src/s_sinf.c */ 2 /* 3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 4 * Optimized by Bruce D. Evans. 5 */ 6 /* 7 * ==================================================== 8 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 9 * 10 * Developed at SunPro, a Sun Microsystems, Inc. business. 11 * Permission to use, copy, modify, and distribute this 12 * software is freely granted, provided that this notice 13 * is preserved. 14 * ==================================================== 15 */ 16 17 #define _GNU_SOURCE 18 #include "libm.h" 19 20 /* Small multiples of pi/2 rounded to double precision. */ 21 static const double 22 s1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */ 23 s2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */ 24 s3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */ 25 s4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */ 26 27 void sincosf(float x, float *sin, float *cos) 28 { 29 double y; 30 float_t s, c; 31 uint32_t ix; 32 unsigned n, sign; 33 34 GET_FLOAT_WORD(ix, x); 35 sign = ix >> 31; 36 ix &= 0x7fffffff; 37 38 /* |x| ~<= pi/4 */ 39 if (ix <= 0x3f490fda) { 40 /* |x| < 2**-12 */ 41 if (ix < 0x39800000) { 42 /* raise inexact if x!=0 and underflow if subnormal */ 43 FORCE_EVAL(ix < 0x00100000 ? x/0x1p120f : x+0x1p120f); 44 *sin = x; 45 *cos = 1.0f; 46 return; 47 } 48 *sin = __sindf(x); 49 *cos = __cosdf(x); 50 return; 51 } 52 53 /* |x| ~<= 5*pi/4 */ 54 if (ix <= 0x407b53d1) { 55 if (ix <= 0x4016cbe3) { /* |x| ~<= 3pi/4 */ 56 if (sign) { 57 *sin = -__cosdf(x + s1pio2); 58 *cos = __sindf(x + s1pio2); 59 } else { 60 *sin = __cosdf(s1pio2 - x); 61 *cos = __sindf(s1pio2 - x); 62 } 63 return; 64 } 65 /* -sin(x+c) is not correct if x+c could be 0: -0 vs +0 */ 66 *sin = -__sindf(sign ? x + s2pio2 : x - s2pio2); 67 *cos = -__cosdf(sign ? x + s2pio2 : x - s2pio2); 68 return; 69 } 70 71 /* |x| ~<= 9*pi/4 */ 72 if (ix <= 0x40e231d5) { 73 if (ix <= 0x40afeddf) { /* |x| ~<= 7*pi/4 */ 74 if (sign) { 75 *sin = __cosdf(x + s3pio2); 76 *cos = -__sindf(x + s3pio2); 77 } else { 78 *sin = -__cosdf(x - s3pio2); 79 *cos = __sindf(x - s3pio2); 80 } 81 return; 82 } 83 *sin = __sindf(sign ? x + s4pio2 : x - s4pio2); 84 *cos = __cosdf(sign ? x + s4pio2 : x - s4pio2); 85 return; 86 } 87 88 /* sin(Inf or NaN) is NaN */ 89 if (ix >= 0x7f800000) { 90 *sin = *cos = x - x; 91 return; 92 } 93 94 /* general argument reduction needed */ 95 n = __rem_pio2f(x, &y); 96 s = __sindf(y); 97 c = __cosdf(y); 98 switch (n&3) { 99 case 0: 100 *sin = s; 101 *cos = c; 102 break; 103 case 1: 104 *sin = c; 105 *cos = -s; 106 break; 107 case 2: 108 *sin = -s; 109 *cos = -c; 110 break; 111 case 3: 112 default: 113 *sin = -c; 114 *cos = s; 115 break; 116 } 117 } 118