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