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/* Modifications for 128-bit long double are 13 Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov> 14 and are incorporated herein by permission of the author. The author 15 reserves the right to distribute this material elsewhere under different 16 copying permissions. These modifications are distributed here under 17 the following terms: 18 19 This library is free software; you can redistribute it and/or 20 modify it under the terms of the GNU Lesser General Public 21 License as published by the Free Software Foundation; either 22 version 2.1 of the License, or (at your option) any later version. 23 24 This library is distributed in the hope that it will be useful, 25 but WITHOUT ANY WARRANTY; without even the implied warranty of 26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 27 Lesser General Public License for more details. 28 29 You should have received a copy of the GNU Lesser General Public 30 License along with this library; if not, see 31 <http://www.gnu.org/licenses/>. */ 32 33/* 34 * __ieee754_jn(n, x), __ieee754_yn(n, x) 35 * floating point Bessel's function of the 1st and 2nd kind 36 * of order n 37 * 38 * Special cases: 39 * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal; 40 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal. 41 * Note 2. About jn(n,x), yn(n,x) 42 * For n=0, j0(x) is called, 43 * for n=1, j1(x) is called, 44 * for n<x, forward recursion us used starting 45 * from values of j0(x) and j1(x). 46 * for n>x, a continued fraction approximation to 47 * j(n,x)/j(n-1,x) is evaluated and then backward 48 * recursion is used starting from a supposed value 49 * for j(n,x). The resulting value of j(0,x) is 50 * compared with the actual value to correct the 51 * supposed value of j(n,x). 52 * 53 * yn(n,x) is similar in all respects, except 54 * that forward recursion is used for all 55 * values of n>1. 56 * 57 */ 58 59#include "quadmath-imp.h" 60 61static const __float128 62 invsqrtpi = 5.6418958354775628694807945156077258584405E-1Q, 63 two = 2, 64 one = 1, 65 zero = 0; 66 67 68__float128 69jnq (int n, __float128 x) 70{ 71 uint32_t se; 72 int32_t i, ix, sgn; 73 __float128 a, b, temp, di, ret; 74 __float128 z, w; 75 ieee854_float128 u; 76 77 78 /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x) 79 * Thus, J(-n,x) = J(n,-x) 80 */ 81 82 u.value = x; 83 se = u.words32.w0; 84 ix = se & 0x7fffffff; 85 86 /* if J(n,NaN) is NaN */ 87 if (ix >= 0x7fff0000) 88 { 89 if ((u.words32.w0 & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3) 90 return x + x; 91 } 92 93 if (n < 0) 94 { 95 n = -n; 96 x = -x; 97 se ^= 0x80000000; 98 } 99 if (n == 0) 100 return (j0q (x)); 101 if (n == 1) 102 return (j1q (x)); 103 sgn = (n & 1) & (se >> 31); /* even n -- 0, odd n -- sign(x) */ 104 x = fabsq (x); 105 106 { 107 SET_RESTORE_ROUNDF128 (FE_TONEAREST); 108 if (x == 0 || ix >= 0x7fff0000) /* if x is 0 or inf */ 109 return sgn == 1 ? -zero : zero; 110 else if ((__float128) n <= x) 111 { 112 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ 113 if (ix >= 0x412D0000) 114 { /* x > 2**302 */ 115 116 /* ??? Could use an expansion for large x here. */ 117 118 /* (x >> n**2) 119 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) 120 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) 121 * Let s=sin(x), c=cos(x), 122 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then 123 * 124 * n sin(xn)*sqt2 cos(xn)*sqt2 125 * ---------------------------------- 126 * 0 s-c c+s 127 * 1 -s-c -c+s 128 * 2 -s+c -c-s 129 * 3 s+c c-s 130 */ 131 __float128 s; 132 __float128 c; 133 sincosq (x, &s, &c); 134 switch (n & 3) 135 { 136 case 0: 137 temp = c + s; 138 break; 139 case 1: 140 temp = -c + s; 141 break; 142 case 2: 143 temp = -c - s; 144 break; 145 case 3: 146 temp = c - s; 147 break; 148 } 149 b = invsqrtpi * temp / sqrtq (x); 150 } 151 else 152 { 153 a = j0q (x); 154 b = j1q (x); 155 for (i = 1; i < n; i++) 156 { 157 temp = b; 158 b = b * ((__float128) (i + i) / x) - a; /* avoid underflow */ 159 a = temp; 160 } 161 } 162 } 163 else 164 { 165 if (ix < 0x3fc60000) 166 { /* x < 2**-57 */ 167 /* x is tiny, return the first Taylor expansion of J(n,x) 168 * J(n,x) = 1/n!*(x/2)^n - ... 169 */ 170 if (n >= 400) /* underflow, result < 10^-4952 */ 171 b = zero; 172 else 173 { 174 temp = x * 0.5; 175 b = temp; 176 for (a = one, i = 2; i <= n; i++) 177 { 178 a *= (__float128) i; /* a = n! */ 179 b *= temp; /* b = (x/2)^n */ 180 } 181 b = b / a; 182 } 183 } 184 else 185 { 186 /* use backward recurrence */ 187 /* x x^2 x^2 188 * J(n,x)/J(n-1,x) = ---- ------ ------ ..... 189 * 2n - 2(n+1) - 2(n+2) 190 * 191 * 1 1 1 192 * (for large x) = ---- ------ ------ ..... 193 * 2n 2(n+1) 2(n+2) 194 * -- - ------ - ------ - 195 * x x x 196 * 197 * Let w = 2n/x and h=2/x, then the above quotient 198 * is equal to the continued fraction: 199 * 1 200 * = ----------------------- 201 * 1 202 * w - ----------------- 203 * 1 204 * w+h - --------- 205 * w+2h - ... 206 * 207 * To determine how many terms needed, let 208 * Q(0) = w, Q(1) = w(w+h) - 1, 209 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2), 210 * When Q(k) > 1e4 good for single 211 * When Q(k) > 1e9 good for double 212 * When Q(k) > 1e17 good for quadruple 213 */ 214 /* determine k */ 215 __float128 t, v; 216 __float128 q0, q1, h, tmp; 217 int32_t k, m; 218 w = (n + n) / (__float128) x; 219 h = 2 / (__float128) x; 220 q0 = w; 221 z = w + h; 222 q1 = w * z - 1; 223 k = 1; 224 while (q1 < 1.0e17Q) 225 { 226 k += 1; 227 z += h; 228 tmp = z * q1 - q0; 229 q0 = q1; 230 q1 = tmp; 231 } 232 m = n + n; 233 for (t = zero, i = 2 * (n + k); i >= m; i -= 2) 234 t = one / (i / x - t); 235 a = t; 236 b = one; 237 /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) 238 * Hence, if n*(log(2n/x)) > ... 239 * single 8.8722839355e+01 240 * double 7.09782712893383973096e+02 241 * long double 1.1356523406294143949491931077970765006170e+04 242 * then recurrent value may overflow and the result is 243 * likely underflow to zero 244 */ 245 tmp = n; 246 v = two / x; 247 tmp = tmp * logq (fabsq (v * tmp)); 248 249 if (tmp < 1.1356523406294143949491931077970765006170e+04Q) 250 { 251 for (i = n - 1, di = (__float128) (i + i); i > 0; i--) 252 { 253 temp = b; 254 b *= di; 255 b = b / x - a; 256 a = temp; 257 di -= two; 258 } 259 } 260 else 261 { 262 for (i = n - 1, di = (__float128) (i + i); i > 0; i--) 263 { 264 temp = b; 265 b *= di; 266 b = b / x - a; 267 a = temp; 268 di -= two; 269 /* scale b to avoid spurious overflow */ 270 if (b > 1e100Q) 271 { 272 a /= b; 273 t /= b; 274 b = one; 275 } 276 } 277 } 278 /* j0() and j1() suffer enormous loss of precision at and 279 * near zero; however, we know that their zero points never 280 * coincide, so just choose the one further away from zero. 281 */ 282 z = j0q (x); 283 w = j1q (x); 284 if (fabsq (z) >= fabsq (w)) 285 b = (t * z / b); 286 else 287 b = (t * w / a); 288 } 289 } 290 if (sgn == 1) 291 ret = -b; 292 else 293 ret = b; 294 } 295 if (ret == 0) 296 { 297 ret = copysignq (FLT128_MIN, ret) * FLT128_MIN; 298 errno = ERANGE; 299 } 300 else 301 math_check_force_underflow (ret); 302 return ret; 303} 304 305 306__float128 307ynq (int n, __float128 x) 308{ 309 uint32_t se; 310 int32_t i, ix; 311 int32_t sign; 312 __float128 a, b, temp, ret; 313 ieee854_float128 u; 314 315 u.value = x; 316 se = u.words32.w0; 317 ix = se & 0x7fffffff; 318 319 /* if Y(n,NaN) is NaN */ 320 if (ix >= 0x7fff0000) 321 { 322 if ((u.words32.w0 & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3) 323 return x + x; 324 } 325 if (x <= 0) 326 { 327 if (x == 0) 328 return ((n < 0 && (n & 1) != 0) ? 1 : -1) / 0.0Q; 329 if (se & 0x80000000) 330 return zero / (zero * x); 331 } 332 sign = 1; 333 if (n < 0) 334 { 335 n = -n; 336 sign = 1 - ((n & 1) << 1); 337 } 338 if (n == 0) 339 return (y0q (x)); 340 { 341 SET_RESTORE_ROUNDF128 (FE_TONEAREST); 342 if (n == 1) 343 { 344 ret = sign * y1q (x); 345 goto out; 346 } 347 if (ix >= 0x7fff0000) 348 return zero; 349 if (ix >= 0x412D0000) 350 { /* x > 2**302 */ 351 352 /* ??? See comment above on the possible futility of this. */ 353 354 /* (x >> n**2) 355 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) 356 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) 357 * Let s=sin(x), c=cos(x), 358 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then 359 * 360 * n sin(xn)*sqt2 cos(xn)*sqt2 361 * ---------------------------------- 362 * 0 s-c c+s 363 * 1 -s-c -c+s 364 * 2 -s+c -c-s 365 * 3 s+c c-s 366 */ 367 __float128 s; 368 __float128 c; 369 sincosq (x, &s, &c); 370 switch (n & 3) 371 { 372 case 0: 373 temp = s - c; 374 break; 375 case 1: 376 temp = -s - c; 377 break; 378 case 2: 379 temp = -s + c; 380 break; 381 case 3: 382 temp = s + c; 383 break; 384 } 385 b = invsqrtpi * temp / sqrtq (x); 386 } 387 else 388 { 389 a = y0q (x); 390 b = y1q (x); 391 /* quit if b is -inf */ 392 u.value = b; 393 se = u.words32.w0 & 0xffff0000; 394 for (i = 1; i < n && se != 0xffff0000; i++) 395 { 396 temp = b; 397 b = ((__float128) (i + i) / x) * b - a; 398 u.value = b; 399 se = u.words32.w0 & 0xffff0000; 400 a = temp; 401 } 402 } 403 /* If B is +-Inf, set up errno accordingly. */ 404 if (! finiteq (b)) 405 errno = ERANGE; 406 if (sign > 0) 407 ret = b; 408 else 409 ret = -b; 410 } 411 out: 412 if (isinfq (ret)) 413 ret = copysignq (FLT128_MAX, ret) * FLT128_MAX; 414 return ret; 415} 416