← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/python-scipy/oneiric-proposed

« back to all changes in this revision

Viewing changes to Lib/special/cephes/psi.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/psi.c
1
 
/*                                                      psi.c
2
 
 *
3
 
 *      Psi (digamma) function
4
 
 *
5
 
 *
6
 
 * SYNOPSIS:
7
 
 *
8
 
 * double x, y, psi();
9
 
 *
10
 
 * y = psi( x );
11
 
 *
12
 
 *
13
 
 * DESCRIPTION:
14
 
 *
15
 
 *              d      -
16
 
 *   psi(x)  =  -- ln | (x)
17
 
 *              dx
18
 
 *
19
 
 * is the logarithmic derivative of the gamma function.
20
 
 * For integer x,
21
 
 *                   n-1
22
 
 *                    -
23
 
 * psi(n) = -EUL  +   >  1/k.
24
 
 *                    -
25
 
 *                   k=1
26
 
 *
27
 
 * This formula is used for 0 < n <= 10.  If x is negative, it
28
 
 * is transformed to a positive argument by the reflection
29
 
 * formula  psi(1-x) = psi(x) + pi cot(pi x).
30
 
 * For general positive x, the argument is made greater than 10
31
 
 * using the recurrence  psi(x+1) = psi(x) + 1/x.
32
 
 * Then the following asymptotic expansion is applied:
33
 
 *
34
 
 *                           inf.   B
35
 
 *                            -      2k
36
 
 * psi(x) = log(x) - 1/2x -   >   -------
37
 
 *                            -        2k
38
 
 *                           k=1   2k x
39
 
 *
40
 
 * where the B2k are Bernoulli numbers.
41
 
 *
42
 
 * ACCURACY:
43
 
 *    Relative error (except absolute when |psi| < 1):
44
 
 * arithmetic   domain     # trials      peak         rms
45
 
 *    DEC       0,30         2500       1.7e-16     2.0e-17
46
 
 *    IEEE      0,30        30000       1.3e-15     1.4e-16
47
 
 *    IEEE      -30,0       40000       1.5e-15     2.2e-16
48
 
 *
49
 
 * ERROR MESSAGES:
50
 
 *     message         condition      value returned
51
 
 * psi singularity    x integer <=0      MAXNUM
52
 
 */
53
 
 
54
 
/*
55
 
Cephes Math Library Release 2.8:  June, 2000
56
 
Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier
57
 
*/
58
 
 
59
 
#include "mconf.h"
60
 
 
61
 
#ifdef UNK
62
 
static double A[] = {
63
 
 8.33333333333333333333E-2,
64
 
-2.10927960927960927961E-2,
65
 
 7.57575757575757575758E-3,
66
 
-4.16666666666666666667E-3,
67
 
 3.96825396825396825397E-3,
68
 
-8.33333333333333333333E-3,
69
 
 8.33333333333333333333E-2
70
 
};
71
 
#endif
72
 
 
73
 
#ifdef DEC
74
 
static unsigned short A[] = {
75
 
0037252,0125252,0125252,0125253,
76
 
0136654,0145314,0126312,0146255,
77
 
0036370,0037017,0101740,0174076,
78
 
0136210,0104210,0104210,0104211,
79
 
0036202,0004040,0101010,0020202,
80
 
0136410,0104210,0104210,0104211,
81
 
0037252,0125252,0125252,0125253
82
 
};
83
 
#endif
84
 
 
85
 
#ifdef IBMPC
86
 
static unsigned short A[] = {
87
 
0x5555,0x5555,0x5555,0x3fb5,
88
 
0x5996,0x9599,0x9959,0xbf95,
89
 
0x1f08,0xf07c,0x07c1,0x3f7f,
90
 
0x1111,0x1111,0x1111,0xbf71,
91
 
0x0410,0x1041,0x4104,0x3f70,
92
 
0x1111,0x1111,0x1111,0xbf81,
93
 
0x5555,0x5555,0x5555,0x3fb5
94
 
};
95
 
#endif
96
 
 
97
 
#ifdef MIEEE
98
 
static unsigned short A[] = {
99
 
0x3fb5,0x5555,0x5555,0x5555,
100
 
0xbf95,0x9959,0x9599,0x5996,
101
 
0x3f7f,0x07c1,0xf07c,0x1f08,
102
 
0xbf71,0x1111,0x1111,0x1111,
103
 
0x3f70,0x4104,0x1041,0x0410,
104
 
0xbf81,0x1111,0x1111,0x1111,
105
 
0x3fb5,0x5555,0x5555,0x5555
106
 
};
107
 
#endif
108
 
 
109
 
#define EUL 0.57721566490153286061
110
 
 
111
 
#ifdef ANSIPROT
112
 
extern double floor ( double );
113
 
extern double log ( double );
114
 
extern double tan ( double );
115
 
extern double polevl ( double, void *, int );
116
 
#else
117
 
double floor(), log(), tan(), polevl();
118
 
#endif
119
 
extern double PI, MAXNUM;
120
 
 
121
 
 
122
 
double psi(x)
123
 
double x;
124
 
{
125
 
double p, q, nz, s, w, y, z;
126
 
int i, n, negative;
127
 
 
128
 
negative = 0;
129
 
nz = 0.0;
130
 
 
131
 
if( x <= 0.0 )
132
 
        {
133
 
        negative = 1;
134
 
        q = x;
135
 
        p = floor(q);
136
 
        if( p == q )
137
 
                {
138
 
                mtherr( "psi", SING );
139
 
                return( MAXNUM );
140
 
                }
141
 
/* Remove the zeros of tan(PI x)
142
 
 * by subtracting the nearest integer from x
143
 
 */
144
 
        nz = q - p;
145
 
        if( nz != 0.5 )
146
 
                {
147
 
                if( nz > 0.5 )
148
 
                        {
149
 
                        p += 1.0;
150
 
                        nz = q - p;
151
 
                        }
152
 
                nz = PI/tan(PI*nz);
153
 
                }
154
 
        else
155
 
                {
156
 
                nz = 0.0;
157
 
                }
158
 
        x = 1.0 - x;
159
 
        }
160
 
 
161
 
/* check for positive integer up to 10 */
162
 
if( (x <= 10.0) && (x == floor(x)) )
163
 
        {
164
 
        y = 0.0;
165
 
        n = x;
166
 
        for( i=1; i<n; i++ )
167
 
                {
168
 
                w = i;
169
 
                y += 1.0/w;
170
 
                }
171
 
        y -= EUL;
172
 
        goto done;
173
 
        }
174
 
 
175
 
s = x;
176
 
w = 0.0;
177
 
while( s < 10.0 )
178
 
        {
179
 
        w += 1.0/s;
180
 
        s += 1.0;
181
 
        }
182
 
 
183
 
if( s < 1.0e17 )
184
 
        {
185
 
        z = 1.0/(s * s);
186
 
        y = z * polevl( z, A, 6 );
187
 
        }
188
 
else
189
 
        y = 0.0;
190
 
 
191
 
y = log(s)  -  (0.5/s)  -  y  -  w;
192
 
 
193
 
done:
194
 
 
195
 
if( negative )
196
 
        {
197
 
        y -= nz;
198
 
        }
199
 
 
200
 
return(y);
201
 
}