Mercurial > hg > octave-max

changeset 3220:3deb1105fbc1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1998-11-19 00:06:30 by jwe]
author jwe
date Thu, 19 Nov 1998 00:06:34 +0000
parents 30770ba4457a
children eba59b8c64dc
files configure.in kpathsea/ChangeLog kpathsea/kdefault.c liboctave/ChangeLog liboctave/lo-specfun.cc liboctave/lo-specfun.h src/ChangeLog src/DLD-FUNCTIONS/besselj.cc src/lex.l src/ov-base-mat.cc src/ov-base-mat.h src/ov-bool-mat.cc src/ov-bool-mat.h src/ov-ch-mat.cc src/ov-ch-mat.h src/ov-complex.h src/ov-cx-mat.cc src/ov-cx-mat.h src/ov-re-mat.cc src/ov-re-mat.h
diffstat 20 files changed, 1136 insertions(+), 442 deletions(-) [+]
line wrap: on
line diff
--- a/configure.in
+++ b/configure.in
@@ -21,7 +21,7 @@
 ### Software Foundation, 59 Temple Place - Suite 330, Boston, MA
 ### 02111-1307, USA. 
 
-AC_REVISION($Revision: 1.310 $)
+AC_REVISION($Revision: 1.311 $)
 AC_PREREQ(2.9)
 AC_INIT(src/octave.cc)
 AC_CONFIG_HEADER(config.h)
@@ -1145,12 +1145,12 @@
   doc/liboctave/Makefile doc/refcard/Makefile emacs/Makefile \
   examples/Makefile liboctave/Makefile src/Makefile \
   libcruft/Makefile libcruft/Makerules \
-  libcruft/blas/Makefile libcruft/dassl/Makefile \
-  libcruft/fftpack/Makefile libcruft/lapack/Makefile \
-  libcruft/linpack/Makefile libcruft/minpack/Makefile \
-  libcruft/misc/Makefile libcruft/odepack/Makefile \
-  libcruft/ordered-qz/Makefile libcruft/quadpack/Makefile \
-  libcruft/ranlib/Makefile libcruft/specfun/Makefile \
+  libcruft/amos/Makefile libcruft/blas/Makefile \
+  libcruft/dassl/Makefile libcruft/fftpack/Makefile \
+  libcruft/lapack/Makefile libcruft/linpack/Makefile \
+  libcruft/minpack/Makefile libcruft/misc/Makefile \
+  libcruft/odepack/Makefile libcruft/ordered-qz/Makefile \
+  libcruft/quadpack/Makefile libcruft/ranlib/Makefile \
   libcruft/slatec-fn/Makefile libcruft/slatec-err/Makefile \
   libcruft/villad/Makefile)
 
--- a/kpathsea/ChangeLog
+++ b/kpathsea/ChangeLog
@@ -1,3 +1,11 @@
+Fri Nov 13 08:24:34 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* kdefault.c (kpse_expand_default): In loop that checks for
+	doubled colon, just break out of the inner loop when a match is
+	found.  Don't check for expansion == path because we've already
+	duplicated path to avoid memory problems (see previous change).
+	From Rafael Laboissiere <rafael@icp.inpg.fr>.
+
 Fri Oct 23 22:05:46 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* kdefault.c (kpse_expand_default): Always return newly allocated
--- a/kpathsea/kdefault.c
+++ b/kpathsea/kdefault.c
@@ -62,7 +62,7 @@
       /* What we'll return if we find none.  */
       expansion = xstrdup (path);
 
-      for (loc = path; *loc && expansion == path; loc++)
+      for (loc = path; *loc; loc++)
         {
           if (IS_ENV_SEP (loc[0]) && IS_ENV_SEP (loc[1]))
             { /* We have a doubled colon.  */
@@ -75,6 +75,8 @@
               /* Copy in FALLBACK, and then the rest of PATH.  */
               strcat (expansion, fallback);
               strcat (expansion, loc + 1);
+
+	      break;
             }
         }
     }
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,8 @@
+Tue Nov 17 23:47:24 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* lo-specfun.cc (besselh, airy, biry): New functions.
+	Update Bessel function support to use library by D. E. Amos.
+
 Thu Nov 12 17:44:15 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* cmd-edit.h (command_editor::readline): Add new variation that
--- a/liboctave/lo-specfun.cc
+++ b/liboctave/lo-specfun.cc
@@ -25,25 +25,46 @@
 #endif
 
 #include "Range.h"
-#include "dColVector.h"
+#include "CColVector.h"
+#include "CMatrix.h"
+#include "dRowVector.h"
 #include "dMatrix.h"
 #include "f77-fcn.h"
 #include "lo-error.h"
+#include "lo-ieee.h"
+#include "lo-specfun.h"
 #include "mx-inlines.cc"
 
 extern "C"
 {
-  int F77_FCN (rjbesl, RJBESL) (const double&, const double&,
-				const int&, double*, int&);
+  int F77_FCN (zbesj, ZBESJ) (const double&, const double&,
+			      const double&, const int&, const int&,
+			      double*, double*, int&, int&);
 
-  int F77_FCN (rybesl, RYBESL) (const double&, const double&,
-				const int&, double*, int&);
+  int F77_FCN (zbesy, ZBESY) (const double&, const double&,
+			      const double&, const int&, const int&,
+			      double*, double*, int&,
+			      double*, double*, int&);
+
+  int F77_FCN (zbesi, ZBESI) (const double&, const double&,
+			      const double&, const int&, const int&,
+			      double*, double*, int&, int&);
 
-  int F77_FCN (ribesl, RIBESL) (const double&, const double&,
-				const int&, const int&, double*, int&);
+  int F77_FCN (zbesk, ZBESK) (const double&, const double&,
+			      const double&, const int&, const int&,
+			      double*, double*, int&, int&);
+
+  int F77_FCN (zbesh, ZBESH) (const double&, const double&,
+			      const double&, const int&, const int&,
+			      const int&, double*, double*, int&, int&);
 
-  int F77_FCN (rkbesl, RKBESL) (const double&, const double&,
-				const int&, const int&, double*, int&);
+  int F77_FCN (zairy, ZAIRY) (const double&, const double&,
+			      const int&, const int&,
+			      double&, double&, int&, int&);
+
+  int F77_FCN (zbiry, ZBIRY) (const double&, const double&,
+			      const int&, const int&,
+			      double&, double&, int&);
 
   int F77_FCN (xdacosh, XDACOSH) (const double&, double&);
 
@@ -132,179 +153,518 @@
   return result;
 }
 
-int
-F77_FCN (ribesl, RIBESL) (const double& x, const double& alpha,
-			  const int& nb, double *result, int& ncalc)
+static inline Complex
+zbesj (const Complex& z, double alpha, int kode, int& ierr);
+
+static inline Complex
+zbesy (const Complex& z, double alpha, int kode, int& ierr);
+
+static inline Complex
+zbesi (const Complex& z, double alpha, int kode, int& ierr);
+
+static inline Complex
+zbesk (const Complex& z, double alpha, int kode, int& ierr);
+
+static inline Complex
+zbesh1 (const Complex& z, double alpha, int kode, int& ierr);
+
+static inline Complex
+zbesh2 (const Complex& z, double alpha, int kode, int& ierr);
+
+static inline Complex
+bessel_return_value (const Complex& val, int ierr)
 {
-  int retval = 0;
-  F77_FCN (ribesl, RIBESL) (x, alpha, nb, 1, result, ncalc);
+  static const Complex inf_val = Complex (octave_Inf, octave_Inf);
+  static const Complex nan_val = Complex (octave_NaN, octave_NaN);
+
+  Complex retval;
+
+  switch (ierr)
+    {
+    case 0:
+    case 3:
+      retval = val;
+      break;
+
+    case 2:
+      retval = inf_val;
+      break;
+
+    default:
+      retval = nan_val;
+      break;
+    }
+
   return retval;
 }
 
