# HG changeset patch
# User Carnë Draug <carandraug@octave.org>
# Date 1430618368 -3600
# Node ID bd565f3e0ecb681c8807a8daba8b8dfc955d63d9
# Parent  1fae49e34a1af4507a545ad7604562bab6c2a42f
psi: refactor to reduce code duplication.

* liboctave/numeric/lo-specfun.cc (psi): reduce code duplication by creating
Lanczos_approximation_psi() to be shared between the functions for complex
and real input.  Also invert the order of the for loops (to sum smaller
numbers first, which will increase accuracy.
* libinterp/corefcn/psi.cc: adjust tolerance for the tests.

diff --git a/libinterp/corefcn/psi.cc b/libinterp/corefcn/psi.cc
--- a/libinterp/corefcn/psi.cc
+++ b/libinterp/corefcn/psi.cc
@@ -127,16 +127,16 @@
 
 ## Interesting identities of the digamma function, in section of 5.1.3
 %!assert (psi (1/3), - em - (3/2) * log(3) - ((sqrt (3) / 6) * pi), eps*10)
-%!assert (psi (1/4), - em -3 * log (2) - pi /2, eps*10)
+%!assert (psi (1/4), - em -3 * log (2) - pi /2)
 %!assert (psi (1/6), - em -2 * log (2) - (3/2) * log (3) - ((sqrt (3) / 2) * pi), eps*10)
 
 ## First 6 zeros of the digamma function, in section of 5.1.5 (and also on
 ## Abramowitz and Stegun, page 258, eq 6.3.19)
 %!assert (psi ( 1.46163214496836234126265954232572132846819620400644), 0, eps)
-%!assert (psi (-0.504083008264455409258269304533302498955385182368579), 0, eps*10)
+%!assert (psi (-0.504083008264455409258269304533302498955385182368579), 0, eps)
 %!assert (psi (-1.573498473162390458778286043690434612655040859116846), 0, eps)
 %!assert (psi (-2.610720868444144650001537715718724207951074010873480), 0, eps*10)
-%!assert (psi (-3.635293366436901097839181566946017713948423861193530), 0, eps)
+%!assert (psi (-3.635293366436901097839181566946017713948423861193530), 0, eps*10)
 %!assert (psi (-4.653237761743142441714598151148207363719069416133868), 0, eps*100)
 
 ## Tests for complex values
@@ -150,10 +150,9 @@
 %!assert (imag (psi (i*z)), 1/2 .* 1./z + 1/2 * pi * coth (pi * z), eps *10)
 
 ## Abramowitz and Stegun, page 259 eq 6.3.12
-%!assert (imag (psi (1/2 + i*z)), 1/2 * pi * tanh (pi * z), eps*10)
+%!assert (imag (psi (1/2 + i*z)), 1/2 * pi * tanh (pi * z), eps)
 
 ## Abramowitz and Stegun, page 259 eq 6.3.13
 %!assert (imag (psi (1 + i*z)), - 1./(2*z) + 1/2 * pi * coth (pi * z), eps*10)
 
 */
-
diff --git a/liboctave/numeric/lo-specfun.cc b/liboctave/numeric/lo-specfun.cc
--- a/liboctave/numeric/lo-specfun.cc
+++ b/liboctave/numeric/lo-specfun.cc
@@ -3728,23 +3728,39 @@
     }
 }
 
-static const double euler_mascheroni = 0.577215664901532860606512090082402431042;
 static const double pi = 3.14159265358979323846;
