1 /* tan (tangent) function with 'long double' argument. 2 3 Copyright (C) 2003-2021 Free Software Foundation, Inc. 4 5 This file is free software: you can redistribute it and/or modify 6 it under the terms of the GNU Lesser General Public License as 7 published by the Free Software Foundation; either version 3 of the 8 License, or (at your option) any later version. 9 10 This file is distributed in the hope that it will be useful, 11 but WITHOUT ANY WARRANTY; without even the implied warranty of 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 GNU Lesser General Public License for more details. 14 15 You should have received a copy of the GNU Lesser General Public License 16 along with this program. If not, see <https://www.gnu.org/licenses/>. */ 17 18 /* s_tanl.c -- long double version of s_tan.c. 19 * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz. 20 */ 21 22 /* @(#)s_tan.c 5.1 93/09/24 */ 23 /* 24 * ==================================================== 25 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 26 * 27 * Developed at SunPro, a Sun Microsystems, Inc. business. 28 * Permission to use, copy, modify, and distribute this 29 * software is freely granted, provided that this notice 30 * is preserved. 31 * ==================================================== 32 */ 33 34 #include <config.h> 35 36 /* Specification. */ 37 #include <math.h> 38 39 #if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE 40 41 long double 42 tanl (long double x) /* */ 43 { 44 return tan (x); 45 } 46 47 #else 48 49 /* Code based on glibc/sysdeps/ieee754/ldbl-128/s_tanl.c 50 and glibc/sysdeps/ieee754/ldbl-128/k_tanl.c. */ 51 52 /* tanl(x) 53 * Return tangent function of x. 54 * 55 * kernel function: 56 * __kernel_tanl ... tangent function on [-pi/4,pi/4] 57 * __ieee754_rem_pio2l ... argument reduction routine 58 * 59 * Method. 60 * Let S,C and T denote the sin, cos and tan respectively on 61 * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 62 * in [-pi/4 , +pi/4], and let n = k mod 4. 63 * We have 64 * 65 * n sin(x) cos(x) tan(x) 66 * ---------------------------------------------------------- 67 * 0 S C T 68 * 1 C -S -1/T 69 * 2 -S -C T 70 * 3 -C S -1/T 71 * ---------------------------------------------------------- 72 * 73 * Special cases: 74 * Let trig be any of sin, cos, or tan. 75 * trig(+-INF) is NaN, with signals; 76 * trig(NaN) is that NaN; 77 * 78 * Accuracy: 79 * TRIG(x) returns trig(x) nearly rounded 80 */ 81 82 # include "trigl.h" 83 84 /* 85 * ==================================================== 86 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 87 * 88 * Developed at SunPro, a Sun Microsystems, Inc. business. 89 * Permission to use, copy, modify, and distribute this 90 * software is freely granted, provided that this notice 91 * is preserved. 92 * ==================================================== 93 */ 94 95 /* 96 Long double expansions contributed by 97 Stephen L. Moshier <moshier@na-net.ornl.gov> 98 */ 99 100 /* __kernel_tanl( x, y, k ) 101 * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854 102 * Input x is assumed to be bounded by ~pi/4 in magnitude. 103 * Input y is the tail of x. 104 * Input k indicates whether tan (if k=1) or 105 * -1/tan (if k= -1) is returned. 106 * 107 * Algorithm 108 * 1. Since tan(-x) = -tan(x), we need only to consider positive x. 109 * 2. if x < 2^-57, return x with inexact if x!=0. 110 * 3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2) 111 * on [0,0.67433]. 112 * 113 * Note: tan(x+y) = tan(x) + tan'(x)*y 114 * ~ tan(x) + (1+x*x)*y 115 * Therefore, for better accuracy in computing tan(x+y), let 116 * r = x^3 * R(x^2) 117 * then 118 * tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y)) 119 * 120 * 4. For x in [0.67433,pi/4], let y = pi/4 - x, then 121 * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y)) 122 * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) 123 */ 124 125 126 static const long double 127 pio4hi = 7.8539816339744830961566084581987569936977E-1L, 128 pio4lo = 2.1679525325309452561992610065108379921906E-35L, 129 130 /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2) 131 0 <= x <= 0.6743316650390625 132 Peak relative error 8.0e-36 */ 133 TH = 3.333333333333333333333333333333333333333E-1L, 134 T0 = -1.813014711743583437742363284336855889393E7L, 135 T1 = 1.320767960008972224312740075083259247618E6L, 136 T2 = -2.626775478255838182468651821863299023956E4L, 137 T3 = 1.764573356488504935415411383687150199315E2L, 138 T4 = -3.333267763822178690794678978979803526092E-1L, 139 140 U0 = -1.359761033807687578306772463253710042010E8L, 141 U1 = 6.494370630656893175666729313065113194784E7L, 142 U2 = -4.180787672237927475505536849168729386782E6L, 143 U3 = 8.031643765106170040139966622980914621521E4L, 144 U4 = -5.323131271912475695157127875560667378597E2L; 145 /* 1.000000000000000000000000000000000000000E0 */ 146 147 148 static long double 149 kernel_tanl (long double x, long double y, int iy) /* */ 150 { 151 long double z, r, v, w, s, u, u1; 152 int invert = 0, sign; 153 154 sign = 1; 155 if (x < 0) 156 { 157 x = -x; 158 y = -y; 159 sign = -1; 160 } 161 162 if (x < 0.000000000000000006938893903907228377647697925567626953125L) /* x < 2**-57 */ 163 { 164 if ((int) x == 0) 165 { /* generate inexact */ 166 if (iy == -1 && x == 0.0) 167 return 1.0L / fabs (x); 168 else 169 return (iy == 1) ? x : -1.0L / x; 170 } 171 } 172 if (x >= 0.6743316650390625) /* |x| >= 0.6743316650390625 */ 173 { 174 invert = 1; 175 176 z = pio4hi - x; 177 w = pio4lo - y; 178 x = z + w; 179 y = 0.0; 180 } 181 z = x * x; 182 r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4))); 183 v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z)))); 184 r = r / v; 185 186 s = z * x; 187 r = y + z * (s * r + y); 188 r += TH * s; 189 w = x + r; 190 if (invert) 191 { 192 v = (long double) iy; 193 w = (v - 2.0 * (x - (w * w / (w + v) - r))); 194 if (sign < 0) 195 w = -w; 196 return w; 197 } 198 if (iy == 1) 199 return w; 200 else 201 { /* if allow error up to 2 ulp, 202 simply return -1.0/(x+r) here */ 203 /* compute -1.0/(x+r) accurately */ 204 u1 = (double) w; 205 v = r - (u1 - x); 206 z = -1.0 / w; 207 u = (double) z; 208 s = 1.0 + u * u1; 209 return u + z * (s + u * v); 210 } 211 } 212 213 long double 214 tanl (long double x) /* */ 215 { 216 long double y[2], z = 0.0L; 217 int n; 218 219 /* tanl(NaN) is NaN */ 220 if (isnanl (x)) 221 return x; 222 223 /* |x| ~< pi/4 */ 224 if (x >= -0.7853981633974483096156608458198757210492 && 225 x <= 0.7853981633974483096156608458198757210492) 226 return kernel_tanl (x, z, 1); 227 228 /* tanl(Inf) is NaN, tanl(0) is 0 */ 229 else if (x + x == x) 230 return x - x; /* NaN */ 231 232 /* argument reduction needed */ 233 else 234 { 235 n = ieee754_rem_pio2l (x, y); 236 /* 1 -- n even, -1 -- n odd */ 237 return kernel_tanl (y[0], y[1], 1 - ((n & 1) << 1)); 238 } 239 } 240 241 #endif 242 243 #if 0 244 int 245 main (void) /* */ 246 { 247 printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492)); 248 printf ("%.16Lg\n", tanl (-0.7853981633974483096156608458198757210492)); 249 printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 *3)); 250 printf ("%.16Lg\n", tanl (-0.7853981633974483096156608458198757210492 *31)); 251 printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 / 2)); 252 printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 * 3/2)); 253 printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 * 5/2)); 254 } 255 #endif