POK
e_jnf.c
1 /*
2  * POK header
3  *
4  * The following file is a part of the POK project. Any modification should
5  * made according to the POK licence. You CANNOT use this file or a part of
6  * this file is this part of a file for your own project
7  *
8  * For more information on the POK licence, please see our LICENCE FILE
9  *
10  * Please follow the coding guidelines described in doc/CODING_GUIDELINES
11  *
12  * Copyright (c) 2007-2009 POK team
13  *
14  * Created by julien on Fri Jan 30 14:41:34 2009
15  */
16 
17 /* e_jnf.c -- float version of e_jn.c.
18  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
19  */
20 
21 /*
22  * ====================================================
23  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
24  *
25  * Developed at SunPro, a Sun Microsystems, Inc. business.
26  * Permission to use, copy, modify, and distribute this
27  * software is freely granted, provided that this notice
28  * is preserved.
29  * ====================================================
30  */
31 
32 #ifdef POK_NEEDS_LIBMATH
33 
34 #include "math_private.h"
35 #include <libm.h>
36 
37 static const float
38 #if 0
39 invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */
40 #endif
41 two = 2.0000000000e+00, /* 0x40000000 */
42 one = 1.0000000000e+00; /* 0x3F800000 */
43 
44 static const float zero = 0.0000000000e+00;
45 
46 float
47 __ieee754_jnf(int n, float x)
48 {
49  int32_t i,hx,ix, sgn;
50  float a, b, temp, di;
51  float z, w;
52 
53  /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
54  * Thus, J(-n,x) = J(n,-x)
55  */
56  GET_FLOAT_WORD(hx,x);
57  ix = 0x7fffffff&hx;
58  /* if J(n,NaN) is NaN */
59  if(ix>0x7f800000) return x+x;
60  if(n<0){
61  n = -n;
62  x = -x;
63  hx ^= 0x80000000;
64  }
65  if(n==0) return(__ieee754_j0f(x));
66  if(n==1) return(__ieee754_j1f(x));
67  sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
68  x = fabsf(x);
69  if(ix==0||ix>=0x7f800000) /* if x is 0 or inf */
70  b = zero;
71  else if((float)n<=x) {
72  /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
73  a = __ieee754_j0f(x);
74  b = __ieee754_j1f(x);
75  for(i=1;i<n;i++){
76  temp = b;
77  b = b*((float)(i+i)/x) - a; /* avoid underflow */
78  a = temp;
79  }
80  } else {
81  if(ix<0x30800000) { /* x < 2**-29 */
82  /* x is tiny, return the first Taylor expansion of J(n,x)
83  * J(n,x) = 1/n!*(x/2)^n - ...
84  */
85  if(n>33) /* underflow */
86  b = zero;
87  else {
88  temp = x*(float)0.5; b = temp;
89  for (a=one,i=2;i<=n;i++) {
90  a *= (float)i; /* a = n! */
91  b *= temp; /* b = (x/2)^n */
92  }
93  b = b/a;
94  }
95  } else {
96  /* use backward recurrence */
97  /* x x^2 x^2
98  * J(n,x)/J(n-1,x) = ---- ------ ------ .....
99  * 2n - 2(n+1) - 2(n+2)
100  *
101  * 1 1 1
102  * (for large x) = ---- ------ ------ .....
103  * 2n 2(n+1) 2(n+2)
104  * -- - ------ - ------ -
105  * x x x
106  *
107  * Let w = 2n/x and h=2/x, then the above quotient
108  * is equal to the continued fraction:
109  * 1
110  * = -----------------------
111  * 1
112  * w - -----------------
113  * 1
114  * w+h - ---------
115  * w+2h - ...
116  *
117  * To determine how many terms needed, let
118  * Q(0) = w, Q(1) = w(w+h) - 1,
119  * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
120  * When Q(k) > 1e4 good for single
121  * When Q(k) > 1e9 good for double
122  * When Q(k) > 1e17 good for quadruple
123  */
124  /* determine k */
125  float t,v;
126  float q0,q1,h,tmp; int32_t k,m;
127  w = (n+n)/(float)x; h = (float)2.0/(float)x;
128  q0 = w; z = w+h; q1 = w*z - (float)1.0; k=1;
129  while(q1<(float)1.0e9) {
130  k += 1; z += h;
131  tmp = z*q1 - q0;
132  q0 = q1;
133  q1 = tmp;
134  }
135  m = n+n;
136  for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
137  a = t;
138  b = one;
139  /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
140  * Hence, if n*(log(2n/x)) > ...
141  * single 8.8722839355e+01
142  * double 7.09782712893383973096e+02
143  * long double 1.1356523406294143949491931077970765006170e+04
144  * then recurrent value may overflow and the result is
145  * likely underflow to zero
146  */
147  tmp = n;
148  v = two/x;
149  tmp = tmp*__ieee754_logf(fabsf(v*tmp));
150  if(tmp<(float)8.8721679688e+01) {
151  for(i=n-1,di=(float)(i+i);i>0;i--){
152  temp = b;
153  b *= di;
154  b = b/x - a;
155  a = temp;
156  di -= two;
157  }
158  } else {
159  for(i=n-1,di=(float)(i+i);i>0;i--){
160  temp = b;
161  b *= di;
162  b = b/x - a;
163  a = temp;
164  di -= two;
165  /* scale b to avoid spurious overflow */
166  if(b>(float)1e10) {
167  a /= b;
168  t /= b;
169  b = one;
170  }
171  }
172  }
173  b = (t*__ieee754_j0f(x)/b);
174  }
175  }
176  if(sgn==1) return -b; else return b;
177 }
178 
179 float
180 __ieee754_ynf(int n, float x)
181 {
182  int32_t i,hx,ix,ib;
183  int32_t sign;
184  float a, b, temp;
185 
186  GET_FLOAT_WORD(hx,x);
187  ix = 0x7fffffff&hx;
188  /* if Y(n,NaN) is NaN */
189  if(ix>0x7f800000) return x+x;
190  if(ix==0) return -one/zero;
191  if(hx<0) return zero/zero;
192  sign = 1;
193  if(n<0){
194  n = -n;
195  sign = 1 - ((n&1)<<1);
196  }
197  if(n==0) return(__ieee754_y0f(x));
198  if(n==1) return(sign*__ieee754_y1f(x));
199  if(ix==0x7f800000) return zero;
200 
201  a = __ieee754_y0f(x);
202  b = __ieee754_y1f(x);
203  /* quit if b is -inf */
204  GET_FLOAT_WORD(ib,b);
205  for(i=1;i<n&&(uint32_t)ib!=0xff800000;i++){
206  temp = b;
207  b = ((float)(i+i)/x)*b - a;
208  GET_FLOAT_WORD(ib,b);
209  a = temp;
210  }
211  if(sign>0) return b; else return -b;
212 }
213 
214 #endif
215