# HG changeset patch
# User John W. Eaton <jwe@octave.org>
# Date 1432390790 14400
# Node ID b5d2f6954c48dc157761da9041e27dd38cfbbf34
# Parent  4e15a4c331e73ccf9899875dd851db33d52003b3# Parent  9866b3202c522da073e19cb58c89a423a0850fb5
maint: Periodic merge of stable to default.

diff -r 9866b3202c52 -r b5d2f6954c48 NEWS
--- a/NEWS	Sat May 23 10:47:03 2015 +0200
+++ b/NEWS	Sat May 23 10:19:50 2015 -0400
@@ -1,3 +1,29 @@
+Summary of important user-visible changes for version 4.2:
+---------------------------------------------------------
+
+ ** Octal ('\NNN') and hex ('\xNN') escape sequences in single quoted
+    strings are now interpreted by the function do_string_escapes().
+    The *printf family of functions now supports octal and hex escape
+    sequences in single-quoted strings for Matlab compatibility.
+
+ ** Special octal and hex escape sequences for the pattern and replacement
+    strings in regular expressions are now interpreted for Matlab compatibility.
+
+    octal: '\oNNN' or '\o{NNN}'
+    hex  : '\xNN'  or '\x{NN}'
+
+ ** mkfifo now interprets the MODE argument as an octal, not decimal, integer.
+    This is consistent with the equivalent shell command. 
+
+ ** The griddata function no longer plots the interpolated mesh if no output
+    argument is requested, instead the vector or array of interpolated values
+    is always returned for Matlab compatibility.
+
+ ** Other new functions added in 4.2:
+
+      psi
+
+
 Summary of important user-visible changes for version 4.0:
 ---------------------------------------------------------
 
diff -r 9866b3202c52 -r b5d2f6954c48 doc/interpreter/arith.txi
--- a/doc/interpreter/arith.txi	Sat May 23 10:47:03 2015 +0200
+++ b/doc/interpreter/arith.txi	Sat May 23 10:19:50 2015 -0400
@@ -315,6 +315,8 @@
 @anchor{XREFgammaln}
 @DOCSTRING(lgamma)
 
+@DOCSTRING(psi)
+
 @node Rational Approximations
 @section Rational Approximations
 
diff -r 9866b3202c52 -r b5d2f6954c48 doc/interpreter/geometry.txi
--- a/doc/interpreter/geometry.txi	Sat May 23 10:47:03 2015 +0200
+++ b/doc/interpreter/geometry.txi	Sat May 23 10:19:50 2015 -0400
@@ -433,7 +433,8 @@
 y = 2*rand (size (x)) - 1;
 z = sin (2*(x.^2+y.^2));
 [xx,yy] = meshgrid (linspace (-1,1,32));
-griddata (x,y,z,xx,yy);
+zz = griddata (x, y, z, xx, yy);
+mesh (xx, yy, zz);
 @end group
 @end example
 
diff -r 9866b3202c52 -r b5d2f6954c48 doc/interpreter/geometryimages.m
--- a/doc/interpreter/geometryimages.m	Sat May 23 10:47:03 2015 +0200
+++ b/doc/interpreter/geometryimages.m	Sat May 23 10:19:50 2015 -0400
@@ -65,7 +65,8 @@
     y = 2 * rand (size (x)) - 1;
     z = sin (2 * (x.^2 + y.^2));
     [xx,yy] = meshgrid (linspace (-1,1,32));
-    griddata (x,y,z,xx,yy);
+    zz = griddata (x, y, z, xx, yy);
+    mesh (xx, yy, zz);
     print ([nm "." typ], d_typ);
   elseif (strcmp (nm, "convhull"))
     x = -3:0.05:3;
diff -r 9866b3202c52 -r b5d2f6954c48 libinterp/corefcn/file-io.cc
--- a/libinterp/corefcn/file-io.cc	Sat May 23 10:47:03 2015 +0200
+++ b/libinterp/corefcn/file-io.cc	Sat May 23 10:19:50 2015 -0400
@@ -2204,6 +2204,8 @@
   return retval;
 }
 
+// FIXME: This routine also exists verbatim in syscalls.cc.
+//        Maybe change to be a general utility routine.
 static int
 convert (int x, int ibase, int obase)
 {
diff -r 9866b3202c52 -r b5d2f6954c48 libinterp/corefcn/module.mk
--- a/libinterp/corefcn/module.mk	Sat May 23 10:47:03 2015 +0200
+++ b/libinterp/corefcn/module.mk	Sat May 23 10:19:50 2015 -0400
@@ -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 \
diff -r 9866b3202c52 -r b5d2f6954c48 libinterp/corefcn/psi.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/psi.cc	Sat May 23 10:19:50 2015 -0400
@@ -0,0 +1,240 @@
+/*
+
+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\
+@deftypefnx {Function File} {} psi (@var{k}, @var{z})\n\
+Compute the psi (polygamma) function.\n\
+\n\
+The polygamma functions are the @var{k}th derivative of the logarithm\n\
+of the gamma function.  If unspecified, @var{k} defaults to zero.  A value\n\
+of zero computes the digamma function, a value of 1, the trigamma function,\n\
+and so on.\n\
+\n\
+The digamma function is defined:\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\
+When computing the digamma function (when @var{k} equals zero), @var{z}\n\
+can have any value real or complex value.  However, for polygamma functions\n\
+(@var{k} higher than 0), @var{z} must be real and non-negative.\n\
+\n\
+@seealso{gamma, gammainc, gammaln}\n\
+@end deftypefn")
+{
+  octave_value retval;
+
+  const octave_idx_type nargin = args.length ();
+  if (nargin < 1 || nargin > 2)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  const octave_value oct_z = (nargin == 1) ? args(0) : args(1);
+  const octave_idx_type k = (nargin == 1) ? 0 : args(0).idx_type_value ();
+  if (error_state || k < 0)
+    {
+      error ("psi: K must be a non-negative integer");
+      return retval;
+    }
+  else if (k == 0)
+    {
+#define FLOAT_BRANCH(T, A, M, E) \
+      if (oct_z.is_ ## T ##_type ()) \
+        { \
+          const A ## NDArray z = oct_z.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++ = psi (*zv++); \
+\
+          retval = psi_z; \
+        }
+
+      if (oct_z.is_complex_type ())
+        {
+          FLOAT_BRANCH(double, Complex, complex_, Complex)
+          else FLOAT_BRANCH(single, FloatComplex, float_complex_, FloatComplex)
+          else
+            {
+              error ("psi: Z must be a floating point");
+            }
+        }
+      else
+        {
+          FLOAT_BRANCH(double, , , double)
+          else FLOAT_BRANCH(single, Float, float_, float)
+          else
+            {
+              error ("psi: Z must be a floating point");
+            }
+        }
+
+#undef FLOAT_BRANCH
+    }
+  else
+    {
+      if (! oct_z.is_real_type ())
+        {
+          error ("psi: Z must be real value for polygamma (K > 0)");
+          return retval;
+        }
+
+#define FLOAT_BRANCH(T, A, M, E) \
+      if (oct_z.is_ ## T ##_type ()) \
+        { \
+          const A ## NDArray z = oct_z.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++) \
+            { \
+              if (*zv < 0) \
+                { \
+                  error ("psi: Z must be non-negative for polygamma (K > 0)"); \
+                  return retval; \
+                } \
+              *psi_zv++ = psi (k, *zv++); \
+            } \
+          retval = psi_z; \
+        }
+
+      FLOAT_BRANCH(double, , , double)
+      else FLOAT_BRANCH(single, Float, float_, float)
+      else
+        {
+          error ("psi: Z must be a floating point for polygamma (K > 0)");
+        }
+
+#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)
+%!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)
+%!assert (psi (-1.573498473162390458778286043690434612655040859116846), 0, eps)
+%!assert (psi (-2.610720868444144650001537715718724207951074010873480), 0, eps*10)
+%!assert (psi (-3.635293366436901097839181566946017713948423861193530), 0, eps*10)
+%!assert (psi (-4.653237761743142441714598151148207363719069416133868), 0, eps*100)
+
+## Tests for complex values
+%!shared z
+%! z = [-10:.1:-.1 .1:.1:20]; # drop the 0
+
+## Abramowitz and Stegun, page 259 eq 6.3.10
+%!assert (real (psi (i*z)), real (psi (1 - i*z)))
+
+## Abramowitz and Stegun, page 259 eq 6.3.11
+%!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)
+
+## 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)
+
+## Abramowitz and Stegun, page 260 eq 6.4.5
+%!test
+%! for z = 0:20
+%!   assert (psi (1, z + 0.5), 0.5 * (pi^2) - 4 * sum ((2*(1:z) -1) .^(-2)), eps*10)
+%! endfor
+
+## Abramowitz and Stegun, page 260 eq 6.4.6
+%!test
+%! z = 0.1:0.1:20;
+%! for n = 0:8
+%!   ## our precision goes down really quick when computing n is too high,
+%!   assert (psi (n, z+1), psi (n, z) + ((-1)^n) * factorial (n) * (z.^(-n-1)), 0.1)
+%! endfor
+
+## Test input validation
+%!error psi ()
+%!error psi (1, 2, 3)
+%!error <Z must be> psi ("non numeric")
+%!error <K must be a non-negative integer> psi (-5, 1)
+%!error <Z must be non-negative for polygamma> psi (5, -1)
+%!error <Z must be a floating point> psi (5, uint8 (-1))
+%!error <Z must be real value for polygamma> psi (5, 5i)
+
+*/
diff -r 9866b3202c52 -r b5d2f6954c48 libinterp/corefcn/regexp.cc
--- a/libinterp/corefcn/regexp.cc	Sat May 23 10:47:03 2015 +0200
+++ b/libinterp/corefcn/regexp.cc	Sat May 23 10:19:50 2015 -0400
@@ -77,13 +77,38 @@
               retval[++i] = 'b';
               break;
 
-#if 0
-// FIXME: To be complete, we need to handle \oN, \o{N}.
-//        The PCRE library already handles \N where N
-//        is an octal number.  New code needs to merely
-//        replace \oN or \o{N} with \N.
-            case 'o': // octal number
-#endif
+            case 'o': // octal input
+            {
+              bool bad_esc_seq = (j+1 >= len);
+
+              bool brace = false;
+              if (! bad_esc_seq && s[++j] == '{')
+                {
+                  brace = true;
+                  j++;
+                }
+
+              int tmpi = 0;
+              size_t k;
+              for (k = j; k < std::min (j+3+brace, len); k++)
+                {
+                  int digit = s[k] - '0';
+                  if (digit < 0 || digit > 7)
+                    break;
+                  tmpi <<= 3;
+                  tmpi += digit;
+                }
+              if (bad_esc_seq || (brace && s[k++] != '}'))
+                {
+                  bad_esc_seq = true;
+                  tmpi = 0;
+                  warning ("malformed octal escape sequence '\\o' --\
+ converting to '\\0'");
+                }
+              retval[i] = tmpi;
+              j = k - 1;
+              break;
+            }
 
             default:  // pass escape sequence through
               retval[i] = '\\';
@@ -150,14 +175,75 @@
               retval[i] = '\v';
               break;
 
-#if 0
-// FIXME: to be complete, we need to handle \oN, \o{N}, \xN, and
-// \x{N}.  Hex digits may be upper or lower case.  Brackets are
-// optional, so \x5Bz is the same as \x{5B}z.
+            case 'o': // octal input
+            {
+              bool bad_esc_seq = (j+1 >= len);
+
+              bool brace = false;
+              if (! bad_esc_seq && s[++j] == '{')
+                {
+                  brace = true;
+                  j++;
+                }
+
+              int tmpi = 0;
+              size_t k;
+              for (k = j; k < std::min (j+3+brace, len); k++)
+                {
+                  int digit = s[k] - '0';
+                  if (digit < 0 || digit > 7)
+                    break;
+                  tmpi <<= 3;
+                  tmpi += digit;
+                }
+              if (bad_esc_seq || (brace && s[k++] != '}'))
+                {
+                  warning ("malformed octal escape sequence '\\o' --\
+ converting to '\\0'");
+                  tmpi = 0;
+                }
+              retval[i] = tmpi;
+              j = k - 1;
+              break;
+            }
 
-            case 'o': // octal number
-            case 'x': // hex number
-#endif
+            case 'x': // hex input
+            {
+              bool bad_esc_seq = (j+1 >= len);
+
+              bool brace = false;
+              if (! bad_esc_seq && s[++j] == '{')
+                {
+                  brace = true;
+                  j++;
+                }
+
+              int tmpi = 0;
+              size_t k;
+              for (k = j; k < std::min (j+2+brace, len); k++)
+                {
+                  if (! isxdigit (s[k]))
+                    break;
+
+                  tmpi <<= 4;
+                  int digit = s[k];
+                  if (digit >= 'a')
+                    tmpi += digit - 'a' + 10;
+                  else if (digit >= 'A')
+                    tmpi += digit - 'A' + 10;
+                  else
+                    tmpi += digit - '0';
+                }
+              if (bad_esc_seq || (brace && s[k++] != '}'))
+                {
+                  warning ("malformed hex escape sequence '\\x' --\
+ converting to '\\0'");
+                  tmpi = 0;
+                }
+              retval[i] = tmpi;
+              j = k - 1;
+              break;
+            }
 
             default:  // pass escape sequence through
               retval[i] = '\\';
diff -r 9866b3202c52 -r b5d2f6954c48 libinterp/corefcn/strfns.cc
--- a/libinterp/corefcn/strfns.cc	Sat May 23 10:47:03 2015 +0200
+++ b/libinterp/corefcn/strfns.cc	Sat May 23 10:19:50 2015 -0400
@@ -210,78 +210,70 @@
   octave_value retval;
 
   int nargin = args.length ();
-
-  if (nargin > 0)
-    {
-      int n_elts = 0;
-
-      size_t max_len = 0;
-
-      std::queue<string_vector> args_as_strings;
+  int n_elts = 0;
+  size_t max_len = 0;
+  std::queue<string_vector> args_as_strings;
 
-      for (int i = 0; i < nargin; i++)
-        {
-          string_vector s = args(i).all_strings ();
-
-          if (error_state)
-            {
-              error ("strvcat: unable to convert some args to strings");
-              return retval;
-            }
+  for (int i = 0; i < nargin; i++)
+    {
+      string_vector s = args(i).all_strings ();
 
-          size_t n = s.length ();
+      if (error_state)
+        {
+          error ("strvcat: unable to convert some args to strings");
+          return retval;
+        }
 
-          // do not count empty strings in calculation of number of elements
-          if (n > 0)
+      size_t n = s.length ();
+
+      // do not count empty strings in calculation of number of elements
+      if (n > 0)
+        {
+          for (size_t j = 0; j < n; j++)
             {
-              for (size_t j = 0; j < n; j++)
-                {
-                  if (s[j].length () > 0)
-                    n_elts++;
-                }
+              if (s[j].length () > 0)
+                n_elts++;
             }
-
-          size_t s_max_len = s.max_length ();
-
-          if (s_max_len > max_len)
-            max_len = s_max_len;
-
-          args_as_strings.push (s);
         }
 
-      string_vector result (n_elts);
+      size_t s_max_len = s.max_length ();
+
+      if (s_max_len > max_len)
+        max_len = s_max_len;
+
+      args_as_strings.push (s);
+    }
+
+  string_vector result (n_elts);
 
-      octave_idx_type k = 0;
+  octave_idx_type k = 0;
 
-      for (int i = 0; i < nargin; i++)
+  for (int i = 0; i < nargin; i++)
+    {
+      string_vector s = args_as_strings.front ();
+      args_as_strings.pop ();
+
+      size_t n = s.length ();
+
+      if (n > 0)
         {
-          string_vector s = args_as_strings.front ();
-          args_as_strings.pop ();
-
-          size_t n = s.length ();
-
-          if (n > 0)
+          for (size_t j = 0; j < n; j++)
             {
-              for (size_t j = 0; j < n; j++)
+              std::string t = s[j];
+              if (t.length () > 0)
                 {
-                  std::string t = s[j];
-                  if (t.length () > 0)
-                    {
-                      size_t t_len = t.length ();
+                  size_t t_len = t.length ();
 
-                      if (max_len > t_len)
-                        t += std::string (max_len - t_len, ' ');
+                  if (max_len > t_len)
+                    t += std::string (max_len - t_len, ' ');
 
-                      result[k++] = t;
-                    }
+                  result[k++] = t;
                 }
             }
         }
+    }
 
-      retval = octave_value (result, '\'');
-    }
-  else
-    print_usage ();
+  retval = octave_value (result, '\'');
 
   return retval;
 }
@@ -298,8 +290,7 @@
 %!assert (strvcat ({100,{100, {""}}}), ["d";"d"])
 %!assert (strvcat (["a";"be"], {"c", 100}), ["a";"be";"c";"d"])
 %!assert (strvcat ("a", "bb", "ccc"), ["a  "; "bb "; "ccc"])
-
-%!error strvcat ()
+%!assert (strvcat (), "")
 */
 
 
diff -r 9866b3202c52 -r b5d2f6954c48 libinterp/corefcn/symtab.cc
--- a/libinterp/corefcn/symtab.cc	Sat May 23 10:47:03 2015 +0200
+++ b/libinterp/corefcn/symtab.cc	Sat May 23 10:19:50 2015 -0400
@@ -44,6 +44,7 @@
 #include "pager.h"
 #include "parse.h"
 #include "pt-arg-list.h"
+#include "pt-pr-code.h"
 #include "symtab.h"
 #include "unwind-prot.h"
 #include "utils.h"
@@ -1810,6 +1811,44 @@
   return retval;
 }
 