-int
-F77_FCN (rkbesl, RKBESL) (const double& x, const double& alpha,
-			  const int& nb, double *result, int& ncalc)
+static inline Complex
+zbesj (const Complex& z, double alpha, int kode, int& ierr)
 {
-  int retval = 0;
-  F77_FCN (rkbesl, RKBESL) (x, alpha, nb, 1, result, ncalc);
+  Complex retval;
+
+  if (alpha >= 0.0)
+    {
+      double yr = 0.0;
+      double yi = 0.0;
+
+      int nz;
+
+      double zr = z.real ();
+      double zi = z.imag ();
+
+      F77_FCN (zbesj, ZBESJ) (zr, zi, alpha, kode, 1, &yr, &yi, nz, ierr);
+
+      if (zi == 0.0 && zr > 0.0)
+	yi = 0.0;
+
+      retval = bessel_return_value (Complex (yr, yi), ierr);
+    }
+  else
+    {
+      alpha = -alpha;
+
+      Complex tmp = cos (M_PI * alpha) * zbesj (z, alpha, kode, ierr);
+
+      if (ierr == 0 || ierr == 3)
+	{
+	  tmp -= sin (M_PI * alpha) * zbesy (z, alpha, kode, ierr);
+
+	  retval = bessel_return_value (tmp, ierr);
+	}
+      else
+	retval = Complex (octave_NaN, octave_NaN);
+    }
+
   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)
+static inline Complex
+zbesy (const Complex& z, double alpha, int kode, int& ierr)
 {
-  Matrix retval;
+  Complex retval;
 
   if (alpha >= 0.0)
     {
-      int nb = 1;
+      double yr = 0.0;
+      double yi = 0.0;
+
+      int nz;
+
+      double wr, wi;
 
-      if (alpha >= 1.0)
+      double zr = z.real ();
+      double zi = z.imag ();
+
+      ierr = 0;
+
+      if (zr == 0.0 && zi == 0.0)
 	{
-	  double integral;
-	  alpha = modf (alpha, &integral);
-	  nb = static_cast<int> (integral) + 1;
+	  yr = -octave_Inf;
+	  yi = 0.0;
+	}
+      else
+	{
+	  F77_FCN (zbesy, ZBESY) (zr, zi, alpha, kode, 1, &yr, &yi, nz,
+				  &wr, &wi, ierr);
+
+	  if (zi == 0.0 && zr > 0.0)
+	    yi = 0.0;
 	}
 
-      Array<double> tmp (nb);
-
-      double *ptmp = tmp.fortran_vec ();
+      return bessel_return_value (Complex (yr, yi), ierr);
+    }
+  else
+    {
+      alpha = -alpha;
 
-      int x_nr = x.rows ();
-      int x_nc = x.cols ();
+      Complex tmp = cos (M_PI * alpha) * zbesy (z, alpha, kode, ierr);
 
-      retval.resize (x_nr, x_nc);
-
-      int ncalc;
-
-      for (int j = 0; j < x_nc; j++)
+      if (ierr == 0 || ierr == 3)
 	{
-	  for (int i = 0; i < x_nr; i++)
-	    {
-	      f (x(i,j), alpha, nb, ptmp, ncalc);
+	  tmp += sin (M_PI * alpha) * zbesj (z, alpha, kode, ierr);
+
+	  retval = bessel_return_value (tmp, ierr);
+	}
+      else
+	retval = Complex (octave_NaN, octave_NaN);
+    }
+
+  return retval;
+}
+
+static inline Complex
+zbesi (const Complex& z, double alpha, int kode, int& ierr)
+{
+  Complex retval;
 
-	      if (ncalc < 0)
-		{
-		  (*current_liboctave_error_handler)
-		    ("%s: error computing function values", fn);
+  if (alpha >= 0.0)
+    {
+      double yr = 0.0;
+      double yi = 0.0;
+
+      int nz;
 
-		  return Matrix ();
-		}
+      double zr = z.real ();
+      double zi = z.imag ();
+
+      F77_FCN (zbesi, ZBESI) (zr, zi, alpha, kode, 1, &yr, &yi, nz, ierr);
 
-	      retval(i,j) = tmp(nb-1);
-	    }
-	}
+      if (zi == 0.0 && zr > 0.0)
+	yi = 0.0;
+
+      retval = bessel_return_value (Complex (yr, yi), ierr);
     }
   else
-    (*current_liboctave_error_handler)
-      ("%s: alpha must be greater than zero", fn);
+    {
+      alpha = -alpha;
+
+      Complex tmp = zbesi (z, alpha, kode, ierr);
+
+      if (ierr == 0 || ierr == 3)
+	{
+	  tmp += (2.0 / M_PI) * sin (M_PI * alpha)
+	    * zbesk (z, alpha, kode, ierr);
+
+	  retval = bessel_return_value (tmp, ierr);
+	}
+      else
+	retval = Complex (octave_NaN, octave_NaN);
+    }
+
+  return retval;
+}
+
+static inline Complex
+zbesk (const Complex& z, double alpha, int kode, int& ierr)
+{
+  Complex retval;
+
+  if (alpha >= 0.0)
+    {
+      double yr = 0.0;
+      double yi = 0.0;
+
+      int nz;
+
+      double zr = z.real ();
+      double zi = z.imag ();
+
+      ierr = 0;
+
+      if (zr == 0.0 && zi == 0.0)
+	{
+	  yr = octave_Inf;
+	  yi = 0.0;
+	}
+      else
+	{
+	  F77_FCN (zbesk, ZBESK) (zr, zi, alpha, kode, 1, &yr, &yi, nz, ierr);
+
+	  if (zi == 0.0 && zr > 0.0)
+	    yi = 0.0;
+	}
+
+      retval = bessel_return_value (Complex (yr, yi), ierr);
+    }
+  else
+    {
+      Complex tmp = zbesk (z, -alpha, kode, ierr);
+
+      retval = bessel_return_value (tmp, ierr);
+    }
 
   return retval;
 }
 
-static Matrix
-do_bessel (fptr f, const char *fn, const Range& alpha_range,
-	   const ColumnVector& x)
+static inline Complex
+zbesh1 (const Complex& z, double alpha, int kode, int& ierr)
 {
-  Matrix retval;
+  Complex retval;
 
-  double alpha_base = alpha_range.base ();
-
-  if (alpha_base >= 0.0)
+  if (alpha >= 0.0)
     {
-      int nb_orig = alpha_range.nelem ();
-      int nb = nb_orig;
-      int offset = 0;
+      double yr = 0.0;
+      double yi = 0.0;
+
+      int nz;
+
+      double zr = z.real ();
+      double zi = z.imag ();
 
-      if (alpha_base >= 1.0)
-	{
-	  double integral;
-	  alpha_base = modf (alpha_base, &integral);
-	  offset = static_cast<int> (integral);
-	  nb += offset;
-	}
+      F77_FCN (zbesh, ZBESH) (zr, zi, alpha, kode, 1, 1, &yr, &yi, nz, ierr);
 
-      Array<double> tmp (nb);
+      retval = bessel_return_value (Complex (yr, yi), ierr);
+    }
+  else
+    {
+      alpha = -alpha;
+
+      static const Complex eye = Complex (0.0, 1.0);
 
-      double *ptmp = tmp.fortran_vec ();
+      Complex tmp = exp (M_PI * alpha * eye) * zbesh1 (z, alpha, kode, ierr);
 
-      int x_len = x.length ();
-
-      retval.resize (x_len, nb_orig);
+      retval = bessel_return_value (tmp, ierr);
+    }
 
-      int ncalc;
+  return retval;
+}
 
