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_tanpi.c for implementation details. 29 */ 30 31 #ifdef __i386__ 32 #include <ieeefp.h> 33 #endif 34 #include <stdint.h> 35 36 #include "fpmath.h" 37 #include "math.h" 38 #include "math_private.h" 39 40 static const double 41 pi_hi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */ 42 pi_lo = -2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */ 43 44 static inline long double 45 __kernel_tanpil(long double x) 46 { 47 long double hi, lo, t; 48 49 if (x < 0.25) { 50 hi = (float)x; 51 lo = x - hi; 52 lo = lo * (pi_lo + pi_hi) + hi * pi_lo; 53 hi *= pi_hi; 54 _2sumF(hi, lo); 55 t = __kernel_tanl(hi, lo, -1); 56 } else if (x > 0.25) { 57 x = 0.5 - x; 58 hi = (float)x; 59 lo = x - hi; 60 lo = lo * (pi_lo + pi_hi) + hi * pi_lo; 61 hi *= pi_hi; 62 _2sumF(hi, lo); 63 t = - __kernel_tanl(hi, lo, 1); 64 } else 65 t = 1; 66 67 return (t); 68 } 69 70 volatile static const double vzero = 0; 71 72 long double 73 tanpil(long double x) 74 { 75 long double ax, hi, lo, t; 76 uint64_t lx, m; 77 uint32_t j0; 78 uint16_t hx, ix; 79 80 EXTRACT_LDBL80_WORDS(hx, lx, x); 81 ix = hx & 0x7fff; 82 INSERT_LDBL80_WORDS(ax, ix, lx); 83 84 ENTERI(); 85 86 if (ix < 0x3fff) { /* |x| < 1 */ 87 if (ix < 0x3ffe) { /* |x| < 0.5 */ 88 if (ix < 0x3fdd) { /* |x| < 0x1p-34 */ 89 if (x == 0) 90 RETURNI(x); 91 INSERT_LDBL80_WORDS(hi, hx, 92 lx & 0xffffffff00000000ull); 93 hi *= 0x1p63L; 94 lo = x * 0x1p63L - hi; 95 t = (pi_lo + pi_hi) * lo + pi_lo * hi + 96 pi_hi * hi; 97 RETURNI(t * 0x1p-63L); 98 } 99 t = __kernel_tanpil(ax); 100 } else if (ax == 0.5) 101 RETURNI((ax - ax) / (ax - ax)); 102 else 103 t = -__kernel_tanpil(1 - ax); 104 RETURNI((hx & 0x8000) ? -t : t); 105 } 106 107 if (ix < 0x403e) { /* 1 <= |x| < 0x1p63 */ 108 /* Determine integer part of ax. */ 109 j0 = ix - 0x3fff + 1; 110 if (j0 < 32) { 111 lx = (lx >> 32) << 32; 112 lx &= ~(((lx << 32)-1) >> j0); 113 } else { 114 m = (uint64_t)-1 >> (j0 + 1); 115 if (lx & m) lx &= ~m; 116 } 117 INSERT_LDBL80_WORDS(x, ix, lx); 118 119 ax -= x; 120 EXTRACT_LDBL80_WORDS(ix, lx, ax); 121 122 if (ix < 0x3ffe) /* |x| < 0.5 */ 123 t = ax == 0 ? 0 : __kernel_tanpil(ax); 124 else if (ax == 0.5) 125 RETURNI((ax - ax) / (ax - ax)); 126 else 127 t = -__kernel_tanpil(1 - ax); 128 RETURNI((hx & 0x8000) ? -t : t); 129 } 130 131 /* x = +-inf or nan. */ 132 if (ix >= 0x7fff) 133 RETURNI(vzero / vzero); 134 135 /* 136 * |x| >= 0x1p63 is always an integer, so return +-0. 137 */ 138 RETURNI(copysignl(0, x)); 139 } 140