+DEFUN (__get_cmdline_fcn_txt__, args, ,
+       "-*- texinfo -*-\n\
+@deftypefn  {Built-in Function} {} __get_cmdline_fcn_txt__ (@var{name})\n\
+Undocumented internal function.\n\
+@end deftypefn")
+{
+  octave_value retval;
+
+  if (args.length () == 1)
+    {
+      std::string name = args(0).string_value ();
+
+      if (! error_state)
+        {
+          octave_value ov = symbol_table::find_cmdline_function (name);
+
+          octave_user_function *f = ov.user_function_value ();
+
+          if (f)
+            {
+              std::ostringstream buf;
+
+              tree_print_code tpc (buf);
+
+              f->accept (tpc);
+
+              retval = buf.str ();
+            }
+        }
+      else
+        error ("__get_cmd_line_function_text__: expecting function name");
+    }
+  else
+    print_usage ();
+
+  return retval;
+}
+
 #if 0
 
 // FIXME: should we have functions like this in Octave?
diff -r 9866b3202c52 -r b5d2f6954c48 libinterp/corefcn/symtab.h
--- a/libinterp/corefcn/symtab.h	Sat May 23 10:47:03 2015 +0200
+++ b/libinterp/corefcn/symtab.h	Sat May 23 10:19:50 2015 -0400
@@ -1548,6 +1548,14 @@
            ? p->second.find_user_function () : octave_value ();
   }
 
+  static octave_value find_cmdline_function (const std::string& name)
+  {
+    fcn_table_iterator p = fcn_table.find (name);
+
+    return (p != fcn_table.end ())
+           ? p->second.find_cmdline_function () : octave_value ();
+  }
+
   static void install_cmdline_function (const std::string& name,
                                         const octave_value& fcn)
   {
diff -r 9866b3202c52 -r b5d2f6954c48 libinterp/corefcn/syscalls.cc
--- a/libinterp/corefcn/syscalls.cc	Sat May 23 10:47:03 2015 +0200
+++ b/libinterp/corefcn/syscalls.cc	Sat May 23 10:19:50 2015 -0400
@@ -819,15 +819,45 @@
   return retval;
 }
 