-      for (int i = 0; i < x_len; i++)
-	{
-	  f (x(i), alpha_base, nb, ptmp, ncalc);
+static inline Complex
+zbesh2 (const Complex& z, double alpha, int kode, int& ierr)
+{
+  Complex retval;
 
-	  if (ncalc < 0)
-	    {
-	      (*current_liboctave_error_handler)
-		("%s: error computing function values", fn);
+  if (alpha >= 0.0)
+    {
+      double yr = 0.0;
+      double yi = 0.0;
+
+      int nz;
 
-	      return Matrix ();
-	    }
+      double zr = z.real ();
+      double zi = z.imag ();
 
-	  for (int j = 0; j < nb_orig; j++)
-	    retval(i,j) = tmp(offset+j);
-	}
+      F77_FCN (zbesh, ZBESH) (zr, zi, alpha, kode, 2, 1, &yr, &yi, nz, ierr);
+
+      retval = bessel_return_value (Complex (yr, yi), ierr);
     }
   else
-    (*current_liboctave_error_handler)
-      ("%s: alpha must be greater than zero", fn);
+    {
+      alpha = -alpha;
+
+      static const Complex eye = Complex (0.0, 1.0);
+
+      Complex tmp = exp (-M_PI * alpha * eye) * zbesh2 (z, alpha, kode, ierr);
+
+      retval = bessel_return_value (tmp, ierr);
+    }
+
+  return retval;
+}
+
+typedef Complex (*fptr) (const Complex&, double, int, int&);
+
+static inline Complex
+do_bessel (fptr f, const char *, double alpha, const Complex& x,
+	   bool scaled, int& ierr)
+{
+  Complex retval;
+
+  retval = f (x, alpha, (scaled ? 2 : 1), ierr);
+
+  return retval;
+}
+
+static inline ComplexMatrix
+do_bessel (fptr f, const char *, double alpha, const ComplexMatrix& x,
+	   bool scaled, Array2<int>& ierr)
+{
+  int nr = x.rows ();
+  int nc = x.cols ();
+
+  ComplexMatrix retval (nr, nc);
+
+  ierr.resize (nr, nc);
+
+  for (int j = 0; j < nc; j++)
+    for (int i = 0; i < nr; i++)
+      retval(i,j) = f (x(i,j), alpha, (scaled ? 2 : 1), ierr(i,j));
+
+  return retval;
+}
+
+static inline ComplexMatrix
+do_bessel (fptr f, const char *, const Matrix& alpha, const Complex& x,
+	   bool scaled, Array2<int>& ierr)
+{
+  int nr = alpha.rows ();
+  int nc = alpha.cols ();
+
+  ComplexMatrix retval (nr, nc);
+
+  ierr.resize (nr, nc);
+
+  for (int j = 0; j < nc; j++)
+    for (int i = 0; i < nr; i++)
+      retval(i,j) = f (x, alpha(i,j), (scaled ? 2 : 1), ierr(i,j));
 
   return retval;
 }
 
-Matrix
-besselj (double alpha, const Matrix& x)
+static inline ComplexMatrix
+do_bessel (fptr f, const char *fn, const Matrix& alpha,
+	   const ComplexMatrix& x, bool scaled, Array2<int>& ierr)
 {
-  return do_bessel (F77_FCN (rjbesl, RJBESL), "besselj", alpha, x);
+  ComplexMatrix retval;
+
+  int x_nr = x.rows ();
+  int x_nc = x.cols ();
+
+  int alpha_nr = alpha.rows ();
+  int alpha_nc = alpha.cols ();
+
+  if (x_nr == alpha_nr && x_nc == alpha_nc)
+    {
+      int nr = x_nr;
+      int nc = x_nc;
+
+      retval.resize (nr, nc);
+
+      ierr.resize (nr, nc);
+
+      for (int j = 0; j < nc; j++)
+	for (int i = 0; i < nr; i++)
+	  retval(i,j) = f (x(i,j), alpha(i,j), (scaled ? 2 : 1), ierr(i,j));
+    }
+  else
+    (*current_liboctave_error_handler)
+      ("%s: the sizes of alpha and x must conform", fn);
+
+  return retval;
 }
 
-Matrix
-bessely (double alpha, const Matrix& x)
+static inline ComplexMatrix
+do_bessel (fptr f, const char *, const RowVector& alpha,
+	   const ComplexColumnVector& x, bool scaled, Array2<int>& ierr)
 {
-  return do_bessel (F77_FCN (rybesl, RYBESL), "bessely", alpha, x);
-}
+  int nr = x.length ();
+  int nc = alpha.length ();
+
+  ComplexMatrix retval (nr, nc);
 
-Matrix
-besseli (double alpha, const Matrix& x)
-{
-  return do_bessel (F77_FCN (ribesl, RIBESL), "besseli", alpha, x);
+  ierr.resize (nr, nc);
+
+  for (int j = 0; j < nc; j++)
+    for (int i = 0; i < nr; i++)
+      retval(i,j) = f (x(i), alpha(j), (scaled ? 2 : 1), ierr(i,j));
+
+  return retval;
 }
 
-Matrix
-besselk (double alpha, const Matrix& x)
+#define SS_BESSEL(name, fcn) \
+  Complex \
+  name (double alpha, const Complex& x, bool scaled, int& ierr) \
+  { \
+    return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
+  }
+
+#define SM_BESSEL(name, fcn) \
+  ComplexMatrix \
+  name (double alpha, const ComplexMatrix& x, bool scaled, \
+	Array2<int>& ierr) \
+  { \
+    return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
+  }
+
+#define MS_BESSEL(name, fcn) \
+  ComplexMatrix \
+  name (const Matrix& alpha, const Complex& x, bool scaled, \
+	Array2<int>& ierr) \
+  { \
+    return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
+  }
+
+#define MM_BESSEL(name, fcn) \
+  ComplexMatrix \
+  name (const Matrix& alpha, const ComplexMatrix& x, bool scaled, \
+	Array2<int>& ierr) \
+  { \
+    return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
+  }
+
+#define RC_BESSEL(name, fcn) \
+  ComplexMatrix \
+  name (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, \
+        Array2<int>& ierr) \
+  { \
+    return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
+  }
+
+#define ALL_BESSEL(name, fcn) \
+  SS_BESSEL (name, fcn) \
+  SM_BESSEL (name, fcn) \
+  MS_BESSEL (name, fcn) \
+  MM_BESSEL (name, fcn) \
+  RC_BESSEL (name, fcn)
+
+ALL_BESSEL (besselj, zbesj)
+ALL_BESSEL (bessely, zbesy)
+ALL_BESSEL (besseli, zbesi)
+ALL_BESSEL (besselk, zbesk)
+ALL_BESSEL (besselh1, zbesh1)
+ALL_BESSEL (besselh2, zbesh2)
+
+Complex
+airy (const Complex& z, bool deriv, bool scaled, int& ierr)
 {
-  return do_bessel (F77_FCN (rkbesl, RKBESL), "besselk", alpha, x);
-}
+  double ar = 0.0;
+  double ai = 0.0;
+
+  int nz;
+
+  double zr = z.real ();
+  double zi = z.imag ();
 
-Matrix
-besselj (const Range& alpha, const ColumnVector& x)
-{
-  return do_bessel (F77_FCN (rjbesl, RJBESL), "besselj", alpha, x);
+  int id = deriv ? 1 : 0;
+
+  int kode = scaled ? 2 : 1;
+
+  F77_FCN (zairy, ZAIRY) (zr, zi, id, kode, ar, ai, nz, ierr);
+
+  return bessel_return_value (Complex (ar, ai), ierr);
 }
 
-Matrix
-bessely (const Range& alpha, const ColumnVector& x)
+Complex
+biry (const Complex& z, bool deriv, bool scaled, int& ierr)
 {
-  return do_bessel (F77_FCN (rybesl, RYBESL), "bessely", alpha, x);
+  double ar = 0.0;
+  double ai = 0.0;
+
+  double zr = z.real ();
+  double zi = z.imag ();
+
+  int id = deriv ? 1 : 0;
+
+  int kode = scaled ? 2 : 1;
+
+  F77_FCN (zbiry, ZBIRY) (zr, zi, id, kode, ar, ai, ierr);
+
+  return bessel_return_value (Complex (ar, ai), ierr);
 }
 
-Matrix
-besseli (const Range& alpha, const ColumnVector& x)
+ComplexMatrix
+airy (const ComplexMatrix& z, bool deriv, bool scaled, Array2<int>& ierr)
 {
-  return do_bessel (F77_FCN (ribesl, RIBESL), "besseli", alpha, x);
+  int nr = z.rows ();
+  int nc = z.cols ();
+
+  ComplexMatrix retval (nr, nc);
+
+  ierr.resize (nr, nc);
+
+  for (int j = 0; j < nc; j++)
+    for (int i = 0; i < nr; i++)
+      retval(i,j) = airy (z(i,j), deriv, scaled, ierr(i,j));
+
+  return retval;
 }
 
