*** empty log message ***

miscmaths.cc

@@ -1761,7 +1761,97 @@ int conjgrad(ColumnVector& x, const Matrix& A, const ColumnVector& b, int maxit)
   return conjgrad(x,A,b,maxit,1e-10);
 }
 
+float csevl(const float x, const ColumnVector& cs, const int n)
+  {
+ 
+    float b0 = 0;
+    float b1 = 0;
+    float b2 = 0;
+    const float twox=2*x;
+    
+    for(int i=1; i<=n; i++)
+      {
+	b2=b1;
+	b1=b0;
+	b0=twox*b1-b2+cs(n+1-i);
+      }
+    
+    return 0.5*(b0-b2);
+  }
 
+  float digamma(const float x)
+  { 
+    int ntapsi(16);
+    int ntpsi(23);
+    ColumnVector psics(ntpsi);
+    ColumnVector apsics(ntapsi);
+
+    psics << -.038057080835217922E0<<
+      .49141539302938713E0<<
+      -.056815747821244730E0<<
+      .008357821225914313E0<<
+      -.001333232857994342E0<<
+      .000220313287069308E0<<
+      -.000037040238178456E0<<
+      .000006283793654854E0<<
+      -.000001071263908506E0<<
+      .000000183128394654E0<<
+      -.000000031353509361E0<<
+      .000000005372808776E0<<
+      -.000000000921168141E0<<
+      .000000000157981265E0<<
+      -.000000000027098646E0<<
+      .000000000004648722E0<<
+      -.000000000000797527E0<<
+      .000000000000136827E0<<
+      -.000000000000023475E0<<
+      .000000000000004027E0<<
+      -.000000000000000691E0<<
+      .000000000000000118E0<<
+      -.000000000000000020E0;
+	  
+    apsics <<-.0204749044678185E0<<
+      -.0101801271534859E0<<
+      .0000559718725387E0<<
+      -.0000012917176570E0<<
+      .0000000572858606E0<<
+      -.0000000038213539E0<<
+      .0000000003397434E0<<
+      -.0000000000374838E0<<
+      .0000000000048990E0<<
+      -.0000000000007344E0<<
+      .0000000000001233E0<<
+      -.0000000000000228E0<<
+      .0000000000000045E0<<
+      -.0000000000000009E0<<
+      .0000000000000002E0<<
+      -.0000000000000000E0;
+
+    float y = fabs(x);
+    float psi;
+
+    if(y<2.0)
+      {
+	// do we need to deal with the following case?
+	// c psi(x) for -2. .lt. x .lt. 2.
+
+	int n = int(floor(x));
+	y = x - n;
+	n = n - 1;
+	psi = csevl(2*y-1, psics, ntpsi);
+	if(n==-1)
+	  {
+	    psi = psi - 1.0/x;
+	  }
+      }
+    else
+      {
+	const float aux = csevl(8/(Sqr(y))-1, apsics, ntapsi);
+	psi = log(x) - 0.5/x + aux;
+      }
+    
+    return psi;
+  }
 
 
 vector<float> ColumnVector2vector(const ColumnVector& col)

miscmaths.h

@@ -214,7 +214,9 @@ namespace MISCMATHS {
   //  (stops when error < reltol * initial error)
   int conjgrad(ColumnVector& x, const Matrix& A, const ColumnVector& b, 
 	       int maxit, float reltol);
-  
+
+  float csevl(const float x, const ColumnVector& cs, const int n);
+  float digamma(const float x);
 
   vector<float> ColumnVector2vector(const ColumnVector& col);