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, write to the Free Software 31 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ 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 <errno.h> 60#include "quadmath-imp.h" 61 62static const __float128 63 invsqrtpi = 5.6418958354775628694807945156077258584405E-1Q, 64 two = 2.0e0Q, 65 one = 1.0e0Q, 66 zero = 0.0Q; 67 68 69__float128 70jnq (int n, __float128 x) 71{ 72 uint32_t se; 73 int32_t i, ix, sgn; 74 __float128 a, b, temp, di; 75 __float128 z, w; 76 ieee854_float128 u; 77 78 79 /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x) 80 * Thus, J(-n,x) = J(n,-x) 81 */ 82 83 u.value = x; 84 se = u.words32.w0; 85 ix = se & 0x7fffffff; 86 87 /* if J(n,NaN) is NaN */ 88 if (ix >= 0x7fff0000) 89 { 90 if ((u.words32.w0 & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3) 91 return x + x; 92 } 93 94 if (n < 0) 95 { 96 n = -n; 97 x = -x; 98 se ^= 0x80000000; 99 } 100 if (n == 0) 101 return (j0q (x)); 102 if (n == 1) 103 return (j1q (x)); 104 sgn = (n & 1) & (se >> 31); /* even n -- 0, odd n -- sign(x) */ 105 x = fabsq (x); 106 107 if (x == 0.0Q || ix >= 0x7fff0000) /* if x is 0 or inf */ 108 b = zero; 109 else if ((__float128) n <= x) 110 { 111 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ 112 if (ix >= 0x412D0000) 113 { /* x > 2**302 */ 114 115 /* ??? Could use an expansion for large x here. */ 116 117 /* (x >> n**2) 118 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) 119 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) 120 * Let s=sin(x), c=cos(x), 121 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then 122 * 123 * n sin(xn)*sqt2 cos(xn)*sqt2 124 * ---------------------------------- 125 * 0 s-c c+s 126 * 1 -s-c -c+s 127 * 2 -s+c -c-s 128 * 3 s+c c-s 129 */ 130 __float128 s; 131 __float128 c; 132 sincosq (x, &s, &c); 133 switch (n & 3) 134 { 135 case 0: 136 temp = c + s; 137 break; 138 case 1: 139 temp = -c + s; 140 break; 141 case 2: 142 temp = -c - s; 143 break; 144 case 3: 145 temp = c - s; 146 break; 147 } 148 b = invsqrtpi * temp / sqrtq (x); 149 } 150 else 151 { 152 a = j0q (x); 153 b = j1q (x); 154 for (i = 1; i < n; i++) 155 { 156 temp = b; 157 b = b * ((__float128) (i + i) / x) - a; /* avoid underflow */ 158 a = temp; 159 } 160 } 161 } 162 else 163 { 164 if (ix < 0x3fc60000) 165 { /* x < 2**-57 */ 166 /* x is tiny, return the first Taylor expansion of J(n,x) 167 * J(n,x) = 1/n!*(x/2)^n - ... 168 */ 169 if (n >= 400) /* underflow, result < 10^-4952 */ 170 b = zero; 171 else 172 { 173 temp = x * 0.5; 174 b = temp; 175 for (a = one, i = 2; i <= n; i++) 176 { 177 a *= (__float128) i; /* a = n! */ 178 b *= temp; /* b = (x/2)^n */ 179 } 180 b = b / a; 181 } 182 } 183 else 184 { 185 /* use backward recurrence */ 186 /* x x^2 x^2 187 * J(n,x)/J(n-1,x) = ---- ------ ------ ..... 188 * 2n - 2(n+1) - 2(n+2) 189 * 190 * 1 1 1 191 * (for large x) = ---- ------ ------ ..... 192 * 2n 2(n+1) 2(n+2) 193 * -- - ------ - ------ - 194 * x x x 195 * 196 * Let w = 2n/x and h=2/x, then the above quotient 197 * is equal to the continued fraction: 198 * 1 199 * = ----------------------- 200 * 1 201 * w - ----------------- 202 * 1 203 * w+h - --------- 204 * w+2h - ... 205 * 206 * To determine how many terms needed, let 207 * Q(0) = w, Q(1) = w(w+h) - 1, 208 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2), 209 * When Q(k) > 1e4 good for single 210 * When Q(k) > 1e9 good for double 211 * When Q(k) > 1e17 good for quadruple 212 */ 213 /* determine k */ 214 __float128 t, v; 215 __float128 q0, q1, h, tmp; 216 int32_t k, m; 217 w = (n + n) / (__float128) x; 218 h = 2.0Q / (__float128) x; 219 q0 = w; 220 z = w + h; 221 q1 = w * z - 1.0Q; 222 k = 1; 223 while (q1 < 1.0e17Q) 224 { 225 k += 1; 226 z += h; 227 tmp = z * q1 - q0; 228 q0 = q1; 229 q1 = tmp; 230 } 231 m = n + n; 232 for (t = zero, i = 2 * (n + k); i >= m; i -= 2) 233 t = one / (i / x - t); 234 a = t; 235 b = one; 236 /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) 237 * Hence, if n*(log(2n/x)) > ... 238 * single 8.8722839355e+01 239 * double 7.09782712893383973096e+02 240 * __float128 1.1356523406294143949491931077970765006170e+04 241 * then recurrent value may overflow and the result is 242 * likely underflow to zero 243 */ 244 tmp = n; 245 v = two / x; 246 tmp = tmp * logq (fabsq (v * tmp)); 247 248 if (tmp < 1.1356523406294143949491931077970765006170e+04Q) 249 { 250 for (i = n - 1, di = (__float128) (i + i); i > 0; i--) 251 { 252 temp = b; 253 b *= di; 254 b = b / x - a; 255 a = temp; 256 di -= two; 257 } 258 } 259 else 260 { 261 for (i = n - 1, di = (__float128) (i + i); i > 0; i--) 262 { 263 temp = b; 264 b *= di; 265 b = b / x - a; 266 a = temp; 267 di -= two; 268 /* scale b to avoid spurious overflow */ 269 if (b > 1e100Q) 270 { 271 a /= b; 272 t /= b; 273 b = one; 274 } 275 } 276 } 277 /* j0() and j1() suffer enormous loss of precision at and 278 * near zero; however, we know that their zero points never 279 * coincide, so just choose the one further away from zero. 280 */ 281 z = j0q (x); 282 w = j1q (x); 283 if (fabsq (z) >= fabsq (w)) 284 b = (t * z / b); 285 else 286 b = (t * w / a); 287 } 288 } 289 if (sgn == 1) 290 return -b; 291 else 292 return b; 293} 294 295__float128 296ynq (int n, __float128 x) 297{ 298 uint32_t se; 299 int32_t i, ix; 300 int32_t sign; 301 __float128 a, b, temp; 302 ieee854_float128 u; 303 304 u.value = x; 305 se = u.words32.w0; 306 ix = se & 0x7fffffff; 307 308 /* if Y(n,NaN) is NaN */ 309 if (ix >= 0x7fff0000) 310 { 311 if ((u.words32.w0 & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3) 312 return x + x; 313 } 314 if (x <= 0.0Q) 315 { 316 if (x == 0.0Q) 317 return -HUGE_VALQ + x; 318 if (se & 0x80000000) 319 return zero / (zero * x); 320 } 321 sign = 1; 322 if (n < 0) 323 { 324 n = -n; 325 sign = 1 - ((n & 1) << 1); 326 } 327 if (n == 0) 328 return (y0q (x)); 329 if (n == 1) 330 return (sign * y1q (x)); 331 if (ix >= 0x7fff0000) 332 return zero; 333 if (ix >= 0x412D0000) 334 { /* x > 2**302 */ 335 336 /* ??? See comment above on the possible futility of this. */ 337 338 /* (x >> n**2) 339 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) 340 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) 341 * Let s=sin(x), c=cos(x), 342 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then 343 * 344 * n sin(xn)*sqt2 cos(xn)*sqt2 345 * ---------------------------------- 346 * 0 s-c c+s 347 * 1 -s-c -c+s 348 * 2 -s+c -c-s 349 * 3 s+c c-s 350 */ 351 __float128 s; 352 __float128 c; 353 sincosq (x, &s, &c); 354 switch (n & 3) 355 { 356 case 0: 357 temp = s - c; 358 break; 359 case 1: 360 temp = -s - c; 361 break; 362 case 2: 363 temp = -s + c; 364 break; 365 case 3: 366 temp = s + c; 367 break; 368 } 369 b = invsqrtpi * temp / sqrtq (x); 370 } 371 else 372 { 373 a = y0q (x); 374 b = y1q (x); 375 /* quit if b is -inf */ 376 u.value = b; 377 se = u.words32.w0 & 0xffff0000; 378 for (i = 1; i < n && se != 0xffff0000; i++) 379 { 380 temp = b; 381 b = ((__float128) (i + i) / x) * b - a; 382 u.value = b; 383 se = u.words32.w0 & 0xffff0000; 384 a = temp; 385 } 386 } 387 /* If B is +-Inf, set up errno accordingly. */ 388 if (! finiteq (b)) 389 errno = ERANGE; 390 if (sign > 0) 391 return b; 392 else 393 return -b; 394} 395