1 /* s_sincosf.c -- float version of s_sincos.c 2 * 3 * Copyright (C) 2013 Elliot Saba 4 * Developed at the University of Washington 5 * 6 * Permission to use, copy, modify, and distribute this 7 * software is freely granted, provided that this notice 8 * is preserved. 9 * ==================================================== 10 */ 11 12 #include "cdefs-compat.h" 13 14 #include <float.h> 15 #include <openlibm_math.h> 16 17 //#define INLINE_KERNEL_COSDF 18 //#define INLINE_KERNEL_SINDF 19 //#define INLINE_REM_PIO2F 20 #include "math_private.h" 21 //#include "e_rem_pio2f.c" 22 //#include "k_cosf.c" 23 //#include "k_sinf.c" 24 25 26 /* Constants used in shortcircuits in sincosf */ 27 static const double 28 sc1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */ 29 sc2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */ 30 sc3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */ 31 sc4pio2 = 4*M_PI_2, /* 0x401921FB, 0x54442D18 */ 32 33 /* Constants used in polynomial approximation of sin/cos */ 34 one = 1.0, 35 S1 = -0x15555554cbac77.0p-55, /* -0.166666666416265235595 */ 36 S2 = 0x111110896efbb2.0p-59, /* 0.0083333293858894631756 */ 37 S3 = -0x1a00f9e2cae774.0p-65, /* -0.000198393348360966317347 */ 38 S4 = 0x16cd878c3b46a7.0p-71, /* 0.0000027183114939898219064 */ 39 C0 = -0x1ffffffd0c5e81.0p-54, /* -0.499999997251031003120 */ 40 C1 = 0x155553e1053a42.0p-57, /* 0.0416666233237390631894 */ 41 C2 = -0x16c087e80f1e27.0p-62, /* -0.00138867637746099294692 */ 42 C3 = 0x199342e0ee5069.0p-68; /* 0.0000243904487962774090654 */ 43 44 static void 45 __kernel_sincosdf( double x, float * s, float * c ) 46 { 47 double r, w, z, v; 48 z = x*x; 49 w = z*z; 50 51 /* cos-specific computation; equivalent to calling 52 __kernel_cos(x,y) and storing in k_c*/ 53 r = C2+z*C3; 54 double k_c = ((one+z*C0) + w*C1) + (w*z)*r; 55 56 /* sin-specific computation; equivalent to calling 57 __kernel_sin(x,y,1) and storing in k_s*/ 58 r = S3+z*S4; 59 v = z*x; 60 double k_s = (x + v*(S1+z*S2)) + v*w*r; 61 62 *c = k_c; 63 *s = k_s; 64 } 65 66 OLM_DLLEXPORT void 67 sincosf(float x, float * s, float * c) { 68 // Worst approximation of sin and cos NA 69 *s = x; 70 *c = x; 71 72 double y; 73 float k_c, k_s; 74 int32_t n, hx, ix; 75 76 GET_FLOAT_WORD(hx,x); 77 ix = hx & 0x7fffffff; 78 79 if(ix <= 0x3f490fda) { /* |x| ~<= pi/4 */ 80 if(ix<0x39800000) { /* |x| < 2**-12 */ 81 /* Check if x is exactly zero */ 82 if(((int)x)==0) { 83 *s = x; 84 *c = 1.0f; 85 return; 86 } 87 } 88 __kernel_sincosdf(x, s, c); 89 return; 90 } 91 /* |x| ~<= 5*pi/4 */ 92 if (ix<=0x407b53d1) { 93 /* |x| ~<= 3pi/4 */ 94 if(ix<=0x4016cbe3) { 95 if(hx>0) { 96 __kernel_sincosdf( sc1pio2 - x, c, s ); 97 } 98 else { 99 __kernel_sincosdf( sc1pio2 + x, c, &k_s ); 100 *s = -k_s; 101 } 102 } else { 103 104 if(hx>0) { 105 __kernel_sincosdf( sc2pio2 - x, s, &k_c ); 106 *c = -k_c; 107 } else { 108 __kernel_sincosdf( -sc2pio2 - x, s, &k_c ); 109 *c = -k_c; 110 } 111 } 112 return; 113 } 114 115 /* |x| ~<= 9*pi/4 */ 116 if(ix<=0x40e231d5) { 117 /* |x| ~> 7*pi/4 */ 118 if(ix<=0x40afeddf) { 119 if(hx>0) { 120 __kernel_sincosdf( x - sc3pio2, c, &k_s ); 121 *s = -k_s; 122 } else { 123 __kernel_sincosdf( x + sc3pio2, &k_c, s ); 124 *c = -k_c; 125 } 126 } 127 else { 128 if( hx > 0 ) { 129 __kernel_sincosdf( x - sc4pio2, s, c ); 130 } else { 131 __kernel_sincosdf( x + sc4pio2, s, c ); 132 } 133 } 134 return; 135 } 136 137 /* cos(Inf or NaN) is NaN */ 138 else if(ix>=0x7f800000) { 139 *c = *s = x-x; 140 } else { 141 /* general argument reduction needed */ 142 n = __ieee754_rem_pio2f(x,&y); 143 144 switch(n&3) { 145 case 0: 146 __kernel_sincosdf( y, s, c ); 147 break; 148 case 1: 149 __kernel_sincosdf( -y, c, s ); 150 break; 151 case 2: 152 __kernel_sincosdf( -y, s, &k_c); 153 *c = -k_c; 154 break; 155 default: 156 __kernel_sincosdf( -y, &k_c, &k_s ); 157 *c = -k_c; 158 *s = -k_s; 159 break; 160 } 161 } 162 163 } 164