URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libgo/] [go/] [math/] [jn.go] - Blame information for rev 747

Details | Compare with Previous | View Log

Line No.	Rev	Author	Line
1	747	jeremybenn	`// Copyright 2010 The Go Authors. All rights reserved.`
2			`// Use of this source code is governed by a BSD-style`
3			`// license that can be found in the LICENSE file.`
4
5			`package math`
6
7			`/*`
8			`Bessel function of the first and second kinds of order n.`
9			`*/`
10
11			`// The original C code and the long comment below are`
12			`// from FreeBSD's /usr/src/lib/msun/src/e_jn.c and`
13			`// came with this notice. The go code is a simplified`
14			`// version of the original C.`
15			`//`
16			`// ====================================================`
17			`// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.`
18			`//`
19			`// Developed at SunPro, a Sun Microsystems, Inc. business.`
20			`// Permission to use, copy, modify, and distribute this`
21			`// software is freely granted, provided that this notice`
22			`// is preserved.`
23			`// ====================================================`
24			`//`
25			`// __ieee754_jn(n, x), __ieee754_yn(n, x)`
26			`// floating point Bessel's function of the 1st and 2nd kind`
27			`// of order n`
28			`//`
29			`// Special cases:`
30			`// y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;`
31			`// y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.`
32			`// Note 2. About jn(n,x), yn(n,x)`
33			`// For n=0, j0(x) is called,`
34			`// for n=1, j1(x) is called,`
35			`// for n`
36			`// from values of j0(x) and j1(x).`
37			`// for n>x, a continued fraction approximation to`
38			`// j(n,x)/j(n-1,x) is evaluated and then backward`
39			`// recursion is used starting from a supposed value`
40			`// for j(n,x). The resulting value of j(0,x) is`
41			`// compared with the actual value to correct the`
42			`// supposed value of j(n,x).`
43			`//`
44			`// yn(n,x) is similar in all respects, except`
45			`// that forward recursion is used for all`
46			`// values of n>1.`
47
48			`// Jn returns the order-n Bessel function of the first kind.`
49			`//`
50			`// Special cases are:`
51			`// Jn(n, ±Inf) = 0`
52			`// Jn(n, NaN) = NaN`
53			`func Jn(n int, x float64) float64 {`
54			`const (`
55			`TwoM29 = 1.0 / (1 << 29) // 2**-29 0x3e10000000000000`
56			`Two302 = 1 << 302 // 2**302 0x52D0000000000000`
57			`)`
58			`// special cases`
59			`switch {`
60			`case IsNaN(x):`
61			`return x`
62			`case IsInf(x, 0):`
63			`return 0`
64			`}`
65			`// J(-n, x) = (-1)*n J(n, x), J(n, -x) = (-1)*n J(n, x)`
66			`// Thus, J(-n, x) = J(n, -x)`
67
68			`if n == 0 {`
69			`return J0(x)`
70			`}`
71			`if x == 0 {`
72			`return 0`
73			`}`
74			`if n < 0 {`
75			`n, x = -n, -x`
76			`}`
77			`if n == 1 {`
78			`return J1(x)`
79			`}`
80			`sign := false`
81			`if x < 0 {`
82			`x = -x`
83			`if n&1 == 1 {`
84			`sign = true // odd n and negative x`
85			`}`
86			`}`
87			`var b float64`
88			`if float64(n) <= x {`
89			`// Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x)`
90			`if x >= Two302 { // x > 2**302`
91
92			`// (x >> n**2)`
93			`// Jn(x) = cos(x-(2n+1)pi/4)sqrt(2/x*pi)`
94			`// Yn(x) = sin(x-(2n+1)pi/4)sqrt(2/x*pi)`
95			`// Let s=sin(x), c=cos(x),`
96			`// xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then`
97			`//`
98			`// n sin(xn)sqt2 cos(xn)sqt2`
99			`// ----------------------------------`
100			`// 0 s-c c+s`
101			`// 1 -s-c -c+s`
102			`// 2 -s+c -c-s`
103			`// 3 s+c c-s`
104
105			`var temp float64`
106			`switch n & 3 {`
107			`case 0:`
108			`temp = Cos(x) + Sin(x)`
109			`case 1:`
110			`temp = -Cos(x) + Sin(x)`
111			`case 2:`
112			`temp = -Cos(x) - Sin(x)`
113			`case 3:`
114			`temp = Cos(x) - Sin(x)`
115			`}`
116			`b = (1 / SqrtPi) * temp / Sqrt(x)`
117			`} else {`
118			`b = J1(x)`
119			`for i, a := 1, J0(x); i < n; i++ {`
120			`a, b = b, b*(float64(i+i)/x)-a // avoid underflow`
121			`}`
122			`}`
123			`} else {`
124			`if x < TwoM29 { // x < 2**-29`
125			`// x is tiny, return the first Taylor expansion of J(n,x)`
126			`// J(n,x) = 1/n!(x/2)*n - ...`
127
128			`if n > 33 { // underflow`
129			`b = 0`
130			`} else {`
131			`temp := x * 0.5`
132			`b = temp`
133			`a := 1.0`
134			`for i := 2; i <= n; i++ {`
135			`a *= float64(i) // a = n!`
136			`b = temp // b = (x/2)*n`
137			`}`
138			`b /= a`
139			`}`
140			`} else {`
141			`// use backward recurrence`
142			`// x x2 x2`
143			`// J(n,x)/J(n-1,x) = ---- ------ ------ .....`
144			`// 2n - 2(n+1) - 2(n+2)`
145			`//`
146			`// 1 1 1`
147			`// (for large x) = ---- ------ ------ .....`
148			`// 2n 2(n+1) 2(n+2)`
149			`// -- - ------ - ------ -`
150			`// x x x`
151			`//`
152			`// Let w = 2n/x and h=2/x, then the above quotient`
153			`// is equal to the continued fraction:`
154			`// 1`
155			`// = -----------------------`
156			`// 1`
157			`// w - -----------------`
158			`// 1`
159			`// w+h - ---------`
160			`// w+2h - ...`
161			`//`
162			`// To determine how many terms needed, let`
163			`// Q(0) = w, Q(1) = w(w+h) - 1,`
164			`// Q(k) = (w+kh)Q(k-1) - Q(k-2),`
165			`// When Q(k) > 1e4 good for single`
166			`// When Q(k) > 1e9 good for double`
167			`// When Q(k) > 1e17 good for quadruple`
168
169			`// determine k`
170			`w := float64(n+n) / x`
171			`h := 2 / x`
172			`q0 := w`
173			`z := w + h`
174			`q1 := w*z - 1`
175			`k := 1`
176			`for q1 < 1e9 {`
177			`k += 1`
178			`z += h`
179			`q0, q1 = q1, z*q1-q0`
180			`}`
181			`m := n + n`
182			`t := 0.0`
183			`for i := 2 * (n + k); i >= m; i -= 2 {`
184			`t = 1 / (float64(i)/x - t)`
185			`}`
186			`a := t`
187			`b = 1`
188			`// estimate log((2/x)*nn!) = nlog(2/x)+nln(n)`
189			`// Hence, if n*(log(2n/x)) > ...`
190			`// single 8.8722839355e+01`
191			`// double 7.09782712893383973096e+02`
192			`// long double 1.1356523406294143949491931077970765006170e+04`
193			`// then recurrent value may overflow and the result is`
194			`// likely underflow to zero`
195
196			`tmp := float64(n)`
197			`v := 2 / x`
198			`tmp = tmp * Log(Abs(v*tmp))`
199			`if tmp < 7.09782712893383973096e+02 {`
200			`for i := n - 1; i > 0; i-- {`
201			`di := float64(i + i)`
202			`a, b = b, b*di/x-a`
203			`di -= 2`
204			`}`
205			`} else {`
206			`for i := n - 1; i > 0; i-- {`
207			`di := float64(i + i)`
208			`a, b = b, b*di/x-a`
209			`di -= 2`
210			`// scale b to avoid spurious overflow`
211			`if b > 1e100 {`
212			`a /= b`
213			`t /= b`
214			`b = 1`
215			`}`
216			`}`
217			`}`
218			`b = t * J0(x) / b`
219			`}`
220			`}`
221			`if sign {`
222			`return -b`
223			`}`
224			`return b`
225			`}`
226
227			`// Yn returns the order-n Bessel function of the second kind.`
228			`//`
229			`// Special cases are:`
230			`// Yn(n, +Inf) = 0`
231			`// Yn(n > 0, 0) = -Inf`
232			`// Yn(n < 0, 0) = +Inf if n is odd, -Inf if n is even`
233			`// Y1(n, x < 0) = NaN`
234			`// Y1(n, NaN) = NaN`
235			`func Yn(n int, x float64) float64 {`
236			`const Two302 = 1 << 302 // 2**302 0x52D0000000000000`
237			`// special cases`
238			`switch {`
239			`case x < 0 \|\| IsNaN(x):`
240			`return NaN()`
241			`case IsInf(x, 1):`
242			`return 0`
243			`}`
244
245			`if n == 0 {`
246			`return Y0(x)`
247			`}`
248			`if x == 0 {`
249			`if n < 0 && n&1 == 1 {`
250			`return Inf(1)`
251			`}`
252			`return Inf(-1)`
253			`}`
254			`sign := false`
255			`if n < 0 {`
256			`n = -n`
257			`if n&1 == 1 {`
258			`sign = true // sign true if n < 0 && \|n\| odd`
259			`}`
260			`}`
261			`if n == 1 {`
262			`if sign {`
263			`return -Y1(x)`
264			`}`
265			`return Y1(x)`
266			`}`
267			`var b float64`
268			`if x >= Two302 { // x > 2**302`
269			`// (x >> n**2)`
270			`// Jn(x) = cos(x-(2n+1)pi/4)sqrt(2/x*pi)`
271			`// Yn(x) = sin(x-(2n+1)pi/4)sqrt(2/x*pi)`
272			`// Let s=sin(x), c=cos(x),`
273			`// xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then`
274			`//`
275			`// n sin(xn)sqt2 cos(xn)sqt2`
276			`// ----------------------------------`
277			`// 0 s-c c+s`
278			`// 1 -s-c -c+s`
279			`// 2 -s+c -c-s`
280			`// 3 s+c c-s`
281
282			`var temp float64`
283			`switch n & 3 {`
284			`case 0:`
285			`temp = Sin(x) - Cos(x)`
286			`case 1:`
287			`temp = -Sin(x) - Cos(x)`
288			`case 2:`
289			`temp = -Sin(x) + Cos(x)`
290			`case 3:`
291			`temp = Sin(x) + Cos(x)`
292			`}`
293			`b = (1 / SqrtPi) * temp / Sqrt(x)`
294			`} else {`
295			`a := Y0(x)`
296			`b = Y1(x)`
297			`// quit if b is -inf`
298			`for i := 1; i < n && !IsInf(b, -1); i++ {`
299			`a, b = b, (float64(i+i)/x)*b-a`
300			`}`
301			`}`
302			`if sign {`
303			`return -b`
304			`}`
305			`return b`
306			`}`

Browse

Tools

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libgo/] [go/] [math/] [jn.go] - Blame information for rev 747