oct-specfun.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 1996-2026 The Octave Project Developers
//
// See the file COPYRIGHT.md in the top-level directory of this
// distribution or <https://octave.org/copyright/>.
//
// 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
// <https://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////
 
#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif
 
#include <cmath>
 
#include <algorithm>
#include <limits>
#include <string>
 
#include "CColVector.h"
#include "CMatrix.h"
#include "CNDArray.h"
#include "Faddeeva.hh"
#include "amos-proto.h"
#include "dMatrix.h"
#include "dNDArray.h"
#include "dRowVector.h"
#include "f77-fcn.h"
#include "fCColVector.h"
#include "fCMatrix.h"
#include "fCNDArray.h"
#include "fMatrix.h"
#include "fNDArray.h"
#include "fRowVector.h"
#include "lo-ieee.h"
#include "mappers.h"
#include "mx-inlines.cc"
#include "oct-error.h"
#include "oct-specfun.h"
#include "slatec-proto.h"
 
OCTAVE_BEGIN_NAMESPACE(octave)
OCTAVE_BEGIN_NAMESPACE(math)
 
static inline Complex
bessel_return_value (const Complex& val, octave_idx_type ierr)
{
  static const Complex inf_val
    = Complex (numeric_limits<double>::Inf (),
               numeric_limits<double>::Inf ());
 
  static const Complex nan_val
    = Complex (numeric_limits<double>::NaN (),
               numeric_limits<double>::NaN ());
 
  Complex retval;
 
  switch (ierr)
    {
    case 0:
    case 3:
    case 4:
      retval = val;
      break;
 
    case 2:
      retval = inf_val;
      break;
 
    default:
      retval = nan_val;
      break;
    }
 
  return retval;
}
 
static inline FloatComplex
bessel_return_value (const FloatComplex& val, octave_idx_type ierr)
{
  static const FloatComplex inf_val
    = FloatComplex (numeric_limits<float>::Inf (),
                    numeric_limits<float>::Inf ());
 
  static const FloatComplex nan_val
    = FloatComplex (numeric_limits<float>::NaN (),
                    numeric_limits<float>::NaN ());
 
  FloatComplex retval;
 
  switch (ierr)
    {
    case 0:
    case 3:
    case 4:
      retval = val;
      break;
 
    case 2:
      retval = inf_val;
      break;
 
    default:
      retval = nan_val;
      break;
    }
 
  return retval;
}
 
Complex
airy (const Complex& z, bool deriv, bool scaled, octave_idx_type& ierr)
{
  double ar = 0.0;
  double ai = 0.0;
 
  double zr = z.real ();
  double zi = z.imag ();
 
  F77_INT id = (deriv ? 1 : 0);
  F77_INT nz, t_ierr;
  F77_INT sc = (scaled ? 2 : 1);
 
  F77_FUNC (zairy, ZAIRY) (zr, zi, id, sc, ar, ai, nz, t_ierr);
 
  ierr = t_ierr;
 
  if (zi == 0.0 && (! scaled || zr >= 0.0))
    ai = 0.0;
 
  return bessel_return_value (Complex (ar, ai), ierr);
}
 
ComplexMatrix
airy (const ComplexMatrix& z, bool deriv, bool scaled,
      Array<octave_idx_type>& ierr)
{
  octave_idx_type nr = z.rows ();
  octave_idx_type nc = z.cols ();
 
  ComplexMatrix retval (nr, nc);
 
  ierr.resize (dim_vector (nr, nc));
 
  for (octave_idx_type j = 0; j < nc; j++)
    for (octave_idx_type i = 0; i < nr; i++)
      retval(i, j) = airy (z(i, j), deriv, scaled, ierr(i, j));
 
  return retval;
}
 
ComplexNDArray
airy (const ComplexNDArray& z, bool deriv, bool scaled,
      Array<octave_idx_type>& ierr)
{
  const dim_vector& dv = z.dims ();
  octave_idx_type nel = dv.numel ();
  ComplexNDArray retval (dv);
 
  ierr.resize (dv);
 
  for (octave_idx_type i = 0; i < nel; i++)
    retval(i) = airy (z(i), deriv, scaled, ierr(i));
 
  return retval;
}
 
FloatComplex
airy (const FloatComplex& z, bool deriv, bool scaled,
      octave_idx_type& ierr)
{
  FloatComplex a;
 
  F77_INT id = (deriv ? 1 : 0);
  F77_INT nz, t_ierr;
  F77_INT sc = (scaled ? 2 : 1);
 
  F77_FUNC (cairy, CAIRY) (F77_CONST_CMPLX_ARG (&z), id, sc,
                           F77_CMPLX_ARG (&a), nz, t_ierr);
 
  ierr = t_ierr;
 
  float ar = a.real ();
  float ai = a.imag ();
 
  if (z.imag () == 0.0 && (! scaled || z.real () >= 0.0))
    ai = 0.0;
 
  return bessel_return_value (FloatComplex (ar, ai), ierr);
}
 
FloatComplexMatrix
airy (const FloatComplexMatrix& z, bool deriv, bool scaled,
      Array<octave_idx_type>& ierr)
{
  octave_idx_type nr = z.rows ();
  octave_idx_type nc = z.cols ();
 
  FloatComplexMatrix retval (nr, nc);
 
  ierr.resize (dim_vector (nr, nc));
 
  for (octave_idx_type j = 0; j < nc; j++)
    for (octave_idx_type i = 0; i < nr; i++)
      retval(i, j) = airy (z(i, j), deriv, scaled, ierr(i, j));
 
  return retval;
}
 
FloatComplexNDArray
airy (const FloatComplexNDArray& z, bool deriv, bool scaled,
      Array<octave_idx_type>& ierr)
{
  const dim_vector& dv = z.dims ();
  octave_idx_type nel = dv.numel ();
  FloatComplexNDArray retval (dv);
 
  ierr.resize (dv);
 
  for (octave_idx_type i = 0; i < nel; i++)
    retval(i) = airy (z(i), deriv, scaled, ierr(i));
 
  return retval;
}
 