+// FIXME: This routine also exists verbatim in file-io.cc.
+//        Maybe change to be a general utility routine.
+static int
+convert (int x, int ibase, int obase)
+{
+  int retval = 0;
+
+  int tmp = x % obase;
+
+  if (tmp > ibase - 1)
+    ::error ("mkfifo: invalid digit");
+  else
+    {
+      retval = tmp;
+      int mult = ibase;
+      while ((x = (x - tmp) / obase))
+        {
+          tmp = x % obase;
+          if (tmp > ibase - 1)
+            {
+              ::error ("mkfifo: invalid digit");
+              break;
+            }
+          retval += mult * tmp;
+          mult *= ibase;
+        }
+    }
+
+  return retval;
+}
+
 DEFUNX ("mkfifo", Fmkfifo, args, ,
         "-*- texinfo -*-\n\
 @deftypefn  {Built-in Function} {@var{err} =} mkfifo (@var{name}, @var{mode})\n\
 @deftypefnx {Built-in Function} {[@var{err}, @var{msg}] =} mkfifo (@var{name}, @var{mode})\n\
 Create a FIFO special file named @var{name} with file mode @var{mode}.\n\
 \n\
-@var{mode} is interpreted as a decimal number (@emph{not} octal) and is\n\
-subject to umask processing.  The final calculated mode is\n\
-@code{@var{mode} - @var{umask}}.\n\
+@var{mode} is interpreted as an octal number and is subject to umask\n\
+processing.  The final calculated mode is @code{@var{mode} - @var{umask}}.\n\
 \n\
 If successful, @var{err} is 0 and @var{msg} is an empty string.\n\
 Otherwise, @var{err} is nonzero and @var{msg} contains a system-dependent\n\
@@ -835,7 +865,7 @@
 @seealso{pipe, umask}\n\
 @end deftypefn")
 {
-  octave_value_list retval;
+  octave_value_list retval (2);
 
   retval(1) = std::string ();
   retval(0) = -1;
@@ -848,23 +878,28 @@
         {
           std::string name = args(0).string_value ();
 
-          if (args(1).is_scalar_type ())
-            {
-              long mode = args(1).long_value ();
+          int octal_mode = args(1).int_value ();
 
-              if (! error_state)
+          if (! error_state)
+            {
+              if (octal_mode < 0)
+                error ("mkfifo: MODE must be a positive integer value");
+              else
                 {
-                  std::string msg;
+                  int mode = convert (octal_mode, 8, 10);
 
-                  int status = octave_mkfifo (name, mode, msg);
+                  if (! error_state)
+                    {
+                      std::string msg;
 
-                  retval(0) = status;
+                      int status = octave_mkfifo (name, mode, msg);
 
-                  if (status < 0)
-                    retval(1) = msg;
+                      retval(0) = status;
+
+                      if (status < 0)
+                        retval(1) = msg;
+                    }
                 }
-              else
-                error ("mkfifo: invalid MODE");
             }
           else
             error ("mkfifo: MODE must be an integer");
@@ -878,6 +913,19 @@
   return retval;
 }
 
+/*
+
+## Test input validation
+%!error mkfifo ()
+%!error mkfifo ("abc")
+%!error mkfifo ("abc", 777, 123)
+%!error <FILE must be a string> mkfifo (123, 456)
+## FIXME: These tests should work, but lasterr is not being set correctly.
+#%!error <MODE must be an integer> mkfifo ("abc", {456})
+#%!error <MODE must be a positive integer value> mkfifo ("abc", -1)
+
+*/
+
 DEFUNX ("pipe", Fpipe, args, ,
         "-*- texinfo -*-\n\
 @deftypefn {Built-in Function} {[@var{read_fd}, @var{write_fd}, @var{err}, @var{msg}] =} pipe ()\n\
diff -r 9866b3202c52 -r b5d2f6954c48 libinterp/corefcn/utils.cc
--- a/libinterp/corefcn/utils.cc	Sat May 23 10:47:03 2015 +0200
+++ b/libinterp/corefcn/utils.cc	Sat May 23 10:19:50 2015 -0400
@@ -631,11 +631,7 @@
         {
           switch (s[++j])
             {
-            case '0':
-              retval[i] = '\0';
-              break;
-
-            case 'a':
+            case 'a': // alarm
               retval[i] = '\a';
               break;
 
@@ -675,6 +671,58 @@
               retval[i] = '"';
               break;
 
+            case '0':
+            case '1':
+            case '2':
+            case '3':
+            case '4':
+            case '5':
+            case '6':
+            case '7': // octal input
+            {
+              size_t k;
+              int tmpi = s[j] - '0';
+              for (k = j+1; k < std::min (j+3, len); k++)
+                {
+                  int digit = s[k] - '0';
+                  if (digit < 0 || digit > 7)
+                    break;
+                  tmpi <<= 3;
+                  tmpi += digit;
+                }
+              retval[i] = tmpi;
+              j = k - 1;
+              break;
+            }
+
+            case 'x': // hex input
+            {
+              size_t k;
+              int tmpi = 0;
+              for (k = j+1; k < std::min (j+3, len); k++)
+                {
+                  if (! isxdigit (s[k]))
+                    break;
+
+                  tmpi <<= 4;
+                  int digit = s[k];
+                  if (digit >= 'a')
+                    tmpi += digit - 'a' + 10;
+                  else if (digit >= 'A')
+                    tmpi += digit - 'A' + 10;
+                  else
+                    tmpi += digit - '0';
+                }
+
+              if (k == j+1)
+                warning ("malformed hex escape sequence '\\x' --\
+ converting to '\\0'");
+
+              retval[i] = tmpi;
+              j = k - 1;
+              break;
+            }
+
             default:
               warning ("unrecognized escape sequence '\\%c' --\
  converting to '%c'", s[j], s[j]);
@@ -745,9 +793,20 @@
 %!assert (do_string_escapes ('\"double-quoted\"'), "\"double-quoted\"")
 %!assert (do_string_escapes ("\\\"double-quoted\\\""), "\"double-quoted\"")
 
+%!assert (do_string_escapes ('A\4B'), ["A" char(4) "B"])
+%!assert (do_string_escapes ('A\45B'), ["A" char(37) "B"])
+%!assert (do_string_escapes ('A\123B'), ["A" char(83) "B"])
+%!assert (sprintf ('\117\143\164\141\166\145'), "Octave")
+
+%!assert (do_string_escapes ('A\x4G'), ["A" char(4) "G"])
+%!assert (do_string_escapes ('A\x4AG'), ["A" char(74) "G"])
+%!assert (sprintf ('\x4f\x63\x74\x61\x76\x65'), "Octave")
+
 %!error do_string_escapes ()
 %!error do_string_escapes ("foo", "bar")
-%!error do_string_escapes (3)
+%!error <STRING argument> do_string_escapes (3)
+%!warning <malformed hex escape sequence> do_string_escapes ('\xG');
+%!warning <unrecognized escape sequence> do_string_escapes ('\G');
 */
 
 const char *
diff -r 9866b3202c52 -r b5d2f6954c48 libinterp/corefcn/variables.cc
--- a/libinterp/corefcn/variables.cc	Sat May 23 10:47:03 2015 +0200
+++ b/libinterp/corefcn/variables.cc	Sat May 23 10:19:50 2015 -0400
@@ -395,6 +395,14 @@
   bool search_dir = type == "dir";
   bool search_file = type == "file";
   bool search_builtin = type == "builtin";
+  bool search_class = type == "class";
+
+  if (! (search_any || search_var || search_dir || search_file ||
+         search_builtin || search_class))
+    {
+      error ("exist: unrecognized type argument \"%s\"", type.c_str ());
+      return 0;
+    }
 
   if (search_any || search_var)
     {
@@ -682,6 +690,7 @@
 %!warning <"class" type argument is not implemented> exist ("a", "class");
 %!error <TYPE must be a string> exist ("a", 1)
 %!error <NAME must be a string> exist (1)
+%!error <unrecognized type argument "foobar"> exist ("a", "foobar")
 
 */
 
diff -r 9866b3202c52 -r b5d2f6954c48 libinterp/octave-value/ov-classdef.cc
--- a/libinterp/octave-value/ov-classdef.cc	Sat May 23 10:47:03 2015 +0200
+++ b/libinterp/octave-value/ov-classdef.cc	Sat May 23 10:19:50 2015 -0400
@@ -1476,17 +1476,26 @@
 
     case '(':
       {
+        const octave_value_list& ival = idx.front ();
+
         refcount++;
-
         cdef_object this_obj (this);
 
-        Array<cdef_object> arr (dim_vector (1, 1), this_obj);
-
-        cdef_object new_obj = cdef_object (new cdef_object_array (arr));
-
-        new_obj.set_class (get_class ());
-
-        retval = new_obj.subsref (type, idx, nargout, skip, cls, auto_add);
+        if (ival.empty ())
+          {
+            skip++;
+            retval(0) = to_ov (this_obj);
+          }
+        else
+          {
+            Array<cdef_object> arr (dim_vector (1, 1), this_obj);
+
+            cdef_object new_obj = cdef_object (new cdef_object_array (arr));
+
+            new_obj.set_class (get_class ());
+
+            retval = new_obj.subsref (type, idx, nargout, skip, cls, auto_add);
+          }
       }
       break;
 
@@ -1626,16 +1635,17 @@
     case '(':
       {
         const octave_value_list& ival = idx.front ();
-        bool is_scalar = true;
-        Array<idx_vector> iv (dim_vector (1, ival.length ()));
 
         if (ival.empty ())
           {
-            ::error ("can't index %s object(s) with empty parentheses",
-                     class_name ().c_str ());
+            refcount++;
+            retval(0) = to_ov (cdef_object (this));
             break;
           }
 
+        bool is_scalar = true;
+        Array<idx_vector> iv (dim_vector (1, ival.length ()));
+
         for (int i = 0; ! error_state && i < ival.length (); i++)
           {
             iv(i) = ival(i).index_vector ();
diff -r 9866b3202c52 -r b5d2f6954c48 libinterp/parse-tree/pt-pr-code.cc
--- a/libinterp/parse-tree/pt-pr-code.cc	Sat May 23 10:47:03 2015 +0200
+++ b/libinterp/parse-tree/pt-pr-code.cc	Sat May 23 10:19:50 2015 -0400
@@ -329,7 +329,7 @@
 
       cmd_list->accept (*this);
 
-      decrement_indent_level ();
+      // endfunction will decrement the indent level.
     }
 
   visit_octave_user_function_trailer (fcn);
@@ -719,6 +719,9 @@
 void
 tree_print_code::visit_no_op_command (tree_no_op_command& cmd)
 {
+  if (cmd.is_end_of_fcn_or_script ())
+    decrement_indent_level ();
+
   indent ();
 
   os << cmd.original_command ();
@@ -873,13 +876,7 @@
     {
       cmd->accept (*this);
 
-      if (! stmt.print_result ())
-        {
-          os << ";";
-          newline (" ");
-        }
-      else
-        newline ();
+      newline ();
     }
   else
     {
diff -r 9866b3202c52 -r b5d2f6954c48 liboctave/cruft/slatec-fn/dpsifn.f
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/cruft/slatec-fn/dpsifn.f	Sat May 23 10:19:50 2015 -0400
@@ -0,0 +1,368 @@
+*DECK DPSIFN
+      SUBROUTINE DPSIFN (X, N, KODE, M, ANS, NZ, IERR)
+C***BEGIN PROLOGUE  DPSIFN
+C***PURPOSE  Compute derivatives of the Psi function.
+C***LIBRARY   SLATEC
+C***CATEGORY  C7C
+C***TYPE      DOUBLE PRECISION (PSIFN-S, DPSIFN-D)
+C***KEYWORDS  DERIVATIVES OF THE GAMMA FUNCTION, POLYGAMMA FUNCTION,
+C             PSI FUNCTION
+C***AUTHOR  Amos, D. E., (SNLA)
+C***DESCRIPTION
+C
+C         The following definitions are used in DPSIFN:
+C
+C      Definition 1
+C         PSI(X) = d/dx (ln(GAMMA(X)), the first derivative of
+C                  the log GAMMA function.
+C      Definition 2
+C                     K   K
+C         PSI(K,X) = d /dx (PSI(X)), the K-th derivative of PSI(X).
+C   ___________________________________________________________________
+C      DPSIFN computes a sequence of SCALED derivatives of
+C      the PSI function; i.e. for fixed X and M it computes
+C      the M-member sequence
+C
+C                    ((-1)**(K+1)/GAMMA(K+1))*PSI(K,X)
+C                       for K = N,...,N+M-1
+C
+C      where PSI(K,X) is as defined above.   For KODE=1, DPSIFN returns
+C      the scaled derivatives as described.  KODE=2 is operative only
+C      when K=0 and in that case DPSIFN returns -PSI(X) + LN(X).  That
+C      is, the logarithmic behavior for large X is removed when KODE=2
+C      and K=0.  When sums or differences of PSI functions are computed
+C      the logarithmic terms can be combined analytically and computed
+C      separately to help retain significant digits.
+C
+C         Note that CALL DPSIFN(X,0,1,1,ANS) results in
+C                   ANS = -PSI(X)
+C
+C     Input      X is DOUBLE PRECISION
+C           X      - Argument, X .gt. 0.0D0
+C           N      - First member of the sequence, 0 .le. N .le. 100
+C                    N=0 gives ANS(1) = -PSI(X)       for KODE=1
+C                                       -PSI(X)+LN(X) for KODE=2
+C           KODE   - Selection parameter
+C                    KODE=1 returns scaled derivatives of the PSI
+C                    function.
+C                    KODE=2 returns scaled derivatives of the PSI
+C                    function EXCEPT when N=0. In this case,
+C                    ANS(1) = -PSI(X) + LN(X) is returned.
+C           M      - Number of members of the sequence, M.ge.1
+C
+C    Output     ANS is DOUBLE PRECISION
+C           ANS    - A vector of length at least M whose first M
+C                    components contain the sequence of derivatives
+C                    scaled according to KODE.
+C           NZ     - Underflow flag
+C                    NZ.eq.0, A normal return
+C                    NZ.ne.0, Underflow, last NZ components of ANS are
+C                             set to zero, ANS(M-K+1)=0.0, K=1,...,NZ
+C           IERR   - Error flag
+C                    IERR=0, A normal return, computation completed
+C                    IERR=1, Input error,     no computation
+C                    IERR=2, Overflow,        X too small or N+M-1 too
+C                            large or both
+C                    IERR=3, Error,           N too large. Dimensioned
+C                            array TRMR(NMAX) is not large enough for N
+C
+C         The nominal computational accuracy is the maximum of unit
+C         roundoff (=D1MACH(4)) and 1.0D-18 since critical constants
+C         are given to only 18 digits.
+C
+C         PSIFN is the single precision version of DPSIFN.
+C
+C *Long Description:
+C
+C         The basic method of evaluation is the asymptotic expansion
+C         for large X.ge.XMIN followed by backward recursion on a two
+C         term recursion relation
+C
+C                  W(X+1) + X**(-N-1) = W(X).
+C
+C         This is supplemented by a series
+C
+C                  SUM( (X+K)**(-N-1) , K=0,1,2,... )
+C
+C         which converges rapidly for large N. Both XMIN and the
+C         number of terms of the series are calculated from the unit
+C         roundoff of the machine environment.
+C
+C***REFERENCES  Handbook of Mathematical Functions, National Bureau
+C                 of Standards Applied Mathematics Series 55, edited
+C                 by M. Abramowitz and I. A. Stegun, equations 6.3.5,
+C                 6.3.18, 6.4.6, 6.4.9 and 6.4.10, pp.258-260, 1964.
+C               D. E. Amos, A portable Fortran subroutine for
+C                 derivatives of the Psi function, Algorithm 610, ACM
+C                 Transactions on Mathematical Software 9, 4 (1983),
+C                 pp. 494-502.
+C***ROUTINES CALLED  D1MACH, I1MACH
+C***REVISION HISTORY  (YYMMDD)
+C   820601  DATE WRITTEN
+C   890531  Changed all specific intrinsics to generic.  (WRB)
+C   890911  Removed unnecessary intrinsics.  (WRB)
+C   891006  Cosmetic changes to prologue.  (WRB)
+C   891006  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   920501  Reformatted the REFERENCES section.  (WRB)
+C***END PROLOGUE  DPSIFN
+      INTEGER I, IERR, J, K, KODE, M, MM, MX, N, NMAX, NN, NP, NX, NZ,
+     *  FN
+      INTEGER I1MACH
+      DOUBLE PRECISION ANS, ARG, B, DEN, ELIM, EPS, FLN,
+     * FX, RLN, RXSQ, R1M4, R1M5, S, SLOPE, T, TA, TK, TOL, TOLS, TRM,
+     * TRMR, TSS, TST, TT, T1, T2, WDTOL, X, XDMLN, XDMY, XINC, XLN,
+     * XM, XMIN, XQ, YINT
+      DOUBLE PRECISION D1MACH
+      DIMENSION B(22), TRM(22), TRMR(100), ANS(*)
+      SAVE NMAX, B
+      DATA NMAX /100/
+C-----------------------------------------------------------------------
+C             BERNOULLI NUMBERS
+C-----------------------------------------------------------------------
+      DATA B(1), B(2), B(3), B(4), B(5), B(6), B(7), B(8), B(9), B(10),
+     * B(11), B(12), B(13), B(14), B(15), B(16), B(17), B(18), B(19),
+     * B(20), B(21), B(22) /1.00000000000000000D+00,
+     * -5.00000000000000000D-01,1.66666666666666667D-01,
+     * -3.33333333333333333D-02,2.38095238095238095D-02,
+     * -3.33333333333333333D-02,7.57575757575757576D-02,
+     * -2.53113553113553114D-01,1.16666666666666667D+00,
+     * -7.09215686274509804D+00,5.49711779448621554D+01,
+     * -5.29124242424242424D+02,6.19212318840579710D+03,
+     * -8.65802531135531136D+04,1.42551716666666667D+06,
+     * -2.72982310678160920D+07,6.01580873900642368D+08,
+     * -1.51163157670921569D+10,4.29614643061166667D+11,
+     * -1.37116552050883328D+13,4.88332318973593167D+14,
+     * -1.92965793419400681D+16/
+C
+C***FIRST EXECUTABLE STATEMENT  DPSIFN
+      IERR = 0
+      NZ=0
+      IF (X.LE.0.0D0) IERR=1
+      IF (N.LT.0) IERR=1
+      IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
+      IF (M.LT.1) IERR=1
+      IF (IERR.NE.0) RETURN
+      MM=M
+      NX = MIN(-I1MACH(15),I1MACH(16))
+      R1M5 = D1MACH(5)
+      R1M4 = D1MACH(4)*0.5D0
+      WDTOL = MAX(R1M4,0.5D-18)
+C-----------------------------------------------------------------------
+C     ELIM = APPROXIMATE EXPONENTIAL OVER AND UNDERFLOW LIMIT
+C-----------------------------------------------------------------------
+      ELIM = 2.302D0*(NX*R1M5-3.0D0)
+      XLN = LOG(X)
+   41 CONTINUE
+      NN = N + MM - 1
+      FN = NN
+      T = (FN+1)*XLN
+C-----------------------------------------------------------------------
+C     OVERFLOW AND UNDERFLOW TEST FOR SMALL AND LARGE X
+C-----------------------------------------------------------------------
+      IF (ABS(T).GT.ELIM) GO TO 290
+      IF (X.LT.WDTOL) GO TO 260
+C-----------------------------------------------------------------------
+C     COMPUTE XMIN AND THE NUMBER OF TERMS OF THE SERIES, FLN+1
+C-----------------------------------------------------------------------
+      RLN = R1M5*I1MACH(14)
+      RLN = MIN(RLN,18.06D0)
+      FLN = MAX(RLN,3.0D0) - 3.0D0
+      YINT = 3.50D0 + 0.40D0*FLN
+      SLOPE = 0.21D0 + FLN*(0.0006038D0*FLN+0.008677D0)
+      XM = YINT + SLOPE*FN
+      MX = INT(XM) + 1
+      XMIN = MX
+      IF (N.EQ.0) GO TO 50
+      XM = -2.302D0*RLN - MIN(0.0D0,XLN)
+      ARG = XM/N
+      ARG = MIN(0.0D0,ARG)
+      EPS = EXP(ARG)
+      XM = 1.0D0 - EPS
+      IF (ABS(ARG).LT.1.0D-3) XM = -ARG
+      FLN = X*XM/EPS
+      XM = XMIN - X
+      IF (XM.GT.7.0D0 .AND. FLN.LT.15.0D0) GO TO 200
+   50 CONTINUE
+      XDMY = X
+      XDMLN = XLN
+      XINC = 0.0D0
+      IF (X.GE.XMIN) GO TO 60
+      NX = INT(X)
+      XINC = XMIN - NX
+      XDMY = X + XINC
+      XDMLN = LOG(XDMY)
+   60 CONTINUE
+C-----------------------------------------------------------------------
+C     GENERATE W(N+MM-1,X) BY THE ASYMPTOTIC EXPANSION
+C-----------------------------------------------------------------------
+      T = FN*XDMLN
+      T1 = XDMLN + XDMLN
+      T2 = T + XDMLN
+      TK = MAX(ABS(T),ABS(T1),ABS(T2))
+      IF (TK.GT.ELIM) GO TO 380
+      TSS = EXP(-T)
+      TT = 0.5D0/XDMY
+      T1 = TT
+      TST = WDTOL*TT
+      IF (NN.NE.0) T1 = TT + 1.0D0/FN
+      RXSQ = 1.0D0/(XDMY*XDMY)
+      TA = 0.5D0*RXSQ
+      T = (FN+1)*TA
+      S = T*B(3)
+      IF (ABS(S).LT.TST) GO TO 80
+      TK = 2.0D0
+      DO 70 K=4,22
+        T = T*((TK+FN+1)/(TK+1.0D0))*((TK+FN)/(TK+2.0D0))*RXSQ
+        TRM(K) = T*B(K)
+        IF (ABS(TRM(K)).LT.TST) GO TO 80
+        S = S + TRM(K)
+        TK = TK + 2.0D0
+   70 CONTINUE
+   80 CONTINUE
+      S = (S+T1)*TSS
+      IF (XINC.EQ.0.0D0) GO TO 100
+C-----------------------------------------------------------------------
+C     BACKWARD RECUR FROM XDMY TO X
+C-----------------------------------------------------------------------
+      NX = INT(XINC)
+      NP = NN + 1
+      IF (NX.GT.NMAX) GO TO 390
+      IF (NN.EQ.0) GO TO 160
+      XM = XINC - 1.0D0
+      FX = X + XM
+C-----------------------------------------------------------------------
+C     THIS LOOP SHOULD NOT BE CHANGED. FX IS ACCURATE WHEN X IS SMALL
+C-----------------------------------------------------------------------
+      DO 90 I=1,NX
+        TRMR(I) = FX**(-NP)
+        S = S + TRMR(I)
+        XM = XM - 1.0D0
+        FX = X + XM
+   90 CONTINUE
+  100 CONTINUE
+      ANS(MM) = S
+      IF (FN.EQ.0) GO TO 180
+C-----------------------------------------------------------------------
+C     GENERATE LOWER DERIVATIVES, J.LT.N+MM-1
+C-----------------------------------------------------------------------
+      IF (MM.EQ.1) RETURN
+      DO 150 J=2,MM
+        FN = FN - 1
+        TSS = TSS*XDMY
+        T1 = TT
+        IF (FN.NE.0) T1 = TT + 1.0D0/FN
+        T = (FN+1)*TA
+        S = T*B(3)
+        IF (ABS(S).LT.TST) GO TO 120
+        TK = 4 + FN
+        DO 110 K=4,22
+          TRM(K) = TRM(K)*(FN+1)/TK
+          IF (ABS(TRM(K)).LT.TST) GO TO 120
+          S = S + TRM(K)
+          TK = TK + 2.0D0
+  110   CONTINUE
+  120   CONTINUE
+        S = (S+T1)*TSS
+        IF (XINC.EQ.0.0D0) GO TO 140
+        IF (FN.EQ.0) GO TO 160
+        XM = XINC - 1.0D0
+        FX = X + XM
+        DO 130 I=1,NX
+          TRMR(I) = TRMR(I)*FX
+          S = S + TRMR(I)
+          XM = XM - 1.0D0
+          FX = X + XM
+  130   CONTINUE
+  140   CONTINUE
+        MX = MM - J + 1
+        ANS(MX) = S
+        IF (FN.EQ.0) GO TO 180
+  150 CONTINUE
+      RETURN
+C-----------------------------------------------------------------------
+C     RECURSION FOR N = 0
+C-----------------------------------------------------------------------
+  160 CONTINUE
+      DO 170 I=1,NX
+        S = S + 1.0D0/(X+NX-I)
+  170 CONTINUE
+  180 CONTINUE
+      IF (KODE.EQ.2) GO TO 190
+      ANS(1) = S - XDMLN
+      RETURN
+  190 CONTINUE
+      IF (XDMY.EQ.X) RETURN
+      XQ = XDMY/X
+      ANS(1) = S - LOG(XQ)
+      RETURN
+C-----------------------------------------------------------------------
+C     COMPUTE BY SERIES (X+K)**(-(N+1)) , K=0,1,2,...
+C-----------------------------------------------------------------------
+  200 CONTINUE
+      NN = INT(FLN) + 1
+      NP = N + 1
+      T1 = (N+1)*XLN
+      T = EXP(-T1)
+      S = T
+      DEN = X
+      DO 210 I=1,NN
+        DEN = DEN + 1.0D0
+        TRM(I) = DEN**(-NP)
+        S = S + TRM(I)
+  210 CONTINUE
+      ANS(1) = S
+      IF (N.NE.0) GO TO 220
+      IF (KODE.EQ.2) ANS(1) = S + XLN
+  220 CONTINUE
+      IF (MM.EQ.1) RETURN
+C-----------------------------------------------------------------------
+C     GENERATE HIGHER DERIVATIVES, J.GT.N
+C-----------------------------------------------------------------------
+      TOL = WDTOL/5.0D0
+      DO 250 J=2,MM
+        T = T/X
+        S = T
+        TOLS = T*TOL
+        DEN = X
+        DO 230 I=1,NN
+          DEN = DEN + 1.0D0
+          TRM(I) = TRM(I)/DEN
+          S = S + TRM(I)
+          IF (TRM(I).LT.TOLS) GO TO 240
+  230   CONTINUE
+  240   CONTINUE
+        ANS(J) = S
+  250 CONTINUE
+      RETURN
+C-----------------------------------------------------------------------
+C     SMALL X.LT.UNIT ROUND OFF
+C-----------------------------------------------------------------------
+  260 CONTINUE
+      ANS(1) = X**(-N-1)
+      IF (MM.EQ.1) GO TO 280
+      K = 1
+      DO 270 I=2,MM
+        ANS(K+1) = ANS(K)/X
+        K = K + 1
+  270 CONTINUE
+  280 CONTINUE
+      IF (N.NE.0) RETURN
+      IF (KODE.EQ.2) ANS(1) = ANS(1) + XLN
+      RETURN
+  290 CONTINUE
+      IF (T.GT.0.0D0) GO TO 380
+      NZ=0
+      IERR=2
+      RETURN
+  380 CONTINUE
+      NZ=NZ+1
+      ANS(MM)=0.0D0
+      MM=MM-1
+      IF (MM.EQ.0) RETURN
+      GO TO 41
+  390 CONTINUE
+      NZ=0
+      IERR=3
+      RETURN
+      END
diff -r 9866b3202c52 -r b5d2f6954c48 liboctave/cruft/slatec-fn/module.mk
--- a/liboctave/cruft/slatec-fn/module.mk	Sat May 23 10:47:03 2015 +0200
+++ b/liboctave/cruft/slatec-fn/module.mk	Sat May 23 10:19:50 2015 -0400
@@ -34,6 +34,7 @@
   cruft/slatec-fn/dlnrel.f \
   cruft/slatec-fn/dpchim.f \
   cruft/slatec-fn/dpchst.f \
+  cruft/slatec-fn/dpsifn.f \
   cruft/slatec-fn/erf.f \
   cruft/slatec-fn/gami.f \
   cruft/slatec-fn/gamit.f \
@@ -44,6 +45,7 @@
   cruft/slatec-fn/inits.f \
   cruft/slatec-fn/pchim.f \
   cruft/slatec-fn/pchst.f \
+  cruft/slatec-fn/psifn.f \
   cruft/slatec-fn/r9lgmc.f \
   cruft/slatec-fn/r9lgit.f \
   cruft/slatec-fn/r9gmit.f \
diff -r 9866b3202c52 -r b5d2f6954c48 liboctave/cruft/slatec-fn/psifn.f
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/cruft/slatec-fn/psifn.f	Sat May 23 10:19:50 2015 -0400
@@ -0,0 +1,368 @@
+*DECK PSIFN
+      SUBROUTINE PSIFN (X, N, KODE, M, ANS, NZ, IERR)
+C***BEGIN PROLOGUE  PSIFN
+C***PURPOSE  Compute derivatives of the Psi function.
+C***LIBRARY   SLATEC
+C***CATEGORY  C7C
+C***TYPE      SINGLE PRECISION (PSIFN-S, DPSIFN-D)
+C***KEYWORDS  DERIVATIVES OF THE GAMMA FUNCTION, POLYGAMMA FUNCTION,
+C             PSI FUNCTION
+C***AUTHOR  Amos, D. E., (SNLA)
+C***DESCRIPTION
+C
+C         The following definitions are used in PSIFN:
+C
+C      Definition 1
+C         PSI(X) = d/dx (ln(GAMMA(X)), the first derivative of
+C                  the LOG GAMMA function.
+C      Definition 2
+C                     K   K
+C         PSI(K,X) = d /dx (PSI(X)), the K-th derivative of PSI(X).
+C   ___________________________________________________________________
+C       PSIFN computes a sequence of SCALED derivatives of
+C       the PSI function; i.e. for fixed X and M it computes
+C       the M-member sequence
+C
+C                  ((-1)**(K+1)/GAMMA(K+1))*PSI(K,X)
+C                    for K = N,...,N+M-1
+C
+C       where PSI(K,X) is as defined above.   For KODE=1, PSIFN returns
+C       the scaled derivatives as described.  KODE=2 is operative only
+C       when K=0 and in that case PSIFN returns -PSI(X) + LN(X).  That
+C       is, the logarithmic behavior for large X is removed when KODE=1
+C       and K=0.  When sums or differences of PSI functions are computed
+C       the logarithmic terms can be combined analytically and computed
+C       separately to help retain significant digits.
+C
+C         Note that CALL PSIFN(X,0,1,1,ANS) results in
+C                   ANS = -PSI(X)
+C
+C     Input
+C           X      - Argument, X .gt. 0.0E0
+C           N      - First member of the sequence, 0 .le. N .le. 100
+C                    N=0 gives ANS(1) = -PSI(X)       for KODE=1
+C                                       -PSI(X)+LN(X) for KODE=2
+C           KODE   - Selection parameter
+C                    KODE=1 returns scaled derivatives of the PSI
+C                    function.
+C                    KODE=2 returns scaled derivatives of the PSI
+C                    function EXCEPT when N=0. In this case,
+C                    ANS(1) = -PSI(X) + LN(X) is returned.
+C           M      - Number of members of the sequence, M .ge. 1
+C
+C    Output
+C           ANS    - A vector of length at least M whose first M
+C                    components contain the sequence of derivatives
+C                    scaled according to KODE.
+C           NZ     - Underflow flag
+C                    NZ.eq.0, A normal return
+C                    NZ.ne.0, Underflow, last NZ components of ANS are
+C                             set to zero, ANS(M-K+1)=0.0, K=1,...,NZ
+C           IERR   - Error flag
+C                    IERR=0, A normal return, computation completed
+C                    IERR=1, Input error,     no computation
+C                    IERR=2, Overflow,        X too small or N+M-1 too
+C                            large or both
+C                    IERR=3, Error,           N too large. Dimensioned
+C                            array TRMR(NMAX) is not large enough for N
+C
+C         The nominal computational accuracy is the maximum of unit
+C         roundoff (=R1MACH(4)) and 1.0E-18 since critical constants
+C         are given to only 18 digits.
+C
+C         DPSIFN is the Double Precision version of PSIFN.
+C
+C *Long Description:
+C
+C         The basic method of evaluation is the asymptotic expansion
+C         for large X.ge.XMIN followed by backward recursion on a two
+C         term recursion relation
+C
+C                  W(X+1) + X**(-N-1) = W(X).
+C
+C         This is supplemented by a series
+C
+C                  SUM( (X+K)**(-N-1) , K=0,1,2,... )
+C
+C         which converges rapidly for large N. Both XMIN and the
+C         number of terms of the series are calculated from the unit
+C         roundoff of the machine environment.
+C
+C***REFERENCES  Handbook of Mathematical Functions, National Bureau
+C                 of Standards Applied Mathematics Series 55, edited
+C                 by M. Abramowitz and I. A. Stegun, equations 6.3.5,
+C                 6.3.18, 6.4.6, 6.4.9 and 6.4.10, pp.258-260, 1964.
+C               D. E. Amos, A portable Fortran subroutine for
+C                 derivatives of the Psi function, Algorithm 610, ACM
+C                 Transactions on Mathematical Software 9, 4 (1983),
+C                 pp. 494-502.
+C***ROUTINES CALLED  I1MACH, R1MACH
+C***REVISION HISTORY  (YYMMDD)
+C   820601  DATE WRITTEN
+C   890531  Changed all specific intrinsics to generic.  (WRB)
+C   890531  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   920501  Reformatted the REFERENCES section.  (WRB)
+C***END PROLOGUE  PSIFN
+      INTEGER I, IERR, J, K, KODE, M, MM, MX, N, NMAX, NN, NP, NX, NZ
+      INTEGER I1MACH
+      REAL ANS, ARG, B, DEN, ELIM, EPS, FLN, FN, FNP, FNS, FX, RLN,
+     * RXSQ, R1M4, R1M5, S, SLOPE, T, TA, TK, TOL, TOLS, TRM, TRMR,
+     * TSS, TST, TT, T1, T2, WDTOL, X, XDMLN, XDMY, XINC, XLN, XM,
+     * XMIN, XQ, YINT
+      REAL R1MACH
+      DIMENSION B(22), TRM(22), TRMR(100), ANS(*)
+      SAVE NMAX, B
+      DATA NMAX /100/
+C-----------------------------------------------------------------------
+C             BERNOULLI NUMBERS
+C-----------------------------------------------------------------------
+      DATA B(1), B(2), B(3), B(4), B(5), B(6), B(7), B(8), B(9), B(10),
+     * B(11), B(12), B(13), B(14), B(15), B(16), B(17), B(18), B(19),
+     * B(20), B(21), B(22) /1.00000000000000000E+00,
+     * -5.00000000000000000E-01,1.66666666666666667E-01,
+     * -3.33333333333333333E-02,2.38095238095238095E-02,
+     * -3.33333333333333333E-02,7.57575757575757576E-02,
+     * -2.53113553113553114E-01,1.16666666666666667E+00,
+     * -7.09215686274509804E+00,5.49711779448621554E+01,
+     * -5.29124242424242424E+02,6.19212318840579710E+03,
+     * -8.65802531135531136E+04,1.42551716666666667E+06,
+     * -2.72982310678160920E+07,6.01580873900642368E+08,
+     * -1.51163157670921569E+10,4.29614643061166667E+11,
+     * -1.37116552050883328E+13,4.88332318973593167E+14,
+     * -1.92965793419400681E+16/
+C
+C***FIRST EXECUTABLE STATEMENT  PSIFN
+      IERR = 0
+      NZ=0
+      IF (X.LE.0.0E0) IERR=1
+      IF (N.LT.0) IERR=1
+      IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
+      IF (M.LT.1) IERR=1
+      IF (IERR.NE.0) RETURN
+      MM=M
+      NX = MIN(-I1MACH(12),I1MACH(13))
+      R1M5 = R1MACH(5)
+      R1M4 = R1MACH(4)*0.5E0
+      WDTOL = MAX(R1M4,0.5E-18)
+C-----------------------------------------------------------------------
+C     ELIM = APPROXIMATE EXPONENTIAL OVER AND UNDERFLOW LIMIT
+C-----------------------------------------------------------------------
+      ELIM = 2.302E0*(NX*R1M5-3.0E0)
+      XLN = LOG(X)
+   41 CONTINUE
+      NN = N + MM - 1
+      FN = NN
+      FNP = FN + 1.0E0
+      T = FNP*XLN
+C-----------------------------------------------------------------------
+C     OVERFLOW AND UNDERFLOW TEST FOR SMALL AND LARGE X
+C-----------------------------------------------------------------------
+      IF (ABS(T).GT.ELIM) GO TO 290
+      IF (X.LT.WDTOL) GO TO 260
+C-----------------------------------------------------------------------
+C     COMPUTE XMIN AND THE NUMBER OF TERMS OF THE SERIES, FLN+1
+C-----------------------------------------------------------------------
+      RLN = R1M5*I1MACH(11)
+      RLN = MIN(RLN,18.06E0)
+      FLN = MAX(RLN,3.0E0) - 3.0E0
+      YINT = 3.50E0 + 0.40E0*FLN
+      SLOPE = 0.21E0 + FLN*(0.0006038E0*FLN+0.008677E0)
+      XM = YINT + SLOPE*FN
+      MX = INT(XM) + 1
+      XMIN = MX
+      IF (N.EQ.0) GO TO 50
+      XM = -2.302E0*RLN - MIN(0.0E0,XLN)
+      FNS = N
+      ARG = XM/FNS
+      ARG = MIN(0.0E0,ARG)
+      EPS = EXP(ARG)
+      XM = 1.0E0 - EPS
+      IF (ABS(ARG).LT.1.0E-3) XM = -ARG
+      FLN = X*XM/EPS
+      XM = XMIN - X
+      IF (XM.GT.7.0E0 .AND. FLN.LT.15.0E0) GO TO 200
+   50 CONTINUE
+      XDMY = X
+      XDMLN = XLN
+      XINC = 0.0E0
+      IF (X.GE.XMIN) GO TO 60
+      NX = INT(X)
+      XINC = XMIN - NX
+      XDMY = X + XINC
+      XDMLN = LOG(XDMY)
+   60 CONTINUE
+C-----------------------------------------------------------------------
+C     GENERATE W(N+MM-1,X) BY THE ASYMPTOTIC EXPANSION
+C-----------------------------------------------------------------------
+      T = FN*XDMLN
+      T1 = XDMLN + XDMLN
+      T2 = T + XDMLN
+      TK = MAX(ABS(T),ABS(T1),ABS(T2))
+      IF (TK.GT.ELIM) GO TO 380
+      TSS = EXP(-T)
+      TT = 0.5E0/XDMY
+      T1 = TT
+      TST = WDTOL*TT
+      IF (NN.NE.0) T1 = TT + 1.0E0/FN
+      RXSQ = 1.0E0/(XDMY*XDMY)
+      TA = 0.5E0*RXSQ
+      T = FNP*TA
+      S = T*B(3)
+      IF (ABS(S).LT.TST) GO TO 80
+      TK = 2.0E0
+      DO 70 K=4,22
+        T = T*((TK+FN+1.0E0)/(TK+1.0E0))*((TK+FN)/(TK+2.0E0))*RXSQ
+        TRM(K) = T*B(K)
+        IF (ABS(TRM(K)).LT.TST) GO TO 80
+        S = S + TRM(K)
+        TK = TK + 2.0E0
+   70 CONTINUE
+   80 CONTINUE
+      S = (S+T1)*TSS
+      IF (XINC.EQ.0.0E0) GO TO 100
+C-----------------------------------------------------------------------
+C     BACKWARD RECUR FROM XDMY TO X
+C-----------------------------------------------------------------------
+      NX = INT(XINC)
+      NP = NN + 1
+      IF (NX.GT.NMAX) GO TO 390
+      IF (NN.EQ.0) GO TO 160
+      XM = XINC - 1.0E0
+      FX = X + XM
+C-----------------------------------------------------------------------
+C     THIS LOOP SHOULD NOT BE CHANGED. FX IS ACCURATE WHEN X IS SMALL
+C-----------------------------------------------------------------------
+      DO 90 I=1,NX
+        TRMR(I) = FX**(-NP)
+        S = S + TRMR(I)
+        XM = XM - 1.0E0
+        FX = X + XM
+   90 CONTINUE
+  100 CONTINUE
+      ANS(MM) = S
+      IF (FN.EQ.0.0E0) GO TO 180
+C-----------------------------------------------------------------------
+C     GENERATE LOWER DERIVATIVES, J.LT.N+MM-1
+C-----------------------------------------------------------------------
+      IF (MM.EQ.1) RETURN
+      DO 150 J=2,MM
+        FNP = FN
+        FN = FN - 1.0E0
+        TSS = TSS*XDMY
+        T1 = TT
+        IF (FN.NE.0.0E0) T1 = TT + 1.0E0/FN
+        T = FNP*TA
+        S = T*B(3)
+        IF (ABS(S).LT.TST) GO TO 120
+        TK = 3.0E0 + FNP
+        DO 110 K=4,22
+          TRM(K) = TRM(K)*FNP/TK
+          IF (ABS(TRM(K)).LT.TST) GO TO 120
+          S = S + TRM(K)
+          TK = TK + 2.0E0
+  110   CONTINUE
+  120   CONTINUE
+        S = (S+T1)*TSS
+        IF (XINC.EQ.0.0E0) GO TO 140
+        IF (FN.EQ.0.0E0) GO TO 160
+        XM = XINC - 1.0E0
+        FX = X + XM
+        DO 130 I=1,NX
+          TRMR(I) = TRMR(I)*FX
+          S = S + TRMR(I)
+          XM = XM - 1.0E0
+          FX = X + XM
+  130   CONTINUE
+  140   CONTINUE
+        MX = MM - J + 1
+        ANS(MX) = S
+        IF (FN.EQ.0.0E0) GO TO 180
+  150 CONTINUE
+      RETURN
+C-----------------------------------------------------------------------
+C     RECURSION FOR N = 0
+C-----------------------------------------------------------------------
+  160 CONTINUE
+      DO 170 I=1,NX
+        S = S + 1.0E0/(X+NX-I)
+  170 CONTINUE
+  180 CONTINUE
+      IF (KODE.EQ.2) GO TO 190
+      ANS(1) = S - XDMLN
+      RETURN
+  190 CONTINUE
+      IF (XDMY.EQ.X) RETURN
+      XQ = XDMY/X
+      ANS(1) = S - LOG(XQ)
+      RETURN
+C-----------------------------------------------------------------------
+C     COMPUTE BY SERIES (X+K)**(-(N+1)) , K=0,1,2,...
+C-----------------------------------------------------------------------
+  200 CONTINUE
+      NN = INT(FLN) + 1
+      NP = N + 1
+      T1 = (FNS+1.0E0)*XLN
+      T = EXP(-T1)
+      S = T
+      DEN = X
+      DO 210 I=1,NN
+        DEN = DEN + 1.0E0
+        TRM(I) = DEN**(-NP)
+        S = S + TRM(I)
+  210 CONTINUE
+      ANS(1) = S
+      IF (N.NE.0) GO TO 220
+      IF (KODE.EQ.2) ANS(1) = S + XLN
+  220 CONTINUE
+      IF (MM.EQ.1) RETURN
+C-----------------------------------------------------------------------
+C     GENERATE HIGHER DERIVATIVES, J.GT.N
+C-----------------------------------------------------------------------
+      TOL = WDTOL/5.0E0
+      DO 250 J=2,MM
+        T = T/X
+        S = T
+        TOLS = T*TOL
+        DEN = X
+        DO 230 I=1,NN
+          DEN = DEN + 1.0E0
+          TRM(I) = TRM(I)/DEN
+          S = S + TRM(I)
+          IF (TRM(I).LT.TOLS) GO TO 240
+  230   CONTINUE
+  240   CONTINUE
+        ANS(J) = S
+  250 CONTINUE
+      RETURN
+C-----------------------------------------------------------------------
+C     SMALL X.LT.UNIT ROUND OFF
+C-----------------------------------------------------------------------
+  260 CONTINUE
+      ANS(1) = X**(-N-1)
+      IF (MM.EQ.1) GO TO 280
+      K = 1
+      DO 270 I=2,MM
+        ANS(K+1) = ANS(K)/X
+        K = K + 1
+  270 CONTINUE
+  280 CONTINUE
+      IF (N.NE.0) RETURN
+      IF (KODE.EQ.2) ANS(1) = ANS(1) + XLN
+      RETURN
+  290 CONTINUE
+      IF (T.GT.0.0E0) GO TO 380
+      NZ=0
+      IERR=2
+      RETURN
+  380 CONTINUE
+      NZ=NZ+1
+      ANS(MM)=0.0E0
+      MM=MM-1
+      IF(MM.EQ.0) RETURN
+      GO TO 41
+  390 CONTINUE
+      IERR=3
+      NZ=0
+      RETURN
+      END
diff -r 9866b3202c52 -r b5d2f6954c48 liboctave/numeric/lo-specfun.cc
--- a/liboctave/numeric/lo-specfun.cc	Sat May 23 10:47:03 2015 +0200
+++ b/liboctave/numeric/lo-specfun.cc	Sat May 23 10:19:50 2015 -0400
@@ -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"
 
@@ -183,6 +186,16 @@
 
   F77_RET_T
   F77_FUNC (algams, ALGAMS) (const float&, float&, float&);
+
+  F77_RET_T
+  F77_FUNC (psifn, PSIFN) (const float*, const octave_idx_type&,
+                           const octave_idx_type&, const octave_idx_type&,
+                           float*, octave_idx_type*, octave_idx_type*);
+
+  F77_RET_T
+  F77_FUNC (dpsifn, DPSIFN) (const double*, const octave_idx_type&,
+                             const octave_idx_type&, const octave_idx_type&,
+                             double*, octave_idx_type*, octave_idx_type*);
 }
 
 #if !defined (HAVE_ACOSH)
@@ -3724,3 +3737,185 @@
       dn = Complex (dd*cc1*dd1/ddd, -m*ss*cc*ss1/ddd);
     }
 }
+
+static const double pi = 3.14159265358979323846;
+
+template<class T>
+static 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;
+  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;
+      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
+      for (octave_idx_type k = z; k > 0; k--)
+        p += 1.0 / (2 * k - 1);
+
+      p = - euler_mascheroni - 2 * log (2) + 2 * (p);
+    }
+  else
+    {
+      // adapted from XLiFE++ gammaFunctions
+
+      T zc = z;
+      // Use formula for derivative of LogGamma(z)
+      if (z < 10)
+        {
+          const signed char n = 10 - z;
+          for (signed char k = n - 1; k >= 0; k--)
+            p -= 1.0 / (k + z);
+          zc += n;
+        }
+      p += Lanczos_approximation_psi (zc);
+    }
+
+  return p;
+}
+
+// explicit instantiations
+template double psi<double> (const double& z);
+template float  psi<float> (const float& z);
+
+template<class T>
+std::complex<T>
+psi (const std::complex<T>& z)
+{
+  // adapted from XLiFE++ gammaFunctions
+
+  typedef typename std::complex<T>::value_type P;
+
+  P z_r  = z.real ();
+  P z_ra = z_r;
+
+  std::complex<T> dgam (0.0, 0.0);
+  if (z.imag () == 0)
+    dgam = std::complex<T> (psi (z_r), 0.0);
+  else if (z_r < 0)
+    dgam = psi (P (1.0) - z)- (P (pi) / tan (P (pi) * z));
+  else
+    {
+      // Use formula for derivative of LogGamma(z)
+      std::complex<T> z_m = z;
+      if (z_ra < 8)
+        {
+          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
+      dgam += Lanczos_approximation_psi (z_m);
+    }
+  return dgam;
+}
+
+// explicit instantiations
+template Complex psi<double> (const Complex& z);
+template FloatComplex psi<float> (const FloatComplex& z);
+
+
+template<typename T>
+static inline void
+fortran_psifn (const T z, const octave_idx_type n, T* ans,
+               octave_idx_type* ierr);
+
+template<>
+inline void
+fortran_psifn<double> (const double z, const octave_idx_type n,
+                       double* ans, octave_idx_type* ierr)
+{
+  octave_idx_type flag = 0;
+  F77_XFCN (dpsifn, DPSIFN, (&z, n, 1, 1, ans, &flag, ierr));
+}
+
+template<>
+inline void
+fortran_psifn<float> (const float z, const octave_idx_type n,
+                      float* ans, octave_idx_type* ierr)
+{
+  octave_idx_type flag = 0;
+  F77_XFCN (psifn, PSIFN, (&z, n, 1, 1, ans, &flag, ierr));
+}
+
+template<class T>
+T
+psi (const octave_idx_type n, const T z)
+{
+  T ans;
+  octave_idx_type ierr = 0;
+  fortran_psifn<T> (z, n, &ans, &ierr);
+  if (ierr == 0)
+    {
+      // Remember that psifn and dpsifn return scales values
+      // When n is 1: do nothing since ((-1)**(n+1)/gamma(n+1)) == 1
+      // When n is 0: change sign since ((-1)**(n+1)/gamma(n+1)) == -1
+      if (n > 1)
+        // FIXME xgamma here is a killer for our precision since it grows
+        //       way too fast
+        ans = ans / (pow (-1.0, n + 1) / xgamma (double (n+1)));
+      else if (n == 0)
+        ans = -ans;
+    }
+  else if (ierr == 2)
+    ans = - octave_Inf;
+  else // we probably never get here
+    ans = octave_NaN;
+
+  return ans;
+}
+
+// explicit instantiations
+template double psi<double> (const octave_idx_type n, const double z);
+template float  psi<float>  (const octave_idx_type n, const float z);
diff -r 9866b3202c52 -r b5d2f6954c48 liboctave/numeric/lo-specfun.h
--- a/liboctave/numeric/lo-specfun.h	Sat May 23 10:47:03 2015 +0200
+++ b/liboctave/numeric/lo-specfun.h	Sat May 23 10:19:50 2015 -0400
@@ -663,4 +663,24 @@
 ellipj (const Complex& u, double m, Complex& sn, Complex& cn, Complex& dn,
         double& err);
 
+//! Digamma function.
+//!
+//! Only defined for double and float.
+template<class T>
+extern OCTAVE_API T psi (const T& z);
+
+//! Digamma function for complex input.
+//!
+//! Only defined for double and float.
+template<class T>
+extern OCTAVE_API std::complex<T> psi (const std::complex<T>& z);
+
+//! Polygamma function.
+//!
+//! Only defined for double and float.
+//! @param n must be non-negative.  If zero, the digamma function is computed.
+//! @param z must be real and non-negative.
+template<class T>
+extern OCTAVE_API T psi (const octave_idx_type n, const T z);
+
 #endif
diff -r 9866b3202c52 -r b5d2f6954c48 m4/acinclude.m4
--- a/m4/acinclude.m4	Sat May 23 10:47:03 2015 +0200
+++ b/m4/acinclude.m4	Sat May 23 10:19:50 2015 -0400
@@ -291,7 +291,7 @@
     [AC_LANG_PUSH(C++)
     ac_octave_save_CPPFLAGS="$CPPFLAGS"
     CPPFLAGS="$QT_CPPFLAGS $CPPFLAGS"
-    AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
+    AC_PREPROC_IFELSE([AC_LANG_PROGRAM([[
         #include <Qsci/qsciglobal.h>
         ]], [[
         #if QSCINTILLA_VERSION < 0x020600
@@ -417,7 +417,7 @@
     [AC_LANG_PUSH(C++)
     ac_octave_save_CPPFLAGS="$CPPFLAGS"
     CPPFLAGS="$QT_CPPFLAGS $CPPFLAGS"
-    AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
+    AC_PREPROC_IFELSE([AC_LANG_PROGRAM([[
         #include <Qt/qglobal.h>
         ]], [[
         #if QT_VERSION < 0x040700
diff -r 9866b3202c52 -r b5d2f6954c48 scripts/general/num2str.m
--- a/scripts/general/num2str.m	Sat May 23 10:47:03 2015 +0200
+++ b/scripts/general/num2str.m	Sat May 23 10:19:50 2015 -0400
@@ -115,7 +115,8 @@
         fmt = "%3d";
       endif
     endif
-    fmt = [deblank(repmat(fmt, 1, columns(x))), "\n"];
+    fmt = do_string_escapes (fmt);  # required now that '\n' is interpreted.
+    fmt = [deblank(repmat (fmt, 1, columns (x))), "\n"];
     nd = ndims (x);
     tmp = sprintf (fmt, permute (x, [2, 1, 3:nd]));
     retval = strtrim (char (ostrsplit (tmp(1:end-1), "\n")));
@@ -204,6 +205,9 @@
 %!xtest
 %! assert (num2str (1e23), "100000000000000000000000");
 
+## Test for bug #44864, extra rows generated from newlines in format
+%!assert (rows (num2str (magic (3), '%3d %3d %3d\n')), 3)
+
 %!error num2str ()
 %!error num2str (1, 2, 3)
 %!error <X must be a numeric> num2str ({1})
diff -r 9866b3202c52 -r b5d2f6954c48 scripts/geometry/griddata.m
--- a/scripts/geometry/griddata.m	Sat May 23 10:47:03 2015 +0200
+++ b/scripts/geometry/griddata.m	Sat May 23 10:19:50 2015 -0400
@@ -133,14 +133,12 @@
     error ("griddata: unknown interpolation METHOD");
   endif
 
-  if (nargout == 3)
+  if (nargout > 1)
     rx = xi;
     ry = yi;
     rz = zi;
-  elseif (nargout == 1)
+  else
     rx = zi;
-  elseif (nargout == 0)
-    mesh (xi, yi, zi);
   endif
 
 endfunction
@@ -153,7 +151,8 @@
 %! y = 2*rand (size (x)) - 1;
 %! z = sin (2*(x.^2 + y.^2));
 %! [xx,yy] = meshgrid (linspace (-1,1,32));
-%! griddata (x,y,z,xx,yy);
+%! zz = griddata (x,y,z,xx,yy);
+%! mesh (xx, yy, zz);
 %! title ("nonuniform grid sampled at 100 points");
 
 %!demo
@@ -163,7 +162,8 @@
 %! y = 2*rand (size (x)) - 1;
 %! z = sin (2*(x.^2 + y.^2));
 %! [xx,yy] = meshgrid (linspace (-1,1,32));
-%! griddata (x,y,z,xx,yy);
+%! zz = griddata (x,y,z,xx,yy);
+%! mesh (xx, yy, zz);
 %! title ("nonuniform grid sampled at 1000 points");
 
 %!demo
@@ -173,7 +173,8 @@
 %! y = 2*rand (size (x)) - 1;
 %! z = sin (2*(x.^2 + y.^2));
 %! [xx,yy] = meshgrid (linspace (-1,1,32));
-%! griddata (x,y,z,xx,yy,"nearest");
+%! zz = griddata (x,y,z,xx,yy,"nearest");
+%! mesh (xx, yy, zz);
 %! title ("nonuniform grid sampled at 1000 points with nearest neighbor");
 
 %!testif HAVE_QHULL
diff -r 9866b3202c52 -r b5d2f6954c48 scripts/gui/private/__fltk_file_filter__.m
--- a/scripts/gui/private/__fltk_file_filter__.m	Sat May 23 10:47:03 2015 +0200
+++ b/scripts/gui/private/__fltk_file_filter__.m	Sat May 23 10:19:50 2015 -0400
@@ -37,7 +37,7 @@
     curr_ext = ostrsplit (curr_ext, ";");
 
     if (length (curr_ext) > 1)
-      curr_ext = regexprep (curr_ext, '\*\.', ',');
+      curr_ext = strrep (curr_ext, '*.', ',');
       curr_ext = strcat (curr_ext{:})(2 : end);
       curr_ext = strcat ("*.{", curr_ext, "}");
     else
@@ -48,8 +48,8 @@
 
     if (c == 2)
       curr_desc = file_filter{idx, 2};
-      curr_desc = regexprep (curr_desc, '\(', '<');
-      curr_desc = regexprep (curr_desc, '\)', '>');
+      curr_desc = strrep (curr_desc, '(', '<');
+      curr_desc = strrep (curr_desc, ')', '>');
     endif
 
     if (length (fltk_str) > 0)
diff -r 9866b3202c52 -r b5d2f6954c48 scripts/help/type.m
--- a/scripts/help/type.m	Sat May 23 10:47:03 2015 +0200
+++ b/scripts/help/type.m	Sat May 23 10:19:50 2015 -0400
@@ -97,6 +97,19 @@
       txt = sprintf ("%s is a dynamically-linked function", name);
     elseif (e == 5)
       txt = sprintf ("%s is a built-in function", name);
+    elseif (e == 103)
+      contents = __get_cmdline_fcn_txt__ (name);
+      if (isempty (contents))
+        txt = sprintf ("%s is a command-line function with no definition",
+                       name);
+      else
+        if (quiet)
+          txt = contents;
+        else
+          txt = sprintf ("%s is the command-line function:\n\n%s",
+                         name, contents);
+        endif
+      endif
     elseif (any (strcmp (__operators__ (), name)))
       txt = sprintf ("%s is an operator", name);
     elseif (any (strcmp (__keywords__ (), name)))
diff -r 9866b3202c52 -r b5d2f6954c48 scripts/image/imformats.m
--- a/scripts/image/imformats.m	Sat May 23 10:47:03 2015 +0200
+++ b/scripts/image/imformats.m	Sat May 23 10:19:50 2015 -0400
@@ -79,7 +79,7 @@
   persistent formats = default_formats ();
 
   if (nargin == 0 && nargout == 0)
-    error ("imformats: pretty print not yet implemented.");
+    pretty_print_formats (formats);
   elseif (nargin >= 1)
     if (isstruct (arg1))
       arrayfun (@is_valid_format, arg1);
@@ -281,6 +281,40 @@
   end_try_catch
 endfunction
 
+function pretty_print_formats (formats)
+  ## define header names (none should be shorter than 3 characters)
+  headers = {"Extension", "isa", "Info", "Read", "Write", "Alpha", "Description"};
+  cols_length = cellfun (@numel, headers);
+
+  ## Adjust the maximal length of the extensions column
+  extensions = cellfun (@strjoin, {formats.ext}, {", "},
+                        "UniformOutput", false);
+  cols_length(1) = max (max (cellfun (@numel, extensions)), cols_length(1));
+  headers{1} = postpad (headers{1}, cols_length(1), " ");
+
+  ## Print the headers
+  disp (strjoin (headers, " | "));
+  under_headers = cellfun (@(x) repmat ("-", 1, numel (x)), headers,
+                           "UniformOutput", false);
+  disp (strjoin (under_headers, "-+-"));
+
+  template = strjoin (arrayfun (@(x) sprintf ("%%-%is", x), cols_length,
+                                "UniformOutput", false), " | ");
+
+  ## Print the function handle for this things won't be a pretty table.  So
+  ## instead we replace them with "yes" or "no", based on the support it has.
+  yes_no_cols = cat (2, {formats.isa}(:), {formats.info}(:), {formats.read}(:),
+                     {formats.write}(:), {formats.alpha}(:));
+  empty = cellfun (@isempty, yes_no_cols);
+  yes_no_cols(empty) = "no";
+  yes_no_cols(! empty) = "yes";
+
+  descriptions = {formats.description};
+  table = cat (2, extensions(:), yes_no_cols, descriptions(:));
+  printf ([template "\n"], table'{:});
+
+endfunction
+
 ## When imread or imfinfo are called, the file must exist or the
 ## function defined by imformats will never be called.  Because
 ## of this, we must create a file for the tests to work.
diff -r 9866b3202c52 -r b5d2f6954c48 scripts/io/strread.m
--- a/scripts/io/strread.m	Sat May 23 10:47:03 2015 +0200
+++ b/scripts/io/strread.m	Sat May 23 10:19:50 2015 -0400
@@ -352,7 +352,7 @@
   ## Remove comments in str
   if (comment_flag)
     ## Expand 'eol_char' here, after option processing which may have set value
-    comment_end = regexprep (comment_end, "eol_char", eol_char);
+    comment_end = strrep (comment_end, "eol_char", eol_char);
     cstart = strfind (str, comment_start);
     cstop  = strfind (str, comment_end);
     ## Treat end of string as additional comment stop
diff -r 9866b3202c52 -r b5d2f6954c48 scripts/plot/draw/private/__contour__.m
--- a/scripts/plot/draw/private/__contour__.m	Sat May 23 10:47:03 2015 +0200
+++ b/scripts/plot/draw/private/__contour__.m	Sat May 23 10:19:50 2015 -0400
@@ -28,40 +28,42 @@
 
   linespec.color = "auto";
   linespec.linestyle = "-";
-  for i = 3:2:nargin
-    arg = varargin{i};
-    if (ischar (arg) || iscellstr (arg))
+  opts = {};
+  i = 3;
+  while (i <= length (varargin))
+    if (ischar (varargin{i}) || iscellstr (varargin{i}))
+      arg = varargin{i};
+      if (i < length (varargin))
+        if (strcmpi (arg, "fill"))
+          filled = varargin{i+1};
+          varargin(i:i+1) = [];
+          continue;
+        elseif (strcmpi (arg, "linecolor"))
+          linespec.color = varargin{i+1};
+          varargin(i:i+1) = [];
+          continue;
+        endif
+      endif
+
       [lspec, valid] = __pltopt__ ("__contour__", arg, false);
       if (valid)
-        have_line_spec = true;
         varargin(i) = [];
-        linespec = lspec;
-        if (isempty (linespec.color))
-          linespec.color = "auto";
+        if (! isempty (lspec.color))
+          linespec.color = lspec.color;
+        endif
+        if (! isempty (lspec.linestyle))
+          linespec.linestyle = lspec.linestyle;
         endif
-        if (isempty (linespec.linestyle))
-          linespec.linestyle = "-";
+      else  # unrecognized option, pass unmodified in opts cell array
+        if (i < length (varargin))
+          opts(end+(1:2)) = varargin(i:i+1);
+          varargin(i:i+1) = [];
+        else
+          error ("__contour__: Uneven number of PROP/VAL pairs");
         endif
-        break;
       endif
-    endif
-  endfor
 
-  opts = {};
-  i = 3;
-  while (i < length (varargin))
-    if (ischar (varargin{i}))
-      if (strcmpi (varargin{i}, "fill"))
-        filled = varargin{i+1};
-        varargin(i:i+1) = [];
-      elseif (strcmpi (varargin{i}, "linecolor"))
-        linespec.color = varargin{i+1};
-        varargin(i:i+1) = [];
-      else
-        opts(end+(1:2)) = varargin(i:i+1);
-        varargin(i:i+1) = [];
-      endif
-    else
+    else  # skip numeric arguments
       i++;
     endif
   endwhile
diff -r 9866b3202c52 -r b5d2f6954c48 scripts/plot/util/ginput.m
--- a/scripts/plot/util/ginput.m	Sat May 23 10:47:03 2015 +0200
+++ b/scripts/plot/util/ginput.m	Sat May 23 10:19:50 2015 -0400
@@ -69,12 +69,15 @@
   ginput_accumulator (0, 0, 0, 0);  # initialize accumulator
 
   orig_windowbuttondownfcn = get (fig, "windowbuttondownfcn");
-  orig_ginput_keypressfcn = get (fig, "keypressfcn");
+  orig_keypressfcn = get (fig, "keypressfcn");
+  orig_closerequestfcn = get (fig, "closerequestfcn");
 
   unwind_protect
 
     set (fig, "windowbuttondownfcn", @ginput_windowbuttondownfcn);
     set (fig, "keypressfcn", @ginput_keypressfcn);
+    set (fig, "closerequestfcn", {@ginput_closerequestfcn,
+                                  orig_closerequestfcn});
 
     do
       if (strcmp (toolkit, "fltk"))
@@ -95,8 +98,11 @@
     endif
 
   unwind_protect_cleanup
-    set (fig, "windowbuttondownfcn", orig_windowbuttondownfcn);
-    set (fig, "keypressfcn", orig_ginput_keypressfcn);
+    if (isfigure (fig))
+      ## Only execute if window still exists
+      set (fig, "windowbuttondownfcn", orig_windowbuttondownfcn);
+      set (fig, "keypressfcn", orig_keypressfcn);
+    endif
   end_unwind_protect
 
   varargout = {x, y, button};
@@ -123,12 +129,12 @@
 
 endfunction
 
-function ginput_windowbuttondownfcn (src, button)
+function ginput_windowbuttondownfcn (~, button)
   point = get (gca (), "currentpoint");
   ginput_accumulator (1, point(1,1), point(1,2), button);
 endfunction
 
-function ginput_keypressfcn (src, evt)
+function ginput_keypressfcn (~, evt)
   point = get (gca (), "currentpoint");
   if (strcmp (evt.Key, "return"))
     ## Enter key stops ginput.
@@ -141,6 +147,11 @@
   endif
 endfunction
 
+function ginput_closerequestfcn (hfig, ~, orig_closerequestfcn)
+  ginput_accumulator (2, NaN, NaN, NaN);  # Stop ginput
+  feval (orig_closerequestfcn);           # Close window with original fcn
+endfunction
+
 
 ## Remove from test statistics.  No real tests possible.
 %!test
diff -r 9866b3202c52 -r b5d2f6954c48 scripts/plot/util/private/__gnuplot_get_var__.m
--- a/scripts/plot/util/private/__gnuplot_get_var__.m	Sat May 23 10:47:03 2015 +0200
+++ b/scripts/plot/util/private/__gnuplot_get_var__.m	Sat May 23 10:19:50 2015 -0400
@@ -52,8 +52,7 @@
   if (use_mkfifo)
     gpin_name = tempname ();
 
-    ## Mode: 0600 == 6*8*8
-    [err, msg] = mkfifo (gpin_name, 6*8*8);
+    [err, msg] = mkfifo (gpin_name, 600);
 
     if (err)
       error ("__gnuplot_get_var__: Can not make FIFO (%s)", msg);
diff -r 9866b3202c52 -r b5d2f6954c48 scripts/plot/util/private/__gnuplot_ginput__.m
--- a/scripts/plot/util/private/__gnuplot_ginput__.m	Sat May 23 10:47:03 2015 +0200
+++ b/scripts/plot/util/private/__gnuplot_ginput__.m	Sat May 23 10:19:50 2015 -0400
@@ -63,8 +63,7 @@
   if (use_mkfifo)
     gpin_name = tempname ();
 
-    ##Mode: 6*8*8 ==  0600
-    [err, msg] = mkfifo (gpin_name, 6*8*8);
+    [err, msg] = mkfifo (gpin_name, 600);
 
     if (err)
       error ("ginput: Can not open fifo (%s)", msg);
diff -r 9866b3202c52 -r b5d2f6954c48 scripts/plot/util/private/__go_draw_axes__.m
--- a/scripts/plot/util/private/__go_draw_axes__.m	Sat May 23 10:47:03 2015 +0200
+++ b/scripts/plot/util/private/__go_draw_axes__.m	Sat May 23 10:19:50 2015 -0400
@@ -2345,12 +2345,12 @@
         ## FIXME: The symbol font doesn't seem to support bold or italic
         ##if (bld)
         ##  if (it)
-        ##    g = regexprep (g, '/Symbol', '/Symbol-bolditalic');
+        ##    g = strrep (g, '/Symbol', '/Symbol-bolditalic');
         ##  else
-        ##    g = regexprep (g, '/Symbol', '/Symbol-bold');
+        ##    g = strrep (g, '/Symbol', '/Symbol-bold');
         ##  endif
         ##elseif (it)
-        ##  g = regexprep (g, '/Symbol', '/Symbol-italic');
+        ##  g = strrep (g, '/Symbol', '/Symbol-italic');
         ##endif
         str = [str(1:s(i) - 1) g str(e(i) + 1:end)];
       elseif (strncmp (f, "rm", 2))
@@ -2406,12 +2406,12 @@
             ## FIXME: The symbol font doesn't seem to support bold or italic
             ##if (bld)
             ##  if (it)
-            ##    g = regexprep (g, '/Symbol', '/Symbol-bolditalic');
+            ##    g = strrep (g, '/Symbol', '/Symbol-bolditalic');
             ##  else
-            ##    g = regexprep (g, '/Symbol', '/Symbol-bold');
+            ##    g = strrep (g, '/Symbol', '/Symbol-bold');
             ##  endif
             ##elseif (it)
-            ##  g = regexprep (g, '/Symbol', '/Symbol-italic');
+            ##  g = strrep (g, '/Symbol', '/Symbol-italic');
             ##endif
             str = [str(1:s(i) - 1) g str(s(i) + length (flds{j}) + 1:end)];
             break;
diff -r 9866b3202c52 -r b5d2f6954c48 scripts/strings/cstrcat.m
--- a/scripts/strings/cstrcat.m	Sat May 23 10:47:03 2015 +0200
+++ b/scripts/strings/cstrcat.m	Sat May 23 10:19:50 2015 -0400
@@ -45,14 +45,16 @@
 
 function st = cstrcat (varargin)
 
-  if (nargin < 1)
-    print_usage ();
-  elseif (! iscellstr (varargin))
+  if (nargin == 0)
+    ## Special because if varargin is empty, iscellstr still returns
+    ## true but then "[varargin{:}]" would be of class double.
+    st = "";
+  elseif (iscellstr (varargin))
+    st = [varargin{:}];
+  else
     error ("cstrcat: expecting arguments to character strings");
   endif
 
-  st = [varargin{:}];
-
 endfunction
 
 
@@ -65,7 +67,8 @@
 %!assert (cstrcat ("foo", "bar"), "foobar")
 %!assert (cstrcat (["a"; "bb"], ["foo"; "bar"]), ["a foo"; "bbbar"])
 
+%!assert (cstrcat (), "")
+
 ## Test input validation
-%!error cstrcat ()
 %!error cstrcat (1, 2)
 
diff -r 9866b3202c52 -r b5d2f6954c48 scripts/strings/strcat.m
--- a/scripts/strings/strcat.m	Sat May 23 10:47:03 2015 +0200
+++ b/scripts/strings/strcat.m	Sat May 23 10:19:50 2015 -0400
@@ -82,10 +82,8 @@
 function st = strcat (varargin)
 
   if (nargin == 0)
-    print_usage ();
-  endif
-
-  if (nargin == 1)
+    st = "";
+  elseif (nargin == 1)
     st = varargin{1};
   else
     ## Convert to cells of strings
@@ -149,5 +147,5 @@
 %!assert (strcat (1, 2), strcat (char (1), char (2)))
 %!assert (strcat ("", 2), strcat ([], char (2)))
 
-%!error strcat ()
+%!assert (strcat (), "")
 
diff -r 9866b3202c52 -r b5d2f6954c48 scripts/testfun/private/html_compare_plot_demos.m
--- a/scripts/testfun/private/html_compare_plot_demos.m	Sat May 23 10:47:03 2015 +0200
+++ b/scripts/testfun/private/html_compare_plot_demos.m	Sat May 23 10:19:50 2015 -0400
@@ -113,7 +113,7 @@
           fprintf (fid, "<img src='%s' style='width: %dpx;'>", ...
                         ffn, in.column_width);
         else
-          err_fn = regexprep(ffn, ".png", ".err");
+          err_fn = strrep (ffn, ".png", ".err");
           if (! exist (err_fn, "file"))
             warning("File %s doesn't exist...", err_fn);
           else
diff -r 9866b3202c52 -r b5d2f6954c48 scripts/testfun/test.m
--- a/scripts/testfun/test.m	Sat May 23 10:47:03 2015 +0200
+++ b/scripts/testfun/test.m	Sat May 23 10:19:50 2015 -0400
@@ -609,14 +609,13 @@
           fflush (__fid);
         endif
         fprintf (__fid, "%s\n", __msg);
-        fflush (__fid);
         ## Show the variable context.
         if (! strcmp (__type, "error") && ! strcmp (__type, "testif")
             && ! all (__shared == " "))
           fputs (__fid, "shared variables ");
           eval (sprintf ("fdisp(__fid,var2struct(%s));", __shared));
-          fflush (__fid);
         endif
+        fflush (__fid);
       endif
       if (! __success && ! __isxtest)
         __all_success = false;
diff -r 9866b3202c52 -r b5d2f6954c48 test/Makefile.am
--- a/test/Makefile.am	Sat May 23 10:47:03 2015 +0200
+++ b/test/Makefile.am	Sat May 23 10:19:50 2015 -0400
@@ -58,6 +58,7 @@
 include bug-36025/module.mk
 include bug-38236/module.mk
 include bug-38691/module.mk
+include bug-44940/module.mk
 include classdef/module.mk
 include classes/module.mk
 include class-concat/module.mk
diff -r 9866b3202c52 -r b5d2f6954c48 test/bug-44940/bug-44940.tst
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/bug-44940/bug-44940.tst	Sat May 23 10:19:50 2015 -0400
@@ -0,0 +1,11 @@
+%!test
+%! a = class_bug44940 ();
+%! b = a;
+%! c = a ();
+%! a.child = 100;
+%! assert (a.child, b.child)
+%! assert (a.child, c.child)
+%! c.child = 500;
+%! assert (a.child, b.child)
+%! assert (a.child, c.child)
+
diff -r 9866b3202c52 -r b5d2f6954c48 test/bug-44940/class_bug44940.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/bug-44940/class_bug44940.m	Sat May 23 10:19:50 2015 -0400
@@ -0,0 +1,5 @@
+classdef class_bug44940 < handle
+  properties
+    child
+  endproperties
+endclassdef
diff -r 9866b3202c52 -r b5d2f6954c48 test/bug-44940/module.mk
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/bug-44940/module.mk	Sat May 23 10:19:50 2015 -0400
@@ -0,0 +1,5 @@
+bug_44940_FCN_FILES = \
+  bug-44940/bug-44940.tst \
+  bug-44940/class_bug44940.m
+
+FCN_FILES += $(bug_44940_FCN_FILES)