← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/special/cephes/jn.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/special/cephes/jn.c
1
 
/*                                                      jn.c
2
 
 *
3
 
 *      Bessel function of integer order
4
 
 *
5
 
 *
6
 
 *
7
 
 * SYNOPSIS:
8
 
 *
9
 
 * int n;
10
 
 * double x, y, jn();
11
 
 *
12
 
 * y = jn( n, x );
13
 
 *
14
 
 *
15
 
 *
16
 
 * DESCRIPTION:
17
 
 *
18
 
 * Returns Bessel function of order n, where n is a
19
 
 * (possibly negative) integer.
20
 
 *
21
 
 * The ratio of jn(x) to j0(x) is computed by backward
22
 
 * recurrence.  First the ratio jn/jn-1 is found by a
23
 
 * continued fraction expansion.  Then the recurrence
24
 
 * relating successive orders is applied until j0 or j1 is
25
 
 * reached.
26
 
 *
27
 
 * If n = 0 or 1 the routine for j0 or j1 is called
28
 
 * directly.
29
 
 *
30
 
 *
31
 
 *
32
 
 * ACCURACY:
33
 
 *
34
 
 *                      Absolute error:
35
 
 * arithmetic   range      # trials      peak         rms
36
 
 *    DEC       0, 30        5500       6.9e-17     9.3e-18
37
 
 *    IEEE      0, 30        5000       4.4e-16     7.9e-17
38
 
 *
39
 
 *
40
 
 * Not suitable for large n or x. Use jv() instead.
41
 
 *
42
 
 */
43
 
 
44
 
/*                                                      jn.c
45
 
Cephes Math Library Release 2.8:  June, 2000
46
 
Copyright 1984, 1987, 2000 by Stephen L. Moshier
47
 
*/
48
 
#include "mconf.h"
49
 
#ifdef ANSIPROT
50
 
extern double fabs ( double );
51
 
extern double j0 ( double );
52
 
extern double j1 ( double );
53
 
#else
54
 
double fabs(), j0(), j1();
55
 
#endif
56
 
extern double MACHEP;
57
 
 
58
 
double jn( n, x )
59
 
int n;
60
 
double x;
61
 
{
62
 
double pkm2, pkm1, pk, xk, r, ans;
63
 
int k, sign;
64
 
 
65
 
if( n < 0 )
66
 
        {
67
 
        n = -n;
68
 
        if( (n & 1) == 0 )      /* -1**n */
69
 
                sign = 1;
70
 
        else
71
 
                sign = -1;
72
 
        }
73
 
else
74
 
        sign = 1;
75
 
 
76
 
if( x < 0.0 )
77
 
        {
78
 
        if( n & 1 )
79
 
                sign = -sign;
80
 
        x = -x;
81
 
        }
82
 
 
83
 
if( n == 0 )
84
 
        return( sign * j0(x) );
85
 
if( n == 1 )
86
 
        return( sign * j1(x) );
87
 
if( n == 2 ) {
88
 
    if (x < 1e-5) {
89
 
        double y = x*x;
90
 
        return sign * 0.125 * y * (1 - y / 12.);
91
 
    } else {
92
 
        return( sign * (2.0 * j1(x) / x  -  j0(x)) );
93
 
    }
94
 
}
95
 
 
96
 
if( x < MACHEP )
97
 
        return( 0.0 );
98
 
 
99
 
/* continued fraction */
100
 
#ifdef DEC
101
 
k = 56;
102
 
#else
103
 
k = 53;
104
 
#endif
105
 
 
106
 
pk = 2 * (n + k);
107
 
ans = pk;
108
 
xk = x * x;
109
 
 
110
 
do
111
 
        {
112
 
        pk -= 2.0;
113
 
        ans = pk - (xk/ans);
114
 
        }
115
 
while( --k > 0 );
116
 
ans = x/ans;
117
 
 
118
 
/* backward recurrence */
119
 
 
120
 
pk = 1.0;
121
 
pkm1 = 1.0/ans;
122
 
k = n-1;
123
 
r = 2 * k;
124
 
 
125
 
do
126
 
        {
127
 
        pkm2 = (pkm1 * r  -  pk * x) / x;
128
 
        pk = pkm1;
129
 
        pkm1 = pkm2;
130
 
        r -= 2.0;
131
 
        }
132
 
while( --k > 0 );
133
 
 
134
 
if( fabs(pk) > fabs(pkm1) )
135
 
        ans = j1(x)/pk;
136
 
else
137
 
        ans = j0(x)/pkm1;
138
 
return( sign * ans );
139
 
}