-Matrix
-besselk (const Range& alpha, const ColumnVector& x)
+ComplexMatrix
+biry (const ComplexMatrix& z, bool deriv, bool scaled, Array2<int>& ierr)
 {
-  return do_bessel (F77_FCN (rkbesl, RKBESL), "besselk", alpha, x);
+  int nr = z.rows ();
+  int nc = z.cols ();
+
+  ComplexMatrix retval (nr, nc);
+
+  ierr.resize (nr, nc);
+
+  for (int j = 0; j < nc; j++)
+    for (int i = 0; i < nr; i++)
+      retval(i,j) = biry (z(i,j), deriv, scaled, ierr(i,j));
+
+  return retval;
 }
 
 static void
--- a/liboctave/lo-specfun.h
+++ b/liboctave/lo-specfun.h
@@ -23,8 +23,13 @@
 #if !defined (octave_liboctave_specfun_h)
 #define octave_liboctave_specfun_h 1
 
-class ColumnVector;
+#include "oct-cmplx.h"
+
+template <class T> class Array2;
 class Matrix;
+class ComplexMatrix;
+class RowVector;
+class ComplexColumnVector;
 class Range;
 
 #if !defined (HAVE_ACOSH)
@@ -50,15 +55,128 @@
 extern double xgamma (double x);
 extern double xlgamma (double x);
 
-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 Complex
+besselj (double alpha, const Complex& x, bool scaled, int& ierr);
+
+extern Complex
+bessely (double alpha, const Complex& x, bool scaled, int& ierr);
+
+extern Complex
+besseli (double alpha, const Complex& x, bool scaled, int& ierr);
+
+extern Complex
+besselk (double alpha, const Complex& x, bool scaled, int& ierr);
+
+extern Complex
+besselh1 (double alpha, const Complex& x, bool scaled, int& ierr);
+
+extern Complex
+besselh2 (double alpha, const Complex& x, bool scaled, int& ierr);
+
+extern ComplexMatrix
+besselj (double alpha, const ComplexMatrix& x, bool scaled,
+	 Array2<int>& ierr);
+
+extern ComplexMatrix
+bessely (double alpha, const ComplexMatrix& x, bool scaled,
+	 Array2<int>& ierr);
+
+extern ComplexMatrix
+besseli (double alpha, const ComplexMatrix& x, bool scaled,
+	 Array2<int>& ierr);
+
+extern ComplexMatrix
+besselk (double alpha, const ComplexMatrix& x, bool scaled,
+	 Array2<int>& ierr);
+
+extern ComplexMatrix
+besselh1 (double alpha, const ComplexMatrix& x, bool scaled,
+	  Array2<int>& ierr);
+
+extern ComplexMatrix
+besselh2 (double alpha, const ComplexMatrix& x, bool scaled,
+	  Array2<int>& ierr);
+
+extern ComplexMatrix
+besselj (const Matrix& alpha, const Complex& x, bool scaled,
+	 Array2<int>& ierr);
+
+extern ComplexMatrix
+bessely (const Matrix& alpha, const Complex& x, bool scaled,
+	 Array2<int>& ierr);
+
+extern ComplexMatrix
+besseli (const Matrix& alpha, const Complex& x, bool scaled,
+	 Array2<int>& ierr);
+
+extern ComplexMatrix
+besselk (const Matrix& alpha, const Complex& x, bool scaled,
+	 Array2<int>& ierr);
 
-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 ComplexMatrix
+besselh1 (const Matrix& alpha, const Complex& x, bool scaled,
+	  Array2<int>& ierr);
+
+extern ComplexMatrix
+besselh2 (const Matrix& alpha, const Complex& x, bool scaled,
+	  Array2<int>& ierr);
+
+extern ComplexMatrix
+besselj (const Matrix& alpha, const ComplexMatrix& x, bool scaled,
+	 Array2<int>& ierr);
+
+extern ComplexMatrix
+bessely (const Matrix& alpha, const ComplexMatrix& x, bool scaled,
+	 Array2<int>& ierr);
+
+extern ComplexMatrix
+besseli (const Matrix& alpha, const ComplexMatrix& x, bool scaled,
+	 Array2<int>& ierr);
+
+extern ComplexMatrix
+besselk (const Matrix& alpha, const ComplexMatrix& x, bool scaled,
+	 Array2<int>& ierr);
+
+extern ComplexMatrix
+besselh1 (const Matrix& alpha, const ComplexMatrix& x, bool scaled,
+	  Array2<int>& ierr);
+
+extern ComplexMatrix
+besselh2 (const Matrix& alpha, const ComplexMatrix& x, bool scaled,
+	  Array2<int>& ierr);
+
+extern ComplexMatrix
+besselj (const RowVector& alpha, const ComplexColumnVector& x, bool scaled,
+	 Array2<int>& ierr);
+
+extern ComplexMatrix
+bessely (const RowVector& alpha, const ComplexColumnVector& x, bool scaled,
+	 Array2<int>& ierr);
+
+extern ComplexMatrix
+besseli (const RowVector& alpha, const ComplexColumnVector& x, bool scaled,
+	 Array2<int>& ierr);
+
+extern ComplexMatrix
+besselk (const RowVector& alpha, const ComplexColumnVector& x, bool scaled,
+	 Array2<int>& ierr);
+
+extern ComplexMatrix
+besselh1 (const RowVector& alpha, const ComplexColumnVector& x, bool scaled,
+	  Array2<int>& ierr);
+
+extern ComplexMatrix
+besselh2 (const RowVector& alpha, const ComplexColumnVector& x, bool scaled,
+	  Array2<int>& ierr);
+
+extern Complex airy (const Complex& z, bool deriv, bool scaled, int& ierr);
+extern Complex biry (const Complex& z, bool deriv, bool scaled, int& ierr);
+
+extern ComplexMatrix
+airy (const ComplexMatrix& z, bool deriv, bool scaled, Array2<int>& ierr);
+
+extern ComplexMatrix
+biry (const ComplexMatrix& z, bool deriv, bool scaled, Array2<int>& ierr);
 
 extern double betainc (double x, double a, double b);
 extern Matrix betainc (double x, double a, const Matrix& b);
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,31 @@
+Wed Nov 18 01:18:46 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* ov-base-mat.cc (octave_base_matrix::is_true): New function.
+	* ov-bool-mat.cc (octave_bool_matrix::is_true): Delete.
+	* ov-ch-mat.cc (octave_char_matrix::is_true): Delete.
+	* ov-cx-mat.cc (octave_complex_matrix::is_true): Delete.
+	* ov-re-mat.cc (octave_matrix::is_true): Delete.
+
+	* ov-base-mat.cc (octave_base_matrix::do_index_op): New function.
+	* ov-bool-mat.cc (octave_bool_matrix::do_index_op): Delete.
+	* ov-cx-mat.cc (octave_complex_matrix::do_index_op): Delete.
+	* ov-re-mat.cc (octave_matrix::do_index_op): Delete.
+
+	* mappers.cc: Don't include lo-specfun.h.
+
+Tue Nov 17 14:35:56 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* besselj.cc (Fbesselh, Fairy); New functions.
+	(Fbesselj, Fbessely, Fbesselk, Fbesseli): Update doc strings.
+	(do_bessel): Handle additional args.
+
+Fri Nov 13 14:47:11 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* lex.l (NUMBER): Allow hexadecimal constants.
+	(looks_like_hex): New function.
+	(handle_number): Check for hexadecimal constants and convert them
+	to unsigned integer values.
+
 Thu Nov 12 11:13:24 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* input.cc (gnu_readline): Check for EOF from command_editor::readline.
--- a/src/DLD-FUNCTIONS/besselj.cc
+++ b/src/DLD-FUNCTIONS/besselj.cc
@@ -32,25 +32,43 @@
 #include "oct-obj.h"
 #include "utils.h"
 
