1 /*- 2 * Copyright (c) 2017 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 9 * notice unmodified, this list of conditions, and the following 10 * disclaimer. 11 * 2. Redistributions in binary form must reproduce the above copyright 12 * notice, this list of conditions and the following disclaimer in the 13 * documentation and/or other materials provided with the distribution. 14 * 15 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 16 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 17 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 18 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, 19 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 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 /* 28 * See ../src/s_sinpi.c for implementation details. 29 */ 30 31 #define INLINE_KERNEL_SINDF 32 #define INLINE_KERNEL_COSDF 33 34 #include "math.h" 35 #include "math_private.h" 36 #include "k_cosf.c" 37 #include "k_sinf.c" 38 39 #define __kernel_cospif(x) (__kernel_cosdf(M_PI * (x))) 40 #define __kernel_sinpif(x) (__kernel_sindf(M_PI * (x))) 41 42 static const float 43 pi_hi = 3.14160156e+00F, /* 0x40491000 */ 44 pi_lo = -8.90890988e-06F; /* 0xb715777a */ 45 46 volatile static const float vzero = 0; 47 48 float 49 sinpif(float x) 50 { 51 float ax, hi, lo, s; 52 uint32_t hx, ix, j0; 53 54 GET_FLOAT_WORD(hx, x); 55 ix = hx & 0x7fffffff; 56 SET_FLOAT_WORD(ax, ix); 57 58 if (ix < 0x3f800000) { /* |x| < 1 */ 59 if (ix < 0x3e800000) { /* |x| < 0.25 */ 60 if (ix < 0x38800000) { /* |x| < 0x1p-14 */ 61 if (x == 0) 62 return (x); 63 SET_FLOAT_WORD(hi, hx & 0xffff0000); 64 hi *= 0x1p23F; 65 lo = x * 0x1p23F - hi; 66 s = (pi_lo + pi_hi) * lo + pi_lo * hi + 67 pi_hi * hi; 68 return (s * 0x1p-23F); 69 } 70 71 s = __kernel_sinpif(ax); 72 return ((hx & 0x80000000) ? -s : s); 73 } 74 75 if (ix < 0x3f000000) /* |x| < 0.5 */ 76 s = __kernel_cospif(0.5F - ax); 77 else if (ix < 0x3f400000) /* |x| < 0.75 */ 78 s = __kernel_cospif(ax - 0.5F); 79 else 80 s = __kernel_sinpif(1 - ax); 81 return ((hx & 0x80000000) ? -s : s); 82 } 83 84 if (ix < 0x4b000000) { /* 1 <= |x| < 0x1p23 */ 85 /* Determine integer part of ax. */ 86 j0 = ((ix >> 23) & 0xff) - 0x7f; 87 ix &= ~(0x007fffff >> j0); 88 SET_FLOAT_WORD(x, ix); 89 90 ax -= x; 91 GET_FLOAT_WORD(ix, ax); 92 93 if (ix == 0) 94 s = 0; 95 else { 96 if (ix < 0x3f000000) { /* |x| < 0.5 */ 97 if (ix < 0x3e800000) /* |x| < 0.25 */ 98 s = __kernel_sinpif(ax); 99 else 100 s = __kernel_cospif(0.5F - ax); 101 } else { 102 if (ix < 0x3f400000) /* |x| < 0.75 */ 103 s = __kernel_cospif(ax - 0.5F); 104 else 105 s = __kernel_sinpif(1 - ax); 106 } 107 108 j0 = (uint32_t)x; 109 s = (j0 & 1) ? -s : s; 110 } 111 return ((hx & 0x80000000) ? -s : s); 112 } 113 114 /* x = +-inf or nan. */ 115 if (ix >= 0x7f800000) 116 return (vzero / vzero); 117 118 /* 119 * |x| >= 0x1p23 is always an integer, so return +-0. 120 */ 121 return (copysignf(0, x)); 122 } 123