--- a/NEWS	Fri May 01 13:14:23 2015 -0700
+++ b/NEWS	Sun Mar 15 03:31:16 2015 +0000
@@ -228,6 +228,7 @@
       open
       ordschur
       pan
+      psi
       qmr
       rotate
       rotate3d

--- a/doc/interpreter/arith.txi	Fri May 01 13:14:23 2015 -0700
+++ b/doc/interpreter/arith.txi	Sun Mar 15 03:31:16 2015 +0000
@@ -315,6 +315,8 @@
 @anchor{XREFgammaln}
 @DOCSTRING(lgamma)
 
+@DOCSTRING(psi)
+
 @node Rational Approximations
 @section Rational Approximations
 

--- a/libinterp/corefcn/module.mk	Fri May 01 13:14:23 2015 -0700
+++ b/libinterp/corefcn/module.mk	Sun Mar 15 03:31:16 2015 +0000
@@ -233,6 +233,7 @@
   corefcn/pr-output.cc \
   corefcn/procstream.cc \
   corefcn/profiler.cc \
+  corefcn/psi.cc \
   corefcn/quad.cc \
   corefcn/quadcc.cc \
   corefcn/qz.cc \

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/psi.cc	Sun Mar 15 03:31:16 2015 +0000
@@ -0,0 +1,130 @@
+/*
+
+Copyright (C) 2015 Carnë Draug
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "ov.h"
+#include "defun.h"
+#include "error.h"
+#include "dNDArray.h"
+#include "fNDArray.h"
+
+#include "lo-specfun.h"
+
+DEFUN (psi, args, ,
+"-*- texinfo -*-\n\
+@deftypefn {Function File} {} psi (@var{z})\n\
+Compute the psi (digamma) function.\n\
+\n\
+@tex\n\
+$$\n\
+\\Psi (z) = {d (log (\\Gamma (z))) \\over dx}\n\
+$$\n\
+@end tex\n\
+@ifnottex\n\
+@example\n\
+@group\n\
+psi (z) = d (log (gamma (z))) / dx\n\
+@end group\n\
+@end example\n\
+@end ifnottex\n\
+\n\
+@seealso{gamma, gammainc, gammaln}\n\
+@end deftypefn")
+{
+  octave_value retval;
+
+  const octave_idx_type nargin = args.length ();
+
+  if (nargin != 1)
+    {
+      print_usage ();
+      return retval;
+    }
+
+#define FLOAT_BRANCH(T, A, M, E) \
+  if (args(0).is_ ## T ##_type ()) \
+    { \
+      const A ## NDArray z = args(0).M ## array_value (); \
+      A ## NDArray psi_z (z.dims ()); \
+\
+      const E* zv = z.data (); \
+      E* psi_zv = psi_z.fortran_vec (); \
+      const octave_idx_type n = z.numel (); \
+      for (octave_idx_type i = 0; i < n; i++) \
+        psi_zv[i] = psi<E> (zv[i]); \
+\
+      retval = psi_z; \
+    } 
+
+  FLOAT_BRANCH(double, , , double)
+  else FLOAT_BRANCH(single, Float, float_, float)
+  else
+    {
+      error ("psi: Z must be a floating point");
+    }
+
+#undef FLOAT_BRANCH
+
+  return retval;
+}
+
+/*
+%!shared em
+%! em = 0.577215664901532860606512090082402431042; # Euler-Mascheroni Constant
+
+%!assert (psi (ones (7, 3, 5)), repmat (-em, [7 3 5]))
+%!assert (psi ([0 1]), [-Inf -em])
+%!assert (psi ([-20:1]), [repmat(-Inf, [1 21]) -em])
+%!assert (psi (single ([0 1])), single ([-Inf -em]))
+
+## Abramowitz and Stegun, page 258, eq 6.3.5
+%!test
+%! z = [-10:.1:-.1 .1:.1:20]; # drop the 0
+%! assert (psi (z + 1), psi (z) + 1 ./ z, eps*1000)
+
+## Abramowitz and Stegun, page 258, eq 6.3.2
+%!assert (psi (1), -em)
+
+## Abramowitz and Stegun, page 258, eq 6.3.3
+%!assert (psi (1/2), -em - 2 * log (2))
+
+## The following tests are from Pascal Sebah and Xavier Gourdon (2002)
+## "Introduction to the Gamma Function"
+
+## 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/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 (-1.573498473162390458778286043690434612655040859116846), 0, eps)
+%!assert (psi (-2.610720868444144650001537715718724207951074010873480), 0, eps*10)
+%!assert (psi (-3.635293366436901097839181566946017713948423861193530), 0, eps)
+%!assert (psi (-4.653237761743142441714598151148207363719069416133868), 0, eps*100)
+*/
+

--- a/liboctave/numeric/lo-specfun.cc	Fri May 01 13:14:23 2015 -0700
+++ b/liboctave/numeric/lo-specfun.cc	Sun Mar 15 03:31:16 2015 +0000
@@ -1,8 +1,10 @@
 /*
 
 Copyright (C) 1996-2015 John W. Eaton
+Copyright (C) 2007-2010 D. Martin
 Copyright (C) 2010 Jaroslav Hajek
 Copyright (C) 2010 VZLU Prague
+Copyright (C) 2015 Carnë Draug
 
 This file is part of Octave.
 
@@ -45,6 +47,7 @@
 #include "lo-specfun.h"
 #include "mx-inlines.cc"
 #include "lo-mappers.h"
+#include "lo-math.h"
 
 #include "Faddeeva.hh"
 
@@ -3724,3 +3727,77 @@
       dn = Complex (dd*cc1*dd1/ddd, -m*ss*cc*ss1/ddd);
     }
 }
+
+template<class T>
+T
+psi (const T& z)
+{
+  const T euler_mascheroni = 0.577215664901532860606512090082402431042;
+  const T pi = 3.14159265358979323846;
+
+  const bool is_int = (xfloor (z) == z);
+
+  T p = 0;
+  if (z <= 0)
+    {
+      // limits - zeros of the gamma function
+      if (is_int)
+        p = -octave_Inf; // Matlab returns -Inf for psi (0)
+      else
+        // Abramowitz and Stegun, page 259, eq 6.3.7
+        p = psi (1 - z) - (pi / tan (pi * z));
+    }
+  else if (is_int)
+    {
+      // 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++)
+        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++)
+        p += 1.0 / (2 * k - 1);
+
+      p = - euler_mascheroni - 2 * log (2) + 2 * (p);
+    }
+  else
+    {
+      // adapted from XLiFE++ gammaFunctions
+
+      // Coefficients for C.Lanczos expansion of DiGamma function
+      // dg_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, ? , ?
+      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 zc = z;
+      // 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++)
+            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 += dg_coeff[k] * overz2k;
+    }
+
+  return p;
+}
+
+// explicit instantiations
+template double psi<double> (const double& z);
+template float  psi<float> (const float& z);
+

--- a/liboctave/numeric/lo-specfun.h	Fri May 01 13:14:23 2015 -0700
+++ b/liboctave/numeric/lo-specfun.h	Sun Mar 15 03:31:16 2015 +0000
@@ -663,4 +663,7 @@
 ellipj (const Complex& u, double m, Complex& sn, Complex& cn, Complex& dn,
         double& err);
 
+template<class T>
+extern OCTAVE_API T psi (const T& z);
+
 #endif