1 /* @(#)s_sincos.c 5.1 13/07/15 */ 2 /* 3 * ==================================================== 4 * Copyright (C) 2013 Elliot Saba. All rights reserved. 5 * 6 * Developed at the University of Washington. 7 * Permission to use, copy, modify, and distribute this 8 * software is freely granted, provided that this notice 9 * is preserved. 10 * ==================================================== 11 */ 12 #include "cdefs-compat.h" 13 14 /* sincos(x, s, c) 15 * Several applications need sine and cosine of the same 16 * angle x. This function computes both at the same time, 17 * and stores the results in *sin and *cos. 18 * 19 * kernel function: 20 * __kernel_sin ... sine function on [-pi/4,pi/4] 21 * __kernel_cos ... cose function on [-pi/4,pi/4] 22 * __ieee754_rem_pio2 ... argument reduction routine 23 * 24 * Method. 25 * Borrow liberally from s_sin.c and s_cos.c, merging 26 * efforts where applicable and returning their values in 27 * appropriate variables, thereby slightly reducing the 28 * amount of work relative to just calling sin/cos(x) 29 * separately 30 * 31 * Special cases: 32 * Let trig be any of sin, cos, or tan. 33 * sincos(+-INF, s, c) is NaN, with signals; 34 * sincos(NaN, s, c) is that NaN; 35 */ 36 37 #include <float.h> 38 #include <openlibm.h> 39 40 //#define INLINE_REM_PIO2 41 #include "math_private.h" 42 //#include "e_rem_pio2.c" 43 44 /* Constants used in polynomial approximation of sin/cos */ 45 static const double 46 one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ 47 half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ 48 S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */ 49 S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */ 50 S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */ 51 S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */ 52 S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */ 53 S6 = 1.58969099521155010221e-10, /* 0x3DE5D93A, 0x5ACFD57C */ 54 C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */ 55 C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */ 56 C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */ 57 C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */ 58 C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */ 59 C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */ 60 61 static void 62 __kernel_sincos( double x, double y, int iy, double * k_s, double * k_c ) 63 { 64 /* Inline calculation of sin/cos, as we can save 65 some work, and we will always need to calculate 66 both values, no matter the result of switch */ 67 double z, w, r, v, hz; 68 z = x*x; 69 w = z*z; 70 71 /* cos-specific computation; equivalent to calling 72 __kernel_cos(x,y) and storing in k_c*/ 73 r = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6)); 74 hz = 0.5*z; 75 v = one-hz; 76 77 *k_c = v + (((one-v)-hz) + (z*r-x*y)); 78 79 /* sin-specific computation; equivalent to calling 80 __kernel_sin(x,y,1) and storing in k_s*/ 81 r = S2+z*(S3+z*S4) + z*w*(S5+z*S6); 82 v = z*x; 83 if(iy == 0) 84 *k_s = x+v*(S1+z*r); 85 else 86 *k_s = x-((z*(half*y-v*r)-y)-v*S1); 87 } 88 89 DLLEXPORT void 90 sincos(double x, double * s, double * c) 91 { 92 double y[2]; 93 int32_t ix; 94 95 /* Store high word of x in ix */ 96 GET_HIGH_WORD(ix,x); 97 98 /* |x| ~< pi/4 */ 99 ix &= 0x7fffffff; 100 if(ix <= 0x3fe921fb) { 101 /* Check for small x for sin and cos */ 102 if(ix<0x3e46a09e) { 103 /* Check for exact zero */ 104 if( (int)x==0 ) { 105 *s = x; 106 *c = 1.0; 107 return; 108 } 109 } 110 /* Call kernel function with 0 extra */ 111 __kernel_sincos(x,0.0,0, s, c); 112 } else if( ix >= 0x7ff00000 ) { 113 /* sincos(Inf or NaN) is NaN */ 114 *s = x-x; 115 *c = x-x; 116 } 117 118 /*argument reduction needed*/ 119 else { 120 double k_c, k_s; 121 122 /* Calculate remainer, then sub out to kernel */ 123 int32_t n = __ieee754_rem_pio2(x,y); 124 __kernel_sincos( y[0], y[1], 1, &k_s, &k_c ); 125 126 /* Figure out permutation of sin/cos outputs to true outputs */ 127 switch(n&3) { 128 case 0: 129 *c = k_c; 130 *s = k_s; 131 break; 132 case 1: 133 *c = -k_s; 134 *s = k_c; 135 break; 136 case 2: 137 *c = -k_c; 138 *s = -k_s; 139 break; 140 default: 141 *c = k_s; 142 *s = -k_c; 143 break; 144 } 145 } 146 } 147 148 #if (LDBL_MANT_DIG == 53) 149 __strong_reference(sincos, sincosl); 150 #endif 151