1 /*- 2 * Copyright (c) 2017, 2023 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, odd, 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 t = 1 / vzero; 102 else 103 t = -__kernel_tanpil(1 - ax); 104 RETURNI((hx & 0x8000) ? -t : t); 105 } 106 107 if (ix < 0x403e) { /* 1 <= |x| < 0x1p63 */ 108 FFLOORL80(x, j0, ix, lx); /* Integer part of ax. */ 109 odd = (uint64_t)x & 1 ? -1 : 1; 110 ax -= x; 111 EXTRACT_LDBL80_WORDS(ix, lx, ax); 112 113 if (ix < 0x3ffe) /* |x| < 0.5 */ 114 t = ix == 0 ? copysignl(0, odd) : __kernel_tanpil(ax); 115 else if (ax == 0.5L) 116 t = odd / vzero; 117 else 118 t = -__kernel_tanpil(1 - ax); 119 RETURNI((hx & 0x8000) ? -t : t); 120 } 121 122 /* x = +-inf or nan. */ 123 if (ix >= 0x7fff) 124 RETURNI(vzero / vzero); 125 126 /* 127 * For 0x1p63 <= |x| < 0x1p64 need to determine if x is an even 128 * or odd integer to set t = +0 or -0. 129 * For |x| >= 0x1p64, it is always an even integer, so t = 0. 130 */ 131 t = ix >= 0x403f ? 0 : (copysignl(0, (lx & 1) ? -1 : 1)); 132 RETURNI((hx & 0x8000) ? -t : t); 133 } 134