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