Revision 477, 1.6 kB
(checked in by mido, 15 years ago)
|
panove, vite, jak jsem peclivej na upravu kodu.. snad se vam bude libit:) konfigurace je v souboru /system/astylerc
|
Rev | Line | |
---|
[329] | 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 | |
---|
[477] | 14 | double psi ( double x ) { |
---|
| 15 | double s, ps, xa, x2; |
---|
| 16 | int n, k; |
---|
| 17 | static double a[] = { |
---|
| 18 | -0.8333333333333e-01, |
---|
| 19 | 0.83333333333333333e-02, |
---|
| 20 | -0.39682539682539683e-02, |
---|
| 21 | 0.41666666666666667e-02, |
---|
| 22 | -0.75757575757575758e-02, |
---|
| 23 | 0.21092796092796093e-01, |
---|
| 24 | -0.83333333333333333e-01, |
---|
| 25 | 0.4432598039215686 |
---|
| 26 | }; |
---|
[329] | 27 | |
---|
[477] | 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 | } else if ( ( xa + 0.5 ) == ( ( int ) ( xa + 0.5 ) ) ) { |
---|
| 41 | n = xa - 0.5; |
---|
| 42 | for ( k = 1; k <= n; k++ ) { |
---|
| 43 | s += 1.0 / ( 2.0 * k - 1.0 ); |
---|
| 44 | } |
---|
| 45 | ps = 2.0 * s - el - 1.386294361119891; |
---|
| 46 | } else { |
---|
| 47 | if ( xa < 10.0 ) { |
---|
| 48 | n = 10 - ( int ) xa; |
---|
| 49 | for ( k = 0; k < n; k++ ) { |
---|
| 50 | s += 1.0 / ( xa + k ); |
---|
| 51 | } |
---|
| 52 | xa += n; |
---|
| 53 | } |
---|
| 54 | x2 = 1.0 / ( xa * xa ); |
---|
| 55 | ps = log ( xa ) - 0.5 / xa + x2 * ( ( ( ( ( ( ( a[7] * x2 + a[6] ) * x2 + a[5] ) * x2 + |
---|
| 56 | a[4] ) * x2 + a[3] ) * x2 + a[2] ) * x2 + a[1] ) * x2 + a[0] ); |
---|
| 57 | ps -= s; |
---|
| 58 | } |
---|
| 59 | if ( x < 0.0 ) |
---|
| 60 | ps = ps - M_PI * cos ( M_PI * x ) / sin ( M_PI * x ) - 1.0 / x; |
---|
| 61 | return ps; |
---|
[329] | 62 | } |
---|