Changeset 328

Timestamp:

04/28/09 16:15:04 (15 years ago)

Author:

dedecius

Message:

Implementation of the digamma function.

Location:

bdm

Files:

• itpp_ext.cpp (modified) (1 diff)
• itpp_ext.h (modified) (1 diff)

bdm/itpp_ext.cpp

@@ -298 +298 @@
                 return std::string(tmp);
         };
+
+// digamma
+// copied from C. Bonds' source
+#include <math.h>
+#define el 0.5772156649015329
+
+double psi(double x) {
+    double s,ps,xa,x2;
+    int n,k;
+    static double a[] = {
+        -0.8333333333333e-01,
+         0.83333333333333333e-02,
+        -0.39682539682539683e-02,
+         0.41666666666666667e-02,
+        -0.75757575757575758e-02,
+         0.21092796092796093e-01,
+        -0.83333333333333333e-01,
+         0.4432598039215686};
+
+    xa = fabs(x);
+    s = 0.0;
+    if ((x == (int)x) && (x <= 0.0)) {
+        ps = 1e308;
+        return ps;
+    }
+    if (xa == (int)xa) {
+        n = xa;
+        for (k=1;k<n;k++) {
+            s += 1.0/k;
+        }
+        ps =  s-el;
+    }
+    else if ((xa+0.5) == ((int)(xa+0.5))) {
+        n = xa-0.5;
+        for (k=1;k<=n;k++) {
+            s += 1.0/(2.0*k-1.0);
+        }
+        ps = 2.0*s-el-1.386294361119891;
+    }
+    else {
+        if (xa < 10.0) {
+            n = 10-(int)xa;
+            for (k=0;k<n;k++) {
+                s += 1.0/(xa+k);
+            }
+            xa += n;
+        }
+        x2 = 1.0/(xa*xa);
+        ps = log(xa)-0.5/xa+x2*(((((((a[7]*x2+a[6])*x2+a[5])*x2+
+            a[4])*x2+a[3])*x2+a[2])*x2+a[1])*x2+a[0]);
+        ps -= s;
+    }
+    if (x < 0.0)
+        ps = ps - M_PI*std::cos(M_PI*x)/std::sin(M_PI*x)-1.0/x;
+        return ps;
 }
+
+}

bdm/itpp_ext.h

@@ -80 +80 @@
         //! reimplementation of matlabs num2str
         std::string num2str(int i);
+
+        //! implementation of digamma (psi) function
+        double psi(double);
 }