/*
Copyright (C) 1996 John W. Eaton
This file is part of Octave.
Octave is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.
Octave is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING. If not, write to the Free
Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
02110-1301, USA.
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include "Range.h"
#include "CColVector.h"
#include "CMatrix.h"
#include "dRowVector.h"
#include "dMatrix.h"
#include "dNDArray.h"
#include "CNDArray.h"
#include "f77-fcn.h"
#include "lo-error.h"
#include "lo-ieee.h"
#include "lo-specfun.h"
#include "mx-inlines.cc"
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
extern "C"
{
F77_RET_T
F77_FUNC (zbesj, ZBESJ) (const double&, const double&, const double&,
const int&, const int&, double*, double*,
int&, int&);
F77_RET_T
F77_FUNC (zbesy, ZBESY) (const double&, const double&, const double&,
const int&, const int&, double*, double*,
int&, double*, double*, int&);
F77_RET_T
F77_FUNC (zbesi, ZBESI) (const double&, const double&, const double&,
const int&, const int&, double*, double*,
int&, int&);
F77_RET_T
F77_FUNC (zbesk, ZBESK) (const double&, const double&, const double&,
const int&, const int&, double*, double*,
int&, int&);
F77_RET_T
F77_FUNC (zbesh, ZBESH) (const double&, const double&, const double&,
const int&, const int&, const int&, double*,
double*, int&, int&);
F77_RET_T
F77_FUNC (zairy, ZAIRY) (const double&, const double&, const int&,
const int&, double&, double&, int&, int&);
F77_RET_T
F77_FUNC (zbiry, ZBIRY) (const double&, const double&, const int&,
const int&, double&, double&, int&);
F77_RET_T
F77_FUNC (xdacosh, XDACOSH) (const double&, double&);
F77_RET_T
F77_FUNC (xdasinh, XDASINH) (const double&, double&);
F77_RET_T
F77_FUNC (xdatanh, XDATANH) (const double&, double&);
F77_RET_T
F77_FUNC (xderf, XDERF) (const double&, double&);
F77_RET_T
F77_FUNC (xderfc, XDERFC) (const double&, double&);
F77_RET_T
F77_FUNC (xdbetai, XDBETAI) (const double&, const double&,
const double&, double&);
F77_RET_T
F77_FUNC (xdgamma, XDGAMMA) (const double&, double&);
F77_RET_T
F77_FUNC (xgammainc, XGAMMAINC) (const double&, const double&, double&);
F77_RET_T
F77_FUNC (dlgams, DLGAMS) (const double&, double&, double&);
}
#if !defined (HAVE_ACOSH)
double
acosh (double x)
{
double retval;
F77_FUNC (xdacosh, XDACOSH) (x, retval);
return retval;
}
#endif
#if !defined (HAVE_ASINH)
double
asinh (double x)
{
double retval;
F77_FUNC (xdasinh, XDASINH) (x, retval);
return retval;
}
#endif
#if !defined (HAVE_ATANH)
double
atanh (double x)
{
double retval;
F77_FUNC (xdatanh, XDATANH) (x, retval);
return retval;
}
#endif
#if !defined (HAVE_ERF)
double
erf (double x)
{
double retval;
F77_FUNC (xderf, XDERF) (x, retval);
return retval;
}
#endif
#if !defined (HAVE_ERFC)
double
erfc (double x)
{
double retval;
F77_FUNC (xderfc, XDERFC) (x, retval);
return retval;
}
#endif
double
xgamma (double x)
{
double result;
F77_FUNC (xdgamma, XDGAMMA) (x, result);
return result;
}
double
xlgamma (double x)
{
double result;
double sgngam;
if (x < 0)
(*current_liboctave_error_handler)
("xlgamma: argument must be nonnegative");
F77_FUNC (dlgams, DLGAMS) (x, result, sgngam);
return result;
}
static inline Complex
zbesj (const Complex& z, double alpha, int kode, int& ierr);
static inline Complex
zbesy (const Complex& z, double alpha, int kode, int& ierr);
static inline Complex
zbesi (const Complex& z, double alpha, int kode, int& ierr);
static inline Complex
zbesk (const Complex& z, double alpha, int kode, int& ierr);
static inline Complex
zbesh1 (const Complex& z, double alpha, int kode, int& ierr);
static inline Complex
zbesh2 (const Complex& z, double alpha, int kode, int& ierr);
static inline Complex
bessel_return_value (const Complex& val, int ierr)
{
static const Complex inf_val = Complex (octave_Inf, octave_Inf);
static const Complex nan_val = Complex (octave_NaN, octave_NaN);
Complex retval;
switch (ierr)
{
case 0:
case 3:
retval = val;
break;
case 2:
retval = inf_val;
break;
default:
retval = nan_val;
break;
}
return retval;
}
static inline bool
is_integer_value (double x)
{
return x == static_cast<double> (static_cast<long> (x));
}
static inline Complex
zbesj (const Complex& z, double alpha, int kode, int& ierr)
{
Complex retval;
if (alpha >= 0.0)
{
double yr = 0.0;
double yi = 0.0;
int nz;
double zr = z.real ();
double zi = z.imag ();
F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr);
if (kode != 2)
{
double expz = exp (std::abs (zi));
yr *= expz;
yi *= expz;
}
if (zi == 0.0 && zr >= 0.0)
yi = 0.0;
retval = bessel_return_value (Complex (yr, yi), ierr);
}
else if (is_integer_value (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 (octave_NaN, octave_NaN);
}
return retval;
}
static inline Complex
zbesy (const Complex& z, double alpha, int kode, int& ierr)
{
Complex retval;
if (alpha >= 0.0)
{
double yr = 0.0;
double yi = 0.0;
int nz;
double wr, wi;
double zr = z.real ();
double zi = z.imag ();
ierr = 0;
if (zr == 0.0 && zi == 0.0)
{
yr = -octave_Inf;
yi = 0.0;
}
else
{
F77_FUNC (zbesy, ZBESY) (zr, zi, alpha, 2, 1, &yr, &yi, nz,
&wr, &wi, ierr);
if (kode != 2)
{
double expz = exp (std::abs (zi));
yr *= expz;
yi *= expz;
}
if (zi == 0.0 && zr >= 0.0)
yi = 0.0;
}
return bessel_return_value (Complex (yr, yi), ierr);
}
else if (is_integer_value (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 (octave_NaN, octave_NaN);
}
return retval;
}
static inline Complex
zbesi (const Complex& z, double alpha, int kode, int& ierr)
{
Complex retval;
if (alpha >= 0.0)
{
double yr = 0.0;
double yi = 0.0;
int nz;
double zr = z.real ();
double zi = z.imag ();
F77_FUNC (zbesi, ZBESI) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr);
if (kode != 2)
{
double expz = exp (std::abs (zr));
yr *= expz;
yi *= expz;
}
if (zi == 0.0 && zr >= 0.0)
yi = 0.0;
retval = bessel_return_value (Complex (yr, yi), ierr);
}
else
{
alpha = -alpha;
Complex tmp = zbesi (z, alpha, kode, ierr);
if (ierr == 0 || ierr == 3)
{
if (! is_integer_value (alpha - 0.5))
tmp += (2.0 / M_PI) * sin (M_PI * alpha)
* zbesk (z, alpha, kode, ierr);
retval = bessel_return_value (tmp, ierr);
}
else
retval = Complex (octave_NaN, octave_NaN);
}
return retval;
}
static inline Complex
zbesk (const Complex& z, double alpha, int kode, int& ierr)
{
Complex retval;
if (alpha >= 0.0)
{
double yr = 0.0;
double yi = 0.0;
int nz;
double zr = z.real ();
double zi = z.imag ();
ierr = 0;
if (zr == 0.0 && zi == 0.0)
{
yr = octave_Inf;
yi = 0.0;
}
else
{
F77_FUNC (zbesk, ZBESK) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr);
if (kode != 2)
{
Complex expz = exp (-z);
double rexpz = real (expz);
double iexpz = imag (expz);
double tmp = yr*rexpz - yi*iexpz;
yi = yr*iexpz + yi*rexpz;
yr = tmp;
}
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, int& ierr)
{
Complex retval;
if (alpha >= 0.0)
{
double yr = 0.0;
double yi = 0.0;
int nz;
double zr = z.real ();
double zi = z.imag ();
F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, 2, 1, 1, &yr, &yi, nz, ierr);
if (kode != 2)
{
Complex expz = exp (Complex (0.0, 1.0) * z);
double rexpz = real (expz);
double iexpz = imag (expz);
double tmp = yr*rexpz - yi*iexpz;
yi = yr*iexpz + yi*rexpz;
yr = tmp;
}
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, int& ierr)
{
Complex retval;
if (alpha >= 0.0)
{
double yr = 0.0;
double yi = 0.0;
int nz;
double zr = z.real ();
double zi = z.imag ();
F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, 2, 2, 1, &yr, &yi, nz, ierr);
if (kode != 2)
{
Complex expz = exp (-Complex (0.0, 1.0) * z);
double rexpz = real (expz);
double iexpz = imag (expz);
double tmp = yr*rexpz - yi*iexpz;
yi = yr*iexpz + yi*rexpz;
yr = tmp;
}
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 (*fptr) (const Complex&, double, int, int&);
static inline Complex
do_bessel (fptr f, const char *, double alpha, const Complex& x,
bool scaled, int& ierr)
{
Complex retval;
retval = f (x, alpha, (scaled ? 2 : 1), ierr);
return retval;
}
static inline ComplexMatrix
do_bessel (fptr f, const char *, double alpha, const ComplexMatrix& x,
bool scaled, Array2<int>& ierr)
{
int nr = x.rows ();
int nc = x.cols ();
ComplexMatrix retval (nr, nc);
ierr.resize (nr, nc);
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
retval(i,j) = f (x(i,j), alpha, (scaled ? 2 : 1), ierr(i,j));
return retval;
}
static inline ComplexMatrix
do_bessel (fptr f, const char *, const Matrix& alpha, const Complex& x,
bool scaled, Array2<int>& ierr)
{
int nr = alpha.rows ();
int nc = alpha.cols ();
ComplexMatrix retval (nr, nc);
ierr.resize (nr, nc);
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
retval(i,j) = f (x, alpha(i,j), (scaled ? 2 : 1), ierr(i,j));
return retval;
}
static inline ComplexMatrix
do_bessel (fptr f, const char *fn, const Matrix& alpha,
const ComplexMatrix& x, bool scaled, Array2<int>& ierr)
{
ComplexMatrix retval;
int x_nr = x.rows ();
int x_nc = x.cols ();
int alpha_nr = alpha.rows ();
int alpha_nc = alpha.cols ();
if (x_nr == alpha_nr && x_nc == alpha_nc)
{
int nr = x_nr;
int nc = x_nc;
retval.resize (nr, nc);
ierr.resize (nr, nc);
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
retval(i,j) = f (x(i,j), alpha(i,j), (scaled ? 2 : 1), ierr(i,j));
}
else
(*current_liboctave_error_handler)
("%s: the sizes of alpha and x must conform", fn);
return retval;
}
static inline ComplexNDArray
do_bessel (fptr f, const char *, double alpha, const ComplexNDArray& x,
bool scaled, ArrayN<int>& ierr)
{
dim_vector dv = x.dims ();
int nel = dv.numel ();
ComplexNDArray retval (dv);
ierr.resize (dv);
for (int i = 0; i < nel; i++)
retval(i) = f (x(i), alpha, (scaled ? 2 : 1), ierr(i));
return retval;
}
static inline ComplexNDArray
do_bessel (fptr f, const char *, const NDArray& alpha, const Complex& x,
bool scaled, ArrayN<int>& ierr)
{
dim_vector dv = alpha.dims ();
int nel = dv.numel ();
ComplexNDArray retval (dv);
ierr.resize (dv);
for (int i = 0; i < nel; i++)
retval(i) = f (x, alpha(i), (scaled ? 2 : 1), ierr(i));
return retval;
}
static inline ComplexNDArray
do_bessel (fptr f, const char *fn, const NDArray& alpha,
const ComplexNDArray& x, bool scaled, ArrayN<int>& ierr)
{
dim_vector dv = x.dims ();
ComplexNDArray retval;
if (dv == alpha.dims ())
{
int nel = dv.numel ();
retval.resize (dv);
ierr.resize (dv);
for (int i = 0; i < nel; i++)
retval(i) = f (x(i), alpha(i), (scaled ? 2 : 1), ierr(i));
}
else
(*current_liboctave_error_handler)
("%s: the sizes of alpha and x must conform", fn);
return retval;
}
static inline ComplexMatrix
do_bessel (fptr f, const char *, const RowVector& alpha,
const ComplexColumnVector& x, bool scaled, Array2<int>& ierr)
{
int nr = x.length ();
int nc = alpha.length ();
ComplexMatrix retval (nr, nc);
ierr.resize (nr, nc);
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
retval(i,j) = f (x(i), alpha(j), (scaled ? 2 : 1), ierr(i,j));
return retval;
}
#define SS_BESSEL(name, fcn) \
Complex \
name (double alpha, const Complex& x, bool scaled, int& ierr) \
{ \
return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
}
#define SM_BESSEL(name, fcn) \
ComplexMatrix \
name (double alpha, const ComplexMatrix& x, bool scaled, \
Array2<int>& ierr) \
{ \
return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
}
#define MS_BESSEL(name, fcn) \
ComplexMatrix \
name (const Matrix& alpha, const Complex& x, bool scaled, \
Array2<int>& ierr) \
{ \
return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
}
#define MM_BESSEL(name, fcn) \
ComplexMatrix \
name (const Matrix& alpha, const ComplexMatrix& x, bool scaled, \
Array2<int>& ierr) \
{ \
return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
}
#define SN_BESSEL(name, fcn) \
ComplexNDArray \
name (double alpha, const ComplexNDArray& x, bool scaled, \
ArrayN<int>& 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, \
ArrayN<int>& 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, \
ArrayN<int>& ierr) \
{ \
return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
}
#define RC_BESSEL(name, fcn) \
ComplexMatrix \
name (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, \
Array2<int>& ierr) \
{ \
return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
}
#define ALL_BESSEL(name, fcn) \
SS_BESSEL (name, fcn) \
SM_BESSEL (name, fcn) \
MS_BESSEL (name, fcn) \
MM_BESSEL (name, fcn) \
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)
Complex
airy (const Complex& z, bool deriv, bool scaled, int& ierr)
{
double ar = 0.0;
double ai = 0.0;
int nz;
double zr = z.real ();
double zi = z.imag ();
int id = deriv ? 1 : 0;
F77_FUNC (zairy, ZAIRY) (zr, zi, id, 2, ar, ai, nz, ierr);
if (! scaled)
{
Complex expz = exp (- 2.0 / 3.0 * z * sqrt(z));
double rexpz = real (expz);
double iexpz = imag (expz);
double tmp = ar*rexpz - ai*iexpz;
ai = ar*iexpz + ai*rexpz;
ar = tmp;
}
if (zi == 0.0 && (! scaled || zr >= 0.0))
ai = 0.0;
return bessel_return_value (Complex (ar, ai), ierr);
}
Complex
biry (const Complex& z, bool deriv, bool scaled, int& ierr)
{
double ar = 0.0;
double ai = 0.0;
double zr = z.real ();
double zi = z.imag ();
int id = deriv ? 1 : 0;
F77_FUNC (zbiry, ZBIRY) (zr, zi, id, 2, ar, ai, ierr);
if (! scaled)
{
Complex expz = exp (std::abs (real (2.0 / 3.0 * z * sqrt (z))));
double rexpz = real (expz);
double iexpz = imag (expz);
double tmp = ar*rexpz - ai*iexpz;
ai = ar*iexpz + ai*rexpz;
ar = tmp;
}
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, Array2<int>& ierr)
{
int nr = z.rows ();
int nc = z.cols ();
ComplexMatrix retval (nr, nc);
ierr.resize (nr, nc);
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
retval(i,j) = airy (z(i,j), deriv, scaled, ierr(i,j));
return retval;
}
ComplexMatrix
biry (const ComplexMatrix& z, bool deriv, bool scaled, Array2<int>& ierr)
{
int nr = z.rows ();
int nc = z.cols ();
ComplexMatrix retval (nr, nc);
ierr.resize (nr, nc);
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
retval(i,j) = biry (z(i,j), deriv, scaled, ierr(i,j));
return retval;
}
ComplexNDArray
airy (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr)
{
dim_vector dv = z.dims ();
int nel = dv.numel ();
ComplexNDArray retval (dv);
ierr.resize (dv);
for (int i = 0; i < nel; i++)
retval (i) = airy (z(i), deriv, scaled, ierr(i));
return retval;
}
ComplexNDArray
biry (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr)
{
dim_vector dv = z.dims ();
int nel = dv.numel ();
ComplexNDArray retval (dv);
ierr.resize (dv);
for (int i = 0; i < nel; i++)
retval (i) = biry (z(i), deriv, scaled, ierr(i));
return retval;
}
static void
gripe_betainc_nonconformant (int r1, int c1, int r2, int c2, int r3,
int c3)
{
(*current_liboctave_error_handler)
("betainc: nonconformant arguments (x is %dx%d, a is %dx%d, b is %dx%d)",
r1, c1, r2, c2, r3, c3);
}
static dim_vector null_dims (0);
static void
gripe_betainc_nonconformant (const dim_vector& d1, const dim_vector& d2,
const dim_vector& d3)
{
std::string d1_str = d1.str ();
std::string d2_str = d2.str ();
std::string d3_str = d3.str ();
(*current_liboctave_error_handler)
("betainc: nonconformant arguments (x is %s, a is %s, b is %s)",
d1_str.c_str (), d2_str.c_str (), d3_str.c_str ());
}
double
betainc (double x, double a, double b)
{
double retval;
F77_FUNC (xdbetai, XDBETAI) (x, a, b, retval);
return retval;
}
Matrix
betainc (double x, double a, const Matrix& b)
{
int nr = b.rows ();
int nc = b.cols ();
Matrix retval (nr, nc);
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
retval(i,j) = betainc (x, a, b(i,j));
return retval;
}
Matrix
betainc (double x, const Matrix& a, double b)
{
int nr = a.rows ();
int nc = a.cols ();
Matrix retval (nr, nc);
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
retval(i,j) = betainc (x, a(i,j), b);
return retval;
}
Matrix
betainc (double x, const Matrix& a, const Matrix& b)
{
Matrix retval;
int a_nr = a.rows ();
int a_nc = a.cols ();
int b_nr = b.rows ();
int b_nc = b.cols ();
if (a_nr == b_nr && a_nc == b_nc)
{
retval.resize (a_nr, a_nc);
for (int j = 0; j < a_nc; j++)
for (int i = 0; i < a_nr; i++)
retval(i,j) = betainc (x, a(i,j), b(i,j));
}
else
gripe_betainc_nonconformant (1, 1, a_nr, a_nc, b_nr, b_nc);
return retval;
}
NDArray
betainc (double x, double a, const NDArray& b)
{
dim_vector dv = b.dims ();
int nel = dv.numel ();
NDArray retval (dv);
for (int i = 0; i < nel; i++)
retval (i) = betainc (x, a, b(i));
return retval;
}
NDArray
betainc (double x, const NDArray& a, double b)
{
dim_vector dv = a.dims ();
int nel = dv.numel ();
NDArray retval (dv);
for (int i = 0; i < nel; i++)
retval (i) = betainc (x, a(i), b);
return retval;
}
NDArray
betainc (double x, const NDArray& a, const NDArray& b)
{
NDArray retval;
dim_vector dv = a.dims ();
if (dv == b.dims ())
{
int nel = dv.numel ();
retval.resize (dv);
for (int i = 0; i < nel; i++)
retval (i) = betainc (x, a(i), b(i));
}
else
gripe_betainc_nonconformant (dim_vector (0), dv, b.dims ());
return retval;
}
Matrix
betainc (const Matrix& x, double a, double b)
{
int nr = x.rows ();
int nc = x.cols ();
Matrix retval (nr, nc);
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
retval(i,j) = betainc (x(i,j), a, b);
return retval;
}
Matrix
betainc (const Matrix& x, double a, const Matrix& b)
{
Matrix retval;
int nr = x.rows ();
int nc = x.cols ();
int b_nr = b.rows ();
int b_nc = b.cols ();
if (nr == b_nr && nc == b_nc)
{
retval.resize (nr, nc);
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
retval(i,j) = betainc (x(i,j), a, b(i,j));
}
else
gripe_betainc_nonconformant (nr, nc, 1, 1, b_nr, b_nc);
return retval;
}
Matrix
betainc (const Matrix& x, const Matrix& a, double b)
{
Matrix retval;
int nr = x.rows ();
int nc = x.cols ();
int a_nr = a.rows ();
int a_nc = a.cols ();
if (nr == a_nr && nc == a_nc)
{
retval.resize (nr, nc);
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
retval(i,j) = betainc (x(i,j), a(i,j), b);
}
else
gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, 1, 1);
return retval;
}
Matrix
betainc (const Matrix& x, const Matrix& a, const Matrix& b)
{
Matrix retval;
int nr = x.rows ();
int nc = x.cols ();
int a_nr = a.rows ();
int a_nc = a.cols ();
int b_nr = b.rows ();
int b_nc = b.cols ();
if (nr == a_nr && nr == b_nr && nc == a_nc && nc == b_nc)
{
retval.resize (nr, nc);
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
retval(i,j) = betainc (x(i,j), a(i,j), b(i,j));
}
else
gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, b_nr, b_nc);
return retval;
}
NDArray
betainc (const NDArray& x, double a, double b)
{
dim_vector dv = x.dims ();
int nel = dv.numel ();
NDArray retval (dv);
for (int i = 0; i < nel; i++)
retval (i) = betainc (x(i), a, b);
return retval;
}
NDArray
betainc (const NDArray& x, double a, const NDArray& b)
{
NDArray retval;
dim_vector dv = x.dims ();
if (dv == b.dims ())
{
int nel = dv.numel ();
retval.resize (dv);
for (int i = 0; i < nel; i++)
retval (i) = betainc (x(i), a, b(i));
}
else
gripe_betainc_nonconformant (dv, dim_vector (0), b.dims ());
return retval;
}
NDArray
betainc (const NDArray& x, const NDArray& a, double b)
{
NDArray retval;
dim_vector dv = x.dims ();
if (dv == a.dims ())
{
int nel = dv.numel ();
retval.resize (dv);
for (int i = 0; i < nel; i++)
retval (i) = betainc (x(i), a(i), b);
}
else
gripe_betainc_nonconformant (dv, a.dims (), dim_vector (0));
return retval;
}
NDArray
betainc (const NDArray& x, const NDArray& a, const NDArray& b)
{
NDArray retval;
dim_vector dv = x.dims ();
if (dv == a.dims () && dv == b.dims ())
{
int nel = dv.numel ();
retval.resize (dv);
for (int i = 0; i < nel; i++)
retval (i) = betainc (x(i), a(i), b(i));
}
else
gripe_betainc_nonconformant (dv, a.dims (), b.dims ());
return retval;
}
// XXX FIXME XXX -- there is still room for improvement here...
double
gammainc (double x, double a, bool& err)
{
double retval;
err = false;
if (a < 0.0 || x < 0.0)
{
(*current_liboctave_error_handler)
("gammainc: A and X must be non-negative");
err = true;
}
else
F77_FUNC (xgammainc, XGAMMAINC) (a, x, retval);
return retval;
}
Matrix
gammainc (double x, const Matrix& a)
{
int nr = a.rows ();
int nc = a.cols ();
Matrix result (nr, nc);
Matrix retval;
bool err;
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
{
result(i,j) = gammainc (x, a(i,j), err);
if (err)
goto done;
}
retval = result;
done:
return retval;
}
Matrix
gammainc (const Matrix& x, double a)
{
int nr = x.rows ();
int nc = x.cols ();
Matrix result (nr, nc);
Matrix retval;
bool err;
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
{
result(i,j) = gammainc (x(i,j), a, err);
if (err)
goto done;
}
retval = result;
done:
return retval;
}
Matrix
gammainc (const Matrix& x, const Matrix& a)
{
Matrix result;
Matrix retval;
int nr = x.rows ();
int nc = x.cols ();
int a_nr = a.rows ();
int a_nc = a.cols ();
if (nr == a_nr && nc == a_nc)
{
result.resize (nr, nc);
bool err;
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
{
result(i,j) = gammainc (x(i,j), a(i,j), err);
if (err)
goto done;
}
retval = result;
}
else
(*current_liboctave_error_handler)
("gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)",
nr, nc, a_nr, a_nc);
done:
return retval;
}
NDArray
gammainc (double x, const NDArray& a)
{
dim_vector dv = a.dims ();
int nel = dv.numel ();
NDArray retval;
NDArray result (dv);
bool err;
for (int i = 0; i < nel; i++)
{
result (i) = gammainc (x, a(i), err);
if (err)
goto done;
}
retval = result;
done:
return retval;
}
NDArray
gammainc (const NDArray& x, double a)
{
dim_vector dv = x.dims ();
int nel = dv.numel ();
NDArray retval;
NDArray result (dv);
bool err;
for (int i = 0; i < nel; i++)
{
result (i) = gammainc (x(i), a, err);
if (err)
goto done;
}
retval = result;
done:
return retval;
}
NDArray
gammainc (const NDArray& x, const NDArray& a)
{
dim_vector dv = x.dims ();
int nel = dv.numel ();
NDArray retval;
NDArray result;
if (dv == a.dims ())
{
result.resize (dv);
bool err;
for (int i = 0; i < nel; i++)
{
result (i) = gammainc (x(i), a(i), err);
if (err)
goto done;
}
retval = result;
}
else
{
std::string x_str = dv.str ();
std::string a_str = a.dims ().str ();
(*current_liboctave_error_handler)
("gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)",
x_str.c_str (), a_str. c_str ());
}
done:
return retval;
}
/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/
syntax highlighted by Code2HTML, v. 0.9.1