1 | /* |
---|---|

2 | * ==================================================== |

3 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |

4 | * |

5 | * Developed at SunPro, a Sun Microsystems, Inc. business. |

6 | * Permission to use, copy, modify, and distribute this |

7 | * software is freely granted, provided that this notice |

8 | * is preserved. |

9 | * ==================================================== |

10 | */ |

11 | |

12 | /* |

13 | Long double expansions are |

14 | Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov> |

15 | and are incorporated herein by permission of the author. The author |

16 | reserves the right to distribute this material elsewhere under different |

17 | copying permissions. These modifications are distributed here under |

18 | the following terms: |

19 | |

20 | This library is free software; you can redistribute it and/or |

21 | modify it under the terms of the GNU Lesser General Public |

22 | License as published by the Free Software Foundation; either |

23 | version 2.1 of the License, or (at your option) any later version. |

24 | |

25 | This library is distributed in the hope that it will be useful, |

26 | but WITHOUT ANY WARRANTY; without even the implied warranty of |

27 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |

28 | Lesser General Public License for more details. |

29 | |

30 | You should have received a copy of the GNU Lesser General Public |

31 | License along with this library; if not, write to the Free Software |

32 | Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ |

33 | |

34 | /* __quadmath_kernel_tanq( x, y, k ) |

35 | * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854 |

36 | * Input x is assumed to be bounded by ~pi/4 in magnitude. |

37 | * Input y is the tail of x. |

38 | * Input k indicates whether tan (if k=1) or |

39 | * -1/tan (if k= -1) is returned. |

40 | * |

41 | * Algorithm |

42 | * 1. Since tan(-x) = -tan(x), we need only to consider positive x. |

43 | * 2. if x < 2^-57, return x with inexact if x!=0. |

44 | * 3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2) |

45 | * on [0,0.67433]. |

46 | * |

47 | * Note: tan(x+y) = tan(x) + tan'(x)*y |

48 | * ~ tan(x) + (1+x*x)*y |

49 | * Therefore, for better accuracy in computing tan(x+y), let |

50 | * r = x^3 * R(x^2) |

51 | * then |

52 | * tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y)) |

53 | * |

54 | * 4. For x in [0.67433,pi/4], let y = pi/4 - x, then |

55 | * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y)) |

56 | * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) |

57 | */ |

58 | |

59 | #include "quadmath-imp.h" |

60 | |

61 | |

62 | |

63 | static const __float128 |

64 | one = 1.0Q, |

65 | pio4hi = 7.8539816339744830961566084581987569936977E-1Q, |

66 | pio4lo = 2.1679525325309452561992610065108379921906E-35Q, |

67 | |

68 | /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2) |

69 | 0 <= x <= 0.6743316650390625 |

70 | Peak relative error 8.0e-36 */ |

71 | TH = 3.333333333333333333333333333333333333333E-1Q, |

72 | T0 = -1.813014711743583437742363284336855889393E7Q, |

73 | T1 = 1.320767960008972224312740075083259247618E6Q, |

74 | T2 = -2.626775478255838182468651821863299023956E4Q, |

75 | T3 = 1.764573356488504935415411383687150199315E2Q, |

76 | T4 = -3.333267763822178690794678978979803526092E-1Q, |

77 | |

78 | U0 = -1.359761033807687578306772463253710042010E8Q, |

79 | U1 = 6.494370630656893175666729313065113194784E7Q, |

80 | U2 = -4.180787672237927475505536849168729386782E6Q, |

81 | U3 = 8.031643765106170040139966622980914621521E4Q, |

82 | U4 = -5.323131271912475695157127875560667378597E2Q; |

83 | /* 1.000000000000000000000000000000000000000E0 */ |

84 | |

85 | |

86 | static __float128 |

87 | __quadmath_kernel_tanq (__float128 x, __float128 y, int iy) |

88 | { |

89 | __float128 z, r, v, w, s; |

90 | int32_t ix, sign = 1; |

91 | ieee854_float128 u, u1; |

92 | |

93 | u.value = x; |

94 | ix = u.words32.w0 & 0x7fffffff; |

95 | if (ix < 0x3fc60000) /* x < 2**-57 */ |

96 | { |

97 | if ((int) x == 0) |

98 | { /* generate inexact */ |

99 | if ((ix | u.words32.w1 | u.words32.w2 | u.words32.w3 |

100 | | (iy + 1)) == 0) |

101 | return one / fabsq (x); |

102 | else if (iy == 1) |

103 | { |

104 | math_check_force_underflow (x); |

105 | return x; |

106 | } |

107 | else |

108 | return -one / x; |

109 | } |

110 | } |

111 | if (ix >= 0x3ffe5942) /* |x| >= 0.6743316650390625 */ |

112 | { |

113 | if ((u.words32.w0 & 0x80000000) != 0) |

114 | { |

115 | x = -x; |

116 | y = -y; |

117 | sign = -1; |

118 | } |

119 | else |

120 | sign = 1; |

121 | z = pio4hi - x; |

122 | w = pio4lo - y; |

123 | x = z + w; |

124 | y = 0.0; |

125 | } |

126 | z = x * x; |

127 | r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4))); |

