Revision 329, 1.7 kB
(checked in by dedecius, 15 years ago)
|
Original source file (digamma /psi/ function).
|
Line | |
---|
1 | // psi.cpp -- psi function for real arguments. |
---|
2 | // Algorithms and coefficient values from "Computation of Special |
---|
3 | // Functions", Zhang and Jin, John Wiley and Sons, 1996. |
---|
4 | // |
---|
5 | // (C) 2003, C. Bond. All rights reserved. |
---|
6 | // |
---|
7 | // Returns psi function for real argument 'x'. |
---|
8 | // NOTE: Returns 1e308 for argument 0 or a negative integer. |
---|
9 | // |
---|
10 | //#include <iostream.h> |
---|
11 | #include <math.h> |
---|
12 | #define el 0.5772156649015329 |
---|
13 | |
---|
14 | double psi(double x) |
---|
15 | { |
---|
16 | double s,ps,xa,x2; |
---|
17 | int n,k; |
---|
18 | static double a[] = { |
---|
19 | -0.8333333333333e-01, |
---|
20 | 0.83333333333333333e-02, |
---|
21 | -0.39682539682539683e-02, |
---|
22 | 0.41666666666666667e-02, |
---|
23 | -0.75757575757575758e-02, |
---|
24 | 0.21092796092796093e-01, |
---|
25 | -0.83333333333333333e-01, |
---|
26 | 0.4432598039215686}; |
---|
27 | |
---|
28 | xa = fabs(x); |
---|
29 | s = 0.0; |
---|
30 | if ((x == (int)x) && (x <= 0.0)) { |
---|
31 | ps = 1e308; |
---|
32 | return ps; |
---|
33 | } |
---|
34 | if (xa == (int)xa) { |
---|
35 | n = xa; |
---|
36 | for (k=1;k<n;k++) { |
---|
37 | s += 1.0/k; |
---|
38 | } |
---|
39 | ps = s-el; |
---|
40 | } |
---|
41 | else if ((xa+0.5) == ((int)(xa+0.5))) { |
---|
42 | n = xa-0.5; |
---|
43 | for (k=1;k<=n;k++) { |
---|
44 | s += 1.0/(2.0*k-1.0); |
---|
45 | } |
---|
46 | ps = 2.0*s-el-1.386294361119891; |
---|
47 | } |
---|
48 | else { |
---|
49 | if (xa < 10.0) { |
---|
50 | n = 10-(int)xa; |
---|
51 | for (k=0;k<n;k++) { |
---|
52 | s += 1.0/(xa+k); |
---|
53 | } |
---|
54 | xa += n; |
---|
55 | } |
---|
56 | x2 = 1.0/(xa*xa); |
---|
57 | ps = log(xa)-0.5/xa+x2*(((((((a[7]*x2+a[6])*x2+a[5])*x2+ |
---|
58 | a[4])*x2+a[3])*x2+a[2])*x2+a[1])*x2+a[0]); |
---|
59 | ps -= s; |
---|
60 | } |
---|
61 | if (x < 0.0) |
---|
62 | ps = ps - M_PI*cos(M_PI*x)/sin(M_PI*x)-1.0/x; |
---|
63 | return ps; |
---|
64 | } |
---|