s_sinl.c (176360) | s_sinl.c (221234) |
---|---|
1/*- 2 * Copyright (c) 2007 Steven G. Kargl 3 * All rights reserved. 4 * 5 * Redistribution and use in source and binary forms, with or without 6 * modification, are permitted provided that the following conditions 7 * are met: 8 * 1. Redistributions of source code must retain the above copyright --- 11 unchanged lines hidden (view full) --- 20 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 21 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 22 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 23 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 24 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 25 */ 26 27#include <sys/cdefs.h> | 1/*- 2 * Copyright (c) 2007 Steven G. Kargl 3 * All rights reserved. 4 * 5 * Redistribution and use in source and binary forms, with or without 6 * modification, are permitted provided that the following conditions 7 * are met: 8 * 1. Redistributions of source code must retain the above copyright --- 11 unchanged lines hidden (view full) --- 20 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 21 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 22 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 23 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 24 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 25 */ 26 27#include <sys/cdefs.h> |
28__FBSDID("$FreeBSD: head/lib/msun/src/s_sinl.c 176360 2008-02-17 07:33:12Z das $"); | 28__FBSDID("$FreeBSD: head/lib/msun/src/s_sinl.c 221234 2011-04-29 23:13:43Z kargl $"); |
29 | 29 |
30/* 31 * Compute sin(x) for x where x is reduced to y = x - k * pi / 2. 32 */ 33 | |
34#include <float.h> 35 36#include "math.h" | 30#include <float.h> 31 32#include "math.h" |
33#define INLINE_REM_PIO2L |
|
37#include "math_private.h" | 34#include "math_private.h" |
38#include "fpmath.h" 39 | |
40#if LDBL_MANT_DIG == 64 | 35#if LDBL_MANT_DIG == 64 |
41#define NX 3 42#define PREC 2 | 36#include "../ld80/e_rem_pio2l.h" |
43#elif LDBL_MANT_DIG == 113 | 37#elif LDBL_MANT_DIG == 113 |
44#define NX 5 45#define PREC 3 | 38#include "../ld128/e_rem_pio2l.h" |
46#else 47#error "Unsupported long double format" 48#endif 49 | 39#else 40#error "Unsupported long double format" 41#endif 42 |
50static const long double two24 = 1.67772160000000000000e+07L; 51 | |
52long double 53sinl(long double x) 54{ 55 union IEEEl2bits z; | 43long double 44sinl(long double x) 45{ 46 union IEEEl2bits z; |
56 int i, e0, s; 57 double xd[NX], yd[PREC]; | 47 int e0, s; 48 long double y[2]; |
58 long double hi, lo; 59 60 z.e = x; 61 s = z.bits.sign; 62 z.bits.sign = 0; 63 64 /* If x = +-0 or x is a subnormal number, then sin(x) = x */ 65 if (z.bits.exp == 0) --- 4 unchanged lines hidden (view full) --- 70 return ((x - x) / (x - x)); 71 72 /* Optimize the case where x is already within range. */ 73 if (z.e < M_PI_4) { 74 hi = __kernel_sinl(z.e, 0, 0); 75 return (s ? -hi : hi); 76 } 77 | 49 long double hi, lo; 50 51 z.e = x; 52 s = z.bits.sign; 53 z.bits.sign = 0; 54 55 /* If x = +-0 or x is a subnormal number, then sin(x) = x */ 56 if (z.bits.exp == 0) --- 4 unchanged lines hidden (view full) --- 61 return ((x - x) / (x - x)); 62 63 /* Optimize the case where x is already within range. */ 64 if (z.e < M_PI_4) { 65 hi = __kernel_sinl(z.e, 0, 0); 66 return (s ? -hi : hi); 67 } 68 |
78 /* Split z.e into a 24-bit representation. */ 79 e0 = ilogbl(z.e) - 23; 80 z.e = scalbnl(z.e, -e0); 81 for (i = 0; i < NX; i++) { 82 xd[i] = (double)((int32_t)z.e); 83 z.e = (z.e - xd[i]) * two24; 84 } | 69 e0 = __ieee754_rem_pio2l(x, y); 70 hi = y[0]; 71 lo = y[1]; |
85 | 72 |
86 /* yd contains the pieces of xd rem pi/2 such that |yd| < pi/4. */ 87 e0 = __kernel_rem_pio2(xd, yd, e0, NX, PREC); 88 89#if PREC == 2 90 hi = (long double)yd[0] + yd[1]; 91 lo = yd[1] - (hi - yd[0]); 92#else /* PREC == 3 */ 93 long double t; 94 t = (long double)yd[2] + yd[1]; 95 hi = t + yd[0]; 96 lo = yd[0] - (hi - t); 97#endif 98 | |
99 switch (e0 & 3) { 100 case 0: 101 hi = __kernel_sinl(hi, lo, 1); 102 break; 103 case 1: 104 hi = __kernel_cosl(hi, lo); 105 break; 106 case 2: 107 hi = - __kernel_sinl(hi, lo, 1); 108 break; 109 case 3: 110 hi = - __kernel_cosl(hi, lo); 111 break; 112 } 113 | 73 switch (e0 & 3) { 74 case 0: 75 hi = __kernel_sinl(hi, lo, 1); 76 break; 77 case 1: 78 hi = __kernel_cosl(hi, lo); 79 break; 80 case 2: 81 hi = - __kernel_sinl(hi, lo, 1); 82 break; 83 case 3: 84 hi = - __kernel_cosl(hi, lo); 85 break; 86 } 87 |
114 return (s ? -hi : hi); | 88 return (hi); |
115} | 89} |