Mercurial > octave
diff liboctave/numeric/lo-specfun.cc @ 31607:aac27ad79be6 stable
maint: Re-indent code after switch to using namespace macros.
* build-env.h, build-env.in.cc, Cell.h, __betainc__.cc, __eigs__.cc,
__ftp__.cc, __ichol__.cc, __ilu__.cc, __isprimelarge__.cc, __magick_read__.cc,
__pchip_deriv__.cc, amd.cc, base-text-renderer.cc, base-text-renderer.h,
besselj.cc, bitfcns.cc, bsxfun.cc, c-file-ptr-stream.h, call-stack.cc,
call-stack.h, ccolamd.cc, cellfun.cc, chol.cc, colamd.cc, dasrt.cc, data.cc,
debug.cc, defaults.cc, defaults.h, det.cc, display.cc, display.h, dlmread.cc,
dynamic-ld.cc, dynamic-ld.h, ellipj.cc, environment.cc, environment.h,
error.cc, error.h, errwarn.h, event-manager.cc, event-manager.h,
event-queue.cc, event-queue.h, fcn-info.cc, fcn-info.h, fft.cc, fft2.cc,
file-io.cc, filter.cc, find.cc, ft-text-renderer.cc, ft-text-renderer.h,
gcd.cc, gl-render.cc, gl-render.h, gl2ps-print.cc, gl2ps-print.h,
graphics-toolkit.cc, graphics-toolkit.h, graphics.cc, gsvd.cc, gtk-manager.cc,
gtk-manager.h, help.cc, help.h, hook-fcn.cc, hook-fcn.h, input.cc, input.h,
interpreter-private.cc, interpreter-private.h, interpreter.cc, interpreter.h,
inv.cc, jsondecode.cc, jsonencode.cc, latex-text-renderer.cc,
latex-text-renderer.h, load-path.cc, load-path.h, load-save.cc, load-save.h,
lookup.cc, ls-hdf5.cc, ls-mat4.cc, ls-mat5.cc, lsode.cc, lu.cc, mappers.cc,
matrix_type.cc, max.cc, mex.cc, mexproto.h, mxarray.h, mxtypes.in.h,
oct-errno.in.cc, oct-hdf5-types.cc, oct-hist.cc, oct-hist.h, oct-map.cc,
oct-map.h, oct-opengl.h, oct-prcstrm.h, oct-process.cc, oct-process.h,
oct-stdstrm.h, oct-stream.cc, oct-stream.h, oct-strstrm.h,
octave-default-image.h, ordqz.cc, ordschur.cc, pager.cc, pager.h, pinv.cc,
pow2.cc, pr-output.cc, psi.cc, qr.cc, quadcc.cc, rand.cc, regexp.cc,
settings.cc, settings.h, sighandlers.cc, sighandlers.h, sparse-xpow.cc,
sqrtm.cc, stack-frame.cc, stack-frame.h, stream-euler.cc, strfns.cc, svd.cc,
syminfo.cc, syminfo.h, symrcm.cc, symrec.cc, symrec.h, symscope.cc, symscope.h,
symtab.cc, symtab.h, sysdep.cc, sysdep.h, text-engine.cc, text-engine.h,
text-renderer.cc, text-renderer.h, time.cc, toplev.cc, typecast.cc,
url-handle-manager.cc, url-handle-manager.h, urlwrite.cc, utils.cc, utils.h,
variables.cc, variables.h, xdiv.cc, __delaunayn__.cc, __init_fltk__.cc,
__init_gnuplot__.cc, __ode15__.cc, __voronoi__.cc, audioread.cc, convhulln.cc,
gzip.cc, cdef-class.cc, cdef-class.h, cdef-fwd.h, cdef-manager.cc,
cdef-manager.h, cdef-method.cc, cdef-method.h, cdef-object.cc, cdef-object.h,
cdef-package.cc, cdef-package.h, cdef-property.cc, cdef-property.h,
cdef-utils.cc, cdef-utils.h, ov-base-diag.cc, ov-base-int.cc, ov-base-mat.cc,
ov-base-mat.h, ov-base-scalar.cc, ov-base.cc, ov-base.h, ov-bool-mat.cc,
ov-bool-mat.h, ov-bool-sparse.cc, ov-bool.cc, ov-builtin.h, ov-cell.cc,
ov-ch-mat.cc, ov-class.cc, ov-class.h, ov-classdef.cc, ov-classdef.h,
ov-complex.cc, ov-cx-diag.cc, ov-cx-mat.cc, ov-cx-sparse.cc, ov-dld-fcn.cc,
ov-dld-fcn.h, ov-fcn-handle.cc, ov-fcn-handle.h, ov-fcn.h, ov-float.cc,
ov-flt-complex.cc, ov-flt-cx-diag.cc, ov-flt-cx-mat.cc, ov-flt-re-diag.cc,
ov-flt-re-mat.cc, ov-flt-re-mat.h, ov-intx.h, ov-java.cc, ov-lazy-idx.cc,
ov-legacy-range.cc, ov-magic-int.cc, ov-mex-fcn.cc, ov-mex-fcn.h,
ov-null-mat.cc, ov-perm.cc, ov-range.cc, ov-re-diag.cc, ov-re-mat.cc,
ov-re-mat.h, ov-re-sparse.cc, ov-scalar.cc, ov-str-mat.cc, ov-struct.cc,
ov-typeinfo.cc, ov-typeinfo.h, ov-usr-fcn.cc, ov-usr-fcn.h, ov.cc, ov.h, ovl.h,
octave.cc, octave.h, op-b-sbm.cc, op-bm-sbm.cc, op-cs-scm.cc, op-fm-fcm.cc,
op-fs-fcm.cc, op-s-scm.cc, op-scm-cs.cc, op-scm-s.cc, op-sm-cs.cc, ops.h,
anon-fcn-validator.cc, anon-fcn-validator.h, bp-table.cc, bp-table.h,
comment-list.cc, comment-list.h, filepos.h, lex.h, oct-lvalue.cc, oct-lvalue.h,
parse.h, profiler.cc, profiler.h, pt-anon-scopes.cc, pt-anon-scopes.h,
pt-arg-list.cc, pt-arg-list.h, pt-args-block.cc, pt-args-block.h,
pt-array-list.cc, pt-array-list.h, pt-assign.cc, pt-assign.h, pt-binop.cc,
pt-binop.h, pt-bp.cc, pt-bp.h, pt-cbinop.cc, pt-cbinop.h, pt-cell.cc,
pt-cell.h, pt-check.cc, pt-check.h, pt-classdef.cc, pt-classdef.h, pt-cmd.h,
pt-colon.cc, pt-colon.h, pt-const.cc, pt-const.h, pt-decl.cc, pt-decl.h,
pt-eval.cc, pt-eval.h, pt-except.cc, pt-except.h, pt-exp.cc, pt-exp.h,
pt-fcn-handle.cc, pt-fcn-handle.h, pt-id.cc, pt-id.h, pt-idx.cc, pt-idx.h,
pt-jump.h, pt-loop.cc, pt-loop.h, pt-mat.cc, pt-mat.h, pt-misc.cc, pt-misc.h,
pt-pr-code.cc, pt-pr-code.h, pt-select.cc, pt-select.h, pt-spmd.cc, pt-spmd.h,
pt-stmt.cc, pt-stmt.h, pt-tm-const.cc, pt-tm-const.h, pt-unop.cc, pt-unop.h,
pt-walk.cc, pt-walk.h, pt.cc, pt.h, token.cc, token.h, Range.cc, Range.h,
idx-vector.cc, idx-vector.h, range-fwd.h, CollocWt.cc, CollocWt.h,
aepbalance.cc, aepbalance.h, chol.cc, chol.h, gepbalance.cc, gepbalance.h,
gsvd.cc, gsvd.h, hess.cc, hess.h, lo-mappers.cc, lo-mappers.h, lo-specfun.cc,
lo-specfun.h, lu.cc, lu.h, oct-convn.cc, oct-convn.h, oct-fftw.cc, oct-fftw.h,
oct-norm.cc, oct-norm.h, oct-rand.cc, oct-rand.h, oct-spparms.cc,
oct-spparms.h, qr.cc, qr.h, qrp.cc, qrp.h, randgamma.cc, randgamma.h,
randmtzig.cc, randmtzig.h, randpoisson.cc, randpoisson.h, schur.cc, schur.h,
sparse-chol.cc, sparse-chol.h, sparse-lu.cc, sparse-lu.h, sparse-qr.cc,
sparse-qr.h, svd.cc, svd.h, child-list.cc, child-list.h, dir-ops.cc, dir-ops.h,
file-ops.cc, file-ops.h, file-stat.cc, file-stat.h, lo-sysdep.cc, lo-sysdep.h,
lo-sysinfo.cc, lo-sysinfo.h, mach-info.cc, mach-info.h, oct-env.cc, oct-env.h,
oct-group.cc, oct-group.h, oct-password.cc, oct-password.h, oct-syscalls.cc,
oct-syscalls.h, oct-time.cc, oct-time.h, oct-uname.cc, oct-uname.h,
action-container.cc, action-container.h, base-list.h, cmd-edit.cc, cmd-edit.h,
cmd-hist.cc, cmd-hist.h, f77-fcn.h, file-info.cc, file-info.h,
lo-array-errwarn.cc, lo-array-errwarn.h, lo-hash.cc, lo-hash.h, lo-ieee.h,
lo-regexp.cc, lo-regexp.h, lo-utils.cc, lo-utils.h, oct-base64.cc,
oct-base64.h, oct-glob.cc, oct-glob.h, oct-inttypes.h, oct-mutex.cc,
oct-mutex.h, oct-refcount.h, oct-shlib.cc, oct-shlib.h, oct-sparse.cc,
oct-sparse.h, oct-string.h, octave-preserve-stream-state.h, pathsearch.cc,
pathsearch.h, quit.cc, quit.h, unwind-prot.cc, unwind-prot.h, url-transfer.cc,
url-transfer.h:
Re-indent code after switch to using namespace macros.
| author | Rik <rik@octave.org> |
|---|---|
| date | Thu, 01 Dec 2022 18:02:15 -0800 |
| parents | e88a07dec498 |
| children | 597f3ee61a48 |
line wrap: on
line diff
--- a/liboctave/numeric/lo-specfun.cc Thu Dec 01 14:23:45 2022 -0800 +++ b/liboctave/numeric/lo-specfun.cc Thu Dec 01 18:02:15 2022 -0800 @@ -59,658 +59,658 @@ OCTAVE_BEGIN_NAMESPACE(math) - static inline Complex - bessel_return_value (const Complex& val, octave_idx_type ierr) +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) { - static const Complex inf_val - = Complex (numeric_limits<double>::Inf (), - numeric_limits<double>::Inf ()); + case 0: + case 3: + case 4: + retval = val; + break; - static const Complex nan_val - = Complex (numeric_limits<double>::NaN (), - numeric_limits<double>::NaN ()); + case 2: + retval = inf_val; + break; - Complex retval; + default: + retval = nan_val; + break; + } - switch (ierr) - { - case 0: - case 3: - case 4: - retval = 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; - case 2: - retval = inf_val; - break; + switch (ierr) + { + case 0: + case 3: + case 4: + retval = val; + break; - default: - retval = nan_val; - break; - } + case 2: + retval = inf_val; + break; - return retval; + default: + retval = nan_val; + break; } - 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 ()); + 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)); - static const FloatComplex nan_val - = FloatComplex (numeric_limits<float>::NaN (), - numeric_limits<float>::NaN ()); + 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) +{ + dim_vector dv = z.dims (); + octave_idx_type nel = dv.numel (); + ComplexNDArray retval (dv); + + ierr.resize (dv); - FloatComplex retval; + 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; - switch (ierr) - { - case 0: - case 3: - case 4: - retval = val; - break; + 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); - case 2: - retval = inf_val; - break; + 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) +{ + dim_vector dv = z.dims (); + octave_idx_type nel = dv.numel (); + FloatComplexNDArray retval (dv); + + ierr.resize (dv); - default: - retval = nan_val; - break; - } + for (octave_idx_type i = 0; i < nel; i++) + retval(i) = airy (z(i), deriv, scaled, ierr(i)); + + 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, 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); - return retval; - } +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); - Complex - airy (const Complex& z, bool deriv, bool scaled, 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 ar = 0.0; - double ai = 0.0; + double yr = 0.0; + double yi = 0.0; + + F77_INT nz, t_ierr; 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) - { - 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); + F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, kode, 1, &yr, &yi, 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; + if (zi == 0.0 && zr >= 0.0) + yi = 0.0; - return bessel_return_value (FloatComplex (ar, ai), ierr); + retval = bessel_return_value (Complex (yr, yi), ierr); } - - FloatComplexMatrix - airy (const FloatComplexMatrix& z, bool deriv, bool scaled, - Array<octave_idx_type>& ierr) + else if (is_integer_value (alpha)) { - octave_idx_type nr = z.rows (); - octave_idx_type nc = z.cols (); - - FloatComplexMatrix retval (nr, nc); + // 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; - ierr.resize (dim_vector (nr, nc)); + 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); - 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; + retval = bessel_return_value (tmp, ierr); + } + else + retval = Complex (numeric_limits<double>::NaN (), + numeric_limits<double>::NaN ()); } - FloatComplexNDArray - airy (const FloatComplexNDArray& z, bool deriv, bool scaled, - Array<octave_idx_type>& ierr) - { - 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; +} - return retval; - } +static inline Complex +zbesy (const Complex& z, double alpha, int kode, octave_idx_type& ierr) +{ + Complex retval; - static inline bool - is_integer_value (double x) + if (alpha >= 0.0) { - return x == static_cast<double> (static_cast<long> (x)); - } + double yr = 0.0; + double yi = 0.0; - 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); + F77_INT nz, t_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); + double wr, wi; - static inline Complex - zbesh1 (const Complex& z, double alpha, int kode, octave_idx_type& ierr); + double zr = z.real (); + double zi = z.imag (); - static inline Complex - zbesh2 (const Complex& z, double alpha, int kode, octave_idx_type& ierr); + ierr = 0; - static inline Complex - zbesj (const Complex& z, double alpha, int kode, octave_idx_type& ierr) - { - Complex retval; - - if (alpha >= 0.0) + if (zr == 0.0 && zi == 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); + 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; + } - retval = bessel_return_value (Complex (yr, yi), ierr); - } - else if (is_integer_value (alpha)) + 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) { - // 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; + 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 (is_integer_value (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 { - 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 (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 (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); + 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 if (is_integer_value (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; + retval = bessel_return_value (Complex (yr, yi), ierr); } - - typedef Complex (*dptr) (const Complex&, double, int, octave_idx_type&); + else + { + Complex tmp = zbesk (z, -alpha, kode, ierr); - 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; + retval = bessel_return_value (tmp, ierr); } - static inline ComplexMatrix - do_bessel (dptr f, const char *, double alpha, const ComplexMatrix& x, - bool scaled, Array<octave_idx_type>& ierr) + return retval; +} + +static inline Complex +zbesh1 (const Complex& z, double alpha, int kode, octave_idx_type& ierr) +{ + Complex retval; + + if (alpha >= 0.0) { - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); + double yr = 0.0; + double yi = 0.0; + + F77_INT nz, t_ierr; - ComplexMatrix retval (nr, nc); + double zr = z.real (); + double zi = z.imag (); - ierr.resize (dim_vector (nr, nc)); + F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, kode, 1, 1, &yr, &yi, nz, + t_ierr); + + ierr = t_ierr; - 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)); + retval = bessel_return_value (Complex (yr, yi), ierr); + } + else + { + alpha = -alpha; - return retval; + 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); } - static inline ComplexMatrix - do_bessel (dptr f, const char *, const Matrix& alpha, const Complex& x, - bool scaled, Array<octave_idx_type>& ierr) + return retval; +} + +static inline Complex +zbesh2 (const Complex& z, double alpha, int kode, octave_idx_type& ierr) +{ + Complex retval; + + if (alpha >= 0.0) { - octave_idx_type nr = alpha.rows (); - octave_idx_type nc = alpha.cols (); + double yr = 0.0; + double yi = 0.0; + + F77_INT nz, t_ierr; - ComplexMatrix retval (nr, nc); + double zr = z.real (); + double zi = z.imag (); - ierr.resize (dim_vector (nr, nc)); + F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, kode, 2, 1, &yr, &yi, nz, + t_ierr); + + ierr = t_ierr; - 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)); + retval = bessel_return_value (Complex (yr, yi), ierr); + } + else + { + alpha = -alpha; - return retval; + 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); } - 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; + return retval; +} + +typedef Complex (*dptr) (const Complex&, double, int, octave_idx_type&); - octave_idx_type x_nr = x.rows (); - octave_idx_type x_nc = x.cols (); +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; +} - octave_idx_type alpha_nr = alpha.rows (); - octave_idx_type alpha_nc = alpha.cols (); +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 (); - if (x_nr != alpha_nr || x_nc != alpha_nc) - (*current_liboctave_error_handler) - ("%s: the sizes of alpha and x must conform", fn); + ComplexMatrix retval (nr, nc); - octave_idx_type nr = x_nr; - octave_idx_type nc = x_nc; + ierr.resize (dim_vector (nr, nc)); - retval.resize (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)); - ierr.resize (dim_vector (nr, nc)); + return retval; +} - 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)); +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); - return retval; - } + 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 ComplexNDArray - do_bessel (dptr f, const char *, double alpha, const ComplexNDArray& x, - bool scaled, Array<octave_idx_type>& ierr) - { - dim_vector dv = x.dims (); - octave_idx_type nel = dv.numel (); - ComplexNDArray retval (dv); +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 (); - ierr.resize (dv); + 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); - for (octave_idx_type i = 0; i < nel; i++) - retval(i) = f (x(i), alpha, (scaled ? 2 : 1), ierr(i)); + octave_idx_type nr = x_nr; + octave_idx_type nc = x_nc; - return retval; - } + retval.resize (nr, nc); + + ierr.resize (dim_vector (nr, nc)); - static inline ComplexNDArray - do_bessel (dptr f, const char *, const NDArray& alpha, const Complex& x, - bool scaled, Array<octave_idx_type>& ierr) - { - dim_vector dv = alpha.dims (); - octave_idx_type nel = dv.numel (); - ComplexNDArray retval (dv); + 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; +} - ierr.resize (dv); +static inline ComplexNDArray +do_bessel (dptr f, const char *, double alpha, const ComplexNDArray& x, + bool scaled, Array<octave_idx_type>& ierr) +{ + dim_vector dv = x.dims (); + octave_idx_type nel = dv.numel (); + ComplexNDArray retval (dv); - for (octave_idx_type i = 0; i < nel; i++) - retval(i) = f (x, alpha(i), (scaled ? 2 : 1), ierr(i)); + ierr.resize (dv); - return retval; - } + 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 *fn, const NDArray& alpha, - const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) - { - dim_vector dv = x.dims (); - ComplexNDArray retval; +static inline ComplexNDArray +do_bessel (dptr f, const char *, const NDArray& alpha, const Complex& x, + bool scaled, Array<octave_idx_type>& ierr) +{ + dim_vector dv = alpha.dims (); + octave_idx_type nel = dv.numel (); + ComplexNDArray retval (dv); - if (dv != alpha.dims ()) - (*current_liboctave_error_handler) - ("%s: the sizes of alpha and x must conform", fn); + 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; +} - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - ierr.resize (dv); +static inline ComplexNDArray +do_bessel (dptr f, const char *fn, const NDArray& alpha, + const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) +{ + dim_vector dv = x.dims (); + ComplexNDArray retval; - for (octave_idx_type i = 0; i < nel; i++) - retval(i) = f (x(i), alpha(i), (scaled ? 2 : 1), ierr(i)); + if (dv != alpha.dims ()) + (*current_liboctave_error_handler) + ("%s: the sizes of alpha and x must conform", fn); + + octave_idx_type nel = dv.numel (); - return retval; - } + 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)); - 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 (); + 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); + ComplexMatrix retval (nr, nc); - ierr.resize (dim_vector (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)); + 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; - } + return retval; +} #define SS_BESSEL(name, fcn) \ Complex \ @@ -785,12 +785,12 @@ 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) +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 @@ -802,458 +802,458 @@ #undef NN_BESSEL #undef RC_BESSEL - static inline FloatComplex - cbesj (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); + +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 - cbesy (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 bool +is_integer_value (float x) +{ + return x == static_cast<float> (static_cast<long> (x)); +} - static inline FloatComplex - cbesi (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; - static inline FloatComplex - cbesk (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr); + F77_FUNC (cbesj, CBESJ) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, + F77_CMPLX_ARG (&y), nz, t_ierr); - static inline FloatComplex - cbesh1 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr); + ierr = t_ierr; + + if (z.imag () == 0.0 && z.real () >= 0.0) + y = FloatComplex (y.real (), 0.0); - static inline FloatComplex - cbesh2 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr); + retval = bessel_return_value (y, ierr); + } + else if (is_integer_value (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; - static inline bool - is_integer_value (float x) - { - return x == static_cast<float> (static_cast<long> (x)); + 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 ()); } - static inline FloatComplex - cbesj (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr) + return retval; +} + +static inline FloatComplex +cbesy (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr) +{ + FloatComplex retval; + + if (alpha >= 0.0) { - FloatComplex retval; + FloatComplex y = 0.0; + + F77_INT nz, t_ierr; - if (alpha >= 0.0) + FloatComplex w; + + ierr = 0; + + if (z.real () == 0.0 && z.imag () == 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); + 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); + } - retval = bessel_return_value (y, ierr); - } - else if (is_integer_value (alpha)) + return bessel_return_value (y, ierr); + } + else if (is_integer_value (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) { - // 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; + 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 { - 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 (is_integer_value (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_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 - { - 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; + retval = bessel_return_value (y, ierr); + } + else + { + FloatComplex tmp = cbesk (z, -alpha, kode, ierr); + + retval = bessel_return_value (tmp, ierr); } - static inline FloatComplex - cbesk (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr) - { - FloatComplex retval; + return retval; +} - if (alpha >= 0.0) - { - FloatComplex y = 0.0; +static inline FloatComplex +cbesh1 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr) +{ + FloatComplex retval; - F77_INT nz, t_ierr; + if (alpha >= 0.0) + { + FloatComplex y = 0.0; - ierr = 0; + F77_INT nz, t_ierr; - 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); + F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, 1, + F77_CMPLX_ARG (&y), nz, t_ierr); - ierr = 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; - retval = bessel_return_value (y, ierr); - } - else - { - FloatComplex tmp = cbesk (z, -alpha, kode, ierr); + static const FloatComplex eye = FloatComplex (0.0, 1.0); - retval = bessel_return_value (tmp, ierr); - } + FloatComplex tmp = exp (static_cast<float> (M_PI) * alpha * eye) + * cbesh1 (z, alpha, kode, ierr); - return retval; + retval = bessel_return_value (tmp, ierr); } - static inline FloatComplex - cbesh1 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr) - { - FloatComplex retval; + return 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; +static inline FloatComplex +cbesh2 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr) +{ + FloatComplex retval; - retval = bessel_return_value (y, ierr); - } - else - { - alpha = -alpha; + if (alpha >= 0.0) + { + FloatComplex y = 0.0;; - static const FloatComplex eye = FloatComplex (0.0, 1.0); + F77_INT nz, t_ierr; - FloatComplex tmp = exp (static_cast<float> (M_PI) * alpha * eye) - * cbesh1 (z, alpha, kode, ierr); + F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 2, 1, + F77_CMPLX_ARG (&y), nz, t_ierr); - retval = bessel_return_value (tmp, ierr); - } + ierr = t_ierr; - return retval; + retval = bessel_return_value (y, ierr); } - - static inline FloatComplex - cbesh2 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr) + else { - 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); + alpha = -alpha; - ierr = t_ierr; - - retval = bessel_return_value (y, ierr); - } - else - { - alpha = -alpha; + static const FloatComplex eye = FloatComplex (0.0, 1.0); - static const FloatComplex eye = FloatComplex (0.0, 1.0); - - FloatComplex tmp = exp (-static_cast<float> (M_PI) * alpha * eye) - * cbesh2 (z, alpha, kode, ierr); + FloatComplex tmp = exp (-static_cast<float> (M_PI) * alpha * eye) + * cbesh2 (z, alpha, kode, ierr); - retval = bessel_return_value (tmp, ierr); - } - - return retval; + retval = bessel_return_value (tmp, ierr); } - typedef FloatComplex (*fptr) (const FloatComplex&, float, int, - octave_idx_type&); + 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; +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); + retval = f (x, alpha, (scaled ? 2 : 1), ierr); - return retval; - } + 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 (); +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); + FloatComplexMatrix retval (nr, nc); - ierr.resize (dim_vector (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)); + 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; - } + 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 (); +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); + FloatComplexMatrix retval (nr, nc); - ierr.resize (dim_vector (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)); + 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; - } + 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; +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 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 (); + 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); + 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; + octave_idx_type nr = x_nr; + octave_idx_type nc = x_nc; - retval.resize (nr, nc); + retval.resize (nr, nc); - ierr.resize (dim_vector (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)); + 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; - } + return retval; +} - static inline FloatComplexNDArray - do_bessel (fptr f, const char *, float alpha, const FloatComplexNDArray& x, - bool scaled, Array<octave_idx_type>& ierr) - { - dim_vector dv = x.dims (); - octave_idx_type nel = dv.numel (); - FloatComplexNDArray retval (dv); +static inline FloatComplexNDArray +do_bessel (fptr f, const char *, float alpha, const FloatComplexNDArray& x, + bool scaled, Array<octave_idx_type>& ierr) +{ + dim_vector dv = x.dims (); + octave_idx_type nel = dv.numel (); + FloatComplexNDArray retval (dv); - ierr.resize (dv); + ierr.resize (dv); - for (octave_idx_type i = 0; i < nel; i++) - retval(i) = f (x(i), alpha, (scaled ? 2 : 1), ierr(i)); + for (octave_idx_type i = 0; i < nel; i++) + retval(i) = f (x(i), alpha, (scaled ? 2 : 1), ierr(i)); - return retval; - } + return retval; +} - static inline FloatComplexNDArray - do_bessel (fptr f, const char *, const FloatNDArray& alpha, - const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr) - { - dim_vector dv = alpha.dims (); - octave_idx_type nel = dv.numel (); - FloatComplexNDArray retval (dv); +static inline FloatComplexNDArray +do_bessel (fptr f, const char *, const FloatNDArray& alpha, + const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr) +{ + dim_vector dv = alpha.dims (); + octave_idx_type nel = dv.numel (); + FloatComplexNDArray retval (dv); - ierr.resize (dv); + ierr.resize (dv); - for (octave_idx_type i = 0; i < nel; i++) - retval(i) = f (x, alpha(i), (scaled ? 2 : 1), ierr(i)); + for (octave_idx_type i = 0; i < nel; i++) + retval(i) = f (x, alpha(i), (scaled ? 2 : 1), ierr(i)); - return retval; - } + 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) - { - dim_vector dv = x.dims (); - FloatComplexNDArray retval; +static inline FloatComplexNDArray +do_bessel (fptr f, const char *fn, const FloatNDArray& alpha, + const FloatComplexNDArray& x, bool scaled, + Array<octave_idx_type>& ierr) +{ + 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); + if (dv != alpha.dims ()) + (*current_liboctave_error_handler) + ("%s: the sizes of alpha and x must conform", fn); - octave_idx_type nel = dv.numel (); + octave_idx_type nel = dv.numel (); - retval.resize (dv); - ierr.resize (dv); + 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)); + for (octave_idx_type i = 0; i < nel; i++) + retval(i) = f (x(i), alpha(i), (scaled ? 2 : 1), ierr(i)); - return retval; - } + 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 (); +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); + FloatComplexMatrix retval (nr, nc); - ierr.resize (dim_vector (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)); + 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; - } + return retval; +} #define SS_BESSEL(name, fcn) \ FloatComplex \ @@ -1330,12 +1330,12 @@ 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) +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 @@ -1347,865 +1347,865 @@ #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; +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)); - double zr = z.real (); - double zi = z.imag (); + 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) +{ + dim_vector dv = z.dims (); + octave_idx_type nel = dv.numel (); + ComplexNDArray retval (dv); + + ierr.resize (dv); - F77_INT id = (deriv ? 1 : 0); - F77_INT t_ierr; - F77_INT sc = (scaled ? 2 : 1); + 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); - F77_FUNC (zbiry, ZBIRY) (zr, zi, id, sc, ar, ai, 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); +} - ierr = t_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); - if (zi == 0.0 && (! scaled || zr >= 0.0)) - ai = 0.0; + 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; +} - return bessel_return_value (Complex (ar, ai), ierr); - } +FloatComplexNDArray +biry (const FloatComplexNDArray& z, bool deriv, bool scaled, + Array<octave_idx_type>& ierr) +{ + 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); } - ComplexMatrix - biry (const ComplexMatrix& z, bool deriv, bool scaled, - Array<octave_idx_type>& ierr) +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) { - 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; - } + (*current_liboctave_warning_with_id_handler) + ("Octave:ellipj-invalid-m", + "ellipj: invalid M value, required value 0 <= M <= 1"); - ComplexNDArray - biry (const ComplexNDArray& z, bool deriv, bool scaled, - Array<octave_idx_type>& ierr) - { - dim_vector dv = z.dims (); - octave_idx_type nel = dv.numel (); - ComplexNDArray retval (dv); + sn = cn = dn = lo_ieee_nan_value (); - ierr.resize (dv); - - for (octave_idx_type i = 0; i < nel; i++) - retval(i) = biry (z(i), deriv, scaled, ierr(i)); - - return retval; + return; } - FloatComplex - biry (const FloatComplex& z, bool deriv, bool scaled, - octave_idx_type& ierr) + double sqrt_eps = std::sqrt (std::numeric_limits<double>::epsilon ()); + if (m < sqrt_eps) { - 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); + // 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; } - - 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) + else if ((1 - m) < sqrt_eps) { - 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); + // 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; } - - 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) + else { - 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) + // 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) { - (*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 (); - + 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; } - - 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 + 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) { - // 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); + 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; - void - ellipj (const Complex& u, double m, Complex& sn, Complex& cn, Complex& dn, - double& err) + ellipj (u.imag (), m1, ss1, cc1, dd1, err); + if (u.real () == 0) { - double m1 = 1 - m, ss1, cc1, dd1; + // 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.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); +} - 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); - } - } +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 error function from the Faddeeva package - Complex - erf (const Complex& x) - { - return Faddeeva::erf (x); - } +// Complex complementary error function from the Faddeeva package +Complex +erfc (const Complex& x) +{ + return Faddeeva::erfc (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 ()); - } +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 ()); +} - // 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. +// 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 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; - 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)))); - // 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; + const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]; - y = yn / yd; + const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0; - 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 (); + y = yn / yd; - 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); - } + 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 (); - return y; - } - - double erfcinv (double x) + if (refine) { - return do_erfcinv (x, true); + // One iteration of Halley's method gives full precision. + double u = (erf (y) - (1-x)) * spi2 * exp (y*y); + y -= u / (1 + y*u); } - float erfcinv (float x) - { - return do_erfcinv (x, false); - } + return y; +} - // 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); - } +double erfcinv (double x) +{ + return do_erfcinv (x, true); +} - 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 ()); - } +float erfcinv (float x) +{ + return do_erfcinv (x, false); +} - // 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); } +// 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 - erfi (const Complex& x) - { - return Faddeeva::erfi (x); - } +Complex +erfcx (const Complex& x) +{ + return Faddeeva::erfcx (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 ()); - } +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 ()); +} - // 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. +// 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); } - 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 - }; +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 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); - } +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 + }; - 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); + static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2. + static const double pbreak = 0.95150; + double ax = fabs (x), y; - return retval; + // 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; } - - FloatComplex - expm1 (const FloatComplex& x) + else if (ax < 1.0) { - FloatComplex retval; + // 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 (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; + 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); } - double - gamma (double x) - { - double result; + return y; +} - // Special cases for (near) compatibility with Matlab instead of tgamma. - // Matlab does not have -0. +double erfinv (double x) +{ + return do_erfinv (x, true); +} - if (x == 0) - result = (math::negative_sign (x) - ? -numeric_limits<double>::Inf () - : numeric_limits<double>::Inf ()); - else if ((x < 0 && math::x_nint (x) == x) - || math::isinf (x)) - result = numeric_limits<double>::Inf (); - else if (math::isnan (x)) - result = numeric_limits<double>::NaN (); - else - result = std::tgamma (x); +float erfinv (float x) +{ + return do_erfinv (x, false); +} - return result; - } +Complex +expm1 (const Complex& x) +{ + Complex retval; - float - gamma (float x) + if (std::abs (x) < 1) { - float result; + 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); - // Special cases for (near) compatibility with Matlab instead of tgamma. - // Matlab does not have -0. + return retval; +} - if (x == 0) - result = (math::negative_sign (x) - ? -numeric_limits<float>::Inf () - : numeric_limits<float>::Inf ()); - else if ((x < 0 && math::x_nint (x) == x) - || math::isinf (x)) - result = numeric_limits<float>::Inf (); - else if (math::isnan (x)) - result = numeric_limits<float>::NaN (); - else - result = std::tgammaf (x); +FloatComplex +expm1 (const FloatComplex& x) +{ + FloatComplex retval; - return result; + 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; +} - Complex - log1p (const Complex& x) - { - Complex retval; +double +gamma (double x) +{ + double result; - double r = x.real (), i = x.imag (); + // Special cases for (near) compatibility with Matlab instead of tgamma. + // Matlab does not have -0. - 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); + if (x == 0) + result = (math::negative_sign (x) + ? -numeric_limits<double>::Inf () + : numeric_limits<double>::Inf ()); + else if ((x < 0 && math::x_nint (x) == 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 retval; - } - - FloatComplex - log1p (const FloatComplex& x) - { - FloatComplex retval; + return result; +} - float r = x.real (), i = x.imag (); +float +gamma (float x) +{ + float result; - 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); + // Special cases for (near) compatibility with Matlab instead of tgamma. + // Matlab does not have -0. - return retval; - } - - static const double pi = 3.14159265358979323846; + if (x == 0) + result = (math::negative_sign (x) + ? -numeric_limits<float>::Inf () + : numeric_limits<float>::Inf ()); + else if ((x < 0 && math::x_nint (x) == x) + || math::isinf (x)) + result = numeric_limits<float>::Inf (); + else if (math::isnan (x)) + result = numeric_limits<float>::NaN (); + else + result = std::tgammaf (x); - template <typename T> - static inline T - xlog (const T& x) - { - return log (x); - } + return result; +} - template <> - inline double - xlog (const double& x) +Complex +log1p (const Complex& x) +{ + Complex retval; + + double r = x.real (), i = x.imag (); + + if (fabs (r) < 0.5 && fabs (i) < 0.5) { - return std::log (x); + 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); - template <> - inline float - xlog (const float& x) - { - return std::log (x); - } + return retval; +} - template <typename T> - static T - lanczos_approximation_psi (const T zc) +FloatComplex +log1p (const FloatComplex& x) +{ + FloatComplex retval; + + float r = x.real (), i = x.imag (); + + if (fabs (r) < 0.5 && fabs (i) < 0.5) { - // 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 - }; + 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); - T overz2 = T (1.0) / (zc * zc); - T overz2k = overz2; + return retval; +} + +static const double pi = 3.14159265358979323846; - 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> +static inline T +xlog (const T& x) +{ + return log (x); +} - template <typename T> - T - xpsi (T z) - { - static const double euler_mascheroni - = 0.577215664901532860606512090082402431042; +template <> +inline double +xlog (const double& x) +{ + return std::log (x); +} - const bool is_int = (std::floor (z) == z); +template <> +inline float +xlog (const float& x) +{ + return std::log (x); +} - 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 (std::floor (z + 0.5) == 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); +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; - p = - euler_mascheroni - 2 * std::log (2) + 2 * (p); - } + 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 - { - // adapted from XLiFE++ gammaFunctions + // 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 (std::floor (z + 0.5) == 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); - 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; + p = - euler_mascheroni - 2 * std::log (2) + 2 * (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) + else { // 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 + T zc = z; + // Use formula for derivative of LogGamma(z) + if (z < 10) { - // 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); + const signed char n = 10 - z; + for (signed char k = n - 1; k >= 0; k--) + p -= 1.0 / (k + z); + zc += n; } - 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; + p += lanczos_approximation_psi (zc); } - template <typename T> - T - xpsi (octave_idx_type n, T z) + 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 { - T ans; - octave_idx_type ierr = 0; - fortran_psifn<T> (z, n, ans, ierr); - if (ierr == 0) + // Use formula for derivative of LogGamma(z) + std::complex<T> z_m = z; + if (z_ra < 8) { - // 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; + 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; } - else if (ierr == 2) - ans = - numeric_limits<T>::Inf (); - else // we probably never get here - ans = numeric_limits<T>::NaN (); + + // 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; +} - return ans; - } +// 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; +} - double psi (octave_idx_type n, double z) { return xpsi (n, z); } - float psi (octave_idx_type n, float z) { return xpsi (n, z); } +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; +} - Complex - rc_lgamma (double x) +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) { - double result; + // 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); + int sgngam; + result = lgamma_r (x, &sgngam); #else - result = std::lgamma (x); - int sgngam = signgam; + result = std::lgamma (x); + int sgngam = signgam; #endif - if (sgngam < 0) - return result + Complex (0., M_PI); - else - return result; - } + if (sgngam < 0) + return result + Complex (0., M_PI); + else + return result; +} - FloatComplex - rc_lgamma (float x) - { - float result; +FloatComplex +rc_lgamma (float x) +{ + float result; #if defined (HAVE_LGAMMAF_R) - int sgngam; - result = lgammaf_r (x, &sgngam); + int sgngam; + result = lgammaf_r (x, &sgngam); #else - result = std::lgammaf (x); - int sgngam = signgam; + result = std::lgammaf (x); + int sgngam = signgam; #endif - if (sgngam < 0) - return result + FloatComplex (0., M_PI); - else - return result; - } + if (sgngam < 0) + return result + FloatComplex (0., 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))); - } +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))); - } +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)