# HG changeset patch
# User jwe
# Date 886712388 0
# Node ID 3d5aefef14e2e14444d382efba7b1d4af5fa4fc0
# Parent  0d640dc625c77f89dc7fe032314dbc0fb9332731
[project @ 1998-02-05 20:58:44 by jwe]

diff -r 0d640dc625c7 -r 3d5aefef14e2 liboctave/ChangeLog
--- a/liboctave/ChangeLog	Thu Feb 05 08:45:07 1998 +0000
+++ b/liboctave/ChangeLog	Thu Feb 05 20:59:48 1998 +0000
@@ -1,5 +1,9 @@
 Thu Feb  5 02:12:38 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
+	* lo-specfun.cc: Don't include dbleBessel.h.
+
+	* Makefile.in (INCLUDES): Delete oct-math.h from the list.
+
 	* dir-ops.h (dir_entry::operator bool ()): Return bool, not void*.
 	* file-stat.h (file_stat::operator bool ()): Likewise.
 	* idx-vector.h (idx_vector::operator bool ()): Likewise.
diff -r 0d640dc625c7 -r 3d5aefef14e2 liboctave/Makefile.in
--- a/liboctave/Makefile.in	Thu Feb 05 08:45:07 1998 +0000
+++ b/liboctave/Makefile.in	Thu Feb 05 20:59:48 1998 +0000
@@ -47,7 +47,7 @@
 	dir-ops.h file-ops.h file-stat.h getopt.h glob-match.h \
 	idx-vector.h lo-ieee.h lo-mappers.h lo-specfun.h lo-sysdep.h \
 	lo-utils.h mach-info.h oct-alloc.h oct-cmplx.h oct-env.h \
-	oct-math.h oct-group.h oct-passwd.h oct-syscalls.h pathlen.h \
+	oct-group.h oct-passwd.h oct-syscalls.h pathlen.h \
 	pathsearch.h prog-args.h statdefs.h str-vec.h sun-utils.h \
 	sysdir.h syswait.h \
 	$(MATRIX_INC) \