128 | v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z)))); |

129 | r = r / v; |

130 | |

131 | s = z * x; |

132 | r = y + z * (s * r + y); |

133 | r += TH * s; |

134 | w = x + r; |

135 | if (ix >= 0x3ffe5942) |

136 | { |

137 | v = (__float128) iy; |

138 | w = (v - 2.0Q * (x - (w * w / (w + v) - r))); |

139 | if (sign < 0) |

140 | w = -w; |

141 | return w; |

142 | } |

143 | if (iy == 1) |

144 | return w; |

145 | else |

146 | { /* if allow error up to 2 ulp, |

147 | simply return -1.0/(x+r) here */ |

148 | /* compute -1.0/(x+r) accurately */ |

149 | u1.value = w; |

150 | u1.words32.w2 = 0; |

151 | u1.words32.w3 = 0; |

152 | v = r - (u1.value - x); /* u1+v = r+x */ |

153 | z = -1.0 / w; |

154 | u.value = z; |

155 | u.words32.w2 = 0; |

156 | u.words32.w3 = 0; |

157 | s = 1.0 + u.value * u1.value; |

158 | return u.value + z * (s + u.value * v); |

159 | } |

160 | } |

161 | |

162 | |

163 | |

164 | |

165 | |

166 | |

167 | |

168 | /* tanq.c -- __float128 version of s_tan.c. |

169 | * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz. |

170 | */ |

171 | |

172 | /* @(#)s_tan.c 5.1 93/09/24 */ |

173 | /* |

174 | * ==================================================== |

175 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |

176 | * |

177 | * Developed at SunPro, a Sun Microsystems, Inc. business. |

178 | * Permission to use, copy, modify, and distribute this |

179 | * software is freely granted, provided that this notice |

180 | * is preserved. |

181 | * ==================================================== |

182 | */ |

183 | |

184 | /* tanl(x) |

185 | * Return tangent function of x. |

186 | * |

187 | * kernel function: |

188 | * __quadmath_kernel_tanq ... tangent function on [-pi/4,pi/4] |

189 | * __quadmath_rem_pio2q ... argument reduction routine |

190 | * |

191 | * Method. |

192 | * Let S,C and T denote the sin, cos and tan respectively on |

193 | * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 |

194 | * in [-pi/4 , +pi/4], and let n = k mod 4. |

195 | * We have |

196 | * |

197 | * n sin(x) cos(x) tan(x) |

198 | * ---------------------------------------------------------- |

199 | * 0 S C T |

200 | * 1 C -S -1/T |

201 | * 2 -S -C T |

202 | * 3 -C S -1/T |

203 | * ---------------------------------------------------------- |

204 | * |

205 | * Special cases: |

206 | * Let trig be any of sin, cos, or tan. |

207 | * trig(+-INF) is NaN, with signals; |

208 | * trig(NaN) is that NaN; |

209 | * |

210 | * Accuracy: |

211 | * TRIG(x) returns trig(x) nearly rounded |

212 | */ |

213 | |

214 | |

215 | __float128 |

216 | tanq (__float128 x) |

217 | { |

218 | __float128 y[2],z=0.0Q; |

219 | int64_t n, ix; |

220 | |

221 | /* High word of x. */ |

222 | GET_FLT128_MSW64(ix,x); |

223 | |

224 | /* |x| ~< pi/4 */ |

225 | ix &= 0x7fffffffffffffffLL; |

226 | if(ix <= 0x3ffe921fb54442d1LL) return __quadmath_kernel_tanq(x,z,1); |

227 | |

228 | /* tanl(Inf or NaN) is NaN */ |

229 | else if (ix>=0x7fff000000000000LL) { |

230 | if (ix == 0x7fff000000000000LL) { |

231 | GET_FLT128_LSW64(n,x); |

232 | } |

233 | return x-x; /* NaN */ |

234 | } |

235 | |

236 | /* argument reduction needed */ |

237 | else { |

238 | n = __quadmath_rem_pio2q(x,y); |

239 | /* 1 -- n even, -1 -- n odd */ |

240 | return __quadmath_kernel_tanq(y[0],y[1],1-((n&1)<<1)); |

241 | } |

242 | } |

243 |