static inline Complex
zbesj (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
 
static inline Complex
zbesy (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
 
static inline Complex
zbesi (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
 
static inline Complex
zbesk (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
 
static inline Complex
zbesh1 (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
 
static inline Complex
zbesh2 (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
 
static inline Complex
zbesj (const Complex& z, double alpha, int kode, octave_idx_type& ierr)
{
  Complex retval;
 
  if (alpha >= 0.0)
    {
      double yr = 0.0;
      double yi = 0.0;
 
      F77_INT nz, t_ierr;
 
      double zr = z.real ();
      double zi = z.imag ();
 
      F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, kode, 1, &yr, &yi, nz, t_ierr);
 
      ierr = t_ierr;
 
      if (zi == 0.0 && zr >= 0.0)
        yi = 0.0;
 
      retval = bessel_return_value (Complex (yr, yi), ierr);
    }
  else if (octave::math::is_integer (alpha))
    {
      // zbesy can overflow as z->0, and cause troubles for generic case below
      alpha = -alpha;
      Complex tmp = zbesj (z, alpha, kode, ierr);
      if ((static_cast<long> (alpha)) & 1)
        tmp = - tmp;
      retval = bessel_return_value (tmp, ierr);
    }
  else
    {
      alpha = -alpha;
 
      Complex tmp = cos (M_PI * alpha) * zbesj (z, alpha, kode, ierr);
 
      if (ierr == 0 || ierr == 3)
        {
          tmp -= sin (M_PI * alpha) * zbesy (z, alpha, kode, ierr);
 
          retval = bessel_return_value (tmp, ierr);
        }
      else
        retval = Complex (numeric_limits<double>::NaN (),
                          numeric_limits<double>::NaN ());
    }
 
  return retval;
}
 
static inline Complex
zbesy (const Complex& z, double alpha, int kode, octave_idx_type& ierr)
{
  Complex retval;
 
  if (alpha >= 0.0)
    {
      double yr = 0.0;
      double yi = 0.0;
 
      F77_INT nz, t_ierr;
 
      double wr, wi;
 
      double zr = z.real ();
      double zi = z.imag ();
 
      ierr = 0;
 
      if (zr == 0.0 && zi == 0.0)
        {
          yr = -numeric_limits<double>::Inf ();
          yi = 0.0;
        }
      else
        {
          F77_FUNC (zbesy, ZBESY) (zr, zi, alpha, kode, 1, &yr, &yi, nz,
                                   &wr, &wi, t_ierr);
 
          ierr = t_ierr;
 
          if (zi == 0.0 && zr >= 0.0)
            yi = 0.0;
        }
 
      return bessel_return_value (Complex (yr, yi), ierr);
    }
  else if (octave::math::is_integer (alpha - 0.5))
    {
      // zbesy can overflow as z->0, and cause troubles for generic case below
      alpha = -alpha;
      Complex tmp = zbesj (z, alpha, kode, ierr);
      if ((static_cast<long> (alpha - 0.5)) & 1)
        tmp = - tmp;
      retval = bessel_return_value (tmp, ierr);
    }
  else
    {
      alpha = -alpha;
 
      Complex tmp = cos (M_PI * alpha) * zbesy (z, alpha, kode, ierr);
 
      if (ierr == 0 || ierr == 3)
        {
          tmp += sin (M_PI * alpha) * zbesj (z, alpha, kode, ierr);
 
          retval = bessel_return_value (tmp, ierr);
        }
      else
        retval = Complex (numeric_limits<double>::NaN (),
                          numeric_limits<double>::NaN ());
    }
 
  return retval;
}
 
static inline Complex
zbesi (const Complex& z, double alpha, int kode, octave_idx_type& ierr)
{
  Complex retval;
 
  if (alpha >= 0.0)
    {
      double yr = 0.0;
      double yi = 0.0;
 
      F77_INT nz, t_ierr;
 
      double zr = z.real ();
      double zi = z.imag ();
 
      F77_FUNC (zbesi, ZBESI) (zr, zi, alpha, kode, 1, &yr, &yi, nz, t_ierr);
 
      ierr = t_ierr;
 
      if (zi == 0.0 && zr >= 0.0)
        yi = 0.0;
 
      retval = bessel_return_value (Complex (yr, yi), ierr);
    }
  else if (octave::math::is_integer (alpha))
    {
      // zbesi can overflow as z->0, and cause troubles for generic case below
      alpha = -alpha;
      Complex tmp = zbesi (z, alpha, kode, ierr);
      retval = bessel_return_value (tmp, ierr);
    }
  else
    {
      alpha = -alpha;
 
      Complex tmp = zbesi (z, alpha, kode, ierr);
 
      if (ierr == 0 || ierr == 3)
        {
          Complex tmp2 = (2.0 / M_PI) * sin (M_PI * alpha)
                         * zbesk (z, alpha, kode, ierr);
 
          if (kode == 2)
            {
              // Compensate for different scaling factor of besk.
              tmp2 *= exp (-z - std::abs (z.real ()));
            }
 
          tmp += tmp2;
 
          retval = bessel_return_value (tmp, ierr);
        }
      else
        retval = Complex (numeric_limits<double>::NaN (),
                          numeric_limits<double>::NaN ());
    }
 
  return retval;
}
 
static inline Complex
zbesk (const Complex& z, double alpha, int kode, octave_idx_type& ierr)
{
  Complex retval;
 
  if (alpha >= 0.0)
    {
      double yr = 0.0;
      double yi = 0.0;
 
      F77_INT nz, t_ierr;
 
      double zr = z.real ();
      double zi = z.imag ();
 
      ierr = 0;
 
      if (zr == 0.0 && zi == 0.0)
        {
          yr = numeric_limits<double>::Inf ();
          yi = 0.0;
        }
      else
        {
          F77_FUNC (zbesk, ZBESK) (zr, zi, alpha, kode, 1, &yr, &yi, nz,
                                   t_ierr);
 
          ierr = t_ierr;
 
          if (zi == 0.0 && zr >= 0.0)
            yi = 0.0;
        }
 
      retval = bessel_return_value (Complex (yr, yi), ierr);
    }
  else
    {
      Complex tmp = zbesk (z, -alpha, kode, ierr);
 
      retval = bessel_return_value (tmp, ierr);
    }
 
  return retval;
}
 
static inline Complex
zbesh1 (const Complex& z, double alpha, int kode, octave_idx_type& ierr)
{
  Complex retval;
 
  if (alpha >= 0.0)
    {
      double yr = 0.0;
      double yi = 0.0;
 
      F77_INT nz, t_ierr;
 
      double zr = z.real ();
      double zi = z.imag ();
 
      F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, kode, 1, 1, &yr, &yi, nz,
                               t_ierr);
 
      ierr = t_ierr;
 
      retval = bessel_return_value (Complex (yr, yi), ierr);
    }
  else
    {
      alpha = -alpha;
 
      static const Complex eye = Complex (0.0, 1.0);
 
      Complex tmp = exp (M_PI * alpha * eye) * zbesh1 (z, alpha, kode, ierr);
 
      retval = bessel_return_value (tmp, ierr);
    }
 
  return retval;
}
 
static inline Complex
zbesh2 (const Complex& z, double alpha, int kode, octave_idx_type& ierr)
{
  Complex retval;
 
  if (alpha >= 0.0)
    {
      double yr = 0.0;
      double yi = 0.0;
 
      F77_INT nz, t_ierr;
 
      double zr = z.real ();
      double zi = z.imag ();
 
      F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, kode, 2, 1, &yr, &yi, nz,
                               t_ierr);
 
      ierr = t_ierr;
 
      retval = bessel_return_value (Complex (yr, yi), ierr);
    }
  else
    {
      alpha = -alpha;
 
      static const Complex eye = Complex (0.0, 1.0);
 
      Complex tmp = exp (-M_PI * alpha * eye) * zbesh2 (z, alpha, kode, ierr);
 
      retval = bessel_return_value (tmp, ierr);
    }
 
  return retval;
}
 
typedef Complex (*dptr) (const Complex&, double, int, octave_idx_type&);
 
static inline Complex
do_bessel (dptr f, const char *, double alpha, const Complex& x,
           bool scaled, octave_idx_type& ierr)
{
  Complex retval;
 
  retval = f (x, alpha, (scaled ? 2 : 1), ierr);
 
  return retval;
}
 
static inline ComplexMatrix
do_bessel (dptr f, const char *, double alpha, const ComplexMatrix& x,
           bool scaled, Array<octave_idx_type>& ierr)
{
  octave_idx_type nr = x.rows ();
  octave_idx_type nc = x.cols ();
 
  ComplexMatrix retval (nr, nc);
 
  ierr.resize (dim_vector (nr, nc));
 
  for (octave_idx_type j = 0; j < nc; j++)
    for (octave_idx_type i = 0; i < nr; i++)
      retval(i, j) = f (x(i, j), alpha, (scaled ? 2 : 1), ierr(i, j));
 
  return retval;
}
 
static inline ComplexMatrix
do_bessel (dptr f, const char *, const Matrix& alpha, const Complex& x,
           bool scaled, Array<octave_idx_type>& ierr)
{
  octave_idx_type nr = alpha.rows ();
  octave_idx_type nc = alpha.cols ();
 
  ComplexMatrix retval (nr, nc);
 
  ierr.resize (dim_vector (nr, nc));
 
  for (octave_idx_type j = 0; j < nc; j++)
    for (octave_idx_type i = 0; i < nr; i++)
      retval(i, j) = f (x, alpha(i, j), (scaled ? 2 : 1), ierr(i, j));
 
  return retval;
}
 
static inline ComplexMatrix
do_bessel (dptr f, const char *fn, const Matrix& alpha,
           const ComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr)
{
  ComplexMatrix retval;
 
  octave_idx_type x_nr = x.rows ();
  octave_idx_type x_nc = x.cols ();
 
  octave_idx_type alpha_nr = alpha.rows ();
  octave_idx_type alpha_nc = alpha.cols ();
 
  if (x_nr != alpha_nr || x_nc != alpha_nc)
    (*current_liboctave_error_handler)
      ("%s: the sizes of alpha and x must conform", fn);
 
  octave_idx_type nr = x_nr;
  octave_idx_type nc = x_nc;
 
  retval.resize (nr, nc);
 
  ierr.resize (dim_vector (nr, nc));
 
  for (octave_idx_type j = 0; j < nc; j++)
    for (octave_idx_type i = 0; i < nr; i++)
      retval(i, j) = f (x(i, j), alpha(i, j), (scaled ? 2 : 1), ierr(i, j));
 
  return retval;
}
 
static inline ComplexNDArray
do_bessel (dptr f, const char *, double alpha, const ComplexNDArray& x,
           bool scaled, Array<octave_idx_type>& ierr)
{
  const dim_vector& dv = x.dims ();
  octave_idx_type nel = dv.numel ();
  ComplexNDArray retval (dv);
 
  ierr.resize (dv);
 
  for (octave_idx_type i = 0; i < nel; i++)
    retval(i) = f (x(i), alpha, (scaled ? 2 : 1), ierr(i));
 
  return retval;
}
 
static inline ComplexNDArray
do_bessel (dptr f, const char *, const NDArray& alpha, const Complex& x,
           bool scaled, Array<octave_idx_type>& ierr)
{
  const dim_vector& dv = alpha.dims ();
  octave_idx_type nel = dv.numel ();
  ComplexNDArray retval (dv);
 
  ierr.resize (dv);
 
  for (octave_idx_type i = 0; i < nel; i++)
    retval(i) = f (x, alpha(i), (scaled ? 2 : 1), ierr(i));
 
  return retval;
}
 
static inline ComplexNDArray
do_bessel (dptr f, const char *fn, const NDArray& alpha,
           const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr)
{
  const dim_vector& dv = x.dims ();
  ComplexNDArray retval;
 
  if (dv != alpha.dims ())
    (*current_liboctave_error_handler)
      ("%s: the sizes of alpha and x must conform", fn);
 
  octave_idx_type nel = dv.numel ();
 
  retval.resize (dv);
  ierr.resize (dv);
 
  for (octave_idx_type i = 0; i < nel; i++)
    retval(i) = f (x(i), alpha(i), (scaled ? 2 : 1), ierr(i));
 
  return retval;
}
 
static inline ComplexMatrix
do_bessel (dptr f, const char *, const RowVector& alpha,
           const ComplexColumnVector& x, bool scaled,
           Array<octave_idx_type>& ierr)
{
  octave_idx_type nr = x.numel ();
  octave_idx_type nc = alpha.numel ();
 
  ComplexMatrix retval (nr, nc);
 
  ierr.resize (dim_vector (nr, nc));
 
  for (octave_idx_type j = 0; j < nc; j++)
    for (octave_idx_type i = 0; i < nr; i++)
      retval(i, j) = f (x(i), alpha(j), (scaled ? 2 : 1), ierr(i, j));
 
  return retval;
}
 
#define SS_BESSEL(name, fcn)                                            \
    Complex                                                             \
    name (double alpha, const Complex& x, bool scaled, octave_idx_type& ierr) \
    {                                                                   \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);            \
    }
 
#define SM_BESSEL(name, fcn)                                    \
    ComplexMatrix                                               \
    name (double alpha, const ComplexMatrix& x, bool scaled,    \
          Array<octave_idx_type>& ierr)                         \
    {                                                           \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);    \
    }
 
#define MS_BESSEL(name, fcn)                                    \
    ComplexMatrix                                               \
    name (const Matrix& alpha, const Complex& x, bool scaled,   \
          Array<octave_idx_type>& ierr)                         \
    {                                                           \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);    \
    }
 