-// Coefficients for C.Lanczos expansion of psi function from XLiFE++ gammaFunctions
-// psi_coef[k] = - (2k+1) * lg_coef[k] (see melina++ gamma functions)
-// -1/12, 3/360,-5/1260, 7/1680,-9/1188, 11*691/360360,-13/156, 15*3617/122400, ? , ?
-static const double psi_coeff[10] = {
-  -0.83333333333333333e-1, 0.83333333333333333e-2,
-  -0.39682539682539683e-2, 0.41666666666666667e-2,
-  -0.75757575757575758e-2, 0.21092796092796093e-1,
-  -0.83333333333333333e-1, 0.4432598039215686,
-  -0.3053954330270122e+1,  0.125318899521531e+2
-};
+
+template<class T>
+T
+Lanczos_approximation_psi (const T zc)
+{
+  // Coefficients for C.Lanczos expansion of psi function from XLiFE++ gammaFunctions
+  // psi_coef[k] = - (2k+1) * lg_coef[k] (see melina++ gamma functions)
+  // -1/12, 3/360,-5/1260, 7/1680,-9/1188, 11*691/360360,-13/156, 15*3617/122400, ? , ?
+  static const T dg_coeff[10] = {
+    -0.83333333333333333e-1, 0.83333333333333333e-2,
+    -0.39682539682539683e-2, 0.41666666666666667e-2,
+    -0.75757575757575758e-2, 0.21092796092796093e-1,
+    -0.83333333333333333e-1, 0.4432598039215686,
+    -0.3053954330270122e+1,  0.125318899521531e+2
+  };
+
+  T overz2  = T (1.0) / (zc * zc);
+  T overz2k = overz2;
+
+  T p = 0;
+  for (octave_idx_type k = 0; k < 10; k++, overz2k *= overz2)
+    p += dg_coeff[k] * overz2k;
+  p += log (zc) - T (0.5) / zc;
+  return p;
+}
 
 template<class T>
 T
 psi (const T& z)
 {
+  static const double euler_mascheroni = 0.577215664901532860606512090082402431042;
+
   const bool is_int = (xfloor (z) == z);
 
   T p = 0;
@@ -3761,15 +3777,13 @@
     {
       // Abramowitz and Stegun, page 258, eq 6.3.2
       p = - euler_mascheroni;
-      const octave_idx_type n = z;
-      for (octave_idx_type k = 1; k < n; k++)
+      for (octave_idx_type k = z - 1; k > 0; k--)
         p += 1.0 / k;
     }
   else if (xfloor (z + 0.5) == z + 0.5)
     {
       // Abramowitz and Stegun, page 258, eq 6.3.3 and 6.3.4
-      const octave_idx_type n = z + 1;
-      for (octave_idx_type k = 1; k < n; k++)
+      for (octave_idx_type k = z; k > 0; k--)
         p += 1.0 / (2 * k - 1);
 
       p = - euler_mascheroni - 2 * log (2) + 2 * (p);
@@ -3782,17 +3796,12 @@
       // Use formula for derivative of LogGamma(z)
       if (z < 10)
         {
-          const octave_idx_type n = 10 - z;
-          for (octave_idx_type k = 0; k < n; k++)
+          const signed char n = 10 - z;
+          for (signed char k = n - 1; k >= 0; k--)
             p -= 1.0 / (k + z);
           zc += n;
         }
-      T overz2  = 1.0 / (zc*zc);
-      T overz2k = overz2;
-
-      p += log (zc) - 0.5 / zc;
-      for (octave_idx_type k = 0; k < 10; k++, overz2k *= overz2)
-        p += psi_coeff[k] * overz2k;
+      p += Lanczos_approximation_psi (zc);
     }
 
   return p;
@@ -3813,46 +3822,35 @@
   P z_r  = z.real ();
   P z_ra = z_r;
 
+  std::complex<T> dgam (0.0, 0.0);
   if (z.imag () == 0)
-    return std::complex<T> (psi (z_r), 0.0);
+    dgam = std::complex<T> (psi (z_r), 0.0);
   else if (z_r < 0)
-    return psi (P (1.0) - z)- (P (pi) / tan (P (pi) * z));
+    dgam = psi (P (1.0) - z)- (P (pi) / tan (P (pi) * z));
   else
     {
       // Use formula for derivative of LogGamma(z)
-
-      std::complex<T> dgam = 0.0;
-      std::complex<T> z_p  = z;
-
-      octave_idx_type n = 0;
-      std::complex<T> z_m = z_p;
+      std::complex<T> z_m = z;
       if (z_ra < 8)
         {
-          n = 8 - octave_idx_type (z_ra);
-          z_m = z_p + std::complex<T> (n, 0.0);
+          unsigned char n = 8 - z_ra;
+          z_m = z + std::complex<T> (n, 0.0);
+
+          // Recurrence formula
+          // for | Re(z) | < 8 , use recursively DiGamma(z) = DiGamma(z+1) - 1/z
+          std::complex<T> z_p = z + P (n - 1);
+          for (unsigned char k = n; k > 0; k--, z_p -= 1.0)
+            dgam -= P (1.0) / z_p;
         }
 
       // for | Re(z) | > 8, use derivative of C.Lanczos expansion for LogGamma
       // psi(z) = log(z) - 1/(2z) - 1/12z^2 + 3/360z^4 - 5/1260z^6 + 7/1680z^8 - 9/1188z^10 + ...
       // (Abramowitz&Stegun, page 259, formula 6.3.18
-      std::complex<T> overz   = P (1.0) / z_m;
-      std::complex<T> overz2  = overz * overz;
-      std::complex<T> overz2k = overz2;
-
-      dgam += log (z_m) - P (0.5) * overz;
-      for (octave_idx_type k = 0; k < 10; k++, overz2k *= overz2)
-        dgam += P (psi_coeff[k]) * overz2k;
-
-      // Recurrence formula
-      // for | Re(z) | < 8 , use recursively DiGamma(z) = DiGamma(z+1) - 1/z
-      for (octave_idx_type k = 0; k < n; k++, z_p += 1.0)
-        dgam -= P (1.0) / z_p;
-
-      return dgam;
+      dgam += Lanczos_approximation_psi (z_m);
     }
+  return dgam;
 }
 
 // explicit instantiations
 template Complex psi<double> (const Complex& z);
 template FloatComplex psi<float> (const FloatComplex& z);
-