-#define DO_BESSEL(type, alpha, x) \
+enum bessel_type
+{
+  BESSEL_J,
+  BESSEL_Y,
+  BESSEL_I,
+  BESSEL_K,
+  BESSEL_H1,
+  BESSEL_H2
+};
+
+#define DO_BESSEL(type, alpha, x, scaled, ierr, result) \
   do \
     { \
       switch (type) \
 	{ \
-	  case 'j': \
-	    retval = besselj (alpha, x); \
+	  case BESSEL_J: \
+	    result = besselj (alpha, x, scaled, ierr); \
+	    break; \
+ \
+	  case BESSEL_Y: \
+	    result = bessely (alpha, x, scaled, ierr); \
 	    break; \
  \
-	  case 'y': \
-	    retval = bessely (alpha, x); \
+	  case BESSEL_I: \
+	    result = besseli (alpha, x, scaled, ierr); \
 	    break; \
  \
-	  case 'i': \
-	    retval = besseli (alpha, x); \
+	  case BESSEL_K: \
+	    result = besselk (alpha, x, scaled, ierr); \
 	    break; \
  \
-	  case 'k': \
-	    retval = besselk (alpha, x); \
+	  case BESSEL_H1: \
+	    result = besselh1 (alpha, x, scaled, ierr); \
+	    break; \
+ \
+	  case BESSEL_H2: \
+	    result = besselh2 (alpha, x, scaled, ierr); \
 	    break; \
  \
 	  default: \
@@ -59,78 +77,155 @@
     } \
   while (0)
 
+static inline Matrix
+int_array2_to_matrix (const Array2<int>& 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) = (double) (a(i,j));
+
+  return retval;
+}
+
 static void
-gripe_bessel_arg_1 (const char *fn)
+gripe_bessel_arg (const char *fn, const char *arg)
 {
-  error ("%s: alpha must be scalar or range with increment equal to 1", fn);
+  error ("%s: expecting scalar or matrix as %s argument", fn, arg);
 }
 
 octave_value_list