#define MM_BESSEL(name, fcn)                                            \
    ComplexMatrix                                                       \
    name (const Matrix& alpha, const ComplexMatrix& x, bool scaled,     \
          Array<octave_idx_type>& ierr)                                 \
    {                                                                   \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);            \
    }
 
#define SN_BESSEL(name, fcn)                                    \
    ComplexNDArray                                              \
    name (double alpha, const ComplexNDArray& x, bool scaled,   \
          Array<octave_idx_type>& ierr)                         \
    {                                                           \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);    \
    }
 
#define NS_BESSEL(name, fcn)                                    \
    ComplexNDArray                                              \
    name (const NDArray& alpha, const Complex& x, bool scaled,  \
          Array<octave_idx_type>& ierr)                         \
    {                                                           \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);    \
    }
 
#define NN_BESSEL(name, fcn)                                            \
    ComplexNDArray                                                      \
    name (const NDArray& alpha, const ComplexNDArray& x, bool scaled,   \
          Array<octave_idx_type>& ierr)                                 \
    {                                                                   \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);            \
    }
 
#define RC_BESSEL(name, fcn)                                            \
    ComplexMatrix                                                       \
    name (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, \
          Array<octave_idx_type>& ierr)                                 \
    {                                                                   \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);            \
    }
 
#define ALL_BESSEL(name, fcn)                   \
    SS_BESSEL (name, fcn)                       \
    SM_BESSEL (name, fcn)                       \
    MS_BESSEL (name, fcn)                       \
    MM_BESSEL (name, fcn)                       \
    SN_BESSEL (name, fcn)                       \
    NS_BESSEL (name, fcn)                       \
    NN_BESSEL (name, fcn)                       \
    RC_BESSEL (name, fcn)
 
ALL_BESSEL (besselj, zbesj)
ALL_BESSEL (bessely, zbesy)
ALL_BESSEL (besseli, zbesi)
ALL_BESSEL (besselk, zbesk)
ALL_BESSEL (besselh1, zbesh1)
ALL_BESSEL (besselh2, zbesh2)
 
#undef ALL_BESSEL
#undef SS_BESSEL
#undef SM_BESSEL
#undef MS_BESSEL
#undef MM_BESSEL
#undef SN_BESSEL
#undef NS_BESSEL
#undef NN_BESSEL
#undef RC_BESSEL
 