diff -r 0d640dc625c7 -r 3d5aefef14e2 liboctave/dbleBessel.cc
--- a/liboctave/dbleBessel.cc	Thu Feb 05 08:45:07 1998 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,221 +0,0 @@
-/*
-
-Copyright (C) 1996 John W. Eaton
-
-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 2, 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, write to the Free
-Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
-
-*/
-
-#if defined (__GNUG__)
-#pragma implementation
-#endif
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "Range.h"
-#include "dColVector.h"
-#include "dMatrix.h"
-#include "dbleBessel.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-#include "mx-inlines.cc"
-
-extern "C"
-{
-  int F77_FCN (rjbesl, RJBESL) (const double&, const double&,
-				const int&, double*, int&);
-
-  int F77_FCN (rybesl, RYBESL) (const double&, const double&,
-				const int&, double*, int&);
-
-  int F77_FCN (ribesl, RIBESL) (const double&, const double&,
-				const int&, const int&, double*, int&);
-
-  int F77_FCN (rkbesl, RKBESL) (const double&, const double&,
-				const int&, const int&, double*, int&);
-}
-
-int
-F77_FCN (ribesl, RIBESL) (const double& x, const double& alpha,
-			  const int& nb, double *result, int& ncalc)
-{
-  int retval = 0;
-  F77_FCN (ribesl, RIBESL) (x, alpha, nb, 1, result, ncalc);
-  return retval;
-}
-
-int
-F77_FCN (rkbesl, RKBESL) (const double& x, const double& alpha,
-			  const int& nb, double *result, int& ncalc)
-{
-  int retval = 0;
-  F77_FCN (rkbesl, RKBESL) (x, alpha, nb, 1, result, ncalc);
-  return retval;
-}
-
-typedef int (*fptr) (const double&, const double&, const int&, double*, int&);
-
-static Matrix
-do_bessel (fptr f, const char *fn, double alpha, const Matrix& x)
-{
-  Matrix retval;
-
-  if (alpha >= 0.0)
-    {
-      int nb = 1;
-
-      if (alpha >= 1.0)
-	{
-	  double integral;
-	  alpha = modf (alpha, &integral);
-	  nb = static_cast<int> (integral) + 1;
-	}
-
-      Array<double> tmp (nb);
-
-      double *ptmp = tmp.fortran_vec ();
-
-      int x_nr = x.rows ();
-      int x_nc = x.cols ();
-
-      retval.resize (x_nr, x_nc);
-
-      int ncalc;
-
-      for (int j = 0; j < x_nc; j++)
-	{
-	  for (int i = 0; i < x_nr; i++)
-	    {
-	      f (x(i,j), alpha, nb, ptmp, ncalc);
-
-	      // XXX FIXME XXX -- check ncalc.
-
-	      retval(i,j) = tmp(nb-1);
-	    }
-	}
-    }
-  else
-    (*current_liboctave_error_handler)
-      ("%s: alpha must be greater than zero", fn);
-
-  return retval;
-}
-
-static Matrix
-do_bessel (fptr f, const char *fn, const Range& alpha_range,
-	   const ColumnVector& x)
-{
-  Matrix retval;
-
-  double alpha_base = alpha_range.base ();
-
-  if (alpha_base >= 0.0)
-    {
-      int nb_orig = alpha_range.nelem ();
-      int nb = nb_orig;
-      int offset = 0;
-
-      if (alpha_base >= 1.0)
-	{
-	  double integral;
-	  alpha_base = modf (alpha_base, &integral);
-	  offset = static_cast<int> (integral);
-	  nb += offset;
-	}
-
-      Array<double> tmp (nb);
-
-      double *ptmp = tmp.fortran_vec ();
-
-      int x_len = x.length ();
-
-      retval.resize (x_len, nb_orig);
-
-      int ncalc;
-
-      for (int i = 0; i < x_len; i++)
-	{
-	  f (x(i), alpha_base, nb, ptmp, ncalc);
-
-	  // XXX FIXME XXX -- check ncalc.
-
-	  for (int j = 0; j < nb_orig; j++)
-	    retval(i,j) = tmp(offset+j);
-	}
-    }
-  else
-    (*current_liboctave_error_handler)
-      ("%s: alpha must be greater than zero", fn);
-
-  return retval;
-}
-
-Matrix
-besselj (double alpha, const Matrix& x)
-{
-  return do_bessel (F77_FCN (rjbesl, RJBESL), "besselj", alpha, x);
-}
-
-Matrix
-bessely (double alpha, const Matrix& x)
-{
-  return do_bessel (F77_FCN (rybesl, RYBESL), "bessely", alpha, x);
-}
-
-Matrix
-besseli (double alpha, const Matrix& x)
-{
-  return do_bessel (F77_FCN (ribesl, RIBESL), "besseli", alpha, x);
-}
-
-Matrix
-besselk (double alpha, const Matrix& x)
-{
-  return do_bessel (F77_FCN (rkbesl, RKBESL), "besselk", alpha, x);
-}
-
-Matrix
-besselj (const Range& alpha, const ColumnVector& x)
-{
-  return do_bessel (F77_FCN (rjbesl, RJBESL), "besselj", alpha, x);
-}
-
-Matrix
-bessely (const Range& alpha, const ColumnVector& x)
-{
-  return do_bessel (F77_FCN (rybesl, RYBESL), "bessely", alpha, x);
-}
-
-Matrix
-besseli (const Range& alpha, const ColumnVector& x)
-{
-  return do_bessel (F77_FCN (ribesl, RIBESL), "besseli", alpha, x);
-}
-
-Matrix
-besselk (const Range& alpha, const ColumnVector& x)
-{
-  return do_bessel (F77_FCN (rkbesl, RKBESL), "besselk", alpha, x);
-}
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
diff -r 0d640dc625c7 -r 3d5aefef14e2 liboctave/dbleBessel.h
--- a/liboctave/dbleBessel.h	Thu Feb 05 08:45:07 1998 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,50 +0,0 @@
-/*
-
-Copyright (C) 1996 John W. Eaton
-
-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 2, 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, write to the Free
-Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
-
-*/
-
-#if !defined (octave_BESSEL_h)
-#define octave_BESSEL_h 1
-
-#if defined (__GNUG__)
-#pragma interface
-#endif
-
-class ColumnVector;
-class Matrix;
-class Range;
-
-extern Matrix besselj (double alpha, const Matrix& x);
-extern Matrix bessely (double alpha, const Matrix& x);
-extern Matrix besseli (double alpha, const Matrix& x);
-extern Matrix besselk (double alpha, const Matrix& x);
-
-extern Matrix besselj (const Range& alpha, const ColumnVector& x);
-extern Matrix bessely (const Range& alpha, const ColumnVector& x);
-extern Matrix besseli (const Range& alpha, const ColumnVector& x);
-extern Matrix besselk (const Range& alpha, const ColumnVector& x);
-
-#endif
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
diff -r 0d640dc625c7 -r 3d5aefef14e2 liboctave/lo-specfun.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/lo-specfun.cc	Thu Feb 05 20:59:48 1998 +0000
@@ -0,0 +1,557 @@
+/*
+
+Copyright (C) 1996 John W. Eaton
+
+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 2, 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, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "Range.h"
+#include "dColVector.h"
+#include "dMatrix.h"
+#include "f77-fcn.h"
+#include "lo-error.h"
+#include "mx-inlines.cc"
+
+extern "C"
+{
+  int F77_FCN (rjbesl, RJBESL) (const double&, const double&,
+				const int&, double*, int&);
+
+  int F77_FCN (rybesl, RYBESL) (const double&, const double&,
+				const int&, double*, int&);
+
+  int F77_FCN (ribesl, RIBESL) (const double&, const double&,
+				const int&, const int&, double*, int&);
+
+  int F77_FCN (rkbesl, RKBESL) (const double&, const double&,
+				const int&, const int&, double*, int&);
+
+  int F77_FCN (xdacosh, XDACOSH) (const double&, double&);
+
+  int F77_FCN (xdasinh, XDASINH) (const double&, double&);
+
+  int F77_FCN (xdatanh, XDATANH) (const double&, double&);
+
+  int F77_FCN (xderf, XDERF) (const double&, double&);
+
+  int F77_FCN (xderfc, XDERFC) (const double&, double&);
+
+  int F77_FCN (xdbetai, XDBETAI) (const double&, const double&,
+				  const double&, double&);
+
+  int F77_FCN (xdgamma, XDGAMMA) (const double&, double&);
+
+  int F77_FCN (xdgami, XDGAMI) (const double&, const double&, double&);
+
+  int F77_FCN (dlgams, DLGAMS) (const double&, double&, double&);
+}
+
+#if !defined (HAVE_ACOSH)
+double
+acosh (double x)
+{
+  double retval;
+  F77_XFCN (dacosh, DACOSH, (x, retval));
+  return retval;
+}
+#endif
+
+#if !defined (HAVE_ASINH)
+double
+asinh (double x)
+{
+  double retval;
+  F77_XFCN (dasinh, DASINH, (x, retval));
+  return retval;
+}
+#endif
+
+#if !defined (HAVE_ATANH)
+double
+atanh (double x)
+{
+  double retval;
+  F77_XFCN (datanh, DATANH, (x, retval));
+  return retval;
+}
+#endif
+
+#if !defined (HAVE_ERF)
+double
+erf (double x)
+{
+  double retval;
+  F77_XFCN (derf, DERF, (x, retval));
+  return retval;
+}
+#endif
+
+#if !defined (HAVE_ERFC)
+double
+erfc (double x)
+{
+  double retval;
+  F77_XFCN (derfc, DERFC, (x, retval));
+  return retval;
+}
+#endif
+
+#if !defined (HAVE_GAMMA)
+double
+gamma (double x)
+{
+  double retval;
+  F77_XFCN (xdgamma, XDGAMMA, (x, retval));
+  return retval;
+}
+#endif
+
+#if !defined (HAVE_LGAMMA)
+// If the system doesn't have lgamma, assume that it doesn't have
+// signgam either.
+
+int signgam;
+
+double
+lgamma (double x)
+{
+  double retval;
+  double sgngam;
+
+  F77_XFCN (dlgams, DLGAMS, (x, retval, sgngam));
+
+  signgam = (int) sgngam;
+
+  return retval;
+}
+#endif
+
+int
+F77_FCN (ribesl, RIBESL) (const double& x, const double& alpha,
+			  const int& nb, double *result, int& ncalc)
+{
+  int retval = 0;
+  F77_FCN (ribesl, RIBESL) (x, alpha, nb, 1, result, ncalc);
+  return retval;
+}
+
+int
+F77_FCN (rkbesl, RKBESL) (const double& x, const double& alpha,
+			  const int& nb, double *result, int& ncalc)
+{
+  int retval = 0;
+  F77_FCN (rkbesl, RKBESL) (x, alpha, nb, 1, result, ncalc);
+  return retval;
+}
+
+typedef int (*fptr) (const double&, const double&, const int&, double*, int&);
+
+static Matrix
+do_bessel (fptr f, const char *fn, double alpha, const Matrix& x)
+{
+  Matrix retval;
+
+  if (alpha >= 0.0)
+    {
+      int nb = 1;
+
+      if (alpha >= 1.0)
+	{
+	  double integral;
+	  alpha = modf (alpha, &integral);
+	  nb = static_cast<int> (integral) + 1;
+	}
+
+      Array<double> tmp (nb);
+
+      double *ptmp = tmp.fortran_vec ();
+
+      int x_nr = x.rows ();
+      int x_nc = x.cols ();
+
+      retval.resize (x_nr, x_nc);
+
+      int ncalc;
+
+      for (int j = 0; j < x_nc; j++)
+	{
+	  for (int i = 0; i < x_nr; i++)
+	    {
+	      f (x(i,j), alpha, nb, ptmp, ncalc);
+
+	      if (ncalc < 0)
+		{
+		  (*current_liboctave_error_handler)
+		    ("%s: error computing function values", fn);
+
+		  return Matrix ();
+		}
+
+	      retval(i,j) = tmp(nb-1);
+	    }
+	}
+    }
+  else
+    (*current_liboctave_error_handler)
+      ("%s: alpha must be greater than zero", fn);
+
+  return retval;
+}
+
+static Matrix
+do_bessel (fptr f, const char *fn, const Range& alpha_range,
+	   const ColumnVector& x)
+{
+  Matrix retval;
+
+  double alpha_base = alpha_range.base ();
+
+  if (alpha_base >= 0.0)
+    {
+      int nb_orig = alpha_range.nelem ();
+      int nb = nb_orig;
+      int offset = 0;
+
+      if (alpha_base >= 1.0)
+	{
+	  double integral;
+	  alpha_base = modf (alpha_base, &integral);
+	  offset = static_cast<int> (integral);
+	  nb += offset;
+	}
+
+      Array<double> tmp (nb);
+
+      double *ptmp = tmp.fortran_vec ();
+
+      int x_len = x.length ();
+
+      retval.resize (x_len, nb_orig);
+
+      int ncalc;
+
+      for (int i = 0; i < x_len; i++)
+	{
+	  f (x(i), alpha_base, nb, ptmp, ncalc);
+
+	  if (ncalc < 0)
+	    {
+	      (*current_liboctave_error_handler)
+		("%s: error computing function values", fn);
+
+	      return Matrix ();
+	    }
+
+	  for (int j = 0; j < nb_orig; j++)
+	    retval(i,j) = tmp(offset+j);
+	}
+    }
+  else
+    (*current_liboctave_error_handler)
+      ("%s: alpha must be greater than zero", fn);
+
+  return retval;
+}
+
+Matrix
+besselj (double alpha, const Matrix& x)
+{
+  return do_bessel (F77_FCN (rjbesl, RJBESL), "besselj", alpha, x);
+}
+
+Matrix
+bessely (double alpha, const Matrix& x)
+{
+  return do_bessel (F77_FCN (rybesl, RYBESL), "bessely", alpha, x);
+}
+
+Matrix
+besseli (double alpha, const Matrix& x)
+{
+  return do_bessel (F77_FCN (ribesl, RIBESL), "besseli", alpha, x);
+}
+
+Matrix
+besselk (double alpha, const Matrix& x)
+{
+  return do_bessel (F77_FCN (rkbesl, RKBESL), "besselk", alpha, x);
+}
+
+Matrix
+besselj (const Range& alpha, const ColumnVector& x)
+{
+  return do_bessel (F77_FCN (rjbesl, RJBESL), "besselj", alpha, x);
+}
+
+Matrix
+bessely (const Range& alpha, const ColumnVector& x)
+{
+  return do_bessel (F77_FCN (rybesl, RYBESL), "bessely", alpha, x);
+}
+
+Matrix
+besseli (const Range& alpha, const ColumnVector& x)
+{
+  return do_bessel (F77_FCN (ribesl, RIBESL), "besseli", alpha, x);
+}
+
+Matrix
+besselk (const Range& alpha, const ColumnVector& x)
+{
+  return do_bessel (F77_FCN (rkbesl, RKBESL), "besselk", alpha, x);
+}
+
+static void
+gripe_betainc_nonconformant (int r1, int c1, int r2, int c2, int r3,
+			     int c3)
+{
+  (*current_liboctave_error_handler)
+   ("betainc: nonconformant arguments (x is %dx%d, a is %dx%d, b is %dx%d)",
+     r1, c1, r2, c2, r3, c3);
+}
+
+double
+betainc (double x, double a, double b)
+{
+  double retval;
+  F77_XFCN (xdbetai, XDBETAI, (x, a, b, retval));
+  return retval;
+}
+
+Matrix
+betainc (double x, double a, const Matrix& b)
+{
+  int nr = b.rows ();
+  int nc = b.cols ();
+
+  Matrix retval (nr, nc);
+
+  for (int j = 0; j < nc; j++)
+    for (int i = 0; i < nr; i++)
+      retval(i,j) = betainc (x, a, b(i,j));
+
+  return retval;
+}
+
+Matrix
+betainc (double x, const Matrix& a, double b)
+{
+  int nr = a.rows ();
+  int nc = a.cols ();
+
+  Matrix retval (nr, nc);
+
+  for (int j = 0; j < nc; j++)
+    for (int i = 0; i < nr; i++)
+      retval(i,j) = betainc (x, a(i,j), b);
+
+  return retval;
+}
+
+Matrix
+betainc (double x, const Matrix& a, const Matrix& b)
+{
+  Matrix retval;
+
+  int a_nr = a.rows ();
+  int a_nc = a.cols ();
+
+  int b_nr = b.rows ();
+  int b_nc = b.cols ();
+
+  if (a_nr == b_nr && a_nc == b_nc)
+    {
+      retval.resize (a_nr, a_nc);
+
+      for (int j = 0; j < a_nc; j++)
+	for (int i = 0; i < a_nr; i++)
+	  retval(i,j) = betainc (x, a(i,j), b(i,j));
+    }
+  else
+    gripe_betainc_nonconformant (1, 1, a_nr, a_nc, b_nr, b_nc);
+
+  return retval;
+}
+
+Matrix
+betainc (const Matrix& x, double a, double b)
+{
+  int nr = x.rows ();
+  int nc = x.cols ();
+
+  Matrix retval (nr, nc);
+
+  for (int j = 0; j < nc; j++)
+    for (int i = 0; i < nr; i++)
+      retval(i,j) = betainc (x(i,j), a, b);
+
+  return retval;
+}
+
+Matrix
+betainc (const Matrix& x, double a, const Matrix& b)
+{
+  Matrix retval;
+
+  int nr = x.rows ();
+  int nc = x.cols ();
+
+  int b_nr = b.rows ();
+  int b_nc = b.cols ();
+
+  if (nr == b_nr && nc == b_nc)
+    {
+      retval.resize (nr, nc);
+
+      for (int j = 0; j < nc; j++)
+	for (int i = 0; i < nr; i++)
+	  retval(i,j) = betainc (x(i,j), a, b(i,j));
+    }
+  else
+    gripe_betainc_nonconformant (nr, nc, 1, 1, b_nr, b_nc);
+
+  return retval;
+}
+
+Matrix
+betainc (const Matrix& x, const Matrix& a, double b)
+{
+  Matrix retval;
+
+  int nr = x.rows ();
+  int nc = x.cols ();
+
+  int a_nr = a.rows ();
+  int a_nc = a.cols ();
+
+  if (nr == a_nr && nc == a_nc)
+    {
+      retval.resize (nr, nc);
+
+      for (int j = 0; j < nc; j++)
+	for (int i = 0; i < nr; i++)
+	  retval(i,j) = betainc (x(i,j), a(i,j), b);
+    }
+  else
+    gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, 1, 1);
+
+  return retval;
+}
+
+Matrix
+betainc (const Matrix& x, const Matrix& a, const Matrix& b)
+{
+  Matrix retval;
+
+  int nr = x.rows ();
+  int nc = x.cols ();
+
+  int a_nr = a.rows ();
+  int a_nc = a.cols ();
+
+  int b_nr = b.rows ();
+  int b_nc = b.cols ();
+
+  if (nr == a_nr && nr == b_nr && nc == a_nc && nc == b_nc)
+    {
+      retval.resize (nr, nc);
+
+      for (int j = 0; j < nc; j++)
+	for (int i = 0; i < nr; i++)
+	  retval(i,j) = betainc (x(i,j), a(i,j), b(i,j));
+    }
+  else
+    gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, b_nr, b_nc);
+
+  return retval;
+}
+
+double
+gammainc (double x, double a)
+{
+  double retval;
+  F77_XFCN (xdgami, XDGAMI, (x, a, retval));
+  return retval;
+}
+
+Matrix
+gammainc (double x, const Matrix& a)
+{
+  int nr = a.rows ();
+  int nc = a.cols ();
+
+  Matrix retval (nr, nc);
+
+  for (int j = 0; j < nc; j++)
+    for (int i = 0; i < nr; i++)
+      retval(i,j) = gammainc (x, a(i,j));
+
+  return retval;
+}
+
+Matrix
+gammainc (const Matrix& x, double a)
+{
+  int nr = x.rows ();
+  int nc = x.cols ();
+
+  Matrix retval (nr, nc);
+
+  for (int j = 0; j < nc; j++)
+    for (int i = 0; i < nr; i++)
+      retval(i,j) = gammainc (x(i,j), a);
+
+  return retval;
+}
+
+Matrix
+gammainc (const Matrix& x, const Matrix& a)
+{
+  Matrix retval;
+
+  int nr = x.rows ();
+  int nc = x.cols ();
+
+  int a_nr = a.rows ();
+  int a_nc = a.cols ();
+
+  if (nr == a_nr && nc == a_nc)
+    {
+      retval.resize (nr, nc);
+
+      for (int j = 0; j < nc; j++)
+	for (int i = 0; i < nr; i++)
+	  retval(i,j) = gammainc (x(i,j), a(i,j));
+    }
+  else
+    (*current_liboctave_error_handler)
+      ("gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)",
+       nr, nc, a_nr, a_nc);
+
+  return retval;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff -r 0d640dc625c7 -r 3d5aefef14e2 liboctave/lo-specfun.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/lo-specfun.h	Thu Feb 05 20:59:48 1998 +0000
@@ -0,0 +1,90 @@
+/*
+
+Copyright (C) 1996, 1997 John W. Eaton
+
+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 2, 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, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if !defined (octave_liboctave_specfun_h)
+#define octave_liboctave_specfun_h 1
+
+class ColumnVector;
+class Matrix;
+class Range;
+
+#if !defined (HAVE_ACOSH)
+extern double acosh (double);
+#endif
+
+#if !defined (HAVE_ASINH)
+extern double asinh (double);
+#endif
+
+#if !defined (HAVE_ATANH)
+extern double atanh (double);
+#endif
+
+#if !defined (HAVE_ERF)
+extern double erf (double);
+#endif
+
+#if !defined (HAVE_ERFC)
+extern double erfc (double);
+#endif
+
+#if !defined (HAVE_GAMMA)
+extern double gamma (double);
+#endif
+
+#if !defined (HAVE_LGAMMA)
+extern double lgamma (double);
+#endif
+
+extern Matrix besselj (double alpha, const Matrix& x);
+extern Matrix bessely (double alpha, const Matrix& x);
+extern Matrix besseli (double alpha, const Matrix& x);
+extern Matrix besselk (double alpha, const Matrix& x);
+
+extern Matrix besselj (const Range& alpha, const ColumnVector& x);
+extern Matrix bessely (const Range& alpha, const ColumnVector& x);
+extern Matrix besseli (const Range& alpha, const ColumnVector& x);
+extern Matrix besselk (const Range& alpha, const ColumnVector& x);
+
+extern double betainc (double x, double a, double b);
+extern Matrix betainc (double x, double a, const Matrix& b);
+extern Matrix betainc (double x, const Matrix& a, double b);
+extern Matrix betainc (double x, const Matrix& a, const Matrix& b);
+
+extern Matrix betainc (const Matrix& x, double a, double b);
+extern Matrix betainc (const Matrix& x, double a, const Matrix& b);
+extern Matrix betainc (const Matrix& x, const Matrix& a, double b);
+extern Matrix betainc (const Matrix& x, const Matrix& a, const Matrix& b);
+
+extern double gammainc (double x, double a);
+extern Matrix gammainc (double x, const Matrix& a);
+extern Matrix gammainc (const Matrix& x, double a);
+extern Matrix gammainc (const Matrix& x, const Matrix& a);
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/