-do_bessel (char type, const char *fn, const octave_value_list& args)
+do_bessel (enum bessel_type type, const char *fn,
+	   const octave_value_list& args, int nargout)
 {
-  octave_value retval;
+  octave_value_list retval;
 
   int nargin = args.length ();
 
-  if (nargin == 2)
+  if (nargin == 2 || nargin == 3)
     {
+      bool scaled = (nargin == 3);
+
       octave_value alpha_arg = args(0);
+      octave_value x_arg = args(1);
 
       if (alpha_arg.is_scalar_type ())
 	{
-	  double alpha = alpha_arg.double_value ();
+	  double alpha = args(0).double_value ();
 
 	  if (! error_state)
 	    {
-	      Matrix x = args(1).matrix_value ();
+	      if (x_arg.is_scalar_type ())
+		{
+		  Complex x = x_arg.complex_value ();
+
+		  if (! error_state)
+		    {
+		      int ierr;
+		      octave_value result;
+
+		      DO_BESSEL (type, alpha, x, scaled, ierr, result);
+
+		      if (nargout > 1)
+			retval(1) = (double) ierr;
 
-	      if (! error_state)
-		DO_BESSEL (type, alpha, x);
+		      retval(0) = result;
+		    }
+		  else
+		    gripe_bessel_arg (fn, "second");
+		}
 	      else
-		error ("%s: expecting matrix as second argument", fn);
+		{
+		  ComplexMatrix x = x_arg.complex_matrix_value ();
+
+		  if (! error_state)
+		    {
+		      Array2<int> ierr;
+		      octave_value result;
+
+		      DO_BESSEL (type, alpha, x, scaled, ierr, result);
+
+		      if (nargout > 1)
+			retval(1) = int_array2_to_matrix (ierr);
+
+		      retval(0) = result;
+		    }
+		  else
+		    gripe_bessel_arg (fn, "second");
+		}
 	    }
 	  else
-	    gripe_bessel_arg_1 (fn);
+	    gripe_bessel_arg (fn, "first");
 	}
       else
 	{
-	  Range alpha;
+	  Matrix alpha = args(0).matrix_value ();
 
-	  if (! alpha_arg.is_range ())
+	  if (! error_state)
 	    {
-	      ColumnVector tmp = alpha_arg.vector_value ();
+	      if (x_arg.is_scalar_type ())
+		{
+		  Complex x = x_arg.complex_value ();
+
+		  if (! error_state)
+		    {
+		      Array2<int> ierr;
+		      octave_value result;
 
-	      if (! error_state)
+		      DO_BESSEL (type, alpha, x, scaled, ierr, result);
+
+		      if (nargout > 1)
+			retval(1) = int_array2_to_matrix (ierr);
+
+		      retval(0) = result;
+		    }
+		  else
+		    gripe_bessel_arg (fn, "second");
+		}
+	      else
 		{
-		  int len = tmp.length ();
+		  ComplexMatrix x = x_arg.complex_matrix_value ();
 
-		  double base = tmp(0);
+		  if (! error_state)
+		    {
+		      if (alpha.rows () == 1 && x.columns () == 1)
+			{
+			  RowVector ralpha = alpha.row (0);
+			  ComplexColumnVector cx = x.column (0);
+
+			  Array2<int> ierr;
+			  octave_value result;
+
+			  DO_BESSEL (type, ralpha, cx, scaled, ierr, result);
 
-		  for (int i = 1; i < len; i++)
-		    {
-		      if (tmp(i) != base + i)
+			  if (nargout > 1)
+			    retval(1) = int_array2_to_matrix (ierr);
+
+			  retval(0) = result;
+			}
+		      else
 			{
-			  gripe_bessel_arg_1 (fn);
-			  break;
+			  Array2<int> ierr;
+			  octave_value result;
+
+			  DO_BESSEL (type, alpha, x, scaled, ierr, result);
+
+			  if (nargout > 1)
+			    retval(1) = int_array2_to_matrix (ierr);
+
+			  retval(0) = result;
 			}
 		    }
-
-		  if (! error_state)
-		    alpha = Range (tmp(0), tmp(len-1));
+		  else
+		    gripe_bessel_arg (fn, "second");
 		}
 	    }
 	  else
-	    alpha = alpha_arg.range_value ();
-
-	  if (! error_state)
-	    {
-	      ColumnVector x = args(1).vector_value ();
-
-	      if (! error_state)
-		DO_BESSEL (type, alpha, x);
-	      else
-		error ("%s: expecting vector as second argument", fn);
-	    }
+	    gripe_bessel_arg (fn, "first");
 	}
     }
   else
@@ -139,72 +234,281 @@
   return retval;
 }
 
-DEFUN_DLD (besselj, args, ,
-  "besselj (alpha, x)\n\
+DEFUN_DLD (besselj, args, nargout,
+  "[J, IERR] = BESSELJ (ALPHA, X [, 1])\n\
 \n\
 Compute Bessel functions of the first kind.\n\
 \n\
-X must be a real matrix, vector or scalar.\n\
+If a third argument is supplied, scale the result by exp(-I*Z) for\n\
+K = 1 or exp(I*Z) for K = 2.\n\
+\n\
+If ALPHA is a scalar, the result is the same size as X.  If X is a\n\
+scalar, the result is the same size as ALPHA.  If ALPHA is a row\n\
+vector and X is a column vector, the result is a matrix with\n\
+length(X) rows and length(ALPHA) columns.  Otherwise, ALPHA and X must\n\
+conform and the result will be the same size.\n\
+\n\
+ALPHA must be real.  X may be complex.\n\
 \n\
-If ALPHA is a scalar, the result is the same size as X.  If ALPHA is a\n\
-range, X must be a vector or scalar, and the result is a matrix with\n\
-length(X) rows and length(ALPHA) columns.\n\
-\n\
-ALPHA must be greater than or equal to zero.  If ALPHA is a range, it\n\
-must have an increment equal to one.")
+If requested, IERR contains the following status information and is\n\
+the same size as the result.\n\
+\n
+  0  normal return\n\
+  1  input error, return NaN\n\
+  2  overflow, return Inf\n\
+  3  loss of significance by argument reduction results in less than\n\
+     half of machine accuracy\n\
+  4  complete loss of significance by argument reduction, return NaN\n\
+  5  error -- no computation, algorithm termination condition not met,\n\
+     return NaN")
 {
-  return do_bessel ('j', "besselj", args);
+  return do_bessel (BESSEL_J, "besselj", args, nargout);
 }
 
-DEFUN_DLD (bessely, args, ,
-  "bessely (alpha, x)\n\
+DEFUN_DLD (bessely, args, nargout,
+  "[Y, IERR] = BESSELY (ALPHA, X [, 1])\n\
 \n\
 Compute Bessel functions of the second kind.\n\
 \n\
-X must be a real matrix, vector or scalar.\n\
+If a third argument is supplied, scale the result by exp(-I*Z) for\n\
+K = 1 or exp(I*Z) for K = 2.\n\
+\n\
+If ALPHA is a scalar, the result is the same size as X.  If X is a\n\
+scalar, the result is the same size as ALPHA.  If ALPHA is a row\n\
+vector and X is a column vector, the result is a matrix with\n\
+length(X) rows and length(ALPHA) columns.  Otherwise, ALPHA and X must\n\
+conform and the result will be the same size.\n\
+\n\
+ALPHA must be real.  X may be complex.\n\
 \n\
-If ALPHA is a scalar, the result is the same size as X.  If ALPHA is a\n\
-range, X must be a vector or scalar, and the result is a matrix with\n\
-length(X) rows and length(ALPHA) columns.\n\
-\n\
-ALPHA must be greater than or equal to zero.  If ALPHA is a range, it\n\
-must have an increment equal to one.")
+If requested, IERR contains the following status information and is\n\
+the same size as the result.\n\
+\n
+  0  normal return\n\
+  1  input error, return NaN\n\
+  2  overflow, return Inf\n\
+  3  loss of significance by argument reduction results in less than\n\
+     half of machine accuracy\n\
+  4  complete loss of significance by argument reduction, return NaN\n\
+  5  error -- no computation, algorithm termination condition not met,\n\
+     return NaN")
 {
-  return do_bessel ('y', "bessely", args);
+  return do_bessel (BESSEL_Y, "bessely", args, nargout);
 }
 
-DEFUN_DLD (besseli, args, ,
-  "besseli (alpha, x)\n\
+DEFUN_DLD (besseli, args, nargout,
+  "[I, IERR] = BESSELI (ALPHA, X [, 1])\n\
 \n\
 Compute modified Bessel functions of the first kind.\n\
 \n\
-X must be a real matrix, vector or scalar.\n\
+If a third argument is supplied, scale the result by exp(-I*Z) for\n\
+K = 1 or exp(I*Z) for K = 2.\n\
+\n\
+If ALPHA is a scalar, the result is the same size as X.  If X is a\n\
+scalar, the result is the same size as ALPHA.  If ALPHA is a row\n\
+vector and X is a column vector, the result is a matrix with\n\
+length(X) rows and length(ALPHA) columns.  Otherwise, ALPHA and X must\n\
+conform and the result will be the same size.\n\
+\n\
+ALPHA must be real.  X may be complex.\n\
 \n\
-If ALPHA is a scalar, the result is the same size as X.  If ALPHA is a\n\
-range, X must be a vector or scalar, and the result is a matrix with\n\
-length(X) rows and length(ALPHA) columns.\n\
-\n\
-ALPHA must be greater than or equal to zero.  If ALPHA is a range, it\n\
-must have an increment equal to one.")
+If requested, IERR contains the following status information and is\n\
+the same size as the result.\n\
+\n
+  0  normal return\n\
+  1  input error, return NaN\n\
+  2  overflow, return Inf\n\
+  3  loss of significance by argument reduction results in less than\n\
+     half of machine accuracy\n\
+  4  complete loss of significance by argument reduction, return NaN\n\
+  5  error -- no computation, algorithm termination condition not met,\n\
+     return NaN")
 {
-  return do_bessel ('i', "besseli", args);
+  return do_bessel (BESSEL_I, "besseli", args, nargout);
 }
 
-DEFUN_DLD (besselk, args, ,
-  "besselk (alpha, x)\n\
+DEFUN_DLD (besselk, args, nargout,
+  "[K, IERR] = BESSELK (ALPHA, X [, 1])\n\
 \n\
 Compute modified Bessel functions of the second kind.\n\
 \n\
-X must be a real matrix, vector or scalar.\n\
+If a third argument is supplied, scale the result by exp(-I*Z) for\n\
+K = 1 or exp(I*Z) for K = 2.\n\
+\n\
+If ALPHA is a scalar, the result is the same size as X.  If X is a\n\
+scalar, the result is the same size as ALPHA.  If ALPHA is a row\n\
+vector and X is a column vector, the result is a matrix with\n\
+length(X) rows and length(ALPHA) columns.  Otherwise, ALPHA and X must\n\
+conform and the result will be the same size.\n\
+\n\
+ALPHA must be real.  X may be complex.\n\
+\n\
+If requested, IERR contains the following status information and is\n\
+the same size as the result.\n\
+\n
+  0  normal return\n\
+  1  input error, return NaN\n\
+  2  overflow, return Inf\n\
+  3  loss of significance by argument reduction results in less than\n\
+     half of machine accuracy\n\
+  4  complete loss of significance by argument reduction, return NaN\n\
+  5  error -- no computation, algorithm termination condition not met,\n\
+     return NaN")
+{
+  return do_bessel (BESSEL_K, "besselk", args, nargout);
+}
+
+DEFUN_DLD (besselh, args, nargout,
+  "[H, IERR] = besselh (ALPHA, K, X [, 1])\n\
+\n\
+Compute Hankel functions of the first (K = 1) or second (K = 2) kind.\n\
+\n\
+If a fourth argument is supplied, scale the result by exp(-I*Z) for\n\
+K = 1 or exp(I*Z) for K = 2.\n\
+\n\
+If ALPHA is a scalar, the result is the same size as X.  If X is a\n\
+scalar, the result is the same size as ALPHA.  If ALPHA is a row\n\
+vector and X is a column vector, the result is a matrix with\n\
+length(X) rows and length(ALPHA) columns.  Otherwise, ALPHA and X must\n\
+conform and the result will be the same size.\n\
+\n\
+ALPHA must be real.  X may be complex.\n\
 \n\
-If ALPHA is a scalar, the result is the same size as X.  If ALPHA is a\n\
-range, X must be a vector or scalar, and the result is a matrix with\n\
-length(X) rows and length(ALPHA) columns.\n\
+If requested, IERR contains the following status information and is\n\
+the same size as the result.\n\
+\n
+  0  normal return\n\
+  1  input error, return NaN\n\
+  2  overflow, return Inf\n\
+  3  loss of significance by argument reduction results in less than\n\
+     half of machine accuracy\n\
+  4  complete loss of significance by argument reduction, return NaN\n\
+  5  error -- no computation, algorithm termination condition not met,\n\
+     return NaN")
+{
+  octave_value_list retval;
+
+  int nargin = args.length ();
+
+  int kind = 1;
+
+  if (nargin == 2)
+    {
+      retval = do_bessel (BESSEL_H1, "besselh", args, nargout);
+    }
+  else if (nargin == 3 || nargin == 4)
+    {
+      double d_kind = args(1).double_value ();
+
+      if (! error_state && D_NINT (d_kind) == d_kind)
+	{
+	  octave_value_list tmp_args;
+
+	  if (nargin == 4)
+	    tmp_args(2) = args(3);
+
+	  tmp_args(1) = args(2);
+	  tmp_args(0) = args(0);
+
+	  if (kind == 1)
+	    retval = do_bessel (BESSEL_H1, "besselh", tmp_args, nargout);
+	  else if (kind == 2)
+	    retval = do_bessel (BESSEL_H2, "besselh", tmp_args, nargout);
+	  else
+	    error ("besselh: expecting K = 1 or 2");
+	}
+      else
+	error ("besselh: invalid value of K");
+    }
+  else
+    print_usage ("besselh");
+
+  return retval;
+}
+
+DEFUN_DLD (airy, args, nargout,
+  "[A, IERR] = airy (K, Z, [, 1])\n\
+\n\
+Compute Airy functions of the first and second kind, and their\n\
+derivatives.\n\
+\n\
+  K   Function   Scale factor (if a third argument is supplied)\n\
+ ---  --------   ----------------------------------------------\n\
+  0   Ai (Z)     exp ((2/3) * Z * sqrt (Z))\n\
+  1   dAi(Z)/dZ  exp ((2/3) * Z * sqrt (Z))\n\
+  2   Bi (Z)     exp (-abs (real ((2/3) * Z *sqrt (Z))))\n\
+  3   dBi(Z)/dZ  exp (-abs (real ((2/3) * Z *sqrt (Z))))\n\
+\n\
+The function call airy (Z) is equivalent to airy (0, Z).\n\
 \n\
-ALPHA must be greater than or equal to zero.  If ALPHA is a range, it\n\
-must have an increment equal to one.")
+The result is the same size as Z.
+\n\
+If requested, IERR contains the following status information and is\n\
+the same size as the result.\n\
+\n
+  0  normal return\n\
+  1  input error, return NaN\n\
+  2  overflow, return Inf\n\
+  3  loss of significance by argument reduction results in less than\n\
+     half of machine accuracy\n\
+  4  complete loss of significance by argument reduction, return NaN\n\
+  5  error -- no computation, algorithm termination condition not met,\n\
+     return NaN")
 {
-  return do_bessel ('k', "besselk", args);
+  octave_value_list retval;
+
+  int nargin = args.length ();
+
+  if (nargin > 0 && nargin < 4)
+    {
+      bool scale = (nargin == 3);
+
+      int kind = 0;
+
+      ComplexMatrix z;
+
+      if (nargin > 1)
+	{
+	  double d_kind = args(0).double_value ();
+
+	  if (! error_state)
+	    {
+	      kind = (int) d_kind;
+
+	      if (kind < 0 || kind > 3)
+		error ("airy: expecting K = 0, 1, 2, or 3");
+	    }	      
+	  else
+	    error ("airy: expecting integer value for K");
+	}
+
+      if (! error_state)
+	{
+	  z = args(nargin == 1 ? 0 : 1).complex_matrix_value ();
+
+	  if (! error_state)
+	    {
+	      Array2<int> ierr;
+	      octave_value result;
+
+	      if (kind > 1)
+		result = biry (z, kind == 3, scale, ierr);
+	      else
+		result = airy (z, kind == 1, scale, ierr);
+
+	      if (nargout > 1)
+		retval(1) = int_array2_to_matrix (ierr);
+
+	      retval(0) = result;
+	    }
+	  else
+	    error ("airy: expecting complex matrix for Z");
+	}
+    }
+  else
+    print_usage ("airy");
+
+  return retval;
 }
 
 /*
--- a/src/lex.l
+++ b/src/lex.l
@@ -201,7 +201,7 @@
 NOT	((\~)|(\!))
 IDENT	([_a-zA-Z][_a-zA-Z0-9]*)
 EXPON	([DdEe][+-]?{D}+)
-NUMBER	(({D}+\.?{D}*{EXPON}?)|(\.{D}+{EXPON}?))
+NUMBER	(({D}+\.?{D}*{EXPON}?)|(\.{D}+{EXPON}?)|(0[xX][0-9a-fA-F]+))
 %%
 
 %{
@@ -1383,6 +1383,12 @@
   return retval;
 }
 
+static inline bool
+looks_like_hex (const char *s, int len)
+{
+  return (len > 2 && s[0] == '0' && (s[1] == 'x' || s[1] == 'X'));
+}
+
 static void
 handle_number (char *yytext)
 {
@@ -1393,8 +1399,17 @@
   if (idx)
     *idx = 'e';
 
-  double value;
-  int nread = sscanf (tmp, "%lf", &value);
+  double value = 0.0;
+  int nread = 0;
+
+  if (looks_like_hex (tmp, strlen (tmp)))
+    {
+      unsigned long ival;
+      nread = sscanf (tmp, "%lx", &ival);
+      value = static_cast<double> (ival);
+    }
+  else
+    nread = sscanf (tmp, "%lf", &value);
 
   delete [] tmp;
 
--- a/src/ov-base-mat.cc
+++ b/src/ov-base-mat.cc
@@ -30,10 +30,76 @@
 
 #include <iostream.h>
 
+#include "oct-obj.h"
 #include "ov-base.h"
 #include "ov-base-mat.h"
 
 template <class MT>
+octave_value
+octave_base_matrix<MT>::do_index_op (const octave_value_list& idx)
+{
+  octave_value retval;
+
+  int len = idx.length ();
+
+  switch (len)
+    {
+    case 2:
+      {
+	idx_vector i = idx (0).index_vector ();
+	idx_vector j = idx (1).index_vector ();
+
+	retval = MT (matrix.index (i, j));
+      }
+      break;
+
+    case 1:
+      {
+	idx_vector i = idx (0).index_vector ();
+
+	retval = MT (matrix.index (i));
+      }
+      break;
+
+    default:
+      {
+	string n = type_name ();
+
+	error ("invalid number of indices (%d) for %s value",
+	       len, n.c_str ());
+      }
+      break;
+    }
+
+  return retval;
+}
+
+template <class MT>
+bool
+octave_base_matrix<MT>::is_true (void) const
+{
+  bool retval = false;
+
+  if (rows () == 0 || columns () == 0)
+    {
+      int flag = Vpropagate_empty_matrices;
+
+      if (flag < 0)
+	warning ("empty matrix used in conditional expression");
+      else if (flag == 0)
+	error ("empty matrix used in conditional expression");
+    }
+  else
+    {
+      boolMatrix m = (matrix.all ()) . all ();
+
+      retval = (m.rows () == 1 && m.columns () == 1 && m (0, 0) != 0.0);
+    }
+
+  return retval;
+}
+
+template <class MT>
 bool
 octave_base_matrix<MT>::print_as_scalar (void) const
 {
--- a/src/ov-base-mat.h
+++ b/src/ov-base-mat.h
@@ -65,6 +65,8 @@
 
   octave_value *clone (void) { return new octave_base_matrix (*this); }
 
+  octave_value do_index_op (const octave_value_list& idx);
+
   int rows (void) const { return matrix.rows (); }
   int columns (void) const { return matrix.columns (); }
 
@@ -80,6 +82,8 @@
 
   bool is_constant (void) const { return true; }
 
+  bool is_true (void) const;
+
   virtual bool print_as_scalar (void) const;
 
   void print (ostream& os, bool pr_as_read_syntax = false) const;
--- a/src/ov-bool-mat.cc
+++ b/src/ov-bool-mat.cc
@@ -78,40 +78,6 @@
   return retval;
 }
 
-octave_value
-octave_bool_matrix::do_index_op (const octave_value_list& idx)
-{
-  octave_value retval;
-
-  int len = idx.length ();
-
-  switch (len)
-    {
-    case 2:
-      {
-	idx_vector i = idx (0).index_vector ();
-	idx_vector j = idx (1).index_vector ();
-
-	retval = boolMatrix (matrix.index (i, j));
-      }
-      break;
-
-    case 1:
-      {
-	idx_vector i = idx (0).index_vector ();
-
-	retval = boolMatrix (matrix.index (i));
-      }
-      break;
-
-    default:
-      error ("invalid number of indices (%d) for matrix value", len);
-      break;
-    }
-
-  return retval;
-}
-
 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
 extern void assign (Array2<bool>&, const Array2<bool>&);
 #endif
@@ -160,30 +126,6 @@
   return false;
 }
 
-bool
-octave_bool_matrix::is_true (void) const
-{
-  bool retval = false;
-
-  if (rows () == 0 || columns () == 0)
-    {
-      int flag = Vpropagate_empty_matrices;
-
-      if (flag < 0)
-	warning ("empty matrix used in conditional expression");
-      else if (flag == 0)
-	error ("empty matrix used in conditional expression");
-    }
-  else
-    {
-      boolMatrix m = (matrix.all ()) . all ();
-
-      retval = (m.rows () == 1 && m.columns () == 1 && m (0, 0));
-    }
-
-  return retval;
-}
-
 double
 octave_bool_matrix::double_value (bool) const
 {
--- a/src/ov-bool-mat.h
+++ b/src/ov-bool-mat.h
@@ -54,13 +54,13 @@
 public:
 
   octave_bool_matrix (void)
-    : octave_base_matrix () { }
+    : octave_base_matrix<boolMatrix> () { }
 
   octave_bool_matrix (const boolMatrix& bm)
-    : octave_base_matrix (bm) { }
+    : octave_base_matrix<boolMatrix> (bm) { }
 
   octave_bool_matrix (const octave_bool_matrix& bm)
-    : octave_base_matrix (bm) { }
+    : octave_base_matrix<boolMatrix> (bm) { }
 
   ~octave_bool_matrix (void) { }
 
@@ -70,8 +70,6 @@
 
   octave_value *try_narrowing_conversion (void);
 
-  octave_value do_index_op (const octave_value_list& idx);
-
   void assign (const octave_value_list& idx, const boolMatrix& rhs);
 
   idx_vector index_vector (void) const { return idx_vector (matrix); }
@@ -93,8 +91,6 @@
 
   bool valid_as_zero_index (void) const { return is_zero_by_zero (); }
 
-  bool is_true (void) const;
-
   double double_value (bool = false) const;
 
   double scalar_value (bool frc_str_conv = false) const
--- a/src/ov-ch-mat.cc
+++ b/src/ov-ch-mat.cc
@@ -62,30 +62,6 @@
   return retval;
 }
 
-bool
-octave_char_matrix::is_true (void) const
-{
-  bool retval = false;
-
-  if (rows () == 0 || columns () == 0)
-    {
-      int flag = Vpropagate_empty_matrices;
-
-      if (flag < 0)
-	warning ("empty matrix used in conditional expression");
-      else if (flag == 0)
-	error ("empty matrix used in conditional expression");
-    }
-  else
-    {
-      Matrix m = (matrix.all ()) . all ();
-
-      retval = (m.rows () == 1 && m.columns () == 1 && m (0, 0) != 0.0);
-    }
-
-  return retval;
-}
-
 double
 octave_char_matrix::double_value (bool) const
 {
--- a/src/ov-ch-mat.h
+++ b/src/ov-ch-mat.h
@@ -94,8 +94,6 @@
   bool valid_as_scalar_index (void) const;
   bool valid_as_zero_index (void) const;
 
-  bool is_true (void) const;
-
   double double_value (bool = false) const;
 
   double scalar_value (bool frc_str_conv = false) const
--- a/src/ov-complex.h
+++ b/src/ov-complex.h
@@ -102,8 +102,6 @@
 
   bool is_true (void) const { return (scalar != 0.0); }
 
-  bool is_empty (void) const { return (rows () == 0 && columns () == 0); }
-
   double double_value (bool = false) const;
 
   double scalar_value (bool frc_str_conv = false) const
--- a/src/ov-cx-mat.cc
+++ b/src/ov-cx-mat.cc
@@ -90,40 +90,6 @@
   return retval;
 }
 
-octave_value
-octave_complex_matrix::do_index_op (const octave_value_list& idx)
-{
-  octave_value retval;
-
-  int len = idx.length ();
-
-  switch (len)
-    {
-    case 2:
-      {
-	idx_vector i = idx (0).index_vector ();
-	idx_vector j = idx (1).index_vector ();
-
-	retval = ComplexMatrix (matrix.index (i, j));
-      }
-      break;
-
-    case 1:
-      {
-	idx_vector i = idx (0).index_vector ();
-
-	retval = ComplexMatrix (matrix.index (i));
-      }
-      break;
-
-    default:
-      error ("invalid number of indices (%d) for complex matrix value", len);
-      break;
-    }
-
-  return retval;
-}
-
 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
 extern void assign (Array2<Complex>&, const Array2<Complex>&);
 #endif
@@ -220,30 +186,6 @@
   return false;
 }
 
-bool
-octave_complex_matrix::is_true (void) const
-{
-  bool retval = false;
-
-  if (rows () == 0 || columns () == 0)
-    {
-      int flag = Vpropagate_empty_matrices;
-
-      if (flag < 0)
-	warning ("empty matrix used in conditional expression");
-      else if (flag == 0)
-	error ("empty matrix used in conditional expression");
-    }
-  else
-    {
-      Matrix m = (matrix.all ()) . all ();
-
-      retval = (m.rows () == 1 && m.columns () == 1 && m (0, 0) != 0.0);
-    }
-
-  return retval;
-}
-
 double
 octave_complex_matrix::double_value (bool force_conversion) const
 {
--- a/src/ov-cx-mat.h
+++ b/src/ov-cx-mat.h
@@ -76,8 +76,6 @@
 
   octave_value *try_narrowing_conversion (void);
 
-  octave_value do_index_op (const octave_value_list& idx);
-
   void assign (const octave_value_list& idx, const ComplexMatrix& rhs);
 
   void assign (const octave_value_list& idx, const Matrix& rhs);
@@ -96,10 +94,6 @@
   bool valid_as_scalar_index (void) const;
   bool valid_as_zero_index (void) const;
 
-  bool is_true (void) const;
-
-  bool is_empty (void) const { return (rows () == 0 && columns () == 0); }
-
   double double_value (bool = false) const;
 
   double scalar_value (bool frc_str_conv = false) const
--- a/src/ov-re-mat.cc
+++ b/src/ov-re-mat.cc
@@ -58,7 +58,7 @@
 
 octave_matrix::octave_matrix (const ColumnVector& v, int pcv)
   : octave_base_matrix<Matrix> ((pcv < 0 && Vprefer_column_vectors) || pcv
-			       ? Matrix (v) : Matrix (v.transpose ())) { }
+				? Matrix (v) : Matrix (v.transpose ())) { }
 
 octave_value *
 octave_matrix::try_narrowing_conversion (void)
@@ -74,40 +74,6 @@
   return retval;
 }
 
-octave_value
-octave_matrix::do_index_op (const octave_value_list& idx)
-{
-  octave_value retval;
-
-  int len = idx.length ();
-
-  switch (len)
-    {
-    case 2:
-      {
-	idx_vector i = idx (0).index_vector ();
-	idx_vector j = idx (1).index_vector ();
-
-	retval = Matrix (matrix.index (i, j));
-      }
-      break;
-
-    case 1:
-      {
-	idx_vector i = idx (0).index_vector ();
-
-	retval = Matrix (matrix.index (i));
-      }
-      break;
-
-    default:
-      error ("invalid number of indices (%d) for matrix value", len);
-      break;
-    }
-
-  return retval;
-}
-
 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
 extern void assign (Array2<double>&, const Array2<double>&);
 #endif
@@ -262,30 +228,6 @@
   return false;
 }
 
-bool
-octave_matrix::is_true (void) const
-{
-  bool retval = false;
-
-  if (rows () == 0 || columns () == 0)
-    {
-      int flag = Vpropagate_empty_matrices;
-
-      if (flag < 0)
-	warning ("empty matrix used in conditional expression");
-      else if (flag == 0)
-	error ("empty matrix used in conditional expression");
-    }
-  else
-    {
-      Matrix m = (matrix.all ()) . all ();
-
-      retval = (m.rows () == 1 && m.columns () == 1 && m (0, 0) != 0.0);
-    }
-
-  return retval;
-}
-
 double
 octave_matrix::double_value (bool) const
 {
--- a/src/ov-re-mat.h
+++ b/src/ov-re-mat.h
@@ -76,8 +76,6 @@
 
   octave_value *try_narrowing_conversion (void);
 
-  octave_value do_index_op (const octave_value_list& idx);
-
   void assign (const octave_value_list& idx, const Matrix& rhs);
 
   void assign_struct_elt (assign_op, const string& elt_nm,
@@ -112,8 +110,6 @@
 
   bool valid_as_zero_index (void) const { return is_zero_by_zero (); }
 
-  bool is_true (void) const;
-
   double double_value (bool = false) const;
 
   double scalar_value (bool frc_str_conv = false) const