static inline FloatComplex
cbesj (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
 
static inline FloatComplex
cbesy (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
 
static inline FloatComplex
cbesi (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
 
static inline FloatComplex
cbesk (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
 
static inline FloatComplex
cbesh1 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
 
static inline FloatComplex
cbesh2 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
 
static inline FloatComplex
cbesj (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr)
{
  FloatComplex retval;
 
  if (alpha >= 0.0)
    {
      FloatComplex y = 0.0;
 
      F77_INT nz, t_ierr;
 
      F77_FUNC (cbesj, CBESJ) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1,
                               F77_CMPLX_ARG (&y), nz, t_ierr);
 
      ierr = t_ierr;
 
      if (z.imag () == 0.0 && z.real () >= 0.0)
        y = FloatComplex (y.real (), 0.0);
 
      retval = bessel_return_value (y, ierr);
    }
  else if (octave::math::is_integer (alpha))
    {
      // zbesy can overflow as z->0, and cause troubles for generic case below
      alpha = -alpha;
      FloatComplex tmp = cbesj (z, alpha, kode, ierr);
      if ((static_cast<long> (alpha)) & 1)
        tmp = - tmp;
      retval = bessel_return_value (tmp, ierr);
    }
  else
    {
      alpha = -alpha;
 
      FloatComplex tmp = cosf (static_cast<float> (M_PI) * alpha)
                         * cbesj (z, alpha, kode, ierr);
 
      if (ierr == 0 || ierr == 3)
        {
          tmp -= sinf (static_cast<float> (M_PI) * alpha)
                 * cbesy (z, alpha, kode, ierr);
 
          retval = bessel_return_value (tmp, ierr);
        }
      else
        retval = FloatComplex (numeric_limits<float>::NaN (),
                               numeric_limits<float>::NaN ());
    }
 
  return retval;
}
 
static inline FloatComplex
cbesy (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr)
{
  FloatComplex retval;
 
  if (alpha >= 0.0)
    {
      FloatComplex y = 0.0;
 
      F77_INT nz, t_ierr;
 
      FloatComplex w;
 
      ierr = 0;
 
      if (z.real () == 0.0 && z.imag () == 0.0)
        {
          y = FloatComplex (-numeric_limits<float>::Inf (), 0.0);
        }
      else
        {
          F77_FUNC (cbesy, CBESY) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1,
                                   F77_CMPLX_ARG (&y), nz,
                                   F77_CMPLX_ARG (&w), t_ierr);
 
          ierr = t_ierr;
 
          if (z.imag () == 0.0 && z.real () >= 0.0)
            y = FloatComplex (y.real (), 0.0);
        }
 
      return bessel_return_value (y, ierr);
    }
  else if (octave::math::is_integer (alpha - 0.5))
    {
      // zbesy can overflow as z->0, and cause troubles for generic case below
      alpha = -alpha;
      FloatComplex tmp = cbesj (z, alpha, kode, ierr);
      if ((static_cast<long> (alpha - 0.5)) & 1)
        tmp = - tmp;
      retval = bessel_return_value (tmp, ierr);
    }
  else
    {
      alpha = -alpha;
 
      FloatComplex tmp = cosf (static_cast<float> (M_PI) * alpha)
                         * cbesy (z, alpha, kode, ierr);
 
      if (ierr == 0 || ierr == 3)
        {
          tmp += sinf (static_cast<float> (M_PI) * alpha)
                 * cbesj (z, alpha, kode, ierr);
 
          retval = bessel_return_value (tmp, ierr);
        }
      else
        retval = FloatComplex (numeric_limits<float>::NaN (),
                               numeric_limits<float>::NaN ());
    }
 
  return retval;
}
 
static inline FloatComplex
cbesi (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr)
{
  FloatComplex retval;
 
  if (alpha >= 0.0)
    {
      FloatComplex y = 0.0;
 
      F77_INT nz, t_ierr;
 
      F77_FUNC (cbesi, CBESI) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1,
                               F77_CMPLX_ARG (&y), nz, t_ierr);
 
      ierr = t_ierr;
 
      if (z.imag () == 0.0 && z.real () >= 0.0)
        y = FloatComplex (y.real (), 0.0);
 
      retval = bessel_return_value (y, ierr);
    }
  else
    {
      alpha = -alpha;
 
      FloatComplex tmp = cbesi (z, alpha, kode, ierr);
 
      if (ierr == 0 || ierr == 3)
        {
          FloatComplex tmp2 = static_cast<float> (2.0 / M_PI)
                              * sinf (static_cast<float> (M_PI) * alpha)
                              * cbesk (z, alpha, kode, ierr);
 
          if (kode == 2)
            {
              // Compensate for different scaling factor of besk.
              tmp2 *= exp (-z - std::abs (z.real ()));
            }
 
          tmp += tmp2;
 
          retval = bessel_return_value (tmp, ierr);
        }
      else
        retval = FloatComplex (numeric_limits<float>::NaN (),
                               numeric_limits<float>::NaN ());
    }
 
  return retval;
}
 
static inline FloatComplex
cbesk (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr)
{
  FloatComplex retval;
 
  if (alpha >= 0.0)
    {
      FloatComplex y = 0.0;
 
      F77_INT nz, t_ierr;
 
      ierr = 0;
 
      if (z.real () == 0.0 && z.imag () == 0.0)
        {
          y = FloatComplex (numeric_limits<float>::Inf (), 0.0);
        }
      else
        {
          F77_FUNC (cbesk, CBESK) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1,
                                   F77_CMPLX_ARG (&y), nz, t_ierr);
 
          ierr = t_ierr;
 
          if (z.imag () == 0.0 && z.real () >= 0.0)
            y = FloatComplex (y.real (), 0.0);
        }
 
      retval = bessel_return_value (y, ierr);
    }
  else
    {
      FloatComplex tmp = cbesk (z, -alpha, kode, ierr);
 
      retval = bessel_return_value (tmp, ierr);
    }
 
  return retval;
}
 
static inline FloatComplex
cbesh1 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr)
{
  FloatComplex retval;
 
  if (alpha >= 0.0)
    {
      FloatComplex y = 0.0;
 
      F77_INT nz, t_ierr;
 
      F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, 1,
                               F77_CMPLX_ARG (&y), nz, t_ierr);
 
      ierr = t_ierr;
 
      retval = bessel_return_value (y, ierr);
    }
  else
    {
      alpha = -alpha;
 
      static const FloatComplex eye = FloatComplex (0.0, 1.0);
 
      FloatComplex tmp = exp (static_cast<float> (M_PI) * alpha * eye)
                         * cbesh1 (z, alpha, kode, ierr);
 
      retval = bessel_return_value (tmp, ierr);
    }
 
  return retval;
}
 
static inline FloatComplex
cbesh2 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr)
{
  FloatComplex retval;
 
  if (alpha >= 0.0)
    {
      FloatComplex y = 0.0;
 
      F77_INT nz, t_ierr;
 
      F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 2, 1,
                               F77_CMPLX_ARG (&y), nz, t_ierr);
 
      ierr = t_ierr;
 
      retval = bessel_return_value (y, ierr);
    }
  else
    {
      alpha = -alpha;
 
      static const FloatComplex eye = FloatComplex (0.0, 1.0);
 
      FloatComplex tmp = exp (-static_cast<float> (M_PI) * alpha * eye)
                         * cbesh2 (z, alpha, kode, ierr);
 
      retval = bessel_return_value (tmp, ierr);
    }
 
  return retval;
}
 
typedef FloatComplex (*fptr) (const FloatComplex&, float, int,
                              octave_idx_type&);
 
static inline FloatComplex
do_bessel (fptr f, const char *, float alpha, const FloatComplex& x,
           bool scaled, octave_idx_type& ierr)
{
  FloatComplex retval;
 
  retval = f (x, alpha, (scaled ? 2 : 1), ierr);
 
  return retval;
}
 
static inline FloatComplexMatrix
do_bessel (fptr f, const char *, float alpha, const FloatComplexMatrix& x,
           bool scaled, Array<octave_idx_type>& ierr)
{
  octave_idx_type nr = x.rows ();
  octave_idx_type nc = x.cols ();
 
  FloatComplexMatrix retval (nr, nc);
 
  ierr.resize (dim_vector (nr, nc));
 
  for (octave_idx_type j = 0; j < nc; j++)
    for (octave_idx_type i = 0; i < nr; i++)
      retval(i, j) = f (x(i, j), alpha, (scaled ? 2 : 1), ierr(i, j));
 
  return retval;
}
 
static inline FloatComplexMatrix
do_bessel (fptr f, const char *, const FloatMatrix& alpha,
           const FloatComplex& x,
           bool scaled, Array<octave_idx_type>& ierr)
{
  octave_idx_type nr = alpha.rows ();
  octave_idx_type nc = alpha.cols ();
 
  FloatComplexMatrix retval (nr, nc);
 
  ierr.resize (dim_vector (nr, nc));
 
  for (octave_idx_type j = 0; j < nc; j++)
    for (octave_idx_type i = 0; i < nr; i++)
      retval(i, j) = f (x, alpha(i, j), (scaled ? 2 : 1), ierr(i, j));
 
  return retval;
}
 
static inline FloatComplexMatrix
do_bessel (fptr f, const char *fn, const FloatMatrix& alpha,
           const FloatComplexMatrix& x, bool scaled,
           Array<octave_idx_type>& ierr)
{
  FloatComplexMatrix retval;
 
  octave_idx_type x_nr = x.rows ();
  octave_idx_type x_nc = x.cols ();
 
  octave_idx_type alpha_nr = alpha.rows ();
  octave_idx_type alpha_nc = alpha.cols ();
 
  if (x_nr != alpha_nr || x_nc != alpha_nc)
    (*current_liboctave_error_handler)
      ("%s: the sizes of alpha and x must conform", fn);
 
  octave_idx_type nr = x_nr;
  octave_idx_type nc = x_nc;
 
  retval.resize (nr, nc);
 
  ierr.resize (dim_vector (nr, nc));
 
  for (octave_idx_type j = 0; j < nc; j++)
    for (octave_idx_type i = 0; i < nr; i++)
      retval(i, j) = f (x(i, j), alpha(i, j), (scaled ? 2 : 1), ierr(i, j));
 
  return retval;
}
 
