--- a/build-aux/mk-opts.pl	Fri Jul 27 17:10:25 2012 -0400
+++ b/build-aux/mk-opts.pl	Fri Jul 27 21:59:10 2012 -0400
@@ -513,7 +513,7 @@
 
 #include "$header"
 
-#include "defun-dld.h"
+#include "defun.h"
 #include "pr-output.h"
 
 #include "oct-obj.h"
@@ -909,11 +909,11 @@
 sub emit_options_function
 {
   print <<"_END_EMIT_OPTIONS_FUNCTION_HDR_";
-DEFUN_DLD ($OPT_FCN_NAME, args, ,
+DEFUN ($OPT_FCN_NAME, args, ,
   "-*- texinfo -*-\\n\\
-\@deftypefn  {Loadable Function} {} $OPT_FCN_NAME ()\\n\\
-\@deftypefnx {Loadable Function} {val =} $OPT_FCN_NAME (\@var{opt})\\n\\
-\@deftypefnx {Loadable Function} {} $OPT_FCN_NAME (\@var{opt}, \@var{val})\\n\\
+\@deftypefn  {Built-in Function} {} $OPT_FCN_NAME ()\\n\\
+\@deftypefnx {Built-in Function} {val =} $OPT_FCN_NAME (\@var{opt})\\n\\
+\@deftypefnx {Built-in Function} {} $OPT_FCN_NAME (\@var{opt}, \@var{val})\\n\\
 $DOC_STRING\\n\\
 \\n\\
 Options include\\n\\

--- a/doc/interpreter/doccheck/aspell-octave.en.pws	Fri Jul 27 17:10:25 2012 -0400
+++ b/doc/interpreter/doccheck/aspell-octave.en.pws	Fri Jul 27 21:59:10 2012 -0400
@@ -450,6 +450,7 @@
 isreal
 issparse
 isvector
+iter
 ith
 iy
 Jacobian
@@ -611,6 +612,7 @@
 nd
 ndgrid
 ne
+Nelder
 neq
 Neudecker
 Neumann

--- a/scripts/help/type.m	Fri Jul 27 17:10:25 2012 -0400
+++ b/scripts/help/type.m	Fri Jul 27 21:59:10 2012 -0400
@@ -118,7 +118,7 @@
 %! typestr = typestr{1}(1:17);
 %! assert (typestr, "var is a variable");
 
-%!assert (type ("dot"){1}, "dot is a dynamically-linked function")
+%!assert (type ("amd"){1}, "amd is a dynamically-linked function")
 %!assert (type ("cat"){1}, "cat is a built-in function")
 %!assert (type ("+"){1}, "+ is an operator")
 %!assert (type ("end"){1}, "end is a keyword")

--- a/scripts/help/which.m	Fri Jul 27 17:10:25 2012 -0400
+++ b/scripts/help/which.m	Fri Jul 27 21:59:10 2012 -0400
@@ -59,8 +59,8 @@
 %! str = which ("ls");
 %! assert (str(end-17:end), strcat ("miscellaneous", filesep (), "ls.m"));
 %!test
-%! str = which ("dot");
-%! assert (str(end-6:end), "dot.oct");
+%! str = which ("amd");
+%! assert (str(end-6:end), "amd.oct");
 
 %!assert (which ("_NO_SUCH_NAME_"), "")
 

--- a/scripts/image/cmunique.m	Fri Jul 27 17:10:25 2012 -0400
+++ b/scripts/image/cmunique.m	Fri Jul 27 21:59:10 2012 -0400
@@ -29,8 +29,8 @@
 ## eliminated.  The output image, @var{Y}, is the original input image with
 ## the indices adjusted to match the new, possibly smaller, colormap.
 ##
-## When the input is an RGB image (an MxNx3 array), the output colormap
-## will contain one entry for every unique color in the original image.  
+## When the input is an RGB image (an @nospell{MxNx3} array), the output
+## colormap will contain one entry for every unique color in the original image.
 ## In the worst case the new map could have as many rows as the number of
 ## pixels in the original image.
 ##

--- a/src/DLD-FUNCTIONS/__contourc__.cc	Fri Jul 27 17:10:25 2012 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,340 +0,0 @@
-/* Contour lines for function evaluated on a grid.
-
-Copyright (C) 2007-2012 Kai Habel
-Copyright (C) 2004, 2007 Shai Ayal
-
-Adapted to an oct file from the stand alone contourl by Victro Munoz
-Copyright (C) 2004 Victor Munoz
-
-Based on contour plot routine (plcont.c) in PLPlot package
-http://plplot.org/
-
-Copyright (C) 1995, 2000, 2001 Maurice LeBrun
-Copyright (C) 2000, 2002 Joao Cardoso
-Copyright (C) 2000, 2001, 2002, 2004  Alan W. Irwin
-Copyright (C) 2004  Andrew Ross
-
-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 <cfloat>
-
-#include "quit.h"
-
-#include "defun-dld.h"
-#include "error.h"
-#include "oct-obj.h"
-
-static Matrix this_contour;
-static Matrix contourc;
-static int elem;
-
-// This is the quanta in which we increase this_contour.
-#define CONTOUR_QUANT 50
-
-// Add a coordinate point (x,y) to this_contour.
-
-static void
-add_point (double x, double y)
-{
-  if (elem % CONTOUR_QUANT == 0)
-    this_contour = this_contour.append (Matrix (2, CONTOUR_QUANT, 0));
-
-  this_contour (0, elem) = x;
-  this_contour (1, elem) = y;
-  elem++;
-}
-
-// Add contents of current contour to contourc.
-// this_contour.cols () - 1;
-
-static void
-end_contour (void)
-{
-  if (elem > 2)
-    {
-      this_contour (1, 0) = elem - 1;
-      contourc = contourc.append (this_contour.extract_n (0, 0, 2, elem));
-    }
-
-  this_contour = Matrix ();
-  elem = 0;
-}
-
-// Start a new contour, and add contents of current one to contourc.
-
-static void
-start_contour (double lvl, double x, double y)
-{
-  end_contour ();
-  this_contour.resize (2, 0);
-  add_point (lvl, 0);
-  add_point (x, y);
-}
-
-static void
-drawcn (const RowVector& X, const RowVector& Y, const Matrix& Z,
-        double lvl, int r, int c, double ct_x, double ct_y,
-        unsigned int start_edge, bool first, charMatrix& mark)
-{
-  double px[4], py[4], pz[4], tmp;
-  unsigned int stop_edge, next_edge, pt[2];
-  int next_r, next_c;
-
-  //get x, y, and z - lvl for current facet
-  px[0] = px[3] = X(c);
-  px[1] = px[2] = X(c+1);
-
-  py[0] = py[1] = Y(r);
-  py[2] = py[3] = Y(r+1);
-
-  pz[3] = Z(r+1, c) - lvl;
-  pz[2] = Z(r+1, c + 1) - lvl;
-  pz[1] = Z(r, c+1) - lvl;
-  pz[0] = Z(r, c) - lvl;
-
-  // Facet edge and point naming assignment.
-  //
-  //  0-----1   .-0-.
-  //  |     |   |   |
-  //  |     |   3   1
-  //  |     |   |   |
-  //  3-----2   .-2-.
-
-  // Get mark value of current facet.
-  char id = static_cast<char> (mark(r, c));
-
-  // Check startedge s.
-  if (start_edge == 255)
-    {
-      // Find start edge.
-      for (unsigned int k = 0; k < 4; k++)
-        if (static_cast<char> (1 << k) & id)
-          start_edge = k;
-    }
-
-  if (start_edge == 255)
-    return;
-
-  // Decrease mark value of current facet for start edge.
-  mark(r, c) -= static_cast<char> (1 << start_edge);
-
-  // Next point (clockwise).
-  pt[0] = start_edge;
-  pt[1] = (pt[0] + 1) % 4;
-
-  // Calculate contour segment start if first of contour.
-  if (first)
-    {
-      tmp = fabs (pz[pt[1]]) / fabs (pz[pt[0]]);
-
-      if (xisnan (tmp))
-        ct_x = ct_y = 0.5;
-      else
-        {
-          ct_x = px[pt[0]] + (px[pt[1]] - px[pt[0]])/(1 + tmp);
-          ct_y = py[pt[0]] + (py[pt[1]] - py[pt[0]])/(1 + tmp);
-        }
-
-      start_contour (lvl, ct_x, ct_y);
-    }
-
-  // Find stop edge.
-  // FIXME -- perhaps this should use a while loop?
-  for (unsigned int k = 1; k <= 4; k++)
-    {
-      if (start_edge == 0 || start_edge == 2)
-        stop_edge = (start_edge + k) % 4;
-      else
-        stop_edge = (start_edge - k) % 4;
-
-      if (static_cast<char> (1 << stop_edge) & id)
-        break;
-    }
-
-  pt[0] = stop_edge;
-  pt[1] = (pt[0] + 1) % 4;
-  tmp = fabs (pz[pt[1]]) / fabs (pz[pt[0]]);
-
-  if (xisnan (tmp))
-    ct_x = ct_y = 0.5;
-  else
-    {
-      ct_x = px[pt[0]] + (px[pt[1]] - px[pt[0]])/(1 + tmp);
-      ct_y = py[pt[0]] + (py[pt[1]] - py[pt[0]])/(1 + tmp);
-    }
-
-  // Add point to contour.
-  add_point (ct_x, ct_y);
-
-  // Decrease id value of current facet for start edge.
-  mark(r, c) -= static_cast<char> (1 << stop_edge);
-
-  // Find next facet.
-  next_c = c;
-  next_r = r;
-
-  if (stop_edge == 0)
-    next_r--;
-  else if (stop_edge == 1)
-    next_c++;
-  else if (stop_edge == 2)
-    next_r++;
-  else if (stop_edge == 3)
-    next_c--;
-
-  // Check if next facet is not done yet.
-  // Go to next facet.
-  if (next_r >= 0 && next_c >= 0 && next_r < mark.rows ()
-      && next_c < mark.cols () && mark(next_r, next_c) > 0)
-    {
-      next_edge = (stop_edge + 2) % 4;
-      drawcn (X, Y, Z, lvl, next_r, next_c, ct_x, ct_y, next_edge, false, mark);
-    }
-}
-
-static void
-mark_facets (const Matrix& Z, charMatrix& mark, double lvl)
-{
-  unsigned int nr = mark.rows ();
-  unsigned int nc = mark.cols ();
-
-  double f[4];
-
-  for (unsigned int c = 0; c < nc; c++)
-    for (unsigned int r = 0; r < nr; r++)
-      {
-        f[0] = Z(r, c) - lvl;
-        f[1] = Z(r, c+1) - lvl;
-        f[3] = Z(r+1, c) - lvl;
-        f[2] = Z(r+1, c+1) - lvl;
-
-        for (unsigned int i = 0; i < 4; i++)
-          if (fabs(f[i]) < DBL_EPSILON)
-            f[i] = DBL_EPSILON;
-
-        if (f[1] * f[2] < 0)
-          mark(r, c) += 2;
-
-        if (f[0] * f[3] < 0)
-          mark(r, c) += 8;
-      }
-
-  for (unsigned int r = 0; r < nr; r++)
-    for (unsigned int c = 0; c < nc; c++)
-      {
-        f[0] = Z(r, c) - lvl;
-        f[1] = Z(r, c+1) - lvl;
-        f[3] = Z(r+1, c) - lvl;
-        f[2] = Z(r+1, c+1) - lvl;
-
-        for (unsigned int i = 0; i < 4; i++)
-          if (fabs(f[i]) < DBL_EPSILON)
-            f[i] = DBL_EPSILON;
-
-        if (f[0] * f[1] < 0)
-          mark(r, c) += 1;
-
-        if (f[2] * f[3] < 0)
-          mark(r, c) += 4;
-      }
-}
-
-static void
-cntr (const RowVector& X, const RowVector& Y, const Matrix& Z, double lvl)
-{
-  unsigned int nr = Z.rows ();
-  unsigned int nc = Z.cols ();
-
-  charMatrix mark (nr - 1, nc - 1, 0);
-
-  mark_facets (Z, mark, lvl);
-
-  // Find contours that start at a domain edge.
-
-  for (unsigned int c = 0; c < nc - 1; c++)
-    {
-      // Top.
-      if (mark(0, c) & 1)
-        drawcn (X, Y, Z, lvl, 0, c, 0.0, 0.0, 0, true, mark);
-
-      // Bottom.
-      if (mark(nr - 2, c) & 4)
-        drawcn (X, Y, Z, lvl, nr - 2, c, 0.0, 0.0, 2, true, mark);
-    }
-
-  for (unsigned int r = 0; r < nr - 1; r++)
-    {
-      // Left.
-      if (mark(r, 0) & 8)
-        drawcn (X, Y, Z, lvl, r, 0, 0.0, 0.0, 3, true, mark);
-
-      // Right.
-      if (mark(r, nc - 2) & 2)
-        drawcn (X, Y, Z, lvl, r, nc - 2, 0.0, 0.0, 1, true, mark);
-    }
-
-  for (unsigned int r = 0; r < nr - 1; r++)
-    for (unsigned int c = 0; c < nc - 1; c++)
-      if (mark (r, c) > 0)
-        drawcn (X, Y, Z, lvl, r, c, 0.0, 0.0, 255, true, mark);
-}
-
-DEFUN_DLD (__contourc__, args, ,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {} __contourc__ (@var{x}, @var{y}, @var{z}, @var{levels})\n\
-Undocumented internal function.\n\
-@end deftypefn")
-{
-  octave_value retval;
-
-  if (args.length () == 4)
-    {
-      RowVector X = args (0).row_vector_value ();
-      RowVector Y = args (1).row_vector_value ();
-      Matrix Z = args (2).matrix_value ();
-      RowVector L = args (3).row_vector_value ();
-
-      if (! error_state)
-        {
-          contourc.resize (2, 0);
-
-          for (int i = 0; i < L.length (); i++)
-            cntr (X, Y, Z, L (i));
-
-          end_contour ();
-
-          retval = contourc;
-        }
-      else
-        error ("__contourc__: invalid argument values");
-    }
-  else
-    print_usage ();
-
-  return retval;
-}
-
-/*
-## No test needed for internal helper function.
-%!assert (1)
-*/

--- a/src/DLD-FUNCTIONS/__dispatch__.cc	Fri Jul 27 17:10:25 2012 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,136 +0,0 @@
-/*
-
-Copyright (C) 2001-2012 John W. Eaton and Paul Kienzle
-
-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 <list>
-#include <map>
-#include <string>
-
-#include "Cell.h"
-#include "oct-map.h"
-#include "defun-dld.h"
-#include "ov.h"
-#include "ov-fcn.h"
-#include "ov-typeinfo.h"
-#include "pager.h"
-#include "parse.h"
-#include "symtab.h"
-#include "variables.h"
-
-DEFUN_DLD (__dispatch__, args, nargout,
-  "Undocumented internal function")
-{
-  octave_value retval;
-
-  int nargin = args.length ();
-
-  std::string f, r, t;
-
-  if (nargin > 0 && nargin < 4)
-    {
-      if (nargin > 0)
-        {
-          f = args(0).string_value ();
-
-          if (error_state)
-            {
-              error ("__dispatch__: first argument must be a function name");
-              return retval;
-            }
-        }
-
-      if (nargin > 1)
-        {
-          r = args(1).string_value ();
-
-          if (error_state)
-            {
-              error ("__dispatch__: second argument must be a function name");
-              return retval;
-            }
-        }
-
-      if (nargin > 2)
-        {
-          t = args(2).string_value ();
-
-          if (error_state)
-            {
-              error ("__dispatch__: third argument must be a type name");
-              return retval;
-            }
-        }
-
-      if (nargin == 1)
-        {
-          if (nargout > 0)
-            {
-              symbol_table::fcn_info::dispatch_map_type dm
-                = symbol_table::get_dispatch (f);
-
-              size_t len = dm.size ();
-
-              Cell type_field (len, 1);
-              Cell name_field (len, 1);
-
-              symbol_table::fcn_info::dispatch_map_type::const_iterator p
-                = dm.begin ();
-
-              for (size_t i = 0; i < len; i++)
-                {
-                  type_field(i) = p->first;
-                  name_field(i) = p->second;
-
-                  p++;
-                }
-
-              octave_scalar_map m;
-
-              m.assign ("type", type_field);
-              m.assign ("name", name_field);
-
-              retval = m;
-            }
-          else
-            symbol_table::print_dispatch (octave_stdout, f);
-        }
-      else if (nargin == 2)
-        {
-          t = r;
-          symbol_table::clear_dispatch (f, t);
-        }
-      else
-        symbol_table::add_dispatch (f, t, r);
-    }
-  else
-    print_usage ();
-
-  return retval;
-}
-
-/*
-## No test needed for internal helper function.
-%!assert (1)
-*/

--- a/src/DLD-FUNCTIONS/__lin_interpn__.cc	Fri Jul 27 17:10:25 2012 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,362 +0,0 @@
-/*
-
-Copyright (C) 2007-2012 Alexander Barth
-
-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 "lo-ieee.h"
-#include "dNDArray.h"
-#include "oct-locbuf.h"
-
-#include "defun-dld.h"
-#include "error.h"
-#include "oct-obj.h"
-
-// equivalent to isvector.m
-
-template <class T>
-bool
-isvector (const T& array)
-{
-  const dim_vector dv = array.dims ();
-  return dv.length () == 2 && (dv(0) == 1 || dv(1) == 1);
-}
-
-// lookup a value in a sorted table (lookup.m)
-template <class T>
-octave_idx_type
-lookup (const T *x, octave_idx_type n, T y)
-{
-  octave_idx_type j;
-
-  if (x[0] < x[n-1])
-    {
-      // increasing x
-
-      if (y > x[n-1] || y < x[0])
-        return -1;
-
-#ifdef EXHAUSTIF
-      for (j = 0; j < n - 1; j++)
-        {
-          if (x[j] <= y && y <= x[j+1])
-            return j;
-        }
-#else
-      octave_idx_type j0 = 0;
-      octave_idx_type j1 = n - 1;
-
-      while (true)
-        {
-          j = (j0+j1)/2;
-
-          if (y <= x[j+1])
-            {
-              if (x[j] <= y)
-                return j;
-
-              j1 = j;
-            }
-
-          if (x[j] <= y)
-            j0 = j;
-        }
-#endif
-    }
-  else
-    {
-      // decreasing x
-      // previous code with x -> -x and y -> -y
-
-      if (y > x[0] || y < x[n-1])
-        return -1;
-
-#ifdef EXHAUSTIF
-      for (j = 0; j < n - 1; j++)
-        {
-          if (x[j+1] <= y && y <= x[j])
-            return j;
-        }
-#else
-      octave_idx_type j0 = 0;
-      octave_idx_type j1 = n - 1;
-
-      while (true)
-        {
-          j = (j0+j1)/2;
-
-          if (y >= x[j+1])
-            {
-              if (x[j] >= y)
-                return j;
-
-              j1 = j;
-            }
-
-          if (x[j] >= y)
-            j0 = j;
-        }
-#endif
-    }
-}
-
-// n-dimensional linear interpolation
-
-template <class T>
-void
-lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale,
-             octave_idx_type Ni, T extrapval, const T **x,
-             const T *v, const T **y, T *vi)
-{
-  bool out = false;
-  int bit;
-
-  OCTAVE_LOCAL_BUFFER (T, coef, 2*n);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, index, n);
-
-  // loop over all points
-  for (octave_idx_type m = 0; m < Ni; m++)
-    {
-      // loop over all dimensions
-      for (int i = 0; i < n; i++)
-        {
-          index[i] = lookup (x[i], size[i], y[i][m]);
-          out = index[i] == -1;
-
-          if (out)
-            break;
-          else
-            {
-              octave_idx_type j = index[i];
-              coef[2*i+1] = (y[i][m] - x[i][j])/(x[i][j+1] - x[i][j]);
-              coef[2*i] = 1 - coef[2*i+1];
-            }
-        }
-
-
-      if (out)
-        vi[m] = extrapval;
-      else
-        {
-          vi[m] = 0;
-
-          // loop over all corners of hypercube (1<<n = 2^n)
-          for (int i = 0; i < (1 << n); i++)
-            {
-              T c = 1;
-              octave_idx_type l = 0;
-
-              // loop over all dimensions
-              for (int j = 0; j < n; j++)
-                {
-                  // test if the jth bit in i is set
-                  bit = i >> j & 1;
-                  l += scale[j] * (index[j] + bit);
-                  c *= coef[2*j+bit];
-                }
-
-              vi[m] += c * v[l];
-            }
-        }
-    }
-}
-
-template <class T, class M>
-octave_value
-lin_interpn (int n, M *X, const M V, M *Y)
-{
-  octave_value retval;
-
-  M Vi = M (Y[0].dims ());
-
-  OCTAVE_LOCAL_BUFFER (const T *, y, n);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, size, n);
-
-  for (int i = 0; i < n; i++)
-    {
-      y[i] = Y[i].data ();
-      size[i] =  V.dims ()(i);
-    }
-
-  OCTAVE_LOCAL_BUFFER (const T *, x, n);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, scale, n);
-
-  const T *v = V.data ();
-  T *vi = Vi.fortran_vec ();
-  octave_idx_type Ni = Vi.numel ();
-
-  T extrapval = octave_NA;
-
-  // offset in memory of each dimension
-
-  scale[0] = 1;
-
-  for (int i = 1; i < n; i++)
-    scale[i] = scale[i-1] * size[i-1];
-
-  // tests if X[0] is a vector, if yes, assume that all elements of X are
-  // in the ndgrid format.
-
-  if (! isvector (X[0]))
-    {
-      for (int i = 0; i < n; i++)
-        {
-          if (X[i].dims () != V.dims ())
-            {
-              error ("interpn: incompatible size of argument number %d", i+1);
-              return retval;
-            }
-          else
-            {
-              M tmp = M (dim_vector (size[i], 1));
-
-              for (octave_idx_type j = 0; j < size[i]; j++)
-                tmp(j) =  X[i](scale[i]*j);
-
-              X[i] = tmp;
-            }
-        }
-    }
-
-  for (int i = 0; i < n; i++)
-    {
-      if (! isvector (X[i]) && X[i].numel () != size[i])
-        {
-          error ("interpn: incompatible size of argument number %d", i+1);
-          return retval;
-        }
-      else
-        x[i] = X[i].data ();
-    }
-
-  lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi);
-
-  retval = Vi;
-
-  return retval;
-}
-
-// Perform @var{n}-dimensional interpolation.  Each element of then
-// @var{n}-dimensional array @var{v} represents a value at a location
-// given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters
-// @var{x1}, @var{x2}, @dots{}, @var{xn} are either @var{n}-dimensional
-// arrays of the same size as the array @var{v} in the \"ndgrid\" format
-// or vectors.  The parameters @var{y1}, @var{y2}, @dots{}, @var{yn} are
-// all @var{n}-dimensional arrays of the same size and represent the
-// points at which the array @var{vi} is interpolated.
-//
-//This function only performs linear interpolation.
-
-DEFUN_DLD (__lin_interpn__, args, ,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {@var{vi} =} __lin_interpn__ (@var{x1}, @var{x2}, @dots{}, @var{xn}, @var{v}, @var{y1}, @var{y2}, @dots{}, @var{yn})\n\
-Undocumented internal function.\n\
-@end deftypefn")
-{
-  octave_value retval;
-
-  int nargin = args.length ();
-
-  if (nargin < 2 ||  nargin % 2 == 0)
-    {
-      print_usage ();
-      return retval;
-    }
-
-  // dimension of the problem
-  int n = (nargin-1)/2;
-
-  if (args(n).is_single_type ())
-    {
-      OCTAVE_LOCAL_BUFFER (FloatNDArray, X, n);
-      OCTAVE_LOCAL_BUFFER (FloatNDArray, Y, n);
-
-      const FloatNDArray V = args(n).float_array_value ();
-
-      if (error_state)
-        {
-          print_usage ();
-          return retval;
-        }
-
-      for (int i = 0; i < n; i++)
-        {
-          X[i] = args(i).float_array_value ();
-          Y[i] = args(n+i+1).float_array_value ();
-
-          if (error_state)
-            {
-              print_usage ();
-              return retval;
-            }
-
-          if (Y[0].dims () != Y[i].dims ())
-            {
-              error ("interpn: incompatible size of argument number %d", n+i+2);
-              return retval;
-            }
-        }
-
-      retval = lin_interpn<float, FloatNDArray> (n, X, V, Y);
-    }
-  else
-    {
-      OCTAVE_LOCAL_BUFFER (NDArray, X, n);
-      OCTAVE_LOCAL_BUFFER (NDArray, Y, n);
-
-      const NDArray V = args(n).array_value ();
-
-      if (error_state)
-        {
-          print_usage ();
-          return retval;
-        }
-
-      for (int i = 0; i < n; i++)
-        {
-          X[i] = args(i).array_value ();
-          Y[i] = args(n+i+1).array_value ();
-
-          if (error_state)
-            {
-              print_usage ();
-              return retval;
-            }
-
-          if (Y[0].dims () != Y[i].dims ())
-            {
-              error ("interpn: incompatible size of argument number %d", n+i+2);
-              return retval;
-            }
-        }
-
-      retval = lin_interpn<double, NDArray> (n, X, V, Y);
-    }
-
-  return retval;
-}
-
-/*
-## No test needed for internal helper function.
-%!assert (1)
-*/

--- a/src/DLD-FUNCTIONS/__pchip_deriv__.cc	Fri Jul 27 17:10:25 2012 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,168 +0,0 @@
-/*
-
-Copyright (C) 2002-2012 Kai Habel
-Copyright (C) 2008-2009 Jaroslav Hajek
-
-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 "defun-dld.h"
-#include "error.h"
-#include "gripes.h"
-#include "oct-obj.h"
-#include "utils.h"
-#include "f77-fcn.h"
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (dpchim, DPCHIM) (const octave_idx_type& n, const double *x,
-                             const double *f, double *d,
-                             const octave_idx_type &incfd,
-                             octave_idx_type *ierr);
-
-  F77_RET_T
-  F77_FUNC (pchim, PCHIM) (const octave_idx_type& n, const float *x,
-                           const float *f, float *d,
-                           const octave_idx_type& incfd,
-                           octave_idx_type *ierr);
-}
-
-// Wrapper for SLATEC/PCHIP function DPCHIM to calculate the derivates
-// for piecewise polynomials.
-
-DEFUN_DLD (__pchip_deriv__, args, ,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {} __pchip_deriv__ (@var{x}, @var{y}, @var{dim})\n\
-Undocumented internal function.\n\
-@end deftypefn")
-{
-  octave_value retval;
-  const int nargin = args.length ();
-
-  bool rows = (nargin == 3 && args (2).uint_value () == 2);
-
-  if (nargin >= 2)
-    {
-      if (args(0).is_single_type () || args(1).is_single_type ())
-        {
-          FloatColumnVector xvec (args(0).float_vector_value ());
-          FloatMatrix ymat (args(1).float_matrix_value ());
-
-          octave_idx_type nx = xvec.length ();
-
-          if (nx < 2)
-            {
-              error ("__pchip_deriv__: X must be at least of length 2");
-              return retval;
-            }
-
-          octave_idx_type nyr = ymat.rows ();
-          octave_idx_type nyc = ymat.columns ();
-
-          if (nx != (rows ? nyc : nyr))
-            {
-              error ("__pchip_deriv__: X and Y dimension mismatch");
-              return retval;
-            }
-
-          const float *yvec = ymat.data ();
-          FloatMatrix dmat (nyr, nyc);
-          float *dvec = dmat.fortran_vec ();
-
-          octave_idx_type ierr;
-          const octave_idx_type incfd = rows ? nyr : 1;
-          const octave_idx_type inc = rows ? 1 : nyr;
-
-          for (octave_idx_type i = (rows ? nyr : nyc); i > 0; i--)
-            {
-              F77_FUNC (pchim, PCHIM) (nx, xvec.data (),
-                                       yvec, dvec, incfd, &ierr);
-
-              yvec += inc;
-              dvec += inc;
-
-              if (ierr < 0)
-                {
-                  error ("PCHIM: error: %i\n", ierr);
-                  return retval;
-                }
-            }
-
-          retval = dmat;
-        }
-      else
-        {
-          ColumnVector xvec (args(0).vector_value ());
-          Matrix ymat (args(1).matrix_value ());
-
-          octave_idx_type nx = xvec.length ();
-
-          if (nx < 2)
-            {
-              error ("__pchip_deriv__: X must be at least of length 2");
-              return retval;
-            }
-
-          octave_idx_type nyr = ymat.rows ();
-          octave_idx_type nyc = ymat.columns ();
-
-          if (nx != (rows ? nyc : nyr))
-            {
-              error ("__pchip_deriv__: X and Y dimension mismatch");
-              return retval;
-            }
-
-          const double *yvec = ymat.data ();
-          Matrix dmat (nyr, nyc);
-          double *dvec = dmat.fortran_vec ();
-
-          octave_idx_type ierr;
-          const octave_idx_type incfd = rows ? nyr : 1;
-          const octave_idx_type inc = rows ? 1 : nyr;
-
-          for (octave_idx_type i = (rows ? nyr : nyc); i > 0; i--)
-            {
-              F77_FUNC (dpchim, DPCHIM) (nx, xvec.data (),
-                                         yvec, dvec, incfd, &ierr);
-
-              yvec += inc;
-              dvec += inc;
-
-              if (ierr < 0)
-                {
-                  error ("DPCHIM: error: %i\n", ierr);
-                  return retval;
-                }
-            }
-
-          retval = dmat;
-        }
-    }
-
-  return retval;
-}
-
-/*
-## No test needed for internal helper function.
-%!assert (1)
-*/

--- a/src/DLD-FUNCTIONS/__qp__.cc	Fri Jul 27 17:10:25 2012 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,535 +0,0 @@
-/*
-
-Copyright (C) 2000-2012 Gabriele Pannocchia
-
-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 <cfloat>
-
-#include "dbleCHOL.h"
-#include "dbleSVD.h"
-#include "mx-m-dm.h"
-#include "EIG.h"
-
-#include "defun-dld.h"
-#include "error.h"
-#include "gripes.h"
-#include "oct-obj.h"
-#include "pr-output.h"
-#include "utils.h"
-
-static Matrix
-null (const Matrix& A, octave_idx_type& rank)
-{
-  Matrix retval;
-
-  rank = 0;
-
-  if (! A.is_empty ())
-    {
-      SVD A_svd (A);
-
-      DiagMatrix S = A_svd.singular_values ();
-
-      ColumnVector s = S.diag ();
-
-      Matrix V = A_svd.right_singular_matrix ();
-
-      octave_idx_type A_nr = A.rows ();
-      octave_idx_type A_nc = A.cols ();
-
-      octave_idx_type tmp = A_nr > A_nc ? A_nr : A_nc;
-
-      double tol = tmp * s(0) * DBL_EPSILON;
-
-      octave_idx_type n = s.length ();
-
-      for (octave_idx_type i = 0; i < n; i++)
-        {
-          if (s(i) > tol)
-            rank++;
-        }
-
-      if (rank < A_nc)
-        retval = V.extract (0, rank, A_nc-1, A_nc-1);
-      else
-        retval.resize (A_nc, 0);
-
-      for (octave_idx_type i = 0; i < retval.numel (); i++)
-        if (std::abs (retval(i)) < DBL_EPSILON)
-          retval(i) = 0;
-    }
-
-  return retval;
-}
-
-static int
-qp (const Matrix& H, const ColumnVector& q,
-    const Matrix& Aeq, const ColumnVector& beq,
-    const Matrix& Ain, const ColumnVector& bin,
-    int maxit,
-    ColumnVector& x, ColumnVector& lambda, int& iter)
-{
-  int info = 0;
-
-  iter = 0;
-
-  double rtol = sqrt (DBL_EPSILON);
-
-  // Problem dimension.
-  octave_idx_type n = x.length ();
-
-  // Dimension of constraints.
-  octave_idx_type n_eq = beq.length ();
-  octave_idx_type n_in = bin.length ();
-
-  // Filling the current active set.
-
-  octave_idx_type n_act = n_eq;
-
-  octave_idx_type n_tot = n_eq + n_in;
-
-  // Equality constraints come first.  We won't check the sign of the
-  // Lagrange multiplier for those.
-
-  Matrix Aact = Aeq;
-  ColumnVector bact = beq;
-  ColumnVector Wact;
-
-  if (n_in > 0)
-    {
-      ColumnVector res = Ain*x - bin;
-
-      for (octave_idx_type i = 0; i < n_in; i++)
-        {
-          res(i) /= (1.0 + std::abs (bin(i)));
-
-          if (res(i) < rtol)
-            {
-              n_act++;
-              Aact = Aact.stack (Ain.row (i));
-              bact.resize (n_act, bin(i));
-              Wact.resize (n_act-n_eq, i);
-            }
-        }
-    }
-
-  // Computing the ???
-
-  EIG eigH (H);
-
-  if (error_state)
-    {
-      error ("qp: failed to compute eigenvalues of H");
-      return -1;
-    }
-
-  ColumnVector eigenvalH = real (eigH.eigenvalues ());
-  Matrix eigenvecH = real (eigH.eigenvectors ());
-  double minReal = eigenvalH.min ();
-  octave_idx_type indminR = 0;
-  for (octave_idx_type i = 0; i < n; i++)
-    {
-      if (minReal == eigenvalH(i))
-        {
-          indminR = i;
-          break;
-        }
-    }
-
-  bool done = false;
-
-  double alpha = 0.0;
-
-  Matrix R;
-  Matrix Y (n, 0, 0.0);
-
-  ColumnVector g (n, 0.0);
-  ColumnVector p (n, 0.0);
-
-  ColumnVector lambda_tmp (n_in, 0.0);
-
-  while (! done)
-    {
-      iter++;
-
-      // Current Gradient
-      // g = q + H * x;
-
-      g = q + H * x;
-
-      if (n_act == 0)
-        {
-          // There are no active constraints.
-
-          if (minReal > 0.0)
-            {
-              // Inverting the Hessian.  Using the Cholesky
-              // factorization since the Hessian is positive
-              // definite.
-
-              CHOL cholH (H);
-
-              R = cholH.chol_matrix ();
-
-              Matrix Hinv = chol2inv (R);
-
-              // Computing the unconstrained step.
-              // p = -Hinv * g;
-
-              p = -Hinv * g;
-
-              info = 0;
-            }
-          else
-            {
-              // Finding the negative curvature of H.
-
-              p = eigenvecH.column (indminR);
-
-              // Following the negative curvature of H.
-
-              if (p.transpose () * g > DBL_EPSILON)
-                p = -p;
-
-              info = 1;
-            }
-
-          // Multipliers are zero.
-          lambda_tmp.fill (0.0);
-        }
-      else
-        {
-          // There are active constraints.
-
-          // Computing the null space.
-
-          octave_idx_type rank;
-
-          Matrix Z = null (Aact, rank);
-
-          octave_idx_type dimZ = n - rank;
-
-          // FIXME -- still remain to handle the case of
-          // non-full rank active set matrix.
-
-          // Computing the Y matrix (orthogonal to Z)
-          Y = Aact.pseudo_inverse ();
-
-          // Reduced Hessian
-          Matrix Zt = Z.transpose ();
-          Matrix rH = Zt * H * Z;
-
-          octave_idx_type pR = 0;
-
-          if (dimZ > 0)
-            {
-              // Computing the Cholesky factorization (pR = 0 means
-              // that the reduced Hessian was positive definite).
-
-              CHOL cholrH (rH, pR);
-              Matrix tR = cholrH.chol_matrix ();
-              if (pR == 0)
-                R = tR;
-            }
-
-          if (pR == 0)
-            {
-              info = 0;
-
-              // Computing the step pz.
-              if (dimZ > 0)
-                {
-                  // Using the Cholesky factorization to invert rH
-
-                  Matrix rHinv = chol2inv (R);
-
-                  ColumnVector pz = -rHinv * Zt * g;
-
-                  // Global step.
-                  p = Z * pz;
-                }
-              else
-                {
-                  // Global step.
-                  p.fill (0.0);
-                }
-            }
-          else
-            {
-              info = 1;
-
-              // Searching for the most negative curvature.
-
-              EIG eigrH (rH);
-
-              if (error_state)
-                {
-                  error ("qp: failed to compute eigenvalues of rH");
-                  return -1;
-                }
-
-              ColumnVector eigenvalrH = real (eigrH.eigenvalues ());
-              Matrix eigenvecrH = real (eigrH.eigenvectors ());
-              double mRrH = eigenvalrH.min ();
-              indminR = 0;
-              for (octave_idx_type i = 0; i < n; i++)
-                {
-                  if (mRrH == eigenvalH(i))
-                    {
-                      indminR = i;
-                      break;
-                    }
-                }
-
-              ColumnVector eVrH = eigenvecrH.column (indminR);
-
-              // Computing the step pz.
-              p = Z * eVrH;
-
-              if (p.transpose () * g > DBL_EPSILON)
-                p = -p;
-            }
-        }
-
-      // Checking the step-size.
-      ColumnVector abs_p (n);
-      for (octave_idx_type i = 0; i < n; i++)
-        abs_p(i) = std::abs (p(i));
-      double max_p = abs_p.max ();
-
-      if (max_p < rtol)
-        {
-          // The step is null.  Checking constraints.
-          if (n_act - n_eq == 0)
-            // Solution is found because no inequality
-            // constraints are active.
-            done = true;
-          else
-            {
-              // Computing the multipliers only for the inequality
-              // constraints that are active.  We do NOT compute
-              // multipliers for the equality constraints.
-              Matrix Yt = Y.transpose ();
-              Yt = Yt.extract_n (n_eq, 0, n_act-n_eq, n);
-              lambda_tmp = Yt * (g + H * p);
-
-              // Checking the multipliers.  We remove the most
-              // negative from the set (if any).
-              double min_lambda = lambda_tmp.min ();
-              if (min_lambda >= 0)
-                {
-                  // Solution is found.
-                  done = true;
-                }
-              else
-                {
-                  octave_idx_type which_eig = 0;
-                  for (octave_idx_type i = 0; i < n_act; i++)
-                    {
-                      if (lambda_tmp(i) == min_lambda)
-                        {
-                          which_eig = i;
-                          break;
-                        }
-                    }
-
-                  // At least one multiplier is negative, we
-                  // remove it from the set.
-
-                  n_act--;
-                  for (octave_idx_type i = which_eig; i < n_act - n_eq; i++)
-                    {
-                      Wact(i) = Wact(i+1);
-                      for (octave_idx_type j = 0; j < n; j++)
-                        Aact(n_eq+i,j) = Aact(n_eq+i+1,j);
-                      bact(n_eq+i) = bact(n_eq+i+1);
-                    }
-
-                  // Resizing the active set.
-                  Wact.resize (n_act-n_eq);
-                  bact.resize (n_act);
-                  Aact.resize (n_act, n);
-                }
-            }
-        }
-      else
-        {
-          // The step is not null.
-          if (n_act - n_eq == n_in)
-            {
-              // All inequality constraints were active.  We can
-              // add the whole step.
-              x += p;
-            }
-          else
-            {
-              // Some constraints were not active.  Checking if
-              // there is a blocking constraint.
-              alpha = 1.0;
-              octave_idx_type is_block = -1;
-
-              for (octave_idx_type i = 0; i < n_in; i++)
-                {
-                  bool found = false;
-
-                  for (octave_idx_type j = 0; j < n_act-n_eq; j++)
-                    {
-                      if (Wact(j) == i)
-                        {
-                          found = true;
-                          break;
-                        }
-                    }
-
-                  if (! found)
-                    {
-                      // The i-th constraint was not in the set.  Is it a
-                      // blocking constraint?
-
-                      RowVector tmp_row = Ain.row (i);
-                      double tmp = tmp_row * p;
-                      double res = tmp_row * x;
-
-                      if (tmp < 0.0)
-                        {
-                          double alpha_tmp = (bin(i) - res) / tmp;
-
-                          if (alpha_tmp < alpha)
-                            {
-                              alpha = alpha_tmp;
-                              is_block = i;
-                            }
-                        }
-                    }
-                }
-
-              // In is_block there is the index of the blocking
-              // constraint (if any).
-              if (is_block >= 0)
-                {
-                  // There is a blocking constraint (index in
-                  // is_block) which is added to the active set.
-                  n_act++;
-                  Aact = Aact.stack (Ain.row (is_block));
-                  bact.resize (n_act, bin(is_block));
-                  Wact.resize (n_act-n_eq, is_block);
-
-                  // Adding the reduced step
-                  x += alpha * p;
-                }
-              else
-                {
-                  // There are no blocking constraints.  Adding the
-                  // whole step.
-                  x += alpha * p;
-                }
-            }
-        }
-
-      if (iter == maxit)
-        {
-          done = true;
-          // warning ("qp_main: maximum number of iteration reached");
-          info = 3;
-        }
-    }
-
-  lambda_tmp = Y.transpose () * (g + H * p);
-
-  // Reordering the Lagrange multipliers.
-
-  lambda.resize (n_tot);
-  lambda.fill (0.0);
-  for (octave_idx_type i = 0; i < n_eq; i++)
-    lambda(i) = lambda_tmp(i);
-
-  for (octave_idx_type i = n_eq; i < n_tot; i++)
-    {
-      for (octave_idx_type j = 0; j < n_act-n_eq; j++)
-        {
-          if (Wact(j) == i - n_eq)
-            {
-              lambda(i) = lambda_tmp(n_eq+j);
-              break;
-            }
-        }
-    }
-
-  return info;
-}
-
-DEFUN_DLD (__qp__, args, ,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {[@var{x}, @var{lambda}, @var{info}, @var{iter}] =} __qp__ (@var{x0}, @var{H}, @var{q}, @var{Aeq}, @var{beq}, @var{Ain}, @var{bin}, @var{maxit})\n\
-Undocumented internal function.\n\
-@end deftypefn")
-{
-  octave_value_list retval;
-
-  if (args.length () == 8)
-    {
-      const ColumnVector x0  (args(0) . vector_value ());
-      const Matrix H         (args(1) . matrix_value ());
-      const ColumnVector q   (args(2) . vector_value ());
-      const Matrix Aeq       (args(3) . matrix_value ());
-      const ColumnVector beq (args(4) . vector_value ());
-      const Matrix Ain       (args(5) . matrix_value ());
-      const ColumnVector bin (args(6) . vector_value ());
-      const int maxit        (args(7) . int_value ());
-
-      if (! error_state)
-        {
-          int iter = 0;
-
-          // Copying the initial guess in the working variable
-          ColumnVector x = x0;
-
-          // Reordering the Lagrange multipliers
-          ColumnVector lambda;
-
-          int info = qp (H, q, Aeq, beq, Ain, bin, maxit, x, lambda, iter);
-
-          if (! error_state)
-            {
-              retval(3) = iter;
-              retval(2) = info;
-              retval(1) = lambda;
-              retval(0) = x;
-            }
-          else
-            error ("qp: internal error");
-        }
-      else
-        error ("__qp__: invalid arguments");
-    }
-  else
-    print_usage ();
-
-  return retval;
-}
-
-/*
-## No test needed for internal helper function.
-%!assert (1)
-*/

--- a/src/DLD-FUNCTIONS/balance.cc	Fri Jul 27 17:10:25 2012 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,389 +0,0 @@
-/*
-
-Copyright (C) 1996-2012 John W. Eaton
-Copyright (C) 2008-2009 Jaroslav Hajek
-
-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/>.
-
-*/
-
-// Author: A. S. Hodel <scotte@eng.auburn.edu>
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <string>
-
-#include "CmplxAEPBAL.h"
-#include "fCmplxAEPBAL.h"
-#include "dbleAEPBAL.h"
-#include "floatAEPBAL.h"
-#include "CmplxGEPBAL.h"
-#include "fCmplxGEPBAL.h"
-#include "dbleGEPBAL.h"
-#include "floatGEPBAL.h"
-#include "quit.h"
-
-#include "defun-dld.h"
-#include "error.h"
-#include "f77-fcn.h"
-#include "gripes.h"
-#include "oct-obj.h"
-#include "utils.h"
-
-DEFUN_DLD (balance, args, nargout,
-  "-*- texinfo -*-\n\
-@deftypefn  {Loadable Function} {@var{AA} =} balance (@var{A})\n\
-@deftypefnx {Loadable Function} {@var{AA} =} balance (@var{A}, @var{opt})\n\
-@deftypefnx {Loadable Function} {[@var{DD}, @var{AA}] =} balance (@var{A}, @var{opt})\n\
-@deftypefnx {Loadable Function} {[@var{D}, @var{P}, @var{AA}] =} balance (@var{A}, @var{opt})\n\
-@deftypefnx {Loadable Function} {[@var{CC}, @var{DD}, @var{AA}, @var{BB}] =} balance (@var{A}, @var{B}, @var{opt})\n\
-\n\
-Compute @code{@var{AA} = @var{DD} \\ @var{A} * @var{DD}} in which @var{AA}\n\
-is a matrix whose row and column norms are roughly equal in magnitude, and\n\
-@code{@var{DD} = @var{P} * @var{D}}, in which @var{P} is a permutation\n\
-matrix and @var{D} is a diagonal matrix of powers of two.  This allows the\n\
-equilibration to be computed without round-off.  Results of eigenvalue\n\
-calculation are typically improved by balancing first.\n\
-\n\
-If two output values are requested, @code{balance} returns\n\
-the diagonal @var{D} and the permutation @var{P} separately as vectors.\n\
-In this case, @code{@var{DD} = eye(n)(:,@var{P}) * diag (@var{D})}, where\n\
-@math{n} is the matrix size.\n\
-\n\
-If four output values are requested, compute @code{@var{AA} =\n\
-@var{CC}*@var{A}*@var{DD}} and @code{@var{BB} = @var{CC}*@var{B}*@var{DD}},\n\
-in which @var{AA} and @var{BB} have non-zero elements of approximately the\n\
-same magnitude and @var{CC} and @var{DD} are permuted diagonal matrices as\n\
-in @var{DD} for the algebraic eigenvalue problem.\n\
-\n\
-The eigenvalue balancing option @var{opt} may be one of:\n\
-\n\
-@table @asis\n\
-@item \"noperm\", \"S\"\n\
-Scale only; do not permute.\n\
-\n\
-@item \"noscal\", \"P\"\n\
-Permute only; do not scale.\n\
-@end table\n\
-\n\
-Algebraic eigenvalue balancing uses standard @sc{lapack} routines.\n\
-\n\
-Generalized eigenvalue problem balancing uses Ward's algorithm\n\
-(SIAM Journal on Scientific and Statistical Computing, 1981).\n\
-@end deftypefn")
-{
-  octave_value_list retval;
-
-  int nargin = args.length ();
-
-  if (nargin < 1 || nargin > 3 || nargout < 0 || nargout > 4)
-    {
-      print_usage ();
-      return retval;
-    }
-
-  // determine if it's AEP or GEP
-  bool AEPcase = nargin == 1 || args(1).is_string ();
-
-  // problem dimension
-  octave_idx_type nn = args(0).rows ();
-
-  if (nn != args(0).columns ())
-    {
-      gripe_square_matrix_required ("balance");
-      return retval;
-    }
-
-  bool isfloat = args(0).is_single_type () ||
-    (! AEPcase && args(1).is_single_type ());
-
-  bool complex_case = (args(0).is_complex_type () ||
-                       (! AEPcase && args(1).is_complex_type ()));
-
-  // Extract argument 1 parameter for both AEP and GEP.
-  Matrix aa;
-  ComplexMatrix caa;
-  FloatMatrix faa;
-  FloatComplexMatrix fcaa;
-
-  if (isfloat)
-    {
-      if (complex_case)
-        fcaa = args(0).float_complex_matrix_value ();
-      else
-        faa = args(0).float_matrix_value ();
-    }
-  else
-    {
-      if (complex_case)
-        caa = args(0).complex_matrix_value ();
-      else
-        aa = args(0).matrix_value ();
-    }
-
-  if (error_state)
-    return retval;
-
-  // Treat AEP/GEP cases.
-  if (AEPcase)
-    {
-      // Algebraic eigenvalue problem.
-      bool noperm = false, noscal = false;
-      if (nargin > 1)
-        {
-          std::string a1s = args(1).string_value ();
-          noperm = a1s == "noperm" || a1s == "S";
-          noscal = a1s == "noscal" || a1s == "P";
-        }
-
-      // balance the AEP
-      if (isfloat)
-        {
-          if (complex_case)
-            {
-              FloatComplexAEPBALANCE result (fcaa, noperm, noscal);
-
-              if (nargout == 0 || nargout == 1)
-                retval(0) = result.balanced_matrix ();
-              else if (nargout == 2)
-                {
-                  retval(1) = result.balanced_matrix ();
-                  retval(0) = result.balancing_matrix ();
-                }
-              else
-                {
-                  retval(2) = result.balanced_matrix ();
-                  retval(1) = result.permuting_vector ();
-                  retval(0) = result.scaling_vector ();
-                }
-
-            }
-          else
-            {
-              FloatAEPBALANCE result (faa, noperm, noscal);
-
-              if (nargout == 0 || nargout == 1)
-                retval(0) = result.balanced_matrix ();
-              else if (nargout == 2)
-                {
-                  retval(1) = result.balanced_matrix ();
-                  retval(0) = result.balancing_matrix ();
-                }
-              else
-                {
-                  retval(2) = result.balanced_matrix ();
-                  retval(1) = result.permuting_vector ();
-                  retval(0) = result.scaling_vector ();
-                }
-            }
-        }
-      else
-        {
-          if (complex_case)
-            {
-              ComplexAEPBALANCE result (caa, noperm, noscal);
-
-              if (nargout == 0 || nargout == 1)
-                retval(0) = result.balanced_matrix ();
-              else if (nargout == 2)
-                {
-                  retval(1) = result.balanced_matrix ();
-                  retval(0) = result.balancing_matrix ();
-                }
-              else
-                {
-                  retval(2) = result.balanced_matrix ();
-                  retval(1) = result.permuting_vector ();
-                  retval(0) = result.scaling_vector ();
-                }
-            }
-          else
-            {
-              AEPBALANCE result (aa, noperm, noscal);
-
-              if (nargout == 0 || nargout == 1)
-                retval(0) = result.balanced_matrix ();
-              else if (nargout == 2)
-                {
-                  retval(1) = result.balanced_matrix ();
-                  retval(0) = result.balancing_matrix ();
-                }
-              else
-                {
-                  retval(2) = result.balanced_matrix ();
-                  retval(1) = result.permuting_vector ();
-                  retval(0) = result.scaling_vector ();
-                }
-            }
-        }
-    }
-  else
-    {
-      std::string bal_job;
-      if (nargout == 1)
-        warning ("balance: used GEP, should have two output arguments");
-
-      // Generalized eigenvalue problem.
-      if (nargin == 2)
-        bal_job = "B";
-      else if (args(2).is_string ())
-        bal_job = args(2).string_value ();
-      else
-        {
-          error ("balance: OPT argument must be a string");
-          return retval;
-        }
-
-      if ((nn != args(1).columns ()) || (nn != args(1).rows ()))
-        {
-          gripe_nonconformant ();
-          return retval;
-        }
-
-      Matrix bb;
-      ComplexMatrix cbb;
-      FloatMatrix fbb;
-      FloatComplexMatrix fcbb;
-
-      if (isfloat)
-        {
-          if (complex_case)
-            fcbb = args(1).float_complex_matrix_value ();
-          else
-            fbb = args(1).float_matrix_value ();
-        }
-      else
-        {
-          if (complex_case)
-            cbb = args(1).complex_matrix_value ();
-          else
-            bb = args(1).matrix_value ();
-        }
-
-      // balance the GEP
-      if (isfloat)
-        {
-          if (complex_case)
-            {
-              FloatComplexGEPBALANCE result (fcaa, fcbb, bal_job);
-
-              switch (nargout)
-                {
-                case 4:
-                  retval(3) = result.balanced_matrix2 ();
-                  // fall through
-                case 3:
-                  retval(2) = result.balanced_matrix ();
-                  retval(1) = result.balancing_matrix2 ();
-                  retval(0) = result.balancing_matrix ();
-                  break;
-                case 2:
-                  retval(1) = result.balancing_matrix2 ();
-                  // fall through
-                case 1:
-                  retval(0) = result.balancing_matrix ();
-                  break;
-                default:
-                  error ("balance: invalid number of output arguments");
-                  break;
-                }
-            }
-          else
-            {
-              FloatGEPBALANCE result (faa, fbb, bal_job);
-
-              switch (nargout)
-                {
-                case 4:
-                  retval(3) = result.balanced_matrix2 ();
-                  // fall through
-                case 3:
-                  retval(2) = result.balanced_matrix ();
-                  retval(1) = result.balancing_matrix2 ();
-                  retval(0) = result.balancing_matrix ();
-                  break;
-                case 2:
-                  retval(1) = result.balancing_matrix2 ();
-                  // fall through
-                case 1:
-                  retval(0) = result.balancing_matrix ();
-                  break;
-                default:
-                  error ("balance: invalid number of output arguments");
-                  break;
-                }
-            }
-        }
-      else
-        {
-          if (complex_case)
-            {
-              ComplexGEPBALANCE result (caa, cbb, bal_job);
-
-              switch (nargout)
-                {
-                case 4:
-                  retval(3) = result.balanced_matrix2 ();
-                  // fall through
-                case 3:
-                  retval(2) = result.balanced_matrix ();
-                  retval(1) = result.balancing_matrix2 ();
-                  retval(0) = result.balancing_matrix ();
-                  break;
-                case 2:
-                  retval(1) = result.balancing_matrix2 ();
-                  // fall through
-                case 1:
-                  retval(0) = result.balancing_matrix ();
-                  break;
-                default:
-                  error ("balance: invalid number of output arguments");
-                  break;
-                }
-            }
-          else
-            {
-              GEPBALANCE result (aa, bb, bal_job);
-
-              switch (nargout)
-                {
-                case 4:
-                  retval(3) = result.balanced_matrix2 ();
-                  // fall through
-                case 3:
-                  retval(2) = result.balanced_matrix ();
-                  retval(1) = result.balancing_matrix2 ();
-                  retval(0) = result.balancing_matrix ();
-                  break;
-                case 2:
-                  retval(1) = result.balancing_matrix2 ();
-                  // fall through
-                case 1:
-                  retval(0) = result.balancing_matrix ();
-                  break;
-                default:
-                  error ("balance: invalid number of output arguments");
-                  break;
-                }
-            }
-        }
-    }
-
-  return retval;
-}

--- a/src/DLD-FUNCTIONS/besselj.cc	Fri Jul 27 17:10:25 2012 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,1246 +0,0 @@
-/*
-
-Copyright (C) 1997-2012 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 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 "lo-specfun.h"
-#include "quit.h"
-
-#include "defun-dld.h"
-#include "error.h"
-#include "gripes.h"
-#include "oct-obj.h"
-#include "utils.h"
-
-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 BESSEL_J: \
-            result = besselj (alpha, x, scaled, ierr); \
-            break; \
- \
-          case BESSEL_Y: \
-            result = bessely (alpha, x, scaled, ierr); \
-            break; \
- \
-          case BESSEL_I: \
-            result = besseli (alpha, x, scaled, ierr); \
-            break; \
- \
-          case BESSEL_K: \
-            result = besselk (alpha, x, scaled, ierr); \
-            break; \
- \
-          case BESSEL_H1: \
-            result = besselh1 (alpha, x, scaled, ierr); \
-            break; \
- \
-          case BESSEL_H2: \
-            result = besselh2 (alpha, x, scaled, ierr); \
-            break; \
- \
-          default: \
-            break; \
-        } \
-    } \
-  while (0)
-
-static void
-gripe_bessel_arg (const char *fn, const char *arg)
-{
-  error ("%s: expecting scalar or matrix as %s argument", fn, arg);
-}
-
-octave_value_list
-do_bessel (enum bessel_type type, const char *fn,
-           const octave_value_list& args, int nargout)
-{
-  octave_value_list retval;
-
-  int nargin = args.length ();
-
-  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_single_type () || x_arg.is_single_type ())
-        {
-          if (alpha_arg.is_scalar_type ())
-            {
-              float alpha = args(0).float_value ();
-
-              if (! error_state)
-                {
-                  if (x_arg.is_scalar_type ())
-                    {
-                      FloatComplex x = x_arg.float_complex_value ();
-
-                      if (! error_state)
-                        {
-                          octave_idx_type ierr;
-                          octave_value result;
-
-                          DO_BESSEL (type, alpha, x, scaled, ierr, result);
-
-                          if (nargout > 1)
-                            retval(1) = static_cast<float> (ierr);
-
-                          retval(0) = result;
-                        }
-                      else
-                        gripe_bessel_arg (fn, "second");
-                    }
-                  else
-                    {
-                      FloatComplexNDArray x = x_arg.float_complex_array_value ();
-
-                      if (! error_state)
-                        {
-                          Array<octave_idx_type> ierr;
-                          octave_value result;
-
-                          DO_BESSEL (type, alpha, x, scaled, ierr, result);
-
-                          if (nargout > 1)
-                            retval(1) = NDArray (ierr);
-
-                          retval(0) = result;
-                        }
-                      else
-                        gripe_bessel_arg (fn, "second");
-                    }
-                }
-              else
-                gripe_bessel_arg (fn, "first");
-            }
-          else
-            {
-              dim_vector dv0 = args(0).dims ();
-              dim_vector dv1 = args(1).dims ();
-
-              bool args0_is_row_vector = (dv0 (1) == dv0.numel ());
-              bool args1_is_col_vector = (dv1 (0) == dv1.numel ());
-
-              if (args0_is_row_vector && args1_is_col_vector)
-                {
-                  FloatRowVector ralpha = args(0).float_row_vector_value ();
-
-                  if (! error_state)
-                    {
-                      FloatComplexColumnVector cx =
-                        x_arg.float_complex_column_vector_value ();
-
-                      if (! error_state)
-                        {
-                          Array<octave_idx_type> ierr;
-                          octave_value result;
-
-                          DO_BESSEL (type, ralpha, cx, scaled, ierr, result);
-
-                          if (nargout > 1)
-                            retval(1) = NDArray (ierr);
-
-                          retval(0) = result;
-                        }
-                      else
-                        gripe_bessel_arg (fn, "second");
-                    }
-                  else
-                    gripe_bessel_arg (fn, "first");
-                }
-              else
-                {
-                  FloatNDArray alpha = args(0).float_array_value ();
-
-                  if (! error_state)
-                    {
-                      if (x_arg.is_scalar_type ())
-                        {
-                          FloatComplex x = x_arg.float_complex_value ();
-
-                          if (! error_state)
-                            {
-                              Array<octave_idx_type> ierr;
-                              octave_value result;
-
-                              DO_BESSEL (type, alpha, x, scaled, ierr, result);
-
-                              if (nargout > 1)
-                                retval(1) = NDArray (ierr);
-
-                              retval(0) = result;
-                            }
-                          else
-                            gripe_bessel_arg (fn, "second");
-                        }
-                      else
-                        {
-                          FloatComplexNDArray x = x_arg.float_complex_array_value ();
-
-                          if (! error_state)
-                            {
-                              Array<octave_idx_type> ierr;
-                              octave_value result;
-
-                              DO_BESSEL (type, alpha, x, scaled, ierr, result);
-
-                              if (nargout > 1)
-                                retval(1) = NDArray (ierr);
-
-                              retval(0) = result;
-                            }
-                          else
-                            gripe_bessel_arg (fn, "second");
-                        }
-                    }
-                  else
-                    gripe_bessel_arg (fn, "first");
-                }
-            }
-        }
-      else
-        {
-          if (alpha_arg.is_scalar_type ())
-            {
-              double alpha = args(0).double_value ();
-
-              if (! error_state)
-                {
-                  if (x_arg.is_scalar_type ())
-                    {
-                      Complex x = x_arg.complex_value ();
-
-                      if (! error_state)
-                        {
-                          octave_idx_type ierr;
-                          octave_value result;
-
-                          DO_BESSEL (type, alpha, x, scaled, ierr, result);
-
-                          if (nargout > 1)
-                            retval(1) = static_cast<double> (ierr);
-
-                          retval(0) = result;
-                        }
-                      else
-                        gripe_bessel_arg (fn, "second");
-                    }
-                  else
-                    {
-                      ComplexNDArray x = x_arg.complex_array_value ();
-
-                      if (! error_state)
-                        {
-                          Array<octave_idx_type> ierr;
-                          octave_value result;
-
-                          DO_BESSEL (type, alpha, x, scaled, ierr, result);
-
-                          if (nargout > 1)
-                            retval(1) = NDArray (ierr);
-
-                          retval(0) = result;
-                        }
-                      else
-                        gripe_bessel_arg (fn, "second");
-                    }
-                }
-              else
-                gripe_bessel_arg (fn, "first");
-            }
-          else
-            {
-              dim_vector dv0 = args(0).dims ();
-              dim_vector dv1 = args(1).dims ();
-
-              bool args0_is_row_vector = (dv0 (1) == dv0.numel ());
-              bool args1_is_col_vector = (dv1 (0) == dv1.numel ());
-
-              if (args0_is_row_vector && args1_is_col_vector)
-                {
-                  RowVector ralpha = args(0).row_vector_value ();
-
-                  if (! error_state)
-                    {
-                      ComplexColumnVector cx =
-                        x_arg.complex_column_vector_value ();
-
-                      if (! error_state)
-                        {
-                          Array<octave_idx_type> ierr;
-                          octave_value result;
-
-                          DO_BESSEL (type, ralpha, cx, scaled, ierr, result);
-
-                          if (nargout > 1)
-                            retval(1) = NDArray (ierr);
-
-                          retval(0) = result;
-                        }
-                      else
-                        gripe_bessel_arg (fn, "second");
-                    }
-                  else
-                    gripe_bessel_arg (fn, "first");
-                }
-              else
-                {
-                  NDArray alpha = args(0).array_value ();
-
-                  if (! error_state)
-                    {
-                      if (x_arg.is_scalar_type ())
-                        {
-                          Complex x = x_arg.complex_value ();
-
-                          if (! error_state)
-                            {
-                              Array<octave_idx_type> ierr;
-                              octave_value result;
-
-                              DO_BESSEL (type, alpha, x, scaled, ierr, result);
-
-                              if (nargout > 1)
-                                retval(1) = NDArray (ierr);
-
-                              retval(0) = result;
-                            }
-                          else
-                            gripe_bessel_arg (fn, "second");
-                        }
-                      else
-                        {
-                          ComplexNDArray x = x_arg.complex_array_value ();
-
-                          if (! error_state)
-                            {
-                              Array<octave_idx_type> ierr;
-                              octave_value result;
-
-                              DO_BESSEL (type, alpha, x, scaled, ierr, result);
-
-                              if (nargout > 1)
-                                retval(1) = NDArray (ierr);
-
-                              retval(0) = result;
-                            }
-                          else
-                            gripe_bessel_arg (fn, "second");
-                        }
-                    }
-                  else
-                    gripe_bessel_arg (fn, "first");
-                }
-            }
-        }
-    }
-  else
-    print_usage ();
-
-  return retval;
-}
-
-DEFUN_DLD (besselj, args, nargout,
-  "-*- texinfo -*-\n\
-@deftypefn  {Loadable Function} {[@var{j}, @var{ierr}] =} besselj (@var{alpha}, @var{x}, @var{opt})\n\
-@deftypefnx {Loadable Function} {[@var{y}, @var{ierr}] =} bessely (@var{alpha}, @var{x}, @var{opt})\n\
-@deftypefnx {Loadable Function} {[@var{i}, @var{ierr}] =} besseli (@var{alpha}, @var{x}, @var{opt})\n\
-@deftypefnx {Loadable Function} {[@var{k}, @var{ierr}] =} besselk (@var{alpha}, @var{x}, @var{opt})\n\
-@deftypefnx {Loadable Function} {[@var{h}, @var{ierr}] =} besselh (@var{alpha}, @var{k}, @var{x}, @var{opt})\n\
-Compute Bessel or Hankel functions of various kinds:\n\
-\n\
-@table @code\n\
-@item besselj\n\
-Bessel functions of the first kind.  If the argument @var{opt} is supplied,\n\
-the result is multiplied by @code{exp (-abs (imag (@var{x})))}.\n\
-\n\
-@item bessely\n\
-Bessel functions of the second kind.  If the argument @var{opt} is supplied,\n\
-the result is multiplied by @code{exp (-abs (imag (@var{x})))}.\n\
-\n\
-@item besseli\n\
-\n\
-Modified Bessel functions of the first kind.  If the argument @var{opt} is\n\
-supplied, the result is multiplied by @code{exp (-abs (real (@var{x})))}.\n\
-\n\
-@item besselk\n\
-\n\
-Modified Bessel functions of the second kind.  If the argument @var{opt} is\n\
-supplied, the result is multiplied by @code{exp (@var{x})}.\n\
-\n\
-@item besselh\n\
-Compute Hankel functions of the first (@var{k} = 1) or second (@var{k}\n\
-= 2) kind.  If the argument @var{opt} is supplied, the result is multiplied\n\
-by @code{exp (-I*@var{x})} for @var{k} = 1 or @code{exp (I*@var{x})} for\n\
-@var{k} = 2.\n\
-@end table\n\
-\n\
-If @var{alpha} is a scalar, the result is the same size as @var{x}.\n\
-If @var{x} is a scalar, the result is the same size as @var{alpha}.\n\
-If @var{alpha} is a row vector and @var{x} is a column vector, the\n\
-result is a matrix with @code{length (@var{x})} rows and\n\
-@code{length (@var{alpha})} columns.  Otherwise, @var{alpha} and\n\
-@var{x} must conform and the result will be the same size.\n\
-\n\
-The value of @var{alpha} must be real.  The value of @var{x} may be\n\
-complex.\n\
-\n\
-If requested, @var{ierr} contains the following status information\n\
-and is the same size as the result.\n\
-\n\
-@enumerate 0\n\
-@item\n\
-Normal return.\n\
-\n\
-@item\n\
-Input error, return @code{NaN}.\n\
-\n\
-@item\n\
-Overflow, return @code{Inf}.\n\
-\n\
-@item\n\
-Loss of significance by argument reduction results in less than\n\
-half of machine accuracy.\n\
-\n\
-@item\n\
-Complete loss of significance by argument reduction, return @code{NaN}.\n\
-\n\
-@item\n\
-Error---no computation, algorithm termination condition not met,\n\
-return @code{NaN}.\n\
-@end enumerate\n\
-@end deftypefn")
-{
-  return do_bessel (BESSEL_J, "besselj", args, nargout);
-}
-
-DEFUN_DLD (bessely, args, nargout,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {[@var{y}, @var{ierr}] =} bessely (@var{alpha}, @var{x}, @var{opt})\n\
-See besselj.\n\
-@end deftypefn")
-{
-  return do_bessel (BESSEL_Y, "bessely", args, nargout);
-}
-
-DEFUN_DLD (besseli, args, nargout,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {[@var{i}, @var{ierr}] =} besseli (@var{alpha}, @var{x}, @var{opt})\n\
-See besselj.\n\
-@end deftypefn")
-{
-  return do_bessel (BESSEL_I, "besseli", args, nargout);
-}
-
-DEFUN_DLD (besselk, args, nargout,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {[@var{k}, @var{ierr}] =} besselk (@var{alpha}, @var{x}, @var{opt})\n\
-See besselj.\n\
-@end deftypefn")
-{
-  return do_bessel (BESSEL_K, "besselk", args, nargout);
-}
-
-DEFUN_DLD (besselh, args, nargout,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {[@var{h}, @var{ierr}] =} besselh (@var{alpha}, @var{k}, @var{x}, @var{opt})\n\
-See besselj.\n\
-@end deftypefn")
-{
-  octave_value_list retval;
-
-  int nargin = args.length ();
-
-  if (nargin == 2)
-    {
-      retval = do_bessel (BESSEL_H1, "besselh", args, nargout);
-    }
-  else if (nargin == 3 || nargin == 4)
-    {
-      octave_idx_type kind = args(1).int_value ();
-
-      if (! error_state)
-        {
-          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 ();
-
-  return retval;
-}
-
-DEFUN_DLD (airy, args, nargout,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {[@var{a}, @var{ierr}] =} airy (@var{k}, @var{z}, @var{opt})\n\
-Compute Airy functions of the first and second kind, and their\n\
-derivatives.\n\
-\n\
-@example\n\
-@group\n\
- K   Function   Scale factor (if 'opt' 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\
-@end group\n\
-@end example\n\
-\n\
-The function call @code{airy (@var{z})} is equivalent to\n\
-@code{airy (0, @var{z})}.\n\
-\n\
-The result is the same size as @var{z}.\n\
-\n\
-If requested, @var{ierr} contains the following status information and\n\
-is the same size as the result.\n\
-\n\
-@enumerate 0\n\
-@item\n\
-Normal return.\n\
-\n\
-@item\n\
-Input error, return @code{NaN}.\n\
-\n\
-@item\n\
-Overflow, return @code{Inf}.\n\
-\n\
-@item\n\
-Loss of significance by argument reduction results in less than half\n\
- of machine accuracy.\n\
-\n\
-@item\n\
-Complete loss of significance by argument reduction, return @code{NaN}.\n\
-\n\
-@item\n\
-Error---no computation, algorithm termination condition not met,\n\
-return @code{NaN}.\n\
-@end enumerate\n\
-@end deftypefn")
-{
-  octave_value_list retval;
-
-  int nargin = args.length ();
-
-  if (nargin > 0 && nargin < 4)
-    {
-      bool scale = (nargin == 3);
-
-      int kind = 0;
-
-      if (nargin > 1)
-        {
-          kind = args(0).int_value ();
-
-          if (! error_state)
-            {
-              if (kind < 0 || kind > 3)
-                error ("airy: expecting K = 0, 1, 2, or 3");
-            }
-          else
-            error ("airy: K must be an integer value");
-        }
-
-      if (! error_state)
-        {
-          int idx = nargin == 1 ? 0 : 1;
-
-          if (args (idx).is_single_type ())
-            {
-              FloatComplexNDArray z = args(idx).float_complex_array_value ();
-
-              if (! error_state)
-                {
-                  Array<octave_idx_type> 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) = NDArray (ierr);
-
-                  retval(0) = result;
-                }
-              else
-                error ("airy: Z must be a complex matrix");
-            }
-          else
-            {
-              ComplexNDArray z = args(idx).complex_array_value ();
-
-              if (! error_state)
-                {
-                  Array<octave_idx_type> 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) = NDArray (ierr);
-
-                  retval(0) = result;
-                }
-              else
-                error ("airy: Z must be a complex matrix");
-            }
-        }
-    }
-  else
-    print_usage ();
-
-  return retval;
-}
-
-/*
-%! # Test values computed with GP/PARI version 2.3.3
-%!
-%!shared alpha, x, jx, yx, ix, kx, nix
-%!
-%! # Bessel functions, even order, positive and negative x
-%! alpha = 2;  x = 1.25;
-%! jx = 0.1710911312405234823613091417;
-%! yx = -1.193199310178553861283790424;
-%! ix = 0.2220184483766341752692212604;
-%! kx = 0.9410016167388185767085460540;
-%!
-%!assert (besselj (alpha,x), jx, 100*eps)
-%!assert (bessely (alpha,x), yx, 100*eps)
-%!assert (besseli (alpha,x), ix, 100*eps)
-%!assert (besselk (alpha,x), kx, 100*eps)
-%!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
-%!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
-%!
-%!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
-%!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
-%!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
-%!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
-%!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
-%!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
-%!
-%!assert (besselj (-alpha,x), jx, 100*eps)
-%!assert (bessely (-alpha,x), yx, 100*eps)
-%!assert (besseli (-alpha,x), ix, 100*eps)
-%!assert (besselk (-alpha,x), kx, 100*eps)
-%!assert (besselh (-alpha,1,x), jx + I*yx, 100*eps)
-%!assert (besselh (-alpha,2,x), jx - I*yx, 100*eps)
-%!
-%!assert (besselj (-alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
-%!assert (bessely (-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
-%!assert (besseli (-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
-%!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
-%!assert (besselh (-alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
-%!assert (besselh (-alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
-%!
-%! x *= -1;
-%! yx = -1.193199310178553861283790424 + 0.3421822624810469647226182835*I;
-%! kx = 0.9410016167388185767085460540 - 0.6974915263814386815610060884*I;
-%!
-%!assert (besselj (alpha,x), jx, 100*eps)
-%!assert (bessely (alpha,x), yx, 100*eps)
-%!assert (besseli (alpha,x), ix, 100*eps)
-%!assert (besselk (alpha,x), kx, 100*eps)
-%!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
-%!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
-%!
-%!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
-%!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
-%!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
-%!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
-%!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
-%!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
-%!
-%! # Bessel functions, odd order, positive and negative x
-%! alpha = 3;  x = 2.5;
-%! jx = 0.2166003910391135247666890035;
-%! yx = -0.7560554967536709968379029772;
-%! ix = 0.4743704087780355895548240179;
-%! kx = 0.2682271463934492027663765197;
-%!
-%!assert (besselj (alpha,x), jx, 100*eps)
-%!assert (bessely (alpha,x), yx, 100*eps)
-%!assert (besseli (alpha,x), ix, 100*eps)
-%!assert (besselk (alpha,x), kx, 100*eps)
-%!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
-%!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
-%!
-%!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
-%!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
-%!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
-%!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
-%!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
-%!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
-%!
-%!assert (besselj (-alpha,x), -jx, 100*eps)
-%!assert (bessely (-alpha,x), -yx, 100*eps)
-%!assert (besseli (-alpha,x), ix, 100*eps)
-%!assert (besselk (-alpha,x), kx, 100*eps)
-%!assert (besselh (-alpha,1,x), -(jx + I*yx), 100*eps)
-%!assert (besselh (-alpha,2,x), -(jx - I*yx), 100*eps)
-%!
-%!assert (besselj (-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
-%!assert (bessely (-alpha,x,1), -yx*exp(-abs(imag(x))), 100*eps)
-%!assert (besseli (-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
-%!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
-%!assert (besselh (-alpha,1,x,1), -(jx + I*yx)*exp(-I*x), 100*eps)
-%!assert (besselh (-alpha,2,x,1), -(jx - I*yx)*exp(I*x), 100*eps)
-%!
-%! x *= -1;
-%! jx = -jx;
-%! yx = 0.7560554967536709968379029772 - 0.4332007820782270495333780070*I;
-%! ix = -ix;
-%! kx = -0.2682271463934492027663765197 - 1.490278591297463775542004240*I;
-%!
-%!assert (besselj (alpha,x), jx, 100*eps)
-%!assert (bessely (alpha,x), yx, 100*eps)
-%!assert (besseli (alpha,x), ix, 100*eps)
-%!assert (besselk (alpha,x), kx, 100*eps)
-%!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
-%!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
-%!
-%!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
-%!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
-%!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
-%!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
-%!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
-%!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
-%!
-%! # Bessel functions, fractional order, positive and negative x
-%!
-%! alpha = 3.5;  x = 2.75;
-%! jx = 0.1691636439842384154644784389;
-%! yx = -0.8301381935499356070267953387;
-%! ix = 0.3930540878794826310979363668;
-%! kx = 0.2844099013460621170288192503;
-%!
-%!assert (besselj (alpha,x), jx, 100*eps)
-%!assert (bessely (alpha,x), yx, 100*eps)
-%!assert (besseli (alpha,x), ix, 100*eps)
-%!assert (besselk (alpha,x), kx, 100*eps)
-%!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
-%!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
-%!
-%!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
-%!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
-%!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
-%!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
-%!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
-%!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
-%!
-%! nix = 0.2119931212254662995364461998;
-%!
-%!assert (besselj (-alpha,x), yx, 100*eps)
-%!assert (bessely (-alpha,x), -jx, 100*eps)
-%!assert (besseli (-alpha,x), nix, 100*eps)
-%!assert (besselk (-alpha,x), kx, 100*eps)
-%!assert (besselh (-alpha,1,x), -I*(jx + I*yx), 100*eps)
-%!assert (besselh (-alpha,2,x), I*(jx - I*yx), 100*eps)
-%!
-%!assert (besselj (-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
-%!assert (bessely (-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
-%!assert (besseli (-alpha,x,1), nix*exp(-abs(real(x))), 100*eps)
-%!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
-%!assert (besselh (-alpha,1,x,1), -I*(jx + I*yx)*exp(-I*x), 100*eps)
-%!assert (besselh (-alpha,2,x,1), I*(jx - I*yx)*exp(I*x), 100*eps)
-%!
-%! x *= -1;
-%! jx *= -I;
-%! yx = -0.8301381935499356070267953387*I;
-%! ix *= -I;
-%! kx = -0.9504059335995575096509874508*I;
-%!
-%!assert (besselj (alpha,x), jx, 100*eps)
-%!assert (bessely (alpha,x), yx, 100*eps)
-%!assert (besseli (alpha,x), ix, 100*eps)
-%!assert (besselk (alpha,x), kx, 100*eps)
-%!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
-%!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
-%!
-%!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
-%!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
-%!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
-%!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
-%!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
-%!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
-%!
-%! # Bessel functions, even order, complex x
-%!
-%! alpha = 2;  x = 1.25 + 3.625 * I;
-%! jx = -1.299533366810794494030065917 + 4.370833116012278943267479589*I;
-%! yx = -4.370357232383223896393056727 - 1.283083391453582032688834041*I;
-%! ix = -0.6717801680341515541002273932 - 0.2314623443930774099910228553*I;
-%! kx = -0.01108009888623253515463783379 + 0.2245218229358191588208084197*I;
-%!
-%!assert (besselj (alpha,x), jx, 100*eps)
-%!assert (bessely (alpha,x), yx, 100*eps)
-%!assert (besseli (alpha,x), ix, 100*eps)
-%!assert (besselk (alpha,x), kx, 100*eps)
-%!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
-%!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
-%!
-%!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
-%!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
-%!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
-%!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
-%!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
-%!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
-%!
-%!assert (besselj (-alpha,x), jx, 100*eps)
-%!assert (bessely (-alpha,x), yx, 100*eps)
-%!assert (besseli (-alpha,x), ix, 100*eps)
-%!assert (besselk (-alpha,x), kx, 100*eps)
-%!assert (besselh (-alpha,1,x), jx + I*yx, 100*eps)
-%!assert (besselh (-alpha,2,x), jx - I*yx, 100*eps)
-%!
-%!assert (besselj (-alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
-%!assert (bessely (-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
-%!assert (besseli (-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
-%!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
-%!assert (besselh (-alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
-%!assert (besselh (-alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
-%!
-%! # Bessel functions, odd order, complex x
-%!
-%! alpha = 3; x = 2.5 + 1.875 * I;
-%! jx = 0.1330721523048277493333458596 + 0.5386295217249660078754395597*I;
-%! yx = -0.6485072392105829901122401551 + 0.2608129289785456797046996987*I;
-%! ix = -0.6182064685486998097516365709 + 0.4677561094683470065767989920*I;
-%! kx = -0.1568585587733540007867882337 - 0.05185853709490846050505141321*I;
-%!
-%!assert (besselj (alpha,x), jx, 100*eps)
-%!assert (bessely (alpha,x), yx, 100*eps)
-%!assert (besseli (alpha,x), ix, 100*eps)
-%!assert (besselk (alpha,x), kx, 100*eps)
-%!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
-%!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
-%!
-%!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
-%!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
-%!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
-%!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
-%!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
-%!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
-%!
-%!assert (besselj (-alpha,x), -jx, 100*eps)
-%!assert (bessely (-alpha,x), -yx, 100*eps)
-%!assert (besseli (-alpha,x), ix, 100*eps)
-%!assert (besselk (-alpha,x), kx, 100*eps)
-%!assert (besselh (-alpha,1,x), -(jx + I*yx), 100*eps)
-%!assert (besselh (-alpha,2,x), -(jx - I*yx), 100*eps)
-%!
-%!assert (besselj (-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
-%!assert (bessely (-alpha,x,1), -yx*exp(-abs(imag(x))), 100*eps)
-%!assert (besseli (-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
-%!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
-%!assert (besselh (-alpha,1,x,1), -(jx + I*yx)*exp(-I*x), 100*eps)
-%!assert (besselh (-alpha,2,x,1), -(jx - I*yx)*exp(I*x), 100*eps)
-%!
-%! # Bessel functions, fractional order, complex x
-%!
-%! alpha = 3.5;  x = 1.75 + 4.125 * I;
-%! jx = -3.018566131370455929707009100 - 0.7585648436793900607704057611*I;
-%! yx = 0.7772278839106298215614791107 - 3.018518722313849782683792010*I;
-%! ix = 0.2100873577220057189038160913 - 0.6551765604618246531254970926*I;
-%! kx = 0.1757147290513239935341488069 + 0.08772348296883849205562558311*I;
-%!
-%!assert (besselj (alpha,x), jx, 100*eps)
-%!assert (bessely (alpha,x), yx, 100*eps)
-%!assert (besseli (alpha,x), ix, 100*eps)
-%!assert (besselk (alpha,x), kx, 100*eps)
-%!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
-%!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
-%!
-%!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
-%!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
-%!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
-%!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
-%!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
-%!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
-%!
-%!  nix = 0.09822388691172060573913739253 - 0.7110230642207380127317227407*I;
-%!
-%!assert (besselj (-alpha,x), yx, 100*eps)
-%!assert (bessely (-alpha,x), -jx, 100*eps)
-%!assert (besseli (-alpha,x), nix, 100*eps)
-%!assert (besselk (-alpha,x), kx, 100*eps)
-%!assert (besselh (-alpha,1,x), -I*(jx + I*yx), 100*eps)
-%!assert (besselh (-alpha,2,x), I*(jx - I*yx), 100*eps)
-%!
-%!assert (besselj (-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
-%!assert (bessely (-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
-%!assert (besseli (-alpha,x,1), nix*exp(-abs(real(x))), 100*eps)
-%!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
-%!assert (besselh (-alpha,1,x,1), -I*(jx + I*yx)*exp(-I*x), 100*eps)
-%!assert (besselh (-alpha,2,x,1), I*(jx - I*yx)*exp(I*x), 100*eps)
-
-
-Tests contributed by Robert T. Short.
-Tests are based on the properties and tables in A&S:
- Abramowitz and Stegun, "Handbook of Mathematical Functions",
- 1972.
-
-For regular Bessel functions, there are 3 tests. These compare octave
-results against Tables 9.1, 9.2, and 9.4 in A&S. Tables 9.1 and 9.2
-are good to only a few decimal places, so any failures should be
-considered a broken implementation. Table 9.4 is an extended table
-for larger orders and arguments. There are some differences between
-Octave and Table 9.4, mostly in the last decimal place but in a very
-few instances the errors are in the last two places. The comparison
-tolerance has been changed to reflect this.
-
-Similarly for modifed Bessel functions, there are 3 tests. These
-compare octave results against Tables 9.8, 9.9, and 9.11 in A&S.
-Tables 9.8 and 9.9 are good to only a few decimal places, so any
-failures should be considered a broken implementation. Table 9.11 is
-an extended table for larger orders and arguments. There are some
-differences between octave and Table 9.11, mostly in the last decimal
-place but in a very few instances the errors are in the last two
-places. The comparison tolerance has been changed to reflect this.
-
-For spherical Bessel functions, there are also three tests, comparing
-octave results to Tables 10.1, 10.2, and 10.4 in A&S. Very similar
-comments may be made here as in the previous lines. At this time,
-modified spherical Bessel function tests are not included.
-
-% Table 9.1 - J and Y for integer orders 0, 1, 2.
-% Compare against excerpts of Table 9.1, Abramowitz and Stegun.
-%!test
-%! n = 0:2;
-%! z = (0:2.5:17.5)';
-%!
-%! Jt = [[ 1.000000000000000,  0.0000000000,  0.0000000000];
-%!       [-0.048383776468198,  0.4970941025,  0.4460590584];
-%!       [-0.177596771314338, -0.3275791376,  0.0465651163];
-%!       [ 0.266339657880378,  0.1352484276, -0.2302734105];
-%!       [-0.245935764451348,  0.0434727462,  0.2546303137];
-%!       [ 0.146884054700421, -0.1654838046, -0.1733614634];
-%!       [-0.014224472826781,  0.2051040386,  0.0415716780];
-%!       [-0.103110398228686, -0.1634199694,  0.0844338303]];
-%!
-%! Yt = [[-Inf,          -Inf,          -Inf        ];
-%!       [ 0.4980703596,  0.1459181380, -0.38133585 ];
-%!       [-0.3085176252,  0.1478631434,  0.36766288 ];
-%!       [ 0.1173132861, -0.2591285105, -0.18641422 ];
-%!       [ 0.0556711673,  0.2490154242, -0.00586808 ];
-%!       [-0.1712143068, -0.1538382565,  0.14660019 ];
-%!       [ 0.2054642960,  0.0210736280, -0.20265448 ];
-%!       [-0.1604111925,  0.0985727987,  0.17167666 ]];
-%!
-%! J = besselj (n,z);
-%! Y = bessely (n,z);
-%! assert (Jt(:,1), J(:,1), 0.5e-10);
-%! assert (Yt(:,1), Y(:,1), 0.5e-10);
-%! assert (Jt(:,2:3), J(:,2:3), 0.5e-10);
-
-Table 9.2 - J and Y for integer orders 3-9.
-
-%!test
-%! n = (3:9);
-%! z = (0:2:20).';
-%!
-%! Jt = [[ 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00];
-%!       [ 1.2894e-01, 3.3996e-02, 7.0396e-03, 1.2024e-03, 1.7494e-04, 2.2180e-05, 2.4923e-06];
-%!       [ 4.3017e-01, 2.8113e-01, 1.3209e-01, 4.9088e-02, 1.5176e-02, 4.0287e-03, 9.3860e-04];
-%!       [ 1.1477e-01, 3.5764e-01, 3.6209e-01, 2.4584e-01, 1.2959e-01, 5.6532e-02, 2.1165e-02];
-%!       [-2.9113e-01,-1.0536e-01, 1.8577e-01, 3.3758e-01, 3.2059e-01, 2.2345e-01, 1.2632e-01];
-%!       [ 5.8379e-02,-2.1960e-01,-2.3406e-01,-1.4459e-02, 2.1671e-01, 3.1785e-01, 2.9186e-01];
-%!       [ 1.9514e-01, 1.8250e-01,-7.3471e-02,-2.4372e-01,-1.7025e-01, 4.5095e-02, 2.3038e-01];
-%!       [-1.7681e-01, 7.6244e-02, 2.2038e-01, 8.1168e-02,-1.5080e-01,-2.3197e-01,-1.1431e-01];
-%!       [-4.3847e-02,-2.0264e-01,-5.7473e-02, 1.6672e-01, 1.8251e-01,-7.0211e-03,-1.8953e-01];
-%!       [ 1.8632e-01, 6.9640e-02,-1.5537e-01,-1.5596e-01, 5.1399e-02, 1.9593e-01, 1.2276e-01];
-%!       [-9.8901e-02, 1.3067e-01, 1.5117e-01,-5.5086e-02,-1.8422e-01,-7.3869e-02, 1.2513e-01]];
-%!
-%! Yt = [[       -Inf,       -Inf,       -Inf,       -Inf,       -Inf,       -Inf,       -Inf];
-%!       [-1.1278e+00,-2.7659e+00,-9.9360e+00,-4.6914e+01,-2.7155e+02,-1.8539e+03,-1.4560e+04];
-%!       [-1.8202e-01,-4.8894e-01,-7.9585e-01,-1.5007e+00,-3.7062e+00,-1.1471e+01,-4.2178e+01];
-%!       [ 3.2825e-01, 9.8391e-02,-1.9706e-01,-4.2683e-01,-6.5659e-01,-1.1052e+00,-2.2907e+00];
-%!       [ 2.6542e-02, 2.8294e-01, 2.5640e-01, 3.7558e-02,-2.0006e-01,-3.8767e-01,-5.7528e-01];
-%!       [-2.5136e-01,-1.4495e-01, 1.3540e-01, 2.8035e-01, 2.0102e-01, 1.0755e-03,-1.9930e-01];
-%!       [ 1.2901e-01,-1.5122e-01,-2.2982e-01,-4.0297e-02, 1.8952e-01, 2.6140e-01, 1.5902e-01];
-%!       [ 1.2350e-01, 2.0393e-01,-6.9717e-03,-2.0891e-01,-1.7209e-01, 3.6816e-02, 2.1417e-01];
-%!       [-1.9637e-01,-7.3222e-05, 1.9633e-01, 1.2278e-01,-1.0425e-01,-2.1399e-01,-1.0975e-01];
-%!       [ 3.3724e-02,-1.7722e-01,-1.1249e-01, 1.1472e-01, 1.8897e-01, 3.2253e-02,-1.6030e-01];
-%!       [ 1.4967e-01, 1.2409e-01,-1.0004e-01,-1.7411e-01,-4.4312e-03, 1.7101e-01, 1.4124e-01]];
-%!
-%! n = (3:9);
-%! z = (0:2:20).';
-%! J = besselj (n,z);
-%! Y = bessely (n,z);
-%!
-%! assert (J(1,:), zeros (1, columns (J)));
-%! assert (J(2:end,:), Jt(2:end,:), -5e-5);
-%! assert (Yt(1,:), Y(1,:));
-%! assert (Y(2:end,:), Yt(2:end,:), -5e-5);
-
-Table 9.4 - J and Y for various integer orders and arguments.
-
-%!test
-%! Jt = [[ 7.651976866e-01,   2.238907791e-01,  -1.775967713e-01,  -2.459357645e-01,  5.581232767e-02,  1.998585030e-02];
-%!       [ 2.497577302e-04,   7.039629756e-03,   2.611405461e-01,  -2.340615282e-01, -8.140024770e-02, -7.419573696e-02];
-%!       [ 2.630615124e-10,   2.515386283e-07,   1.467802647e-03,   2.074861066e-01, -1.138478491e-01, -5.473217694e-02];
-%!       [ 2.297531532e-17,   7.183016356e-13,   4.796743278e-07,   4.507973144e-03, -1.082255990e-01,  1.519812122e-02];
-%!       [ 3.873503009e-25,   3.918972805e-19,   2.770330052e-11,   1.151336925e-05, -1.167043528e-01,  6.221745850e-02];
-%!       [ 3.482869794e-42,   3.650256266e-33,   2.671177278e-21,   1.551096078e-12,  4.843425725e-02,  8.146012958e-02];
-%!       [ 1.107915851e-60,   1.196077458e-48,   8.702241617e-33,   6.030895312e-21, -1.381762812e-01,  7.270175482e-02];
-%!       [ 2.906004948e-80,   3.224095839e-65,   2.294247616e-45,   1.784513608e-30,  1.214090219e-01, -3.869833973e-02];
-%!       [ 8.431828790e-189,  1.060953112e-158,  6.267789396e-119,  6.597316064e-89,  1.115927368e-21,  9.636667330e-02]];
-%!
-%! Yt = [[ 8.825696420e-02,   5.103756726e-01,  -3.085176252e-01,   5.567116730e-02, -9.806499547e-02, -7.724431337e-02]
-%!       [-2.604058666e+02,  -9.935989128e+00,  -4.536948225e-01,   1.354030477e-01, -7.854841391e-02, -2.948019628e-02]
-%!       [-1.216180143e+08,  -1.291845422e+05,  -2.512911010e+01,  -3.598141522e-01,  5.723897182e-03,  5.833157424e-02]
-%!       [-9.256973276e+14,  -2.981023646e+10,  -4.694049564e+04,  -6.364745877e+00,  4.041280205e-02,  7.879068695e-02]
-%!       [-4.113970315e+22,  -4.081651389e+16,  -5.933965297e+08,  -1.597483848e+03,  1.644263395e-02,  5.124797308e-02]
-%!       [-3.048128783e+39,  -2.913223848e+30,  -4.028568418e+18,  -7.256142316e+09, -1.164572349e-01,  6.138839212e-03]
-%!       [-7.184874797e+57,  -6.661541235e+45,  -9.216816571e+29,  -1.362803297e+18, -4.530801120e-02,  4.074685217e-02]
-%!       [-2.191142813e+77,  -1.976150576e+62,  -2.788837017e+42,  -3.641066502e+27, -2.103165546e-01,  7.650526394e-02]
-%!       [-3.775287810e+185, -3.000826049e+155, -5.084863915e+115, -4.849148271e+85, -3.293800188e+18, -1.669214114e-01]];
-%!
-%! n = [(0:5:20).';30;40;50;100];
-%! z = [1,2,5,10,50,100];
-%! J = besselj (n.', z.').';
-%! Y = bessely (n.', z.').';
-%! assert (J, Jt, -1e-9);
-%! assert (Y, Yt, -1e-9);
-
-Table 9.8 - I and K for integer orders 0, 1, 2.
-
-%!test
-%! n  = 0:2;
-%! z1 = [0.1;2.5;5.0];
-%! z2 = [7.5;10.0;15.0;20.0];
-%! rtbl = [[ 0.9071009258   0.0452984468   0.1251041992   2.6823261023  10.890182683    1.995039646  ];
-%!         [ 0.2700464416   0.2065846495   0.2042345837   0.7595486903   0.9001744239   0.759126289  ];
-%!         [ 0.1835408126   0.1639722669   0.7002245988   0.5478075643   0.6002738588   0.132723593  ];
-%!         [ 0.1483158301   0.1380412115   0.111504840    0.4505236991   0.4796689336   0.57843541   ];
-%!         [ 0.1278333372   0.1212626814   0.103580801    0.3916319344   0.4107665704   0.47378525   ];
-%!         [ 0.1038995314   0.1003741751   0.090516308    0.3210023535   0.3315348950   0.36520701   ];
-%!         [ 0.0897803119   0.0875062222   0.081029690    0.2785448768   0.2854254970   0.30708743   ]];
-%!
-%! tbl = [besseli(n,z1,1), besselk(n,z1,1)];
-%! tbl(:,3) = tbl(:,3) .* (exp (z1) .* z1.^(-2));
-%! tbl(:,6) = tbl(:,6) .* (exp (-z1) .* z1.^(2));
-%! tbl = [tbl;[besseli(n,z2,1),besselk(n,z2,1)]];
-%!
-%! assert (tbl, rtbl, -2e-8);
-
-Table 9.9 - I and K for orders 3-9.
-
-%!test
-%! It = [[  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00];
-%!       [  2.8791e-02  6.8654e-03  1.3298e-03  2.1656e-04  3.0402e-05  3.7487e-06  4.1199e-07];
-%!       [  6.1124e-02  2.5940e-02  9.2443e-03  2.8291e-03  7.5698e-04  1.7968e-04  3.8284e-05];
-%!       [  7.4736e-02  4.1238e-02  1.9752e-02  8.3181e-03  3.1156e-03  1.0484e-03  3.1978e-04];
-%!       [  7.9194e-02  5.0500e-02  2.8694e-02  1.4633e-02  6.7449e-03  2.8292e-03  1.0866e-03];
-%!       [  7.9830e-02  5.5683e-02  3.5284e-02  2.0398e-02  1.0806e-02  5.2694e-03  2.3753e-03];
-%!       [  7.8848e-02  5.8425e-02  3.9898e-02  2.5176e-02  1.4722e-02  8.0010e-03  4.0537e-03];
-%!       [  7.7183e-02  5.9723e-02  4.3056e-02  2.8969e-02  1.8225e-02  1.0744e-02  5.9469e-03];
-%!       [  7.5256e-02  6.0155e-02  4.5179e-02  3.1918e-02  2.1240e-02  1.3333e-02  7.9071e-03];
-%!       [  7.3263e-02  6.0059e-02  4.6571e-02  3.4186e-02  2.3780e-02  1.5691e-02  9.8324e-03];
-%!       [  7.1300e-02  5.9640e-02  4.7444e-02  3.5917e-02  2.5894e-02  1.7792e-02  1.1661e-02]];
-%!
-%! Kt = [[ Inf         Inf         Inf         Inf         Inf         Inf         Inf];
-%!      [  4.7836e+00  1.6226e+01  6.9687e+01  3.6466e+02  2.2576e+03  1.6168e+04  1.3160e+05];
-%!      [  1.6317e+00  3.3976e+00  8.4268e+00  2.4465e+01  8.1821e+01  3.1084e+02  1.3252e+03];
-%!      [  9.9723e-01  1.6798e+00  3.2370e+00  7.0748e+00  1.7387e+01  4.7644e+01  1.4444e+02];
-%!      [  7.3935e-01  1.1069e+00  1.8463e+00  3.4148e+00  6.9684e+00  1.5610e+01  3.8188e+01];
-%!      [  6.0028e-01  8.3395e-01  1.2674e+00  2.1014e+00  3.7891e+00  7.4062e+00  1.5639e+01];
-%!      [  5.1294e-01  6.7680e-01  9.6415e-01  1.4803e+00  2.4444e+00  4.3321e+00  8.2205e+00];
-%!      [  4.5266e-01  5.7519e-01  7.8133e-01  1.1333e+00  1.7527e+00  2.8860e+00  5.0510e+00];
-%!      [  4.0829e-01  5.0414e-01  6.6036e-01  9.1686e-01  1.3480e+00  2.0964e+00  3.4444e+00];
-%!      [  3.7411e-01  4.5162e-01  5.7483e-01  7.7097e-01  1.0888e+00  1.6178e+00  2.5269e+00];
-%!      [  3.4684e-01  4.1114e-01  5.1130e-01  6.6679e-01  9.1137e-01  1.3048e+00  1.9552e+00]];
-%!
-%! n = (3:9);
-%! z = (0:2:20).';
-%! I = besseli (n,z,1);
-%! K = besselk (n,z,1);
-%!
-%! assert (abs (I(1,:)), zeros (1, columns (I)));
-%! assert (I(2:end,:), It(2:end,:), -5e-5);
-%! assert (Kt(1,:), K(1,:));
-%! assert (K(2:end,:), Kt(2:end,:), -5e-5);
-
-Table 9.11 - I and K for various integer orders and arguments.
-
-%!test
-%! It = [[   1.266065878e+00    2.279585302e+00    2.723987182e+01    2.815716628e+03     2.93255378e+20     1.07375171e+42 ];
-%!       [   2.714631560e-04    9.825679323e-03    2.157974547e+00    7.771882864e+02     2.27854831e+20     9.47009387e+41 ];
-%!       [   2.752948040e-10    3.016963879e-07    4.580044419e-03    2.189170616e+01     1.07159716e+20     6.49897552e+41 ];
-%!       [   2.370463051e-17    8.139432531e-13    1.047977675e-06    1.043714907e-01     3.07376455e+19     3.47368638e+41 ];
-%!       [   3.966835986e-25    4.310560576e-19    5.024239358e-11    1.250799736e-04     5.44200840e+18     1.44834613e+41 ];
-%!       [   3.539500588e-42    3.893519664e-33    3.997844971e-21    7.787569783e-12     4.27499365e+16     1.20615487e+40 ];
-%!       [   1.121509741e-60    1.255869192e-48    1.180426980e-32    2.042123274e-20     6.00717897e+13     3.84170550e+38 ];
-%!       [   2.934635309e-80    3.353042830e-65    2.931469647e-45    4.756894561e-30     1.76508024e+10     4.82195809e+36 ];
-%!       [   8.473674008e-189   1.082171475e-158   7.093551489e-119   1.082344202e-88     2.72788795e-16     4.64153494e+21 ]];
-%!
-%! Kt = [[   4.210244382e-01    1.138938727e-01    3.691098334e-03    1.778006232e-05     3.41016774e-23     4.65662823e-45 ];
-%!       [   3.609605896e+02    9.431049101e+00    3.270627371e-02    5.754184999e-05     4.36718224e-23     5.27325611e-45 ];
-%!       [   1.807132899e+08    1.624824040e+05    9.758562829e+00    1.614255300e-03     9.15098819e-23     7.65542797e-45 ];
-%!       [   1.403066801e+15    4.059213332e+10    3.016976630e+04    2.656563849e-01     3.11621117e-22     1.42348325e-44 ];
-%!       [   6.294369360e+22    5.770856853e+16    4.827000521e+08    1.787442782e+02     1.70614838e-21     3.38520541e-44 ];
-%!       [   4.706145527e+39    4.271125755e+30    4.112132063e+18    2.030247813e+09     2.00581681e-19     3.97060205e-43 ];
-%!       [   1.114220651e+58    9.940839886e+45    1.050756722e+30    5.938224681e+17     1.29986971e-16     1.20842080e-41 ];
-%!       [   3.406896854e+77    2.979981740e+62    3.394322243e+42    2.061373775e+27     4.00601349e-13     9.27452265e-40 ];
-%!       [   5.900333184e+185   4.619415978e+155   7.039860193e+115   4.596674084e+85     1.63940352e+13     7.61712963e-25 ]];
-%!
-%! n = [(0:5:20).';30;40;50;100];
-%! z = [1,2,5,10,50,100];
-%! I = besseli (n.', z.').';
-%! K = besselk (n.', z.').';
-%! assert (I, It, -5e-9);
-%! assert (K, Kt, -5e-9);
-
-The next section checks that negative integer orders and positive
-integer orders are appropriately related.
-
-%!test
-%! n = (0:2:20);
-%! assert (besselj (n,1), besselj (-n,1), 1e-8);
-%! assert (-besselj (n+1,1), besselj (-n-1,1), 1e-8);
-
-besseli (n,z) = besseli (-n,z);
-
-%!test
-%! n = (0:2:20);
-%! assert (besseli (n,1), besseli (-n,1), 1e-8);
-
-Table 10.1 - j and y for integer orders 0, 1, 2.
-Compare against excerpts of Table 10.1, Abramowitz and Stegun.
-
-%!test
-%! n = (0:2);
-%! z = [0.1;(2.5:2.5:10.0).'];
-%!
-%! jt = [[ 9.9833417e-01  3.33000119e-02  6.6619061e-04 ];
-%!       [ 2.3938886e-01  4.16212989e-01  2.6006673e-01 ];
-%!       [-1.9178485e-01 -9.50894081e-02  1.3473121e-01 ];
-%!       [    1.2507e-01     -2.9542e-02    -1.3688e-01 ];
-%!       [   -5.4402e-02      7.8467e-02     7.7942e-02 ]];
-%!
-%! yt = [[-9.9500417e+00  -1.0049875e+02 -3.0050125e+03 ];
-%!       [ 3.2045745e-01  -1.1120588e-01 -4.5390450e-01 ];
-%!       [-5.6732437e-02   1.8043837e-01  1.6499546e-01 ];
-%!       [   -4.6218e-02     -1.3123e-01    -6.2736e-03 ];
-%!       [    8.3907e-02      6.2793e-02    -6.5069e-02 ]];
-%!
-%! j = sqrt ((pi/2)./z) .* besselj (n+1/2,z);
-%! y = sqrt ((pi/2)./z) .* bessely (n+1/2,z);
-%! assert (jt, j, -5e-5);
-%! assert (yt, y, -5e-5);
-
-Table 10.2 - j and y for orders 3-8.
-Compare against excerpts of Table 10.2, Abramowitzh and Stegun.
-
- Important note: In A&S, y_4(0.1) = -1.0507e+7, but Octave returns
- y_4(0.1) = -1.0508e+07 (-10507503.75).  If I compute the same term using
- a series, the difference is in the eighth significant digit so I left
- the Octave results in place.
-
-%!test
-%! n = (3:8);
-%! z = (0:2.5:10).';  z(1) = 0.1;
-%!
-%! jt = [[ 9.5185e-06  1.0577e-07  9.6163e-10  7.3975e-12  4.9319e-14  2.9012e-16];
-%!       [ 1.0392e-01  3.0911e-02  7.3576e-03  1.4630e-03  2.5009e-04  3.7516e-05];
-%!       [ 2.2982e-01  1.8702e-01  1.0681e-01  4.7967e-02  1.7903e-02  5.7414e-03];
-%!       [-6.1713e-02  7.9285e-02  1.5685e-01  1.5077e-01  1.0448e-01  5.8188e-02];
-%!       [-3.9496e-02 -1.0559e-01 -5.5535e-02  4.4501e-02  1.1339e-01  1.2558e-01]];
-%!
-%! yt = [[-1.5015e+05 -1.0508e+07 -9.4553e+08 -1.0400e+11 -1.3519e+13 -2.0277e+15];
-%!       [-7.9660e-01 -1.7766e+00 -5.5991e+00 -2.2859e+01 -1.1327e+02 -6.5676e+02];
-%!       [-1.5443e-02 -1.8662e-01 -3.2047e-01 -5.1841e-01 -1.0274e+00 -2.5638e+00];
-%!       [ 1.2705e-01  1.2485e-01  2.2774e-02 -9.1449e-02 -1.8129e-01 -2.7112e-01];
-%!       [-9.5327e-02 -1.6599e-03  9.3834e-02  1.0488e-01  4.2506e-02 -4.1117e-02]];
-%!
-%! j = sqrt ((pi/2)./z) .* besselj (n+1/2,z);
-%! y = sqrt ((pi/2)./z) .* bessely (n+1/2,z);
-%!
-%! assert (jt, j, -5e-5);
-%! assert (yt, y, -5e-5);
-
-Table 10.4 - j and y for various integer orders and arguments.
-
-%!test
-%! jt = [[ 8.414709848e-01    4.546487134e-01   -1.917848549e-01   -5.440211109e-02   -5.247497074e-03   -5.063656411e-03];
-%!       [ 9.256115861e-05    2.635169770e-03    1.068111615e-01   -5.553451162e-02   -2.004830056e-02   -9.290148935e-03];
-%!       [ 7.116552640e-11    6.825300865e-08    4.073442442e-04    6.460515449e-02   -1.503922146e-02   -1.956578597e-04];
-%!       [ 5.132686115e-18    1.606982166e-13    1.084280182e-07    1.063542715e-03   -1.129084539e-02    7.877261748e-03];
-%!       [ 7.537795722e-26    7.632641101e-20    5.427726761e-12    2.308371961e-06   -1.578502990e-02    1.010767128e-02];
-%!       [ 5.566831267e-43    5.836617888e-34    4.282730217e-22    2.512057385e-13   -1.494673454e-03    8.700628514e-03];
-%!       [ 1.538210374e-61    1.660978779e-49    1.210347583e-33    8.435671634e-22   -2.606336952e-02    1.043410851e-02];
-%!       [ 3.615274717e-81    4.011575290e-66    2.857479350e-46    2.230696023e-31    1.882910737e-02    5.797140882e-04];
-%!       [7.444727742e-190   9.367832591e-160   5.535650303e-120    5.832040182e-90    1.019012263e-22    1.088047701e-02]];
-%!
-%! yt = [[ -5.403023059e-01    2.080734183e-01   -5.673243709e-02    8.390715291e-02   -1.929932057e-02   -8.623188723e-03]
-%!       [ -9.994403434e+02   -1.859144531e+01   -3.204650467e-01    9.383354168e-02   -6.971131965e-04    3.720678486e-03]
-%!       [ -6.722150083e+08   -3.554147201e+05   -2.665611441e+01   -1.724536721e-01    1.352468751e-02    1.002577737e-02]
-%!       [ -6.298007233e+15   -1.012182944e+11   -6.288146513e+04   -3.992071745e+00    1.712319725e-02    6.258641510e-03]
-%!       [ -3.239592219e+23   -1.605436493e+17   -9.267951403e+08   -1.211210605e+03    1.375953130e-02    5.631729379e-05]
-%!       [ -2.946428547e+40   -1.407393871e+31   -7.760717570e+18   -6.908318646e+09   -2.241226812e-02   -5.412929349e-03]
-%!       [ -8.028450851e+58   -3.720929322e+46   -2.055758716e+30   -1.510304919e+18    4.978797221e-05   -7.048420407e-04]
-%!       [ -2.739192285e+78   -1.235021944e+63   -6.964109188e+42   -4.528227272e+27   -4.190000150e-02    1.074782297e-02]
-%!       [-6.683079463e+186  -2.655955830e+156  -1.799713983e+116   -8.573226309e+85   -1.125692891e+18   -2.298385049e-02]];
-%!
-%! n = [(0:5:20).';30;40;50;100];
-%! z = [1,2,5,10,50,100];
-%! j = sqrt ((pi/2)./z) .* besselj ((n+1/2).', z.').';
-%! y = sqrt ((pi/2)./z) .* bessely ((n+1/2).', z.').';
-%! assert (j, jt, -1e-9);
-%! assert (y, yt, -1e-9);
-*/

--- a/src/DLD-FUNCTIONS/betainc.cc	Fri Jul 27 17:10:25 2012 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,491 +0,0 @@
-/*
-
-Copyright (C) 1997-2012 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 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 "lo-specfun.h"
-
-#include "defun-dld.h"
-#include "error.h"
-#include "gripes.h"
-#include "oct-obj.h"
-#include "utils.h"
-
-// FIXME: These functions do not need to be dynamically loaded.  They should
-//        be placed elsewhere in the Octave code hierarchy.
-
-DEFUN_DLD (betainc, args, ,
-  "-*- texinfo -*-\n\
-@deftypefn {Mapping Function} {} betainc (@var{x}, @var{a}, @var{b})\n\
-Return the regularized incomplete Beta function,\n\
-@tex\n\
-$$\n\
- I (x, a, b) = {1 \\over {B (a, b)}} \\int_0^x t^{(a-z)} (1-t)^{(b-1)} dt.\n\
-$$\n\
-@end tex\n\
-@ifnottex\n\
-@c Set example in small font to prevent overfull line\n\
-\n\
-@smallexample\n\
-@group\n\
-                                   x\n\
-                          1       /\n\
-betainc (x, a, b) = -----------   | t^(a-1) (1-t)^(b-1) dt.\n\
-                    beta (a, b)   /\n\
-                               t=0\n\
-@end group\n\
-@end smallexample\n\
-\n\
-@end ifnottex\n\
-\n\
-If @var{x} has more than one component, both @var{a} and @var{b} must be\n\
-scalars.  If @var{x} is a scalar, @var{a} and @var{b} must be of\n\
-compatible dimensions.\n\
-@seealso{betaincinv, beta, betaln}\n\
-@end deftypefn")
-{
-  octave_value retval;
-
-  int nargin = args.length ();
-
-  if (nargin == 3)
-    {
-      octave_value x_arg = args(0);
-      octave_value a_arg = args(1);
-      octave_value b_arg = args(2);
-
-      // FIXME Can we make a template version of the duplicated code below
-      if (x_arg.is_single_type () || a_arg.is_single_type () ||
-          b_arg.is_single_type ())
-        {
-          if (x_arg.is_scalar_type ())
-            {
-              float x = x_arg.float_value ();
-
-              if (a_arg.is_scalar_type ())
-                {
-                  float a = a_arg.float_value ();
-
-                  if (! error_state)
-                    {
-                      if (b_arg.is_scalar_type ())
-                        {
-                          float b = b_arg.float_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                      else
-                        {
-                          Array<float> b = b_arg.float_array_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                    }
-                }
-              else
-                {
-                  Array<float> a = a_arg.float_array_value ();
-
-                  if (! error_state)
-                    {
-                      if (b_arg.is_scalar_type ())
-                        {
-                          float b = b_arg.float_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                      else
-                        {
-                          Array<float> b = b_arg.float_array_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                    }
-                }
-            }
-          else
-            {
-              Array<float> x = x_arg.float_array_value ();
-
-              if (a_arg.is_scalar_type ())
-                {
-                  float a = a_arg.float_value ();
-
-                  if (! error_state)
-                    {
-                      if (b_arg.is_scalar_type ())
-                        {
-                          float b = b_arg.float_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                      else
-                        {
-                          Array<float> b = b_arg.float_array_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                    }
-                }
-              else
-                {
-                  Array<float> a = a_arg.float_array_value ();
-
-                  if (! error_state)
-                    {
-                      if (b_arg.is_scalar_type ())
-                        {
-                          float b = b_arg.float_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                      else
-                        {
-                          Array<float> b = b_arg.float_array_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                    }
-                }
-            }
-        }
-      else
-        {
-          if (x_arg.is_scalar_type ())
-            {
-              double x = x_arg.double_value ();
-
-              if (a_arg.is_scalar_type ())
-                {
-                  double a = a_arg.double_value ();
-
-                  if (! error_state)
-                    {
-                      if (b_arg.is_scalar_type ())
-                        {
-                          double b = b_arg.double_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                      else
-                        {
-                          Array<double> b = b_arg.array_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                    }
-                }
-              else
-                {
-                  Array<double> a = a_arg.array_value ();
-
-                  if (! error_state)
-                    {
-                      if (b_arg.is_scalar_type ())
-                        {
-                          double b = b_arg.double_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                      else
-                        {
-                          Array<double> b = b_arg.array_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                    }
-                }
-            }
-          else
-            {
-              Array<double> x = x_arg.array_value ();
-
-              if (a_arg.is_scalar_type ())
-                {
-                  double a = a_arg.double_value ();
-
-                  if (! error_state)
-                    {
-                      if (b_arg.is_scalar_type ())
-                        {
-                          double b = b_arg.double_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                      else
-                        {
-                          Array<double> b = b_arg.array_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                    }
-                }
-              else
-                {
-                  Array<double> a = a_arg.array_value ();
-
-                  if (! error_state)
-                    {
-                      if (b_arg.is_scalar_type ())
-                        {
-                          double b = b_arg.double_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                      else
-                        {
-                          Array<double> b = b_arg.array_value ();
-
-                          if (! error_state)
-                            retval = betainc (x, a, b);
-                        }
-                    }
-                }
-            }
-        }
-    }
-  else
-    print_usage ();
-
-  return retval;
-}
-
-/*
-## Double precision
-%!test
-%! a = [1, 1.5, 2, 3];
-%! b = [4, 3, 2, 1];
-%! v1 = betainc (1,a,b);
-%! v2 = [1,1,1,1];
-%! x = [.2, .4, .6, .8];
-%! v3 = betainc (x, a, b);
-%! v4 = 1 - betainc (1.-x, b, a);
-%! assert (v1, v2, sqrt (eps));
-%! assert (v3, v4, sqrt (eps));
-
-## Single precision
-%!test
-%! a = single ([1, 1.5, 2, 3]);
-%! b = single ([4, 3, 2, 1]);
-%! v1 = betainc (1,a,b);
-%! v2 = single ([1,1,1,1]);
-%! x = single ([.2, .4, .6, .8]);
-%! v3 = betainc (x, a, b);
-%! v4 = 1 - betainc (1.-x, b, a);
-%! assert (v1, v2, sqrt (eps ("single")));
-%! assert (v3, v4, sqrt (eps ("single")));
-
-## Mixed double/single precision
-%!test
-%! a = single ([1, 1.5, 2, 3]);
-%! b = [4, 3, 2, 1];
-%! v1 = betainc (1,a,b);
-%! v2 = single ([1,1,1,1]);
-%! x = [.2, .4, .6, .8];
-%! v3 = betainc (x, a, b);
-%! v4 = 1-betainc (1.-x, b, a);
-%! assert (v1, v2, sqrt (eps ("single")));
-%! assert (v3, v4, sqrt (eps ("single")));
-
-%!error betainc ()
-%!error betainc (1)
-%!error betainc (1,2)
-%!error betainc (1,2,3,4)
-*/
-
-DEFUN_DLD (betaincinv, args, ,
-    "-*- texinfo -*-\n\
-@deftypefn {Mapping Function} {} betaincinv (@var{y}, @var{a}, @var{b})\n\
-Compute the inverse of the incomplete Beta function, i.e., @var{x} such that\n\
-\n\
-@example\n\
-@var{y} == betainc (@var{x}, @var{a}, @var{b}) \n\
-@end example\n\
-@seealso{betainc, beta, betaln}\n\
-@end deftypefn")
-{
-  octave_value retval;
-
-  int nargin = args.length ();
-
-  if (nargin == 3)
-    {
-      octave_value x_arg = args(0);
-      octave_value a_arg = args(1);
-      octave_value b_arg = args(2);
-
-      if (x_arg.is_scalar_type ())
-        {
-          double x = x_arg.double_value ();
-
-          if (a_arg.is_scalar_type ())
-            {
-              double a = a_arg.double_value ();
-
-              if (! error_state)
-                {
-                  if (b_arg.is_scalar_type ())
-                    {
-                      double b = b_arg.double_value ();
-
-                      if (! error_state)
-                        retval = betaincinv (x, a, b);
-                    }
-                  else
-                    {
-                      Array<double> b = b_arg.array_value ();
-
-                      if (! error_state)
-                        retval = betaincinv (x, a, b);
-                    }
-                }
-            }
-          else
-            {
-              Array<double> a = a_arg.array_value ();
-
-              if (! error_state)
-                {
-                  if (b_arg.is_scalar_type ())
-                    {
-                      double b = b_arg.double_value ();
-
-                      if (! error_state)
-                        retval = betaincinv (x, a, b);
-                    }
-                  else
-                    {
-                      Array<double> b = b_arg.array_value ();
-
-                      if (! error_state)
-                        retval = betaincinv (x, a, b);
-                    }
-                }
-            }
-        }
-      else
-        {
-          Array<double> x = x_arg.array_value ();
-
-          if (a_arg.is_scalar_type ())
-            {
-              double a = a_arg.double_value ();
-
-              if (! error_state)
-                {
-                  if (b_arg.is_scalar_type ())
-                    {
-                      double b = b_arg.double_value ();
-
-                      if (! error_state)
-                        retval = betaincinv (x, a, b);
-                    }
-                  else
-                    {
-                      Array<double> b = b_arg.array_value ();
-
-                      if (! error_state)
-                        retval = betaincinv (x, a, b);
-                    }
-                }
-            }
-          else
-            {
-              Array<double> a = a_arg.array_value ();
-
-              if (! error_state)
-                {
-                  if (b_arg.is_scalar_type ())
-                    {
-                      double b = b_arg.double_value ();
-
-                      if (! error_state)
-                        retval = betaincinv (x, a, b);
-                    }
-                  else
-                    {
-                      Array<double> b = b_arg.array_value ();
-
-                      if (! error_state)
-                        retval = betaincinv (x, a, b);
-                    }
-                }
-            }
-        }
-
-      // FIXME: It would be better to have an algorithm for betaincinv which
-      // accepted float inputs and returned float outputs.  As it is, we do
-      // extra work to calculate betaincinv to double precision and then throw
-      // that precision away.
-      if (x_arg.is_single_type () || a_arg.is_single_type () ||
-          b_arg.is_single_type ())
-        {
-          retval = Array<float> (retval.array_value ());
-        }
-    }
-  else
-    print_usage ();
-
-  return retval;
-}
-
-/*
-%!assert (betaincinv ([0.875 0.6875], [1 2], 3), [0.5 0.5], sqrt (eps))
-%!assert (betaincinv (0.5, 3, 3), 0.5, sqrt (eps))
-%!assert (betaincinv (0.34375, 4, 3), 0.5, sqrt (eps))
-%!assert (betaincinv (0.2265625, 5, 3), 0.5, sqrt (eps))
-%!assert (betaincinv (0.14453125, 6, 3), 0.5, sqrt (eps))
-%!assert (betaincinv (0.08984375, 7, 3), 0.5, sqrt (eps))
-%!assert (betaincinv (0.0546875, 8, 3), 0.5, sqrt (eps))
-%!assert (betaincinv (0.03271484375, 9, 3), 0.5, sqrt (eps))
-%!assert (betaincinv (0.019287109375, 10, 3), 0.5, sqrt (eps))
-
-## Test class single as well
-%!assert (betaincinv ([0.875 0.6875], [1 2], single (3)), [0.5 0.5], sqrt (eps ("single")))
-%!assert (betaincinv (0.5, 3, single (3)), 0.5, sqrt (eps ("single")))
-%!assert (betaincinv (0.34375, 4, single (3)), 0.5, sqrt (eps ("single")))
-
-## Extreme values
-%!assert (betaincinv (0, 42, 42), 0, sqrt (eps))
-%!assert (betaincinv (1, 42, 42), 1, sqrt (eps))
-
-%!error betaincinv ()
-%!error betaincinv (1, 2)
-*/
-