static inline FloatComplexNDArray
do_bessel (fptr f, const char *, float alpha, const FloatComplexNDArray& x,
           bool scaled, Array<octave_idx_type>& ierr)
{
  const dim_vector& dv = x.dims ();
  octave_idx_type nel = dv.numel ();
  FloatComplexNDArray retval (dv);
 
  ierr.resize (dv);
 
  for (octave_idx_type i = 0; i < nel; i++)
    retval(i) = f (x(i), alpha, (scaled ? 2 : 1), ierr(i));
 
  return retval;
}
 
static inline FloatComplexNDArray
do_bessel (fptr f, const char *, const FloatNDArray& alpha,
           const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr)
{
  const dim_vector& dv = alpha.dims ();
  octave_idx_type nel = dv.numel ();
  FloatComplexNDArray retval (dv);
 
  ierr.resize (dv);
 
  for (octave_idx_type i = 0; i < nel; i++)
    retval(i) = f (x, alpha(i), (scaled ? 2 : 1), ierr(i));
 
  return retval;
}
 
static inline FloatComplexNDArray
do_bessel (fptr f, const char *fn, const FloatNDArray& alpha,
           const FloatComplexNDArray& x, bool scaled,
           Array<octave_idx_type>& ierr)
{
  const dim_vector& dv = x.dims ();
  FloatComplexNDArray retval;
 
  if (dv != alpha.dims ())
    (*current_liboctave_error_handler)
      ("%s: the sizes of alpha and x must conform", fn);
 
  octave_idx_type nel = dv.numel ();
 
  retval.resize (dv);
  ierr.resize (dv);
 
  for (octave_idx_type i = 0; i < nel; i++)
    retval(i) = f (x(i), alpha(i), (scaled ? 2 : 1), ierr(i));
 
  return retval;
}
 
static inline FloatComplexMatrix
do_bessel (fptr f, const char *, const FloatRowVector& alpha,
           const FloatComplexColumnVector& x, bool scaled,
           Array<octave_idx_type>& ierr)
{
  octave_idx_type nr = x.numel ();
  octave_idx_type nc = alpha.numel ();
 
  FloatComplexMatrix retval (nr, nc);
 
  ierr.resize (dim_vector (nr, nc));
 
  for (octave_idx_type j = 0; j < nc; j++)
    for (octave_idx_type i = 0; i < nr; i++)
      retval(i, j) = f (x(i), alpha(j), (scaled ? 2 : 1), ierr(i, j));
 
  return retval;
}
 
#define SS_BESSEL(name, fcn)                                    \
    FloatComplex                                                \
    name (float alpha, const FloatComplex& x, bool scaled,      \
          octave_idx_type& ierr)                                \
    {                                                           \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);    \
    }
 
#define SM_BESSEL(name, fcn)                                            \
    FloatComplexMatrix                                                  \
    name (float alpha, const FloatComplexMatrix& x, bool scaled,        \
          Array<octave_idx_type>& ierr)                                 \
    {                                                                   \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);            \
    }
 
#define MS_BESSEL(name, fcn)                                            \
    FloatComplexMatrix                                                  \
    name (const FloatMatrix& alpha, const FloatComplex& x, bool scaled, \
          Array<octave_idx_type>& ierr)                                 \
    {                                                                   \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);            \
    }
 
#define MM_BESSEL(name, fcn)                                            \
    FloatComplexMatrix                                                  \
    name (const FloatMatrix& alpha, const FloatComplexMatrix& x,        \
          bool scaled, Array<octave_idx_type>& ierr)                    \
    {                                                                   \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);            \
    }
 
#define SN_BESSEL(name, fcn)                                            \
    FloatComplexNDArray                                                 \
    name (float alpha, const FloatComplexNDArray& x, bool scaled,       \
          Array<octave_idx_type>& ierr)                                 \
    {                                                                   \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);            \
    }
 
#define NS_BESSEL(name, fcn)                                    \
    FloatComplexNDArray                                         \
    name (const FloatNDArray& alpha, const FloatComplex& x,     \
          bool scaled, Array<octave_idx_type>& ierr)            \
    {                                                           \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);    \
    }
 
#define NN_BESSEL(name, fcn)                                            \
    FloatComplexNDArray                                                 \
    name (const FloatNDArray& alpha, const FloatComplexNDArray& x,      \
          bool scaled, Array<octave_idx_type>& ierr)                    \
    {                                                                   \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);            \
    }
 
#define RC_BESSEL(name, fcn)                                    \
    FloatComplexMatrix                                          \
    name (const FloatRowVector& alpha,                          \
          const FloatComplexColumnVector& x, bool scaled,       \
          Array<octave_idx_type>& ierr)                         \
    {                                                           \
      return do_bessel (fcn, #name, alpha, x, scaled, ierr);    \
    }
 
#define ALL_BESSEL(name, fcn)                   \
    SS_BESSEL (name, fcn)                       \
    SM_BESSEL (name, fcn)                       \
    MS_BESSEL (name, fcn)                       \
    MM_BESSEL (name, fcn)                       \
    SN_BESSEL (name, fcn)                       \
    NS_BESSEL (name, fcn)                       \
    NN_BESSEL (name, fcn)                       \
    RC_BESSEL (name, fcn)
 
ALL_BESSEL (besselj, cbesj)
ALL_BESSEL (bessely, cbesy)
ALL_BESSEL (besseli, cbesi)
ALL_BESSEL (besselk, cbesk)
ALL_BESSEL (besselh1, cbesh1)
ALL_BESSEL (besselh2, cbesh2)
 
#undef ALL_BESSEL
#undef SS_BESSEL
#undef SM_BESSEL
#undef MS_BESSEL
#undef MM_BESSEL
#undef SN_BESSEL
#undef NS_BESSEL
#undef NN_BESSEL
#undef RC_BESSEL
 
Complex
biry (const Complex& z, bool deriv, bool scaled, octave_idx_type& ierr)
{
  double ar = 0.0;
  double ai = 0.0;
 
  double zr = z.real ();
  double zi = z.imag ();
 
  F77_INT id = (deriv ? 1 : 0);
  F77_INT t_ierr;
  F77_INT sc = (scaled ? 2 : 1);
 
  F77_FUNC (zbiry, ZBIRY) (zr, zi, id, sc, ar, ai, t_ierr);
 
  ierr = t_ierr;
 
  if (zi == 0.0 && (! scaled || zr >= 0.0))
    ai = 0.0;
 
  return bessel_return_value (Complex (ar, ai), ierr);
}
 
ComplexMatrix
biry (const ComplexMatrix& z, bool deriv, bool scaled,
      Array<octave_idx_type>& ierr)
{
  octave_idx_type nr = z.rows ();
  octave_idx_type nc = z.cols ();
 
  ComplexMatrix retval (nr, nc);
 
  ierr.resize (dim_vector (nr, nc));
 
  for (octave_idx_type j = 0; j < nc; j++)
    for (octave_idx_type i = 0; i < nr; i++)
      retval(i, j) = biry (z(i, j), deriv, scaled, ierr(i, j));
 
  return retval;
}
 
ComplexNDArray
biry (const ComplexNDArray& z, bool deriv, bool scaled,
      Array<octave_idx_type>& ierr)
{
  const dim_vector& dv = z.dims ();
  octave_idx_type nel = dv.numel ();
  ComplexNDArray retval (dv);
 
  ierr.resize (dv);
 
  for (octave_idx_type i = 0; i < nel; i++)
    retval(i) = biry (z(i), deriv, scaled, ierr(i));
 
  return retval;
}
 
FloatComplex
biry (const FloatComplex& z, bool deriv, bool scaled,
      octave_idx_type& ierr)
{
  FloatComplex a;
 
  F77_INT id = (deriv ? 1 : 0);
  F77_INT t_ierr;
  F77_INT sc = (scaled ? 2 : 1);
 
  F77_FUNC (cbiry, CBIRY) (F77_CONST_CMPLX_ARG (&z), id, sc,
                           F77_CMPLX_ARG (&a), t_ierr);
 
  ierr = t_ierr;
 
  float ar = a.real ();
  float ai = a.imag ();
 
  if (z.imag () == 0.0 && (! scaled || z.real () >= 0.0))
    ai = 0.0;
 
  return bessel_return_value (FloatComplex (ar, ai), ierr);
}
 
FloatComplexMatrix
biry (const FloatComplexMatrix& z, bool deriv, bool scaled,
      Array<octave_idx_type>& ierr)
{
  octave_idx_type nr = z.rows ();
  octave_idx_type nc = z.cols ();
 
  FloatComplexMatrix retval (nr, nc);
 
  ierr.resize (dim_vector (nr, nc));
 
  for (octave_idx_type j = 0; j < nc; j++)
    for (octave_idx_type i = 0; i < nr; i++)
      retval(i, j) = biry (z(i, j), deriv, scaled, ierr(i, j));
 
  return retval;
}
 
FloatComplexNDArray
biry (const FloatComplexNDArray& z, bool deriv, bool scaled,
      Array<octave_idx_type>& ierr)
{
  const dim_vector& dv = z.dims ();
  octave_idx_type nel = dv.numel ();
  FloatComplexNDArray retval (dv);
 
  ierr.resize (dv);
 
  for (octave_idx_type i = 0; i < nel; i++)
    retval(i) = biry (z(i), deriv, scaled, ierr(i));
 
  return retval;
}
 
// Real and complex Dawson function (= scaled erfi) from Faddeeva package
double
dawson (double x) { return Faddeeva::Dawson (x); }
float
dawson (float x) { return Faddeeva::Dawson (x); }
 
Complex
dawson (const Complex& x)
{
  return Faddeeva::Dawson (x);
}
 
FloatComplex
dawson (const FloatComplex& x)
{
  Complex xd (x.real (), x.imag ());
  Complex ret;
  ret = Faddeeva::Dawson (xd, std::numeric_limits<float>::epsilon ());
  return FloatComplex (ret.real (), ret.imag ());
}
 
void
ellipj (double u, double m, double& sn, double& cn, double& dn, double& err)
{
  static const int Nmax = 16;
  double m1, t=0, si_u, co_u, se_u, ta_u, b, c[Nmax], a[Nmax], phi;
  int n, Nn, ii;
 
  if (m < 0 || m > 1)
    {
      (*current_liboctave_warning_with_id_handler)
        ("Octave:ellipj-invalid-m",
         "ellipj: invalid M value, required value 0 <= M <= 1");
 
      sn = cn = dn = lo_ieee_nan_value ();
 
      return;
    }
 
  double sqrt_eps = std::sqrt (std::numeric_limits<double>::epsilon ());
  if (m < sqrt_eps)
    {
      // For small m, (Abramowitz and Stegun, Section 16.13)
      si_u = sin (u);
      co_u = cos (u);
      t = 0.25*m*(u - si_u*co_u);
      sn = si_u - t * co_u;
      cn = co_u + t * si_u;
      dn = 1 - 0.5*m*si_u*si_u;
    }
  else if ((1 - m) < sqrt_eps)
    {
      // For m1 = (1-m) small (Abramowitz and Stegun, Section 16.15)
      m1 = 1 - m;
      si_u = sinh (u);
      co_u = cosh (u);
      ta_u = tanh (u);
      se_u = 1/co_u;
      sn = ta_u + 0.25*m1*(si_u*co_u - u)*se_u*se_u;
      cn = se_u - 0.25*m1*(si_u*co_u - u)*ta_u*se_u;
      dn = se_u + 0.25*m1*(si_u*co_u + u)*ta_u*se_u;
    }
  else
    {
      // Arithmetic-Geometric Mean (AGM) algorithm
      //   (Abramowitz and Stegun, Section 16.4)
      a[0] = 1;
      b    = std::sqrt (1 - m);
      c[0] = std::sqrt (m);
      for (n = 1; n < Nmax; ++n)
        {
          a[n] = (a[n - 1] + b)/2;
          c[n] = (a[n - 1] - b)/2;
          b = std::sqrt (a[n - 1]*b);
          if (c[n]/a[n] < std::numeric_limits<double>::epsilon ()) break;
        }
      if (n >= Nmax - 1)
        {
          err = 1;
          return;
        }
      Nn = n;
      for (ii = 1; n > 0; ii *= 2, --n) {}  // ii = pow(2,Nn)
      phi = ii*a[Nn]*u;
      for (n = Nn; n > 0; --n)
        {
          phi = (std::asin ((c[n]/a[n])* sin (phi)) + phi)/2;
        }
      sn = sin (phi);
      cn = cos (phi);
      dn = std::sqrt (1 - m*sn*sn);
    }
}
 
void
ellipj (const Complex& u, double m, Complex& sn, Complex& cn, Complex& dn,
        double& err)
{
  double m1 = 1 - m, ss1, cc1, dd1;
 
  ellipj (u.imag (), m1, ss1, cc1, dd1, err);
  if (u.real () == 0)
    {
      // u is pure imag: Jacoby imag. transf.
      sn = Complex (0, ss1/cc1);
      cn = 1/cc1;         //    cn.imag = 0;
      dn = dd1/cc1;       //    dn.imag = 0;
    }
  else
    {
      // u is generic complex
      double ss, cc, dd, ddd;
 
      ellipj (u.real (), m, ss, cc, dd, err);
      ddd = cc1*cc1 + m*ss*ss*ss1*ss1;
      sn = Complex (ss*dd1/ddd, cc*dd*ss1*cc1/ddd);
      cn = Complex (cc*cc1/ddd, -ss*dd*ss1*dd1/ddd);
      dn = Complex (dd*cc1*dd1/ddd, -m*ss*cc*ss1/ddd);
    }
}
 
// Complex error function from the Faddeeva package
Complex
erf (const Complex& x)
{
  return Faddeeva::erf (x);
}
 
FloatComplex
erf (const FloatComplex& x)
{
  Complex xd (x.real (), x.imag ());
  Complex ret = Faddeeva::erf (xd, std::numeric_limits<float>::epsilon ());
  return FloatComplex (ret.real (), ret.imag ());
}
 
// Complex complementary error function from the Faddeeva package
Complex
erfc (const Complex& x)
{
  return Faddeeva::erfc (x);
}
 
FloatComplex
erfc (const FloatComplex& x)
{
  Complex xd (x.real (), x.imag ());
  Complex ret = Faddeeva::erfc (xd, std::numeric_limits<float>::epsilon ());
  return FloatComplex (ret.real (), ret.imag ());
}
 
// The algorithm for erfcinv is an adaptation of the erfinv algorithm
// above from P. J. Acklam.  It has been modified to run over the
// different input domain of erfcinv.  See the notes for erfinv for an
// explanation.
 
static double
do_erfcinv (double x, bool refine)
{
  // Coefficients of rational approximation.
  static const double a[] =
  {
    -2.806989788730439e+01,  1.562324844726888e+02,
      -1.951109208597547e+02,  9.783370457507161e+01,
      -2.168328665628878e+01,  1.772453852905383e+00
    };
  static const double b[] =
  {
    -5.447609879822406e+01,  1.615858368580409e+02,
      -1.556989798598866e+02,  6.680131188771972e+01,
      -1.328068155288572e+01
    };
  static const double c[] =
  {
    -5.504751339936943e-03, -2.279687217114118e-01,
      -1.697592457770869e+00, -1.802933168781950e+00,
      3.093354679843505e+00,  2.077595676404383e+00
    };
  static const double d[] =
  {
    7.784695709041462e-03,  3.224671290700398e-01,
    2.445134137142996e+00,  3.754408661907416e+00
  };
 
  static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2.
  static const double pbreak_lo = 0.04850;  // 1-pbreak
  static const double pbreak_hi = 1.95150;  // 1+pbreak
  double y;
 
  // Select case.
  if (x >= pbreak_lo && x <= pbreak_hi)
    {
      // Middle region.
      const double q = 0.5*(1-x), r = q*q;
      const double yn = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q;
      const double yd = ((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0;
      y = yn / yd;
    }
  else if (x > 0.0 && x < 2.0)
    {
      // Tail region.
      const double q = (x < 1
                        ? std::sqrt (-2*std::log (0.5*x))
                        : std::sqrt (-2*std::log (0.5*(2-x))));
 
      const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5];
 
      const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0;
 
      y = yn / yd;
 
      if (x < pbreak_lo)
        y = -y;
    }
  else if (x == 0.0)
    return numeric_limits<double>::Inf ();
  else if (x == 2.0)
    return -numeric_limits<double>::Inf ();
  else
    return numeric_limits<double>::NaN ();
 
  if (refine)
    {
      // One iteration of Halley's method gives full precision.
      double u = (erf (y) - (1-x)) * spi2 * exp (y*y);
      y -= u / (1 + y*u);
    }
 
  return y;
}
 
double
erfcinv (double x)
{
  return do_erfcinv (x, true);
}
 
float
erfcinv (float x)
{
  return do_erfcinv (x, false);
}
 
// Real and complex scaled complementary error function from Faddeeva pkg.
double
erfcx (double x) { return Faddeeva::erfcx (x); }
float
erfcx (float x) { return Faddeeva::erfcx (x); }
 
Complex
erfcx (const Complex& x)
{
  return Faddeeva::erfcx (x);
}
 
FloatComplex
erfcx (const FloatComplex& x)
{
  Complex xd (x.real (), x.imag ());
  Complex ret;
  ret = Faddeeva::erfcx (xd, std::numeric_limits<float>::epsilon ());
  return FloatComplex (ret.real (), ret.imag ());
}
 
// Real and complex imaginary error function from Faddeeva package
double
erfi (double x) { return Faddeeva::erfi (x); }
float
erfi (float x) { return Faddeeva::erfi (x); }
 
Complex
erfi (const Complex& x)
{
  return Faddeeva::erfi (x);
}
 
FloatComplex
erfi (const FloatComplex& x)
{
  Complex xd (x.real (), x.imag ());
  Complex ret = Faddeeva::erfi (xd, std::numeric_limits<float>::epsilon ());
  return FloatComplex (ret.real (), ret.imag ());
}
 
// This algorithm is due to P. J. Acklam.
//
// See http://home.online.no/~pjacklam/notes/invnorm/
//
// The rational approximation has relative accuracy 1.15e-9 in the whole
// region.  For doubles, it is refined by a single step of Halley's 3rd
// order method.  For single precision, the accuracy is already OK, so
// we skip it to get faster evaluation.
 
static double
do_erfinv (double x, bool refine)
{
  // Coefficients of rational approximation.
  static const double a[] =
  {
    -2.806989788730439e+01,  1.562324844726888e+02,
      -1.951109208597547e+02,  9.783370457507161e+01,
      -2.168328665628878e+01,  1.772453852905383e+00
    };
  static const double b[] =
  {
    -5.447609879822406e+01,  1.615858368580409e+02,
      -1.556989798598866e+02,  6.680131188771972e+01,
      -1.328068155288572e+01
    };
  static const double c[] =
  {
    -5.504751339936943e-03, -2.279687217114118e-01,
      -1.697592457770869e+00, -1.802933168781950e+00,
      3.093354679843505e+00,  2.077595676404383e+00
    };
  static const double d[] =
  {
    7.784695709041462e-03,  3.224671290700398e-01,
    2.445134137142996e+00,  3.754408661907416e+00
  };
 
  static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2.
  static const double pbreak = 0.95150;
  double ax = fabs (x), y;
 
  // Select case.
  if (ax <= pbreak)
    {
      // Middle region.
      const double q = 0.5 * x, r = q*q;
      const double yn = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q;
      const double yd = ((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0;
      y = yn / yd;
    }
  else if (ax < 1.0)
    {
      // Tail region.
      const double q = std::sqrt (-2*std::log (0.5*(1-ax)));
      const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5];
      const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0;
      y = yn / yd * math::signum (-x);
    }
  else if (ax == 1.0)
    return numeric_limits<double>::Inf () * math::signum (x);
  else
    return numeric_limits<double>::NaN ();
 
  if (refine)
    {
      // One iteration of Halley's method gives full precision.
      double u = (erf (y) - x) * spi2 * exp (y*y);
      y -= u / (1 + y*u);
    }
 
  return y;
}
 
double
erfinv (double x)
{
  return do_erfinv (x, true);
}
 
float
erfinv (float x)
{
  return do_erfinv (x, false);
}
 
Complex
expm1 (const Complex& x)
{
  Complex retval;
 
  if (std::abs (x) < 1)
    {
      double im = x.imag ();
      double u = expm1 (x.real ());
      double v = sin (im/2);
      v = -2*v*v;
      retval = Complex (u*v + u + v, (u+1) * sin (im));
    }
  else
    retval = std::exp (x) - Complex (1);
 
  return retval;
}
 
FloatComplex
expm1 (const FloatComplex& x)
{
  FloatComplex retval;
 
  if (std::abs (x) < 1)
    {
      float im = x.imag ();
      float u = expm1 (x.real ());
      float v = sin (im/2);
      v = -2*v*v;
      retval = FloatComplex (u*v + u + v, (u+1) * sin (im));
    }
  else
    retval = std::exp (x) - FloatComplex (1);
 
  return retval;
}
 
double
gamma (double x)
{
  double result;
 
  // Special cases for (near) compatibility with Matlab instead of tgamma.
  // Matlab does not have -0.
 
  if (x == 0)
    result = (math::negative_sign (x)
              ? -numeric_limits<double>::Inf ()
              : numeric_limits<double>::Inf ());
  else if ((x < 0 && math::is_integer (x))
           || math::isinf (x))
    result = numeric_limits<double>::Inf ();
  else if (math::isnan (x))
    result = numeric_limits<double>::NaN ();
  else
    result = std::tgamma (x);
 
  return result;
}
 
float
gamma (float x)
{
  float result;
 
  // Special cases for (near) compatibility with Matlab instead of tgamma.
  // Matlab does not have -0.
 
  if (x == 0)
    result = (math::negative_sign (x)
              ? -numeric_limits<float>::Inf ()
              : numeric_limits<float>::Inf ());
  else if ((x < 0 && math::is_integer (x))
           || math::isinf (x))
    result = numeric_limits<float>::Inf ();
  else if (math::isnan (x))
    result = numeric_limits<float>::NaN ();
  else
    result = std::tgammaf (x);
 
  return result;
}
 
Complex
log1p (const Complex& x)
{
  Complex retval;
 
  double r = x.real (), i = x.imag ();
 
  if (fabs (r) < 0.5 && fabs (i) < 0.5)
    {
      double u = 2*r + r*r + i*i;
      retval = Complex (log1p (u / (1+std::sqrt (u+1))),
                        atan2 (i, 1 + r));
    }
  else
    retval = std::log (Complex (1) + x);
 
  return retval;
}
 
FloatComplex
log1p (const FloatComplex& x)
{
  FloatComplex retval;
 
  float r = x.real (), i = x.imag ();
 
  if (fabs (r) < 0.5 && fabs (i) < 0.5)
    {
      float u = 2*r + r*r + i*i;
      retval = FloatComplex (log1p (u / (1+std::sqrt (u+1))),
                             atan2 (i, 1 + r));
    }
  else
    retval = std::log (FloatComplex (1) + x);
 
  return retval;
}
 
static const double pi = 3.14159265358979323846;
 
template <typename T>
static inline T
xlog (const T& x)
{
  return log (x);
}
 
template <>
inline double
xlog (const double& x)
{
  return std::log (x);
}
 
template <>
inline float
xlog (const float& x)
{
  return std::log (x);
}
 
template <typename T>
static T
lanczos_approximation_psi (const T zc)
{
  // Coefficients for C.Lanczos expansion of psi function from XLiFE++
  // gammaFunctions psi_coef[k] = - (2k+1) * lg_coef[k] (see melina++
  // gamma functions -1/12, 3/360,-5/1260, 7/1680,-9/1188,
  // 11*691/360360,-13/156, 15*3617/122400, ? , ?
  static const T dg_coeff[10] =
  {
    -0.83333333333333333e-1, 0.83333333333333333e-2,
      -0.39682539682539683e-2, 0.41666666666666667e-2,
      -0.75757575757575758e-2, 0.21092796092796093e-1,
      -0.83333333333333333e-1, 0.4432598039215686,
      -0.3053954330270122e+1,  0.125318899521531e+2
    };
 
  T overz2  = T (1.0) / (zc * zc);
  T overz2k = overz2;
 
  T p = 0;
  for (octave_idx_type k = 0; k < 10; k++, overz2k *= overz2)
    p += dg_coeff[k] * overz2k;
  p += xlog (zc) - T (0.5) / zc;
  return p;
}
 
template <typename T>
T
xpsi (T z)
{
  static const double euler_mascheroni
    = 0.577215664901532860606512090082402431042;
 
  const bool is_int = (std::floor (z) == z);
 
  T p = 0;
  if (z <= 0)
    {
      // limits - zeros of the gamma function
      if (is_int)
        p = -numeric_limits<T>::Inf (); // Matlab returns -Inf for psi (0)
      else
        // Abramowitz and Stegun, page 259, eq 6.3.7
        p = psi (1 - z) - (pi / tan (pi * z));
    }
  else if (is_int)
    {
      // Abramowitz and Stegun, page 258, eq 6.3.2
      p = - euler_mascheroni;
      for (octave_idx_type k = z - 1; k > 0; k--)
        p += 1.0 / k;
    }
  else if (z - std::floor (z) == 0.5)
    {
      // Abramowitz and Stegun, page 258, eq 6.3.3 and 6.3.4
      for (octave_idx_type k = z; k > 0; k--)
        p += 1.0 / (2 * k - 1);
 
      p = - euler_mascheroni - 2 * std::log (2) + 2 * (p);
    }
  else
    {
      // adapted from XLiFE++ gammaFunctions
 
      T zc = z;
      // Use formula for derivative of LogGamma(z)
      if (z < 10)
        {
          const signed char n = 10 - z;
          for (signed char k = n - 1; k >= 0; k--)
            p -= 1.0 / (k + z);
          zc += n;
        }
      p += lanczos_approximation_psi (zc);
    }
 
  return p;
}
 
// explicit instantiations
double
psi (double z) { return xpsi (z); }
float
psi (float z) { return xpsi (z); }
 
template <typename T>
std::complex<T>
xpsi (const std::complex<T>& z)
{
  // adapted from XLiFE++ gammaFunctions
 
  typedef typename std::complex<T>::value_type P;
 
  P z_r  = z.real ();
  P z_ra = z_r;
 
  std::complex<T> dgam (0.0, 0.0);
  if (z.imag () == 0)
    dgam = std::complex<T> (psi (z_r), 0.0);
  else if (z_r < 0)
    dgam = psi (P (1.0) - z)- (P (pi) / tan (P (pi) * z));
  else
    {
      // Use formula for derivative of LogGamma(z)
      std::complex<T> z_m = z;
      if (z_ra < 8)
        {
          unsigned char n = 8 - z_ra;
          z_m = z + std::complex<T> (n, 0.0);
 
          // Recurrence formula.  For | Re(z) | < 8, use recursively
          //
          //   DiGamma(z) = DiGamma(z+1) - 1/z
          std::complex<T> z_p = z + P (n - 1);
          for (unsigned char k = n; k > 0; k--, z_p -= 1.0)
            dgam -= P (1.0) / z_p;
        }
 
      // for | Re(z) | > 8, use derivative of C.Lanczos expansion for
      // LogGamma
      //
      //   psi(z) = log(z) - 1/(2z) - 1/12z^2 + 3/360z^4 - 5/1260z^6
      //     + 7/1680z^8 - 9/1188z^10 + ...
      //
      // (Abramowitz&Stegun, page 259, formula 6.3.18
      dgam += lanczos_approximation_psi (z_m);
    }
  return dgam;
}
 
// explicit instantiations
Complex
psi (const Complex& z) { return xpsi (z); }
FloatComplex
psi (const FloatComplex& z) { return xpsi (z); }
 
template <typename T>
static inline void
fortran_psifn (T z, octave_idx_type n, T& ans, octave_idx_type& ierr);
 
template <>
inline void
fortran_psifn<double> (double z, octave_idx_type n_arg,
                       double& ans, octave_idx_type& ierr)
{
  F77_INT n = to_f77_int (n_arg);
  F77_INT flag = 0;
  F77_INT t_ierr;
  F77_XFCN (dpsifn, DPSIFN, (z, n, 1, 1, ans, flag, t_ierr));
  ierr = t_ierr;
}
 
template <>
inline void
fortran_psifn<float> (float z, octave_idx_type n_arg,
                      float& ans, octave_idx_type& ierr)
{
  F77_INT n = to_f77_int (n_arg);
  F77_INT flag = 0;
  F77_INT t_ierr;
  F77_XFCN (psifn, PSIFN, (z, n, 1, 1, ans, flag, t_ierr));
  ierr = t_ierr;
}
 
template <typename T>
T
xpsi (octave_idx_type n, T z)
{
  T ans;
  octave_idx_type ierr = 0;
  fortran_psifn<T> (z, n, ans, ierr);
  if (ierr == 0)
    {
      // Remember that psifn and dpsifn return scales values
      // When n is 1: do nothing since ((-1)**(n+1)/gamma(n+1)) == 1
      // When n is 0: change sign since ((-1)**(n+1)/gamma(n+1)) == -1
      if (n > 1)
        // FIXME: xgamma here is a killer for our precision since it grows
        //        way too fast.
        ans = ans / (std::pow (-1.0, n + 1) / gamma (double (n+1)));
      else if (n == 0)
        ans = -ans;
    }
  else if (ierr == 2)
    ans = - numeric_limits<T>::Inf ();
  else // we probably never get here
    ans = numeric_limits<T>::NaN ();
 
  return ans;
}
 
double
psi (octave_idx_type n, double z) { return xpsi (n, z); }
float
psi (octave_idx_type n, float z) { return xpsi (n, z); }
 
Complex
rc_lgamma (double x)
{
  double result;
 
#if defined (HAVE_LGAMMA_R)
  int sgngam;
  result = lgamma_r (x, &sgngam);
#else
  result = std::lgamma (x);
  int sgngam = signgam;
#endif
 
  if (sgngam < 0 && std::isfinite (result))
    return Complex (result, M_PI);
  else
    return result;
}
 
FloatComplex
rc_lgamma (float x)
{
  float result;
 
#if defined (HAVE_LGAMMAF_R)
  int sgngam;
  result = lgammaf_r (x, &sgngam);
#else
  result = std::lgammaf (x);
  int sgngam = signgam;
#endif
 
  if (sgngam < 0 && std::isfinite (result))
    return FloatComplex (result, M_PI);
  else
    return result;
}
 
Complex
rc_log1p (double x)
{
  return (x < -1.0
          ? Complex (std::log (-(1.0 + x)), M_PI)
          : Complex (log1p (x)));
}
 
FloatComplex
rc_log1p (float x)
{
  return (x < -1.0f
          ? FloatComplex (std::log (-(1.0f + x)), M_PI)
          : FloatComplex (log1p (x)));
}
 
OCTAVE_END_NAMESPACE(math)
OCTAVE_END_NAMESPACE(octave)