Mercurial > octave
changeset 24928:f38ac278ff7d
maint: merge stable to default.
| author | Rik <rik@octave.org> |
|---|---|
| date | Mon, 19 Mar 2018 13:50:10 -0700 |
| parents | 2f6698dd7dad (current diff) c280560d9c96 (diff) |
| children | ecb5688b875f |
| files | configure.ac libinterp/corefcn/betainc.cc liboctave/external/slatec-fn/betai.f liboctave/external/slatec-fn/dbetai.f liboctave/external/slatec-fn/dgami.f liboctave/external/slatec-fn/dgamit.f liboctave/external/slatec-fn/gami.f liboctave/external/slatec-fn/gamit.f liboctave/external/slatec-fn/xbetai.f liboctave/external/slatec-fn/xdbetai.f liboctave/external/slatec-fn/xdgami.f liboctave/external/slatec-fn/xdgamit.f liboctave/external/slatec-fn/xgmainc.f liboctave/external/slatec-fn/xsgmainc.f |
| diffstat | 43 files changed, 2594 insertions(+), 2661 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Sun Mar 18 18:43:04 2018 -0700 +++ b/NEWS Mon Mar 19 13:50:10 2018 -0700 @@ -44,6 +44,15 @@ hex2num would accept a cell array of strings of arbitrary dimension but would always return a column vector. + ** New special functions cosint, sinint, and gammaincinv have been added. + + ** Special functions in Octave have been rewritten for larger input + domains, better accuracy, and additional options. + * gammainc now accepts negative real values for X. + * improved accuracy for gammainc, betainc, betaincinv, expint. + * gammainc has new options "scaledlower" and "scaledupper". + * betainc, betaincinv have new option "upper". + ** The "names" option used in regular expressions now returns a struct array, rather than a struct with a cell array for each field. This change was made for Matlab compatibility. @@ -252,7 +261,9 @@ camva camzoom corrcoef + cosint erase + gammaincinv getframe groot gsvd @@ -269,6 +280,7 @@ repelem rgb2gray rticks + sinint tfqmr thetaticks vecnorm
--- a/configure.ac Sun Mar 18 18:43:04 2018 -0700 +++ b/configure.ac Mon Mar 19 13:50:10 2018 -0700 @@ -172,11 +172,9 @@ ## Where Octave will search for fallback font files shipped with distribution. OCTAVE_SET_DEFAULT([octfontsdir], '${datadir}/octave/${version}/fonts') - ## Where Octave will look for startup files. OCTAVE_SET_DEFAULT([startupfiledir], '${fcnfiledir}/startup') OCTAVE_SET_DEFAULT([localstartupfiledir], '${localfcnfiledir}/startup') - ## Where Octave will look for man and info files. OCTAVE_SET_DEFAULT([man1dir], '${mandir}/man1') OCTAVE_SET_DEFAULT([man1ext], '.1') @@ -2908,7 +2906,6 @@ AM_CONDITIONAL([AMCOND_BUILD_DOCS], [test $ENABLE_DOCS = yes]) ### Determine whether Mercurial ID should be embedded in library binaries. - ENABLE_HG_ID=yes AC_ARG_ENABLE([hg-id], [AS_HELP_STRING([--disable-hg-id],
--- a/doc/interpreter/arith.txi Sun Mar 18 18:43:04 2018 -0700 +++ b/doc/interpreter/arith.txi Mon Mar 19 13:50:10 2018 -0700 @@ -295,6 +295,8 @@ @DOCSTRING(commutation_matrix) +@DOCSTRING(cosint) + @DOCSTRING(duplication_matrix) @DOCSTRING(dawson) @@ -321,6 +323,8 @@ @DOCSTRING(gammainc) +@DOCSTRING(gammaincinv) + @DOCSTRING(legendre) @anchor{XREFgammaln} @@ -328,6 +332,8 @@ @DOCSTRING(psi) +@DOCSTRING(sinint) + @node Rational Approximations @section Rational Approximations
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/corefcn/__betainc__.cc Mon Mar 19 13:50:10 2018 -0700 @@ -0,0 +1,202 @@ +/* + +Copyright (C) 2018 Michele Ginesi +Copyright (C) 2018 Stefan Schlögl + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 3 of the License, or +(at your option) any later version. + +Octave is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#if defined (HAVE_CONFIG_H) +# include "config.h" +#endif + +#include "defun.h" +#include "dNDArray.h" +#include "fNDArray.h" + +DEFUN (__betainc__, args, , + doc: /* -*- texinfo -*- +@deftypefn {} {@var{y} =} __betainc__ (@var{x}, @var{a}, @var{b}) +Continued fraction for incomplete beta function. +@end deftypefn */) +{ + int nargin = args.length (); + + if (nargin != 3) + print_usage (); + + bool is_single = (args(0).is_single_type () || args(1).is_single_type () + || args(2).is_single_type ()); + + // Total number of scenarios: get maximum of length of all vectors + int numel_x = args(0).numel (); + int numel_a = args(1).numel (); + int numel_b = args(2).numel (); + int len = std::max (std::max (numel_x, numel_a), numel_b); + + octave_value_list retval; + // Initialize output dimension vector + dim_vector output_dv (len, 1); + + // Lentz's algorithm in two cases: single and double precision + if (is_single) + { + // Initialize output and inputs + FloatColumnVector output (output_dv); + FloatNDArray x, a, b; + + if (numel_x == 1) + x = FloatNDArray (output_dv, args(0).float_scalar_value ()); + else + x = args(0).float_array_value (); + + + if (numel_a == 1) + a = FloatNDArray (output_dv, args(1).float_scalar_value ()); + else + a = args(1).float_array_value (); + + if (numel_b == 1) + b = FloatNDArray (output_dv, args(2).float_scalar_value ()); + else + b = args(2).float_array_value (); + + // Initialize variables used in algorithm + static const float tiny = pow (2, -50); + static const float eps = std::numeric_limits<float>::epsilon (); + float xj, x2, y, Cj, Dj, aj, bj, Deltaj, alpha_j, beta_j; + int j, maxit; + maxit = 200; + + // Loop over all elements + for (octave_idx_type i = 0; i < len; ++i) + { + // Catch Ctrl+C + OCTAVE_QUIT; + + // Variable initialization for the current element + xj = x(i); + y = tiny; + Cj = y; + Dj = 0; + aj = a(i); + bj = b(i); + Deltaj = 0; + alpha_j = 1; + beta_j = aj - (aj * (aj + bj)) / (aj + 1) * xj; + x2 = xj * xj; + j = 1; + + // Lentz's algorithm + while ((std::abs ((Deltaj - 1)) > eps) && (j < maxit)) + { + Dj = beta_j + alpha_j * Dj; + if (Dj == 0) + Dj = tiny; + Cj = beta_j + alpha_j / Cj; + if (Cj == 0) + Cj = tiny; + Dj = 1 / Dj; + Deltaj = Cj * Dj; + y *= Deltaj; + alpha_j = ((aj + j - 1) * (aj + bj + j -1) * (bj - j) * j) + / ((aj + 2 * j - 1) * (aj + 2 * j - 1)) * x2; + beta_j = aj + 2 * j + ((j * (bj - j)) / (aj + 2 * j - 1) + - ((aj + j) * (aj + bj + j)) / (aj + 2 * j + 1)) * xj; + j++; + } + + output(i) = y; + } + + retval(0) = output; + } + else + { + // Initialize output and inputs + ColumnVector output (output_dv); + NDArray x, a, b; + + if (numel_x == 1) + x = NDArray (output_dv, args(0).scalar_value ()); + else + x = args(0).array_value (); + + if (numel_a == 1) + a = NDArray (output_dv, args(1).scalar_value ()); + else + a = args(1).array_value (); + + if (numel_b == 1) + b = NDArray (output_dv, args(2).scalar_value ()); + else + b = args(2).array_value (); + + // Initialize variables used in algorithm + static const double tiny = pow (2, -100); + static const double eps = std::numeric_limits<double>::epsilon (); + double xj, x2, y, Cj, Dj, aj, bj, Deltaj, alpha_j, beta_j; + int j, maxit; + maxit = 200; + + // Loop over all elements + for (octave_idx_type i = 0; i < len; ++i) + { + // Catch Ctrl+C + OCTAVE_QUIT; + + // Variable initialization for the current element + xj = x(i); + y = tiny; + Cj = y; + Dj = 0; + aj = a(i); + bj = b(i); + Deltaj = 0; + alpha_j = 1; + beta_j = aj - (aj * (aj + bj)) / (aj + 1) * xj; + x2 = xj * xj; + j = 1; + + // Lentz's algorithm + while ((std::abs ((Deltaj - 1)) > eps) && (j < maxit)) + { + Dj = beta_j + alpha_j * Dj; + if (Dj == 0) + Dj = tiny; + Cj = beta_j + alpha_j / Cj; + if (Cj == 0) + Cj = tiny; + Dj = 1 / Dj; + Deltaj = Cj * Dj; + y *= Deltaj; + alpha_j = ((aj + j - 1) * (aj + bj + j - 1) * (bj - j) * j) + / ((aj + 2 * j - 1) * (aj + 2 * j - 1)) * x2; + beta_j = aj + 2 * j + ((j * (bj - j)) / (aj + 2 * j - 1) + - ((aj + j) * (aj + bj + j)) / (aj + 2 * j + 1)) * xj; + j++; + } + + output(i) = y; + } + + retval(0) = output; + } + + return retval; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/corefcn/__expint__.cc Mon Mar 19 13:50:10 2018 -0700 @@ -0,0 +1,187 @@ +/* + +Copyright (C) 2018 Michele Ginesi + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 3 of the License, or +(at your option) any later version. + +Octave is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#if defined (HAVE_CONFIG_H) +# include "config.h" +#endif + +#include "CNDArray.h" +#include "defun.h" +#include "fCNDArray.h" + +DEFUN (__expint__, args, , + doc: /* -*- texinfo -*- +@deftypefn {} {@var{y} =} __expint__ (@var{x}) +Continued fraction expansion for the exponential integral. +@end deftypefn */) +{ + int nargin = args.length (); + + if (nargin != 1) + print_usage (); + + octave_value_list retval; + + bool is_single = args(0).is_single_type (); + + int numel_x = args(0).numel (); + + // Initialize output dimension vector + dim_vector output_dv (numel_x, 1); + + // Lentz's algorithm in two cases: single and double precision + if (is_single) + { + // Initialize output and inputs + FloatComplexColumnVector output (output_dv); + FloatComplexNDArray x; + + if (numel_x == 1) + x = FloatComplexNDArray (output_dv, args(0).float_complex_value ()); + else + x = args(0).float_complex_array_value (); + + // Initialize variables used in algorithm + static const FloatComplex tiny = pow (2, -50); + static const float eps = std::numeric_limits<float>::epsilon (); + FloatComplex cone (1.0, 0.0); + FloatComplex czero (0.0, 0.0); + FloatComplex xj = x(0); + FloatComplex y = tiny; + FloatComplex Cj = y; + FloatComplex Dj = czero; + FloatComplex alpha_j = cone; + FloatComplex beta_j = czero; + FloatComplex Deltaj = czero; + int j = 1; + int maxit = 100; + + // Loop over all elements + for (octave_idx_type i = 0; i < numel_x; ++i) + { + // Catch Ctrl+C + OCTAVE_QUIT; + + // Variable initialization for the current element + xj = x(i); + y = tiny; + Cj = y; + Dj = czero; + alpha_j = cone; + beta_j = xj; + Deltaj = czero; + j = 1; + + // Lentz's algorithm + while ((std::abs (Deltaj - cone) > eps) && (j < maxit)) + { + Dj = beta_j + alpha_j * Dj; + if (Dj == czero) + Dj = tiny; + Cj = beta_j + alpha_j / Cj; + if (Cj == czero) + Cj = tiny; + Dj = cone / Dj; + Deltaj = Cj * Dj; + y *= Deltaj; + alpha_j = floor ((j + 1) / 2); + if ((j % 2) == 0) + beta_j = xj; + else + beta_j = cone; + j++; + } + + output(i) = y; + } + retval(0) = output; + } + else + { + // Initialize output and inputs + ComplexColumnVector output (output_dv); + ComplexNDArray x; + + if (numel_x == 1) + x = ComplexNDArray (output_dv, args(0).complex_value ()); + else + x = args(0).complex_array_value (); + + // Initialize variables used in algorithm + static const Complex tiny = pow (2, -100); + static const double eps = std::numeric_limits<double>::epsilon (); + Complex cone (1.0, 0.0); + Complex czero (0.0, 0.0); + Complex xj = x(0); + Complex y = tiny; + Complex Cj = y; + Complex Dj = czero; + Complex alpha_j = cone; + Complex beta_j = xj; + Complex Deltaj = czero; + int j = 1; + int maxit = 200; + + // Loop over all scenarios + for (octave_idx_type i = 0; i < numel_x; ++i) + { + // Catch Ctrl+C + OCTAVE_QUIT; + + // Variable initialization for the current element + xj = x(i); + y = tiny; + Cj = y; + Dj = czero; + alpha_j = cone; + beta_j = xj; + Deltaj = czero; + j = 1; + + // Lentz's algorithm + while ((std::abs (Deltaj - cone) > eps) && (j < maxit)) + { + Dj = beta_j + alpha_j * Dj; + if (Dj == czero) + Dj = tiny; + Cj = beta_j + alpha_j / Cj; + if (Cj == czero) + Cj = tiny; + Dj = cone / Dj; + Deltaj = Cj * Dj; + y *= Deltaj; + alpha_j = floor ((j + 1) / 2); + if ((j % 2) == 0) + beta_j = xj; + else + beta_j = cone; + j++; + } + + output(i) = y; + } + + retval(0) = output; + } + + return retval; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/corefcn/__gammainc__.cc Mon Mar 19 13:50:10 2018 -0700 @@ -0,0 +1,167 @@ +/* + +Copyright (C) 2017 Nir Krakauer +Copyright (C) 2018 Michele Ginesi +Copyright (C) 2018 Stefan Schlögl + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 3 of the License, or +(at your option) any later version. + +Octave is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#if defined (HAVE_CONFIG_H) +# include "config.h" +#endif + +#include "defun.h" +#include "fNDArray.h" + +DEFUN (__gammainc__, args, , + doc: /* -*- texinfo -*- +@deftypefn {} {@var{y} =} __gammainc__ (@var{x}, @var{a}) +Continued fraction for incomplete gamma function. +@end deftypefn */) +{ + int nargin = args.length (); + + if (nargin != 2) + print_usage (); + + bool is_single = args(0).is_single_type () || args(1).is_single_type (); + + // Total number of scenarios: get maximum of length of all vectors + int numel_x = args(0).numel (); + int numel_a = args(1).numel (); + int len = std::max (numel_x, numel_a); + + octave_value_list retval; + // Initialize output dimension vector + dim_vector output_dv (len, 1); + + // Lentz's algorithm in two cases: single and double precision + if (is_single) + { + // Initialize output and inputs + FloatColumnVector output (output_dv); + FloatNDArray x, a; + + if (numel_x == 1) + x = FloatNDArray (output_dv, args(0).float_scalar_value ()); + else + x = args(0).float_array_value (); + + if (numel_a == 1) + a = FloatNDArray (output_dv, args(1).float_scalar_value ()); + else + a = args(1).float_array_value (); + + // Initialize variables used in algorithm + static const float tiny = pow (2, -50); + static const float eps = std::numeric_limits<float>::epsilon(); + float y, Cj, Dj, bj, aj, Deltaj; + int j, maxit; + maxit = 200; + + // Loop over all elements + for (octave_idx_type i = 0; i < len; ++i) + { + // Catch Ctrl+C + OCTAVE_QUIT; + + // Variable initialization for the current element + y = tiny; + Cj = y; + Dj = 0; + bj = x(i) - a(i) + 1; + aj = a(i); + Deltaj = 0; + j = 1; + + // Lentz's algorithm + while ((std::abs ((Deltaj - 1) / y) > eps) && (j < maxit)) + { + Cj = bj + aj/Cj; + Dj = 1 / (bj + aj*Dj); + Deltaj = Cj * Dj; + y *= Deltaj; + bj += 2; + aj = j * (a(i) - j); + j++; + } + + output(i) = y; + } + + retval(0) = output; + } + else + { + // Initialize output and inputs + ColumnVector output (output_dv); + NDArray x, a; + + if (numel_x == 1) + x = NDArray (output_dv, args(0).scalar_value ()); + else + x = args(0).array_value (); + + if (numel_a == 1) + a = NDArray (output_dv, args(1).scalar_value ()); + else + a = args(1).array_value (); + + // Initialize variables used in algorithm + static const double tiny = pow (2, -100); + static const double eps = std::numeric_limits<double>::epsilon(); + double y, Cj, Dj, bj, aj, Deltaj; + int j, maxit; + maxit = 200; + + // Loop over all elements + for (octave_idx_type i = 0; i < len; ++i) + { + // Catch Ctrl+C + OCTAVE_QUIT; + + // Variable initialization for the current element + y = tiny; + Cj = y; + Dj = 0; + bj = x(i) - a(i) + 1; + aj = a(i); + Deltaj = 0; + j = 1; + + // Lentz's algorithm + while ((std::abs ((Deltaj - 1) / y) > eps) && (j < maxit)) + { + Cj = bj + aj/Cj; + Dj = 1 / (bj + aj*Dj); + Deltaj = Cj * Dj; + y *= Deltaj; + bj += 2; + aj = j * (a(i) - j); + j++; + } + + output(i) = y; + } + + retval(0) = output; + } + + return retval; +}
--- a/libinterp/corefcn/besselj.cc Sun Mar 18 18:43:04 2018 -0700 +++ b/libinterp/corefcn/besselj.cc Mon Mar 19 13:50:10 2018 -0700 @@ -335,7 +335,7 @@ accuracy. @item -Complete loss of significance by argument reduction, return @code{NaN}. +Loss of significance by argument reduction, output may be inaccurate. @item Error---no computation, algorithm termination condition not met, return @@ -648,7 +648,7 @@ of machine accuracy. @item -Complete loss of significance by argument reduction, return @code{NaN}. +Loss of significance by argument reduction, output may be inaccurate. @item Error---no computation, algorithm termination condition not met,
--- a/libinterp/corefcn/betainc.cc Sun Mar 18 18:43:04 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,424 +0,0 @@ -/* - -Copyright (C) 1997-2017 John W. Eaton - -This file is part of Octave. - -Octave is free software: you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by -the Free Software Foundation, either version 3 of the License, or -(at your option) any later version. - -Octave is distributed in the hope that it will be useful, but -WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU General Public License for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<https://www.gnu.org/licenses/>. - -*/ - -#if defined (HAVE_CONFIG_H) -# include "config.h" -#endif - -#include "lo-specfun.h" - -#include "defun.h" -#include "error.h" -#include "errwarn.h" -#include "ovl.h" -#include "utils.h" - - -DEFUN (betainc, args, , - doc: /* -*- texinfo -*- -@deftypefn {} {} betainc (@var{x}, @var{a}, @var{b}) -Compute the regularized incomplete Beta function. - -The regularized incomplete Beta function is defined by -@tex -$$ - I (x, a, b) = {1 \over {B (a, b)}} \int_0^x t^{(a-z)} (1-t)^{(b-1)} dt. -$$ -@end tex -@ifnottex -@c Set example in small font to prevent overfull line - -@smallexample -@group - x - 1 / -betainc (x, a, b) = ----------- | t^(a-1) (1-t)^(b-1) dt. - beta (a, b) / - t=0 -@end group -@end smallexample - -@end ifnottex - -If @var{x} has more than one component, both @var{a} and @var{b} must be -scalars. If @var{x} is a scalar, @var{a} and @var{b} must be of -compatible dimensions. -@seealso{betaincinv, beta, betaln} -@end deftypefn */) -{ - if (args.length () != 3) - print_usage (); - - octave_value retval; - - octave_value x_arg = args(0); - octave_value a_arg = args(1); - octave_value b_arg = args(2); - - // FIXME: Can we make a template version of the duplicated code below - if (x_arg.is_single_type () || a_arg.is_single_type () - || b_arg.is_single_type ()) - { - if (x_arg.is_scalar_type ()) - { - float x = x_arg.float_value (); - - if (a_arg.is_scalar_type ()) - { - float a = a_arg.float_value (); - - if (b_arg.is_scalar_type ()) - { - float b = b_arg.float_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<float> b = b_arg.float_array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - else - { - Array<float> a = a_arg.float_array_value (); - - if (b_arg.is_scalar_type ()) - { - float b = b_arg.float_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<float> b = b_arg.float_array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - } - else - { - Array<float> x = x_arg.float_array_value (); - - if (a_arg.is_scalar_type ()) - { - float a = a_arg.float_value (); - - if (b_arg.is_scalar_type ()) - { - float b = b_arg.float_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<float> b = b_arg.float_array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - else - { - Array<float> a = a_arg.float_array_value (); - - if (b_arg.is_scalar_type ()) - { - float b = b_arg.float_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<float> b = b_arg.float_array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - } - } - else - { - if (x_arg.is_scalar_type ()) - { - double x = x_arg.double_value (); - - if (a_arg.is_scalar_type ()) - { - double a = a_arg.double_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - else - { - Array<double> a = a_arg.array_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - } - else - { - Array<double> x = x_arg.array_value (); - - if (a_arg.is_scalar_type ()) - { - double a = a_arg.double_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - else - { - Array<double> a = a_arg.array_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betainc (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betainc (x, a, b); - } - } - } - } - - return retval; -} - -/* -## Double precision -%!test -%! a = [1, 1.5, 2, 3]; -%! b = [4, 3, 2, 1]; -%! v1 = betainc (1,a,b); -%! v2 = [1,1,1,1]; -%! x = [.2, .4, .6, .8]; -%! v3 = betainc (x, a, b); -%! v4 = 1 - betainc (1.-x, b, a); -%! assert (v1, v2, sqrt (eps)); -%! assert (v3, v4, sqrt (eps)); - -## Single precision -%!test -%! a = single ([1, 1.5, 2, 3]); -%! b = single ([4, 3, 2, 1]); -%! v1 = betainc (1,a,b); -%! v2 = single ([1,1,1,1]); -%! x = single ([.2, .4, .6, .8]); -%! v3 = betainc (x, a, b); -%! v4 = 1 - betainc (1.-x, b, a); -%! assert (v1, v2, sqrt (eps ("single"))); -%! assert (v3, v4, sqrt (eps ("single"))); - -## Mixed double/single precision -%!test -%! a = single ([1, 1.5, 2, 3]); -%! b = [4, 3, 2, 1]; -%! v1 = betainc (1,a,b); -%! v2 = single ([1,1,1,1]); -%! x = [.2, .4, .6, .8]; -%! v3 = betainc (x, a, b); -%! v4 = 1-betainc (1.-x, b, a); -%! assert (v1, v2, sqrt (eps ("single"))); -%! assert (v3, v4, sqrt (eps ("single"))); - -%!error betainc () -%!error betainc (1) -%!error betainc (1,2) -%!error betainc (1,2,3,4) -*/ - -DEFUN (betaincinv, args, , - doc: /* -*- texinfo -*- -@deftypefn {} {} betaincinv (@var{y}, @var{a}, @var{b}) -Compute the inverse of the incomplete Beta function. - -The inverse is the value @var{x} such that - -@example -@var{y} == betainc (@var{x}, @var{a}, @var{b}) -@end example -@seealso{betainc, beta, betaln} -@end deftypefn */) -{ - if (args.length () != 3) - print_usage (); - - octave_value retval; - - octave_value x_arg = args(0); - octave_value a_arg = args(1); - octave_value b_arg = args(2); - - if (x_arg.is_scalar_type ()) - { - double x = x_arg.double_value (); - - if (a_arg.is_scalar_type ()) - { - double a = a_arg.double_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betaincinv (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betaincinv (x, a, b); - } - } - else - { - Array<double> a = a_arg.array_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betaincinv (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betaincinv (x, a, b); - } - } - } - else - { - Array<double> x = x_arg.array_value (); - - if (a_arg.is_scalar_type ()) - { - double a = a_arg.double_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betaincinv (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betaincinv (x, a, b); - } - } - else - { - Array<double> a = a_arg.array_value (); - - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - retval = octave::math::betaincinv (x, a, b); - } - else - { - Array<double> b = b_arg.array_value (); - - retval = octave::math::betaincinv (x, a, b); - } - } - } - - // FIXME: It would be better to have an algorithm for betaincinv which - // accepted float inputs and returned float outputs. As it is, we do - // extra work to calculate betaincinv to double precision and then throw - // that precision away. - if (x_arg.is_single_type () || a_arg.is_single_type () - || b_arg.is_single_type ()) - { - retval = Array<float> (retval.array_value ()); - } - - return retval; -} - -/* -%!assert (betaincinv ([0.875 0.6875], [1 2], 3), [0.5 0.5], sqrt (eps)) -%!assert (betaincinv (0.5, 3, 3), 0.5, sqrt (eps)) -%!assert (betaincinv (0.34375, 4, 3), 0.5, sqrt (eps)) -%!assert (betaincinv (0.2265625, 5, 3), 0.5, sqrt (eps)) -%!assert (betaincinv (0.14453125, 6, 3), 0.5, sqrt (eps)) -%!assert (betaincinv (0.08984375, 7, 3), 0.5, sqrt (eps)) -%!assert (betaincinv (0.0546875, 8, 3), 0.5, sqrt (eps)) -%!assert (betaincinv (0.03271484375, 9, 3), 0.5, sqrt (eps)) -%!assert (betaincinv (0.019287109375, 10, 3), 0.5, sqrt (eps)) - -## Test class single as well -%!assert (betaincinv ([0.875 0.6875], [1 2], single (3)), [0.5 0.5], sqrt (eps ("single"))) -%!assert (betaincinv (0.5, 3, single (3)), 0.5, sqrt (eps ("single"))) -%!assert (betaincinv (0.34375, 4, single (3)), 0.5, sqrt (eps ("single"))) - -## Extreme values -%!assert (betaincinv (0, 42, 42), 0, sqrt (eps)) -%!assert (betaincinv (1, 42, 42), 1, sqrt (eps)) - -%!error betaincinv () -%!error betaincinv (1, 2) -*/
--- a/libinterp/corefcn/module.mk Sun Mar 18 18:43:04 2018 -0700 +++ b/libinterp/corefcn/module.mk Mon Mar 19 13:50:10 2018 -0700 @@ -107,8 +107,11 @@ COREFCN_SRC = \ %reldir%/Cell.cc \ + %reldir%/__betainc__.cc \ %reldir%/__contourc__.cc \ %reldir%/__dsearchn__.cc \ + %reldir%/__expint__.cc \ + %reldir%/__gammainc__.cc \ %reldir%/__ichol__.cc \ %reldir%/__ilu__.cc \ %reldir%/__lin_interpn__.cc \ @@ -118,7 +121,6 @@ %reldir%/__qp__.cc \ %reldir%/balance.cc \ %reldir%/besselj.cc \ - %reldir%/betainc.cc \ %reldir%/bitfcns.cc \ %reldir%/bsxfun.cc \ %reldir%/c-file-ptr-stream.cc \ @@ -155,7 +157,6 @@ %reldir%/filter.cc \ %reldir%/find.cc \ %reldir%/ft-text-renderer.cc \ - %reldir%/gammainc.cc \ %reldir%/gcd.cc \ %reldir%/getgrent.cc \ %reldir%/getpwent.cc \
--- a/liboctave/external/amos/README Sun Mar 18 18:43:04 2018 -0700 +++ b/liboctave/external/amos/README Mon Mar 19 13:50:10 2018 -0700 @@ -13,3 +13,21 @@ jwe@octave.org Wed Nov 11 17:29:50 1998 + +More files have been modified to recover Matlab compatibility for +Bessel functions. Now the output is always computed, independently +from the magnitude of the input + + cbesh.f + cbesi.f + cbesj.f + cbesk.f + zbesh.f + zbesi.f + zbesj.f + zbesk.f + +Michele Ginesi +michele.ginesi@gmail.com + +Sat Jul 22 11:43:40 2017
--- a/liboctave/external/amos/cbesh.f Sun Mar 18 18:43:04 2018 -0700 +++ b/liboctave/external/amos/cbesh.f Mon Mar 19 13:50:10 2018 -0700 @@ -219,6 +219,7 @@ C----------------------------------------------------------------------- C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE C----------------------------------------------------------------------- + 35 CONTINUE UFL = R1MACH(1)*1.0E+3 IF (AZ.LT.UFL) GO TO 220 IF (FNU.GT.FNUL) GO TO 90 @@ -325,7 +326,6 @@ IERR=5 RETURN 240 CONTINUE - NZ=0 IERR=4 - RETURN + GO TO 35 END
--- a/liboctave/external/amos/cbesi.f Sun Mar 18 18:43:04 2018 -0700 +++ b/liboctave/external/amos/cbesi.f Mon Mar 19 13:50:10 2018 -0700 @@ -201,6 +201,7 @@ AA=SQRT(AA) IF(AZ.GT.AA) IERR=3 IF(FN.GT.AA) IERR=3 + 35 CONTINUE ZN = Z CSGN = CONE IF (XX.GE.0.0E0) GO TO 40 @@ -252,7 +253,6 @@ IERR=5 RETURN 140 CONTINUE - NZ=0 IERR=4 - RETURN + GO TO 35 END
--- a/liboctave/external/amos/cbesj.f Sun Mar 18 18:43:04 2018 -0700 +++ b/liboctave/external/amos/cbesj.f Mon Mar 19 13:50:10 2018 -0700 @@ -199,6 +199,7 @@ C CALCULATE CSGN=EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE C WHEN FNU IS LARGE C----------------------------------------------------------------------- + 35 CONTINUE INU = INT(FNU) INUH = INU/2 IR = INU - 2*INUH @@ -247,7 +248,6 @@ IERR=5 RETURN 140 CONTINUE - NZ=0 IERR=4 - RETURN + GO TO 35 END
--- a/liboctave/external/amos/cbesk.f Sun Mar 18 18:43:04 2018 -0700 +++ b/liboctave/external/amos/cbesk.f Mon Mar 19 13:50:10 2018 -0700 @@ -202,6 +202,7 @@ C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE C----------------------------------------------------------------------- C UFL = EXP(-ELIM) + 35 CONTINUE UFL = R1MACH(1)*1.0E+3 IF (AZ.LT.UFL) GO TO 180 IF (FNU.GT.FNUL) GO TO 80 @@ -270,7 +271,6 @@ IERR=5 RETURN 210 CONTINUE - NZ=0 IERR=4 - RETURN + GO TO 35 END
--- a/liboctave/external/amos/zbesh.f Sun Mar 18 18:43:04 2018 -0700 +++ b/liboctave/external/amos/zbesh.f Mon Mar 19 13:50:10 2018 -0700 @@ -220,6 +220,7 @@ C----------------------------------------------------------------------- C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE C----------------------------------------------------------------------- + 35 CONTINUE UFL = D1MACH(1)*1.0D+3 IF (AZ.LT.UFL) GO TO 230 IF (FNU.GT.FNUL) GO TO 90 @@ -342,7 +343,6 @@ IERR=5 RETURN 260 CONTINUE - NZ=0 IERR=4 - RETURN + GO TO 35 END
--- a/liboctave/external/amos/zbesi.f Sun Mar 18 18:43:04 2018 -0700 +++ b/liboctave/external/amos/zbesi.f Mon Mar 19 13:50:10 2018 -0700 @@ -202,6 +202,7 @@ AA = DSQRT(AA) IF (AZ.GT.AA) IERR=3 IF (FN.GT.AA) IERR=3 + 35 CONTINUE ZNR = ZR ZNI = ZI CSGNR = CONER @@ -263,7 +264,6 @@ IERR=5 RETURN 260 CONTINUE - NZ=0 IERR=4 - RETURN + GO TO 35 END
--- a/liboctave/external/amos/zbesj.f Sun Mar 18 18:43:04 2018 -0700 +++ b/liboctave/external/amos/zbesj.f Mon Mar 19 13:50:10 2018 -0700 @@ -200,6 +200,7 @@ C CALCULATE CSGN=EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE C WHEN FNU IS LARGE C----------------------------------------------------------------------- + 35 CONTINUE CII = 1.0D0 INU = INT(SNGL(FNU)) INUH = INU/2 @@ -260,7 +261,7 @@ IERR=5 RETURN 260 CONTINUE - NZ=0 IERR=4 + GO TO 35 RETURN END
--- a/liboctave/external/amos/zbesk.f Sun Mar 18 18:43:04 2018 -0700 +++ b/liboctave/external/amos/zbesk.f Mon Mar 19 13:50:10 2018 -0700 @@ -205,6 +205,7 @@ C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE C----------------------------------------------------------------------- C UFL = DEXP(-ELIM) + 35 CONTINUE UFL = D1MACH(1)*1.0D+3 IF (AZ.LT.UFL) GO TO 180 IF (FNU.GT.FNUL) GO TO 80 @@ -275,7 +276,6 @@ IERR=5 RETURN 260 CONTINUE - NZ=0 IERR=4 - RETURN + GO TO 35 END
--- a/liboctave/external/slatec-fn/betai.f Sun Mar 18 18:43:04 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,118 +0,0 @@ -*DECK BETAI - REAL FUNCTION BETAI (X, PIN, QIN) -C***BEGIN PROLOGUE BETAI -C***PURPOSE Calculate the incomplete Beta function. -C***LIBRARY SLATEC (FNLIB) -C***CATEGORY C7F -C***TYPE SINGLE PRECISION (BETAI-S, DBETAI-D) -C***KEYWORDS FNLIB, INCOMPLETE BETA FUNCTION, SPECIAL FUNCTIONS -C***AUTHOR Fullerton, W., (LANL) -C***DESCRIPTION -C -C BETAI calculates the REAL incomplete beta function. -C -C The incomplete beta function ratio is the probability that a -C random variable from a beta distribution having parameters PIN and -C QIN will be less than or equal to X. -C -C -- Input Arguments -- All arguments are REAL. -C X upper limit of integration. X must be in (0,1) inclusive. -C PIN first beta distribution parameter. PIN must be .GT. 0.0. -C QIN second beta distribution parameter. QIN must be .GT. 0.0. -C -C***REFERENCES Nancy E. Bosten and E. L. Battiste, Remark on Algorithm -C 179, Communications of the ACM 17, 3 (March 1974), -C pp. 156. -C***ROUTINES CALLED ALBETA, R1MACH, XERMSG -C***REVISION HISTORY (YYMMDD) -C 770401 DATE WRITTEN -C 890531 Changed all specific intrinsics to generic. (WRB) -C 890531 REVISION DATE from Version 3.2 -C 891214 Prologue converted to Version 4.0 format. (BAB) -C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) -C 900326 Removed duplicate information from DESCRIPTION section. -C (WRB) -C 920528 DESCRIPTION and REFERENCES sections revised. (WRB) -C***END PROLOGUE BETAI - LOGICAL FIRST - SAVE EPS, ALNEPS, SML, ALNSML, FIRST - DATA FIRST /.TRUE./ -C***FIRST EXECUTABLE STATEMENT BETAI - IF (FIRST) THEN - EPS = R1MACH(3) - ALNEPS = LOG(EPS) - SML = R1MACH(1) - ALNSML = LOG(SML) - ENDIF - FIRST = .FALSE. -C - IF (X .LT. 0. .OR. X .GT. 1.0) CALL XERMSG ('SLATEC', 'BETAI', - + 'X IS NOT IN THE RANGE (0,1)', 1, 2) - IF (PIN .LE. 0. .OR. QIN .LE. 0.) CALL XERMSG ('SLATEC', 'BETAI', - + 'P AND/OR Q IS LE ZERO', 2, 2) -C - Y = X - P = PIN - Q = QIN - IF (Q.LE.P .AND. X.LT.0.8) GO TO 20 - IF (X.LT.0.2) GO TO 20 - Y = 1.0 - Y - P = QIN - Q = PIN -C - 20 IF ((P+Q)*Y/(P+1.).LT.EPS) GO TO 80 -C -C EVALUATE THE INFINITE SUM FIRST. -C TERM WILL EQUAL Y**P/BETA(PS,P) * (1.-PS)I * Y**I / FAC(I) -C - PS = Q - AINT(Q) - IF (PS.EQ.0.) PS = 1.0 - XB = P*LOG(Y) - ALBETA(PS, P) - LOG(P) - BETAI = 0.0 - IF (XB.LT.ALNSML) GO TO 40 -C - BETAI = EXP (XB) - TERM = BETAI*P - IF (PS.EQ.1.0) GO TO 40 -C - N = MAX (ALNEPS/LOG(Y), 4.0E0) - DO 30 I=1,N - TERM = TERM*(I-PS)*Y/I - BETAI = BETAI + TERM/(P+I) - 30 CONTINUE -C -C NOW EVALUATE THE FINITE SUM, MAYBE. -C - 40 IF (Q.LE.1.0) GO TO 70 -C - XB = P*LOG(Y) + Q*LOG(1.0-Y) - ALBETA(P,Q) - LOG(Q) - IB = MAX (XB/ALNSML, 0.0E0) - TERM = EXP (XB - IB*ALNSML) - C = 1.0/(1.0-Y) - P1 = Q*C/(P+Q-1.) -C - FINSUM = 0.0 - N = Q - IF (Q.EQ.REAL(N)) N = N - 1 - DO 50 I=1,N - IF (P1.LE.1.0 .AND. TERM/EPS.LE.FINSUM) GO TO 60 - TERM = (Q-I+1)*C*TERM/(P+Q-I) -C - IF (TERM.GT.1.0) IB = IB - 1 - IF (TERM.GT.1.0) TERM = TERM*SML -C - IF (IB.EQ.0) FINSUM = FINSUM + TERM - 50 CONTINUE -C - 60 BETAI = BETAI + FINSUM - 70 IF (Y.NE.X .OR. P.NE.PIN) BETAI = 1.0 - BETAI - BETAI = MAX (MIN (BETAI, 1.0), 0.0) - RETURN -C - 80 BETAI = 0.0 - XB = P*LOG(MAX(Y,SML)) - LOG(P) - ALBETA(P,Q) - IF (XB.GT.ALNSML .AND. Y.NE.0.) BETAI = EXP (XB) - IF (Y.NE.X .OR. P.NE.PIN) BETAI = 1.0 - BETAI - RETURN -C - END
--- a/liboctave/external/slatec-fn/dbetai.f Sun Mar 18 18:43:04 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,121 +0,0 @@ - -*DECK DBETAI - DOUBLE PRECISION FUNCTION DBETAI (X, PIN, QIN) -C***BEGIN PROLOGUE DBETAI -C***PURPOSE Calculate the incomplete Beta function. -C***LIBRARY SLATEC (FNLIB) -C***CATEGORY C7F -C***TYPE DOUBLE PRECISION (BETAI-S, DBETAI-D) -C***KEYWORDS FNLIB, INCOMPLETE BETA FUNCTION, SPECIAL FUNCTIONS -C***AUTHOR Fullerton, W., (LANL) -C***DESCRIPTION -C -C DBETAI calculates the DOUBLE PRECISION incomplete beta function. -C -C The incomplete beta function ratio is the probability that a -C random variable from a beta distribution having parameters PIN and -C QIN will be less than or equal to X. -C -C -- Input Arguments -- All arguments are DOUBLE PRECISION. -C X upper limit of integration. X must be in (0,1) inclusive. -C PIN first beta distribution parameter. PIN must be .GT. 0.0. -C QIN second beta distribution parameter. QIN must be .GT. 0.0. -C -C***REFERENCES Nancy E. Bosten and E. L. Battiste, Remark on Algorithm -C 179, Communications of the ACM 17, 3 (March 1974), -C pp. 156. -C***ROUTINES CALLED D1MACH, DLBETA, XERMSG -C***REVISION HISTORY (YYMMDD) -C 770701 DATE WRITTEN -C 890531 Changed all specific intrinsics to generic. (WRB) -C 890911 Removed unnecessary intrinsics. (WRB) -C 890911 REVISION DATE from Version 3.2 -C 891214 Prologue converted to Version 4.0 format. (BAB) -C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) -C 920528 DESCRIPTION and REFERENCES sections revised. (WRB) -C***END PROLOGUE DBETAI - DOUBLE PRECISION X, PIN, QIN, ALNEPS, ALNSML, C, EPS, FINSUM, P, - 1 PS, Q, SML, TERM, XB, XI, Y, D1MACH, DLBETA, P1 - LOGICAL FIRST - SAVE EPS, ALNEPS, SML, ALNSML, FIRST - DATA FIRST /.TRUE./ -C***FIRST EXECUTABLE STATEMENT DBETAI - IF (FIRST) THEN - EPS = D1MACH(3) - ALNEPS = LOG (EPS) - SML = D1MACH(1) - ALNSML = LOG (SML) - ENDIF - FIRST = .FALSE. -C - IF (X .LT. 0.D0 .OR. X .GT. 1.D0) CALL XERMSG ('SLATEC', 'DBETAI', - + 'X IS NOT IN THE RANGE (0,1)', 1, 2) - IF (PIN .LE. 0.D0 .OR. QIN .LE. 0.D0) CALL XERMSG ('SLATEC', - + 'DBETAI', 'P AND/OR Q IS LE ZERO', 2, 2) -C - Y = X - P = PIN - Q = QIN - IF (Q.LE.P .AND. X.LT.0.8D0) GO TO 20 - IF (X.LT.0.2D0) GO TO 20 - Y = 1.0D0 - Y - P = QIN - Q = PIN -C - 20 IF ((P+Q)*Y/(P+1.D0).LT.EPS) GO TO 80 -C -C EVALUATE THE INFINITE SUM FIRST. TERM WILL EQUAL -C Y**P/BETA(PS,P) * (1.-PS)-SUB-I * Y**I / FAC(I) . -C - PS = Q - AINT(Q) - IF (PS.EQ.0.D0) PS = 1.0D0 - XB = P*LOG(Y) - DLBETA(PS,P) - LOG(P) - DBETAI = 0.0D0 - IF (XB.LT.ALNSML) GO TO 40 -C - DBETAI = EXP (XB) - TERM = DBETAI*P - IF (PS.EQ.1.0D0) GO TO 40 - N = MAX (ALNEPS/LOG(Y), 4.0D0) - DO 30 I=1,N - XI = I - TERM = TERM * (XI-PS)*Y/XI - DBETAI = DBETAI + TERM/(P+XI) - 30 CONTINUE -C -C NOW EVALUATE THE FINITE SUM, MAYBE. -C - 40 IF (Q.LE.1.0D0) GO TO 70 -C - XB = P*LOG(Y) + Q*LOG(1.0D0-Y) - DLBETA(P,Q) - LOG(Q) - IB = MAX (XB/ALNSML, 0.0D0) - TERM = EXP(XB - IB*ALNSML) - C = 1.0D0/(1.D0-Y) - P1 = Q*C/(P+Q-1.D0) -C - FINSUM = 0.0D0 - N = Q - IF (Q.EQ.DBLE(N)) N = N - 1 - DO 50 I=1,N - IF (P1.LE.1.0D0 .AND. TERM/EPS.LE.FINSUM) GO TO 60 - XI = I - TERM = (Q-XI+1.0D0)*C*TERM/(P+Q-XI) -C - IF (TERM.GT.1.0D0) IB = IB - 1 - IF (TERM.GT.1.0D0) TERM = TERM*SML -C - IF (IB.EQ.0) FINSUM = FINSUM + TERM - 50 CONTINUE -C - 60 DBETAI = DBETAI + FINSUM - 70 IF (Y.NE.X .OR. P.NE.PIN) DBETAI = 1.0D0 - DBETAI - DBETAI = MAX (MIN (DBETAI, 1.0D0), 0.0D0) - RETURN -C - 80 DBETAI = 0.0D0 - XB = P*LOG(MAX(Y,SML)) - LOG(P) - DLBETA(P,Q) - IF (XB.GT.ALNSML .AND. Y.NE.0.0D0) DBETAI = EXP(XB) - IF (Y.NE.X .OR. P.NE.PIN) DBETAI = 1.0D0 - DBETAI -C - RETURN - END
--- a/liboctave/external/slatec-fn/dgami.f Sun Mar 18 18:43:04 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +0,0 @@ - -*DECK DGAMI - DOUBLE PRECISION FUNCTION DGAMI (A, X) -C***BEGIN PROLOGUE DGAMI -C***PURPOSE Evaluate the incomplete Gamma function. -C***LIBRARY SLATEC (FNLIB) -C***CATEGORY C7E -C***TYPE DOUBLE PRECISION (GAMI-S, DGAMI-D) -C***KEYWORDS FNLIB, INCOMPLETE GAMMA FUNCTION, SPECIAL FUNCTIONS -C***AUTHOR Fullerton, W., (LANL) -C***DESCRIPTION -C -C Evaluate the incomplete gamma function defined by -C -C DGAMI = integral from T = 0 to X of EXP(-T) * T**(A-1.0) . -C -C DGAMI is evaluated for positive values of A and non-negative values -C of X. A slight deterioration of 2 or 3 digits accuracy will occur -C when DGAMI is very large or very small, because logarithmic variables -C are used. The function and both arguments are double precision. -C -C***REFERENCES (NONE) -C***ROUTINES CALLED DGAMIT, DLNGAM, XERMSG -C***REVISION HISTORY (YYMMDD) -C 770701 DATE WRITTEN -C 890531 Changed all specific intrinsics to generic. (WRB) -C 890531 REVISION DATE from Version 3.2 -C 891214 Prologue converted to Version 4.0 format. (BAB) -C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) -C***END PROLOGUE DGAMI - DOUBLE PRECISION A, X, FACTOR, DLNGAM, DGAMIT -C***FIRST EXECUTABLE STATEMENT DGAMI - IF (A .LE. 0.D0) CALL XERMSG ('SLATEC', 'DGAMI', - + 'A MUST BE GT ZERO', 1, 2) - IF (X .LT. 0.D0) CALL XERMSG ('SLATEC', 'DGAMI', - + 'X MUST BE GE ZERO', 2, 2) -C - DGAMI = 0.D0 - IF (X.EQ.0.0D0) RETURN -C -C THE ONLY ERROR POSSIBLE IN THE EXPRESSION BELOW IS A FATAL OVERFLOW. - FACTOR = EXP (DLNGAM(A) + A*LOG(X)) -C - DGAMI = FACTOR * DGAMIT (A, X) -C - RETURN - END
--- a/liboctave/external/slatec-fn/dgamit.f Sun Mar 18 18:43:04 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,119 +0,0 @@ -*DECK DGAMIT - DOUBLE PRECISION FUNCTION DGAMIT (A, X) -C***BEGIN PROLOGUE DGAMIT -C***PURPOSE Calculate Tricomi's form of the incomplete Gamma function. -C***LIBRARY SLATEC (FNLIB) -C***CATEGORY C7E -C***TYPE DOUBLE PRECISION (GAMIT-S, DGAMIT-D) -C***KEYWORDS COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB, -C SPECIAL FUNCTIONS, TRICOMI -C***AUTHOR Fullerton, W., (LANL) -C***DESCRIPTION -C -C Evaluate Tricomi's incomplete Gamma function defined by -C -C DGAMIT = X**(-A)/GAMMA(A) * integral from 0 to X of EXP(-T) * -C T**(A-1.) -C -C for A .GT. 0.0 and by analytic continuation for A .LE. 0.0. -C GAMMA(X) is the complete gamma function of X. -C -C DGAMIT is evaluated for arbitrary real values of A and for non- -C negative values of X (even though DGAMIT is defined for X .LT. -C 0.0), except that for X = 0 and A .LE. 0.0, DGAMIT is infinite, -C which is a fatal error. -C -C The function and both arguments are DOUBLE PRECISION. -C -C A slight deterioration of 2 or 3 digits accuracy will occur when -C DGAMIT is very large or very small in absolute value, because log- -C arithmic variables are used. Also, if the parameter A is very -C close to a negative integer (but not a negative integer), there is -C a loss of accuracy, which is reported if the result is less than -C half machine precision. -C -C***REFERENCES W. Gautschi, A computational procedure for incomplete -C gamma functions, ACM Transactions on Mathematical -C Software 5, 4 (December 1979), pp. 466-481. -C W. Gautschi, Incomplete gamma functions, Algorithm 542, -C ACM Transactions on Mathematical Software 5, 4 -C (December 1979), pp. 482-489. -C***ROUTINES CALLED D1MACH, D9GMIT, D9LGIC, D9LGIT, DGAMR, DLGAMS, -C DLNGAM, XERCLR, XERMSG -C***REVISION HISTORY (YYMMDD) -C 770701 DATE WRITTEN -C 890531 Changed all specific intrinsics to generic. (WRB) -C 890531 REVISION DATE from Version 3.2 -C 891214 Prologue converted to Version 4.0 format. (BAB) -C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) -C 920528 DESCRIPTION and REFERENCES sections revised. (WRB) -C***END PROLOGUE DGAMIT - DOUBLE PRECISION A, X, AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX, - 1 BOT, H, SGA, SGNGAM, SQEPS, T, D1MACH, DGAMR, D9GMIT, D9LGIT, - 2 DLNGAM, D9LGIC - LOGICAL FIRST - SAVE ALNEPS, SQEPS, BOT, FIRST - DATA FIRST /.TRUE./ -C***FIRST EXECUTABLE STATEMENT DGAMIT - IF (FIRST) THEN - ALNEPS = -LOG (D1MACH(3)) - SQEPS = SQRT(D1MACH(4)) - BOT = LOG (D1MACH(1)) - ENDIF - FIRST = .FALSE. -C - IF (X .LT. 0.D0) CALL XERMSG ('SLATEC', 'DGAMIT', 'X IS NEGATIVE' - + , 2, 2) -C - IF (X.NE.0.D0) ALX = LOG (X) - SGA = 1.0D0 - IF (A.NE.0.D0) SGA = SIGN (1.0D0, A) - AINTA = AINT (A + 0.5D0*SGA) - AEPS = A - AINTA -C - IF (X.GT.0.D0) GO TO 20 - DGAMIT = 0.0D0 - IF (AINTA.GT.0.D0 .OR. AEPS.NE.0.D0) DGAMIT = DGAMR(A+1.0D0) - RETURN -C - 20 IF (X.GT.1.D0) GO TO 30 - IF (A.GE.(-0.5D0) .OR. AEPS.NE.0.D0) CALL DLGAMS (A+1.0D0, ALGAP1, - 1 SGNGAM) - DGAMIT = D9GMIT (A, X, ALGAP1, SGNGAM, ALX) - RETURN -C - 30 IF (A.LT.X) GO TO 40 - T = D9LGIT (A, X, DLNGAM(A+1.0D0)) - IF (T.LT.BOT) CALL XERCLR - DGAMIT = EXP (T) - RETURN -C - 40 ALNG = D9LGIC (A, X, ALX) -C -C EVALUATE DGAMIT IN TERMS OF LOG (DGAMIC (A, X)) -C - H = 1.0D0 - IF (AEPS.EQ.0.D0 .AND. AINTA.LE.0.D0) GO TO 50 -C - CALL DLGAMS (A+1.0D0, ALGAP1, SGNGAM) - T = LOG (ABS(A)) + ALNG - ALGAP1 - IF (T.GT.ALNEPS) GO TO 60 -C - IF (T.GT.(-ALNEPS)) H = 1.0D0 - SGA * SGNGAM * EXP(T) - IF (ABS(H).GT.SQEPS) GO TO 50 -C - CALL XERCLR - CALL XERMSG ('SLATEC', 'DGAMIT', 'RESULT LT HALF PRECISION', 1, - + 1) -C - 50 T = -A*ALX + LOG(ABS(H)) - IF (T.LT.BOT) CALL XERCLR - DGAMIT = SIGN (EXP(T), H) - RETURN -C - 60 T = T - A*ALX - IF (T.LT.BOT) CALL XERCLR - DGAMIT = -SGA * SGNGAM * EXP(T) - RETURN -C - END
--- a/liboctave/external/slatec-fn/gami.f Sun Mar 18 18:43:04 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,45 +0,0 @@ -*DECK GAMI - FUNCTION GAMI (A, X) -C***BEGIN PROLOGUE GAMI -C***PURPOSE Evaluate the incomplete Gamma function. -C***LIBRARY SLATEC (FNLIB) -C***CATEGORY C7E -C***TYPE SINGLE PRECISION (GAMI-S, DGAMI-D) -C***KEYWORDS FNLIB, INCOMPLETE GAMMA FUNCTION, SPECIAL FUNCTIONS -C***AUTHOR Fullerton, W., (LANL) -C***DESCRIPTION -C -C Evaluate the incomplete gamma function defined by -C -C GAMI = integral from T = 0 to X of EXP(-T) * T**(A-1.0) . -C -C GAMI is evaluated for positive values of A and non-negative values -C of X. A slight deterioration of 2 or 3 digits accuracy will occur -C when GAMI is very large or very small, because logarithmic variables -C are used. GAMI, A, and X are single precision. -C -C***REFERENCES (NONE) -C***ROUTINES CALLED ALNGAM, GAMIT, XERMSG -C***REVISION HISTORY (YYMMDD) -C 770701 DATE WRITTEN -C 890531 Changed all specific intrinsics to generic. (WRB) -C 890531 REVISION DATE from Version 3.2 -C 891214 Prologue converted to Version 4.0 format. (BAB) -C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) -C***END PROLOGUE GAMI -C***FIRST EXECUTABLE STATEMENT GAMI - IF (A .LE. 0.0) CALL XERMSG ('SLATEC', 'GAMI', - + 'A MUST BE GT ZERO', 1, 2) - IF (X .LT. 0.0) CALL XERMSG ('SLATEC', 'GAMI', - + 'X MUST BE GE ZERO', 2, 2) -C - GAMI = 0.0 - IF (X.EQ.0.0) RETURN -C -C THE ONLY ERROR POSSIBLE IN THE EXPRESSION BELOW IS A FATAL OVERFLOW. - FACTOR = EXP (ALNGAM(A) + A*LOG(X) ) -C - GAMI = FACTOR * GAMIT(A, X) -C - RETURN - END
--- a/liboctave/external/slatec-fn/gamit.f Sun Mar 18 18:43:04 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,112 +0,0 @@ -*DECK GAMIT - REAL FUNCTION GAMIT (A, X) -C***BEGIN PROLOGUE GAMIT -C***PURPOSE Calculate Tricomi's form of the incomplete Gamma function. -C***LIBRARY SLATEC (FNLIB) -C***CATEGORY C7E -C***TYPE SINGLE PRECISION (GAMIT-S, DGAMIT-D) -C***KEYWORDS COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB, -C SPECIAL FUNCTIONS, TRICOMI -C***AUTHOR Fullerton, W., (LANL) -C***DESCRIPTION -C -C Evaluate Tricomi's incomplete gamma function defined by -C -C GAMIT = X**(-A)/GAMMA(A) * integral from 0 to X of EXP(-T) * -C T**(A-1.) -C -C for A .GT. 0.0 and by analytic continuation for A .LE. 0.0. -C GAMMA(X) is the complete gamma function of X. -C -C GAMIT is evaluated for arbitrary real values of A and for non- -C negative values of X (even though GAMIT is defined for X .LT. -C 0.0), except that for X = 0 and A .LE. 0.0, GAMIT is infinite, -C which is a fatal error. -C -C The function and both arguments are REAL. -C -C A slight deterioration of 2 or 3 digits accuracy will occur when -C GAMIT is very large or very small in absolute value, because log- -C arithmic variables are used. Also, if the parameter A is very -C close to a negative integer (but not a negative integer), there is -C a loss of accuracy, which is reported if the result is less than -C half machine precision. -C -C***REFERENCES W. Gautschi, A computational procedure for incomplete -C gamma functions, ACM Transactions on Mathematical -C Software 5, 4 (December 1979), pp. 466-481. -C W. Gautschi, Incomplete gamma functions, Algorithm 542, -C ACM Transactions on Mathematical Software 5, 4 -C (December 1979), pp. 482-489. -C***ROUTINES CALLED ALGAMS, ALNGAM, GAMR, R1MACH, R9GMIT, R9LGIC, -C R9LGIT, XERCLR, XERMSG -C***REVISION HISTORY (YYMMDD) -C 770701 DATE WRITTEN -C 890531 Changed all specific intrinsics to generic. (WRB) -C 890531 REVISION DATE from Version 3.2 -C 891214 Prologue converted to Version 4.0 format. (BAB) -C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) -C 920528 DESCRIPTION and REFERENCES sections revised. (WRB) -C***END PROLOGUE GAMIT - LOGICAL FIRST - SAVE ALNEPS, SQEPS, BOT, FIRST - DATA FIRST /.TRUE./ -C***FIRST EXECUTABLE STATEMENT GAMIT - IF (FIRST) THEN - ALNEPS = -LOG(R1MACH(3)) - SQEPS = SQRT(R1MACH(4)) - BOT = LOG(R1MACH(1)) - ENDIF - FIRST = .FALSE. -C - IF (X .LT. 0.0) CALL XERMSG ('SLATEC', 'GAMIT', 'X IS NEGATIVE', - + 2, 2) -C - IF (X.NE.0.0) ALX = LOG(X) - SGA = 1.0 - IF (A.NE.0.0) SGA = SIGN (1.0, A) - AINTA = AINT (A+0.5*SGA) - AEPS = A - AINTA -C - IF (X.GT.0.0) GO TO 20 - GAMIT = 0.0 - IF (AINTA.GT.0.0 .OR. AEPS.NE.0.0) GAMIT = GAMR(A+1.0) - RETURN -C - 20 IF (X.GT.1.0) GO TO 40 - IF (A.GE.(-0.5) .OR. AEPS.NE.0.0) CALL ALGAMS (A+1.0, ALGAP1, - 1 SGNGAM) - GAMIT = R9GMIT (A, X, ALGAP1, SGNGAM, ALX) - RETURN -C - 40 IF (A.LT.X) GO TO 50 - T = R9LGIT (A, X, ALNGAM(A+1.0)) - IF (T.LT.BOT) CALL XERCLR - GAMIT = EXP(T) - RETURN -C - 50 ALNG = R9LGIC (A, X, ALX) -C -C EVALUATE GAMIT IN TERMS OF LOG(GAMIC(A,X)) -C - H = 1.0 - IF (AEPS.EQ.0.0 .AND. AINTA.LE.0.0) GO TO 60 - CALL ALGAMS (A+1.0, ALGAP1, SGNGAM) - T = LOG(ABS(A)) + ALNG - ALGAP1 - IF (T.GT.ALNEPS) GO TO 70 - IF (T.GT.(-ALNEPS)) H = 1.0 - SGA*SGNGAM*EXP(T) - IF (ABS(H).GT.SQEPS) GO TO 60 - CALL XERCLR - CALL XERMSG ('SLATEC', 'GAMIT', 'RESULT LT HALF PRECISION', 1, 1) -C - 60 T = -A*ALX + LOG(ABS(H)) - IF (T.LT.BOT) CALL XERCLR - GAMIT = SIGN (EXP(T), H) - RETURN -C - 70 T = T - A*ALX - IF (T.LT.BOT) CALL XERCLR - GAMIT = -SGA*SGNGAM*EXP(T) - RETURN -C - END
--- a/liboctave/external/slatec-fn/module.mk Sun Mar 18 18:43:04 2018 -0700 +++ b/liboctave/external/slatec-fn/module.mk Mon Mar 19 13:50:10 2018 -0700 @@ -3,16 +3,12 @@ %reldir%/alngam.f \ %reldir%/alnrel.f \ %reldir%/algams.f \ - %reldir%/betai.f \ %reldir%/csevl.f \ %reldir%/d9gmit.f \ %reldir%/d9lgic.f \ %reldir%/d9lgit.f \ %reldir%/d9lgmc.f \ - %reldir%/dbetai.f \ %reldir%/dcsevl.f \ - %reldir%/dgami.f \ - %reldir%/dgamit.f \ %reldir%/dgamlm.f \ %reldir%/dgamma.f \ %reldir%/dgamr.f \ @@ -23,8 +19,6 @@ %reldir%/dpchim.f \ %reldir%/dpchst.f \ %reldir%/dpsifn.f \ - %reldir%/gami.f \ - %reldir%/gamit.f \ %reldir%/gamlim.f \ %reldir%/gamma.f \ %reldir%/gamr.f \ @@ -36,12 +30,6 @@ %reldir%/r9lgmc.f \ %reldir%/r9lgit.f \ %reldir%/r9gmit.f \ - %reldir%/r9lgic.f \ - %reldir%/xdbetai.f \ - %reldir%/xdgami.f \ - %reldir%/xdgamit.f \ - %reldir%/xgmainc.f \ - %reldir%/xsgmainc.f \ - %reldir%/xbetai.f + %reldir%/r9lgic.f DIRSTAMP_FILES += %reldir%/$(octave_dirstamp)
--- a/liboctave/external/slatec-fn/xbetai.f Sun Mar 18 18:43:04 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ - subroutine xbetai (x, a, b, result) - external betai - real x, a, b, result, betai - result = betai (x, a, b) - return - end
--- a/liboctave/external/slatec-fn/xdbetai.f Sun Mar 18 18:43:04 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ - subroutine xdbetai (x, a, b, result) - external dbetai - double precision x, a, b, result, dbetai - result = dbetai (x, a, b) - return - end
--- a/liboctave/external/slatec-fn/xdgami.f Sun Mar 18 18:43:04 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ - subroutine xdgami (a, x, result) - external dgami - double precision a, x, result, dgami - result = dgami (a, x) - return - end
--- a/liboctave/external/slatec-fn/xdgamit.f Sun Mar 18 18:43:04 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ - subroutine xdgamit (a, x, result) - external dgamit - double precision a, x, result, dgamit - result = dgamit (a, x) - return - end
--- a/liboctave/external/slatec-fn/xgmainc.f Sun Mar 18 18:43:04 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,100 +0,0 @@ - subroutine xgammainc (a, x, result) - -c -- jwe, based on DGAMIT. -c -c -- Do a better job than dgami for large values of x. - - double precision a, x, result - intrinsic exp, log, sqrt, sign, aint - external dgami, dlngam, d9lgit, d9lgic, d9gmit - -C external dgamr -C DOUBLE PRECISION DGAMR - - DOUBLE PRECISION AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX, - $ BOT, H, SGA, SGNGAM, SQEPS, T, D1MACH, D9GMIT, - $ D9LGIC, D9LGIT, DLNGAM, DGAMI - - LOGICAL FIRST - - SAVE ALNEPS, SQEPS, BOT, FIRST - - DATA FIRST /.TRUE./ - - if (x .eq. 0.0d0) then - - if (a .eq. 0.0d0) then - result = 1.0d0 - else - result = 0.0d0 - endif - - else - - IF (FIRST) THEN - ALNEPS = -LOG (D1MACH(3)) - SQEPS = SQRT(D1MACH(4)) - BOT = LOG (D1MACH(1)) - ENDIF - FIRST = .FALSE. -C - IF (X .LT. 0.D0) CALL XERMSG ('SLATEC', 'XGMAINC', 'X IS NEGATIVE' - + , 2, 2) -C - IF (X.NE.0.D0) ALX = LOG (X) - SGA = 1.0D0 - IF (A.NE.0.D0) SGA = SIGN (1.0D0, A) - AINTA = AINT (A + 0.5D0*SGA) - AEPS = A - AINTA -C -C IF (X.GT.0.D0) GO TO 20 -C DGAMIT = 0.0D0 -C IF (AINTA.GT.0.D0 .OR. AEPS.NE.0.D0) DGAMIT = DGAMR(A+1.0D0) -C RETURN -C - 20 IF (X.GT.1.D0) GO TO 30 - IF (A.GE.(-0.5D0) .OR. AEPS.NE.0.D0) CALL DLGAMS (A+1.0D0, ALGAP1, - 1 SGNGAM) -C DGAMIT = D9GMIT (A, X, ALGAP1, SGNGAM, ALX) - result = exp (a*alx + log (D9GMIT (A, X, ALGAP1, SGNGAM, ALX))) - RETURN -C - 30 IF (A.LT.X) GO TO 40 - T = D9LGIT (A, X, DLNGAM(A+1.0D0)) - IF (T.LT.BOT) CALL XERCLR -C DGAMIT = EXP (T) - result = EXP (a*alx + T) - RETURN -C - 40 ALNG = D9LGIC (A, X, ALX) -C -C EVALUATE DGAMIT IN TERMS OF LOG (DGAMIC (A, X)) -C - H = 1.0D0 - IF (AEPS.EQ.0.D0 .AND. AINTA.LE.0.D0) GO TO 50 -C - CALL DLGAMS (A+1.0D0, ALGAP1, SGNGAM) - T = LOG (ABS(A)) + ALNG - ALGAP1 - IF (T.GT.ALNEPS) GO TO 60 -C - IF (T.GT.(-ALNEPS)) H = 1.0D0 - SGA * SGNGAM * EXP(T) - IF (ABS(H).GT.SQEPS) GO TO 50 -C - CALL XERCLR - CALL XERMSG ('SLATEC', 'XGMAINC', 'RESULT LT HALF PRECISION', 1, - + 1) -C -C 50 T = -A*ALX + LOG(ABS(H)) -C IF (T.LT.BOT) CALL XERCLR -C DGAMIT = SIGN (EXP(T), H) - 50 result = H - RETURN -C -C 60 T = T - A*ALX - 60 IF (T.LT.BOT) CALL XERCLR - result = -SGA * SGNGAM * EXP(T) - RETURN - - endif - return - end
--- a/liboctave/external/slatec-fn/xsgmainc.f Sun Mar 18 18:43:04 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,100 +0,0 @@ - subroutine xsgammainc (a, x, result) - -c -- jwe, based on GAMIT. -c -c -- Do a better job than gami for large values of x. - - real a, x, result - intrinsic exp, log, sqrt, sign, aint - external gami, alngam, r9lgit, r9lgic, r9gmit - -C external gamr -C real GAMR - - REAL AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX, - $ BOT, H, SGA, SGNGAM, SQEPS, T, R1MACH, R9GMIT, - $ R9LGIC, R9LGIT, ALNGAM, GAMI - - LOGICAL FIRST - - SAVE ALNEPS, SQEPS, BOT, FIRST - - DATA FIRST /.TRUE./ - - if (x .eq. 0.0e0) then - - if (a .eq. 0.0e0) then - result = 1.0e0 - else - result = 0.0e0 - endif - - else - - IF (FIRST) THEN - ALNEPS = -LOG (R1MACH(3)) - SQEPS = SQRT(R1MACH(4)) - BOT = LOG (R1MACH(1)) - ENDIF - FIRST = .FALSE. -C - IF (X .LT. 0.E0) CALL XERMSG ('SLATEC', 'XGMAINC', 'X IS NEGATIVE' - + , 2, 2) -C - IF (X.NE.0.E0) ALX = LOG (X) - SGA = 1.0E0 - IF (A.NE.0.E0) SGA = SIGN (1.0E0, A) - AINTA = AINT (A + 0.5E0*SGA) - AEPS = A - AINTA -C -C IF (X.GT.0.E0) GO TO 20 -C GAMIT = 0.0E0 -C IF (AINTA.GT.0.E0 .OR. AEPS.NE.0.E0) GAMIT = GAMR(A+1.0E0) -C RETURN -C - 20 IF (X.GT.1.E0) GO TO 30 - IF (A.GE.(-0.5E0) .OR. AEPS.NE.0.E0) CALL ALGAMS (A+1.0E0, ALGAP1, - 1 SGNGAM) -C GAMIT = R9GMIT (A, X, ALGAP1, SGNGAM, ALX) - result = exp (a*alx + log (R9GMIT (A, X, ALGAP1, SGNGAM, ALX))) - RETURN -C - 30 IF (A.LT.X) GO TO 40 - T = R9LGIT (A, X, ALNGAM(A+1.0E0)) - IF (T.LT.BOT) CALL XERCLR -C GAMIT = EXP (T) - result = EXP (a*alx + T) - RETURN -C - 40 ALNG = R9LGIC (A, X, ALX) -C -C EVALUATE GAMIT IN TERMS OF LOG (DGAMIC (A, X)) -C - H = 1.0E0 - IF (AEPS.EQ.0.E0 .AND. AINTA.LE.0.E0) GO TO 50 -C - CALL ALGAMS (A+1.0E0, ALGAP1, SGNGAM) - T = LOG (ABS(A)) + ALNG - ALGAP1 - IF (T.GT.ALNEPS) GO TO 60 -C - IF (T.GT.(-ALNEPS)) H = 1.0E0 - SGA * SGNGAM * EXP(T) - IF (ABS(H).GT.SQEPS) GO TO 50 -C - CALL XERCLR - CALL XERMSG ('SLATEC', 'XGMAINC', 'RESULT LT HALF PRECISION', 1, - + 1) -C -C 50 T = -A*ALX + LOG(ABS(H)) -C IF (T.LT.BOT) CALL XERCLR -C GAMIT = SIGN (EXP(T), H) - 50 result = H - RETURN -C -C 60 T = T - A*ALX - 60 IF (T.LT.BOT) CALL XERCLR - result = -SGA * SGNGAM * EXP(T) - RETURN - - endif - return - end
--- a/liboctave/numeric/lo-specfun.cc Sun Mar 18 18:43:04 2018 -0700 +++ b/liboctave/numeric/lo-specfun.cc Mon Mar 19 13:50:10 2018 -0700 @@ -77,6 +77,7 @@ { case 0: case 3: + case 4: retval = val; break; @@ -109,6 +110,7 @@ { case 0: case 3: + case 4: retval = val; break; @@ -304,17 +306,10 @@ double zr = z.real (); double zi = z.imag (); - F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, 2, 1, &yr, &yi, nz, t_ierr); + F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, kode, 1, &yr, &yi, nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - double expz = exp (std::abs (zi)); - yr *= expz; - yi *= expz; - } - if (zi == 0.0 && zr >= 0.0) yi = 0.0; @@ -375,18 +370,11 @@ } else { - F77_FUNC (zbesy, ZBESY) (zr, zi, alpha, 2, 1, &yr, &yi, nz, + F77_FUNC (zbesy, ZBESY) (zr, zi, alpha, kode, 1, &yr, &yi, nz, &wr, &wi, t_ierr); ierr = t_ierr; - if (kode != 2) - { - double expz = exp (std::abs (zi)); - yr *= expz; - yi *= expz; - } - if (zi == 0.0 && zr >= 0.0) yi = 0.0; } @@ -437,17 +425,10 @@ double zr = z.real (); double zi = z.imag (); - F77_FUNC (zbesi, ZBESI) (zr, zi, alpha, 2, 1, &yr, &yi, nz, t_ierr); + F77_FUNC (zbesi, ZBESI) (zr, zi, alpha, kode, 1, &yr, &yi, nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - double expz = exp (std::abs (zr)); - yr *= expz; - yi *= expz; - } - if (zi == 0.0 && zr >= 0.0) yi = 0.0; @@ -513,24 +494,11 @@ } else { - F77_FUNC (zbesk, ZBESK) (zr, zi, alpha, 2, 1, &yr, &yi, nz, + F77_FUNC (zbesk, ZBESK) (zr, zi, alpha, kode, 1, &yr, &yi, nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - Complex expz = exp (-z); - - double rexpz = expz.real (); - double iexpz = expz.imag (); - - double tmp = yr*rexpz - yi*iexpz; - - yi = yr*iexpz + yi*rexpz; - yr = tmp; - } - if (zi == 0.0 && zr >= 0.0) yi = 0.0; } @@ -562,24 +530,11 @@ double zr = z.real (); double zi = z.imag (); - F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, 2, 1, 1, &yr, &yi, nz, + F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, kode, 1, 1, &yr, &yi, nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - Complex expz = exp (Complex (0.0, 1.0) * z); - - double rexpz = expz.real (); - double iexpz = expz.imag (); - - double tmp = yr*rexpz - yi*iexpz; - - yi = yr*iexpz + yi*rexpz; - yr = tmp; - } - retval = bessel_return_value (Complex (yr, yi), ierr); } else @@ -611,24 +566,11 @@ double zr = z.real (); double zi = z.imag (); - F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, 2, 2, 1, &yr, &yi, nz, + F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, kode, 2, 1, &yr, &yi, nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - Complex expz = exp (-Complex (0.0, 1.0) * z); - - double rexpz = expz.real (); - double iexpz = expz.imag (); - - double tmp = yr*rexpz - yi*iexpz; - - yi = yr*iexpz + yi*rexpz; - yr = tmp; - } - retval = bessel_return_value (Complex (yr, yi), ierr); } else @@ -922,17 +864,11 @@ F77_INT nz, t_ierr; - F77_FUNC (cbesj, CBESJ) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 1, + F77_FUNC (cbesj, CBESJ) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, F77_CMPLX_ARG (&y), nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - float expz = exp (std::abs (z.imag ())); - y *= expz; - } - if (z.imag () == 0.0 && z.real () >= 0.0) y = FloatComplex (y.real (), 0.0); @@ -990,18 +926,12 @@ } else { - F77_FUNC (cbesy, CBESY) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 1, + 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 (kode != 2) - { - float expz = exp (std::abs (z.imag ())); - y *= expz; - } - if (z.imag () == 0.0 && z.real () >= 0.0) y = FloatComplex (y.real (), 0.0); } @@ -1050,17 +980,11 @@ F77_INT nz, t_ierr; - F77_FUNC (cbesi, CBESI) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 1, + F77_FUNC (cbesi, CBESI) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, F77_CMPLX_ARG (&y), nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - float expz = exp (std::abs (z.real ())); - y *= expz; - } - if (z.imag () == 0.0 && z.real () >= 0.0) y = FloatComplex (y.real (), 0.0); @@ -1115,24 +1039,11 @@ } else { - F77_FUNC (cbesk, CBESK) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 1, + F77_FUNC (cbesk, CBESK) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, F77_CMPLX_ARG (&y), nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - FloatComplex expz = exp (-z); - - float rexpz = expz.real (); - float iexpz = expz.imag (); - - float tmp_r = y.real () * rexpz - y.imag () * iexpz; - float tmp_i = y.real () * iexpz + y.imag () * rexpz; - - y = FloatComplex (tmp_r, tmp_i); - } - if (z.imag () == 0.0 && z.real () >= 0.0) y = FloatComplex (y.real (), 0.0); } @@ -1160,24 +1071,11 @@ F77_INT nz, t_ierr; - F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 1, 1, + F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, 1, F77_CMPLX_ARG (&y), nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - FloatComplex expz = exp (FloatComplex (0.0, 1.0) * z); - - float rexpz = expz.real (); - float iexpz = expz.imag (); - - float tmp_r = y.real () * rexpz - y.imag () * iexpz; - float tmp_i = y.real () * iexpz + y.imag () * rexpz; - - y = FloatComplex (tmp_r, tmp_i); - } - retval = bessel_return_value (y, ierr); } else @@ -1206,24 +1104,11 @@ F77_INT nz, t_ierr; - F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 2, 1, + F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 2, 1, F77_CMPLX_ARG (&y), nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - FloatComplex expz = exp (-FloatComplex (0.0, 1.0) * z); - - float rexpz = expz.real (); - float iexpz = expz.imag (); - - float tmp_r = y.real () * rexpz - y.imag () * iexpz; - float tmp_i = y.real () * iexpz + y.imag () * rexpz; - - y = FloatComplex (tmp_r, tmp_i); - } - retval = bessel_return_value (y, ierr); } else @@ -1489,765 +1374,6 @@ #undef NN_BESSEL #undef RC_BESSEL - OCTAVE_NORETURN static void - err_betainc_nonconformant (const dim_vector& d1, const dim_vector& d2, - const dim_vector& d3) - { - std::string d1_str = d1.str (); - std::string d2_str = d2.str (); - std::string d3_str = d3.str (); - - (*current_liboctave_error_handler) - ("betainc: nonconformant arguments (x is %s, a is %s, b is %s)", - d1_str.c_str (), d2_str.c_str (), d3_str.c_str ()); - } - - OCTAVE_NORETURN static void - err_betaincinv_nonconformant (const dim_vector& d1, const dim_vector& d2, - const dim_vector& d3) - { - std::string d1_str = d1.str (); - std::string d2_str = d2.str (); - std::string d3_str = d3.str (); - - (*current_liboctave_error_handler) - ("betaincinv: nonconformant arguments (x is %s, a is %s, b is %s)", - d1_str.c_str (), d2_str.c_str (), d3_str.c_str ()); - } - - double - betainc (double x, double a, double b) - { - double retval; - F77_XFCN (xdbetai, XDBETAI, (x, a, b, retval)); - return retval; - } - - Array<double> - betainc (double x, double a, const Array<double>& b) - { - dim_vector dv = b.dims (); - octave_idx_type nel = dv.numel (); - - Array<double> retval (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x, a, b(i)); - - return retval; - } - - Array<double> - betainc (double x, const Array<double>& a, double b) - { - dim_vector dv = a.dims (); - octave_idx_type nel = dv.numel (); - - Array<double> retval (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x, a(i), b); - - return retval; - } - - Array<double> - betainc (double x, const Array<double>& a, const Array<double>& b) - { - Array<double> retval; - dim_vector dv = a.dims (); - - if (dv != b.dims ()) - err_betainc_nonconformant (dim_vector (0, 0), dv, b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x, a(i), b(i)); - - return retval; - } - - Array<double> - betainc (const Array<double>& x, double a, double b) - { - dim_vector dv = x.dims (); - octave_idx_type nel = dv.numel (); - - Array<double> retval (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a, b); - - return retval; - } - - Array<double> - betainc (const Array<double>& x, double a, const Array<double>& b) - { - Array<double> retval; - dim_vector dv = x.dims (); - - if (dv != b.dims ()) - err_betainc_nonconformant (dv, dim_vector (0, 0), b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a, b(i)); - - return retval; - } - - Array<double> - betainc (const Array<double>& x, const Array<double>& a, double b) - { - Array<double> retval; - dim_vector dv = x.dims (); - - if (dv != a.dims ()) - err_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0)); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a(i), b); - - return retval; - } - - Array<double> - betainc (const Array<double>& x, const Array<double>& a, - const Array<double>& b) - { - Array<double> retval; - dim_vector dv = x.dims (); - - if (dv != a.dims () || dv != b.dims ()) - err_betainc_nonconformant (dv, a.dims (), b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a(i), b(i)); - - return retval; - } - - float - betainc (float x, float a, float b) - { - float retval; - F77_XFCN (xbetai, XBETAI, (x, a, b, retval)); - return retval; - } - - Array<float> - betainc (float x, float a, const Array<float>& b) - { - dim_vector dv = b.dims (); - octave_idx_type nel = dv.numel (); - - Array<float> retval (dv); - - float *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x, a, b(i)); - - return retval; - } - - Array<float> - betainc (float x, const Array<float>& a, float b) - { - dim_vector dv = a.dims (); - octave_idx_type nel = dv.numel (); - - Array<float> retval (dv); - - float *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x, a(i), b); - - return retval; - } - - Array<float> - betainc (float x, const Array<float>& a, const Array<float>& b) - { - Array<float> retval; - dim_vector dv = a.dims (); - - if (dv != b.dims ()) - err_betainc_nonconformant (dim_vector (0, 0), dv, b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - float *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x, a(i), b(i)); - - return retval; - } - - Array<float> - betainc (const Array<float>& x, float a, float b) - { - dim_vector dv = x.dims (); - octave_idx_type nel = dv.numel (); - - Array<float> retval (dv); - - float *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a, b); - - return retval; - } - - Array<float> - betainc (const Array<float>& x, float a, const Array<float>& b) - { - Array<float> retval; - dim_vector dv = x.dims (); - - if (dv != b.dims ()) - err_betainc_nonconformant (dv, dim_vector (0, 0), b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - float *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a, b(i)); - - return retval; - } - - Array<float> - betainc (const Array<float>& x, const Array<float>& a, float b) - { - Array<float> retval; - dim_vector dv = x.dims (); - - if (dv != a.dims ()) - err_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0)); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - float *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a(i), b); - - return retval; - } - - Array<float> - betainc (const Array<float>& x, const Array<float>& a, - const Array<float>& b) - { - Array<float> retval; - dim_vector dv = x.dims (); - - if (dv != a.dims () || dv != b.dims ()) - err_betainc_nonconformant (dv, a.dims (), b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - float *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betainc (x(i), a(i), b(i)); - - return retval; - } - - // - // Incomplete Beta function ratio - // - // Algorithm based on the one by John Burkardt. - // See http://people.sc.fsu.edu/~jburkardt/cpp_src/asa109/asa109.html - // - // The original code is distributed under the GNU LGPL v3 license. - // - // Reference: - // - // KL Majumder, GP Bhattacharjee, - // Algorithm AS 63: - // The incomplete Beta Integral, - // Applied Statistics, - // Volume 22, Number 3, 1973, pages 409-411. - // - double - betain (double x, double p, double q, double beta, bool& err) - { - double acu = 0.1E-14, ai, cx; - bool indx; - int ns; - double pp, psq, qq, rx, temp, term, value, xx; - - value = x; - err = false; - - // Check the input arguments. - - if ((p <= 0.0 || q <= 0.0) || (x < 0.0 || 1.0 < x)) - { - err = true; - return value; - } - - // Special cases. - - if (x == 0.0 || x == 1.0) - { - return value; - } - - // Change tail if necessary and determine S. - - psq = p + q; - cx = 1.0 - x; - - if (p < psq * x) - { - xx = cx; - cx = x; - pp = q; - qq = p; - indx = true; - } - else - { - xx = x; - pp = p; - qq = q; - indx = false; - } - - term = 1.0; - ai = 1.0; - value = 1.0; - ns = static_cast<int> (qq + cx * psq); - - // Use the Soper reduction formula. - - rx = xx / cx; - temp = qq - ai; - if (ns == 0) - { - rx = xx; - } - - for ( ; ; ) - { - term *= temp * rx / (pp + ai); - value += term; - temp = fabs (term); - - if (temp <= acu && temp <= acu * value) - { - value *= exp (pp * std::log (xx) - + (qq - 1.0) * std::log (cx) - beta) / pp; - - if (indx) - { - value = 1.0 - value; - } - break; - } - - ai += 1.0; - ns -= 1; - - if (0 <= ns) - { - temp = qq - ai; - if (ns == 0) - { - rx = xx; - } - } - else - { - temp = psq; - psq += 1.0; - } - } - - return value; - } - - // - // Inverse of the incomplete Beta function - // - // Algorithm based on the one by John Burkardt. - // See http://people.sc.fsu.edu/~jburkardt/cpp_src/asa109/asa109.html - // - // The original code is distributed under the GNU LGPL v3 license. - // - // Reference: - // - // GW Cran, KJ Martin, GE Thomas, - // Remark AS R19 and Algorithm AS 109: - // A Remark on Algorithms AS 63: The Incomplete Beta Integral - // and AS 64: Inverse of the Incomplete Beta Integeral, - // Applied Statistics, - // Volume 26, Number 1, 1977, pages 111-114. - // - double - betaincinv (double y, double p, double q) - { - double a, acu, adj, fpu, g, h; - int iex; - bool indx; - double pp, prev, qq, r, s, sae = -37.0, sq, t, tx, value, w, xin, ycur, yprev; - - double beta = lgamma (p) + lgamma (q) - lgamma (p + q); - bool err = false; - fpu = std::pow (10.0, sae); - value = y; - - // Test for admissibility of parameters. - - if (p <= 0.0 || q <= 0.0) - (*current_liboctave_error_handler) ("betaincinv: wrong parameters"); - if (y < 0.0 || 1.0 < y) - (*current_liboctave_error_handler) ("betaincinv: wrong parameter Y"); - - if (y == 0.0 || y == 1.0) - return value; - - // Change tail if necessary. - - if (0.5 < y) - { - a = 1.0 - y; - pp = q; - qq = p; - indx = true; - } - else - { - a = y; - pp = p; - qq = q; - indx = false; - } - - // Calculate the initial approximation. - - r = std::sqrt (- std::log (a * a)); - - ycur = r - (2.30753 + 0.27061 * r) / (1.0 + (0.99229 + 0.04481 * r) * r); - - if (1.0 < pp && 1.0 < qq) - { - r = (ycur * ycur - 3.0) / 6.0; - s = 1.0 / (pp + pp - 1.0); - t = 1.0 / (qq + qq - 1.0); - h = 2.0 / (s + t); - w = ycur * std::sqrt (h + r) / h - (t - s) * (r + 5.0 / 6.0 - 2.0 / (3.0 * h)); - value = pp / (pp + qq * exp (w + w)); - } - else - { - r = qq + qq; - t = 1.0 / (9.0 * qq); - t = r * std::pow (1.0 - t + ycur * std::sqrt (t), 3); - - if (t <= 0.0) - { - value = 1.0 - exp ((std::log ((1.0 - a) * qq) + beta) / qq); - } - else - { - t = (4.0 * pp + r - 2.0) / t; - - if (t <= 1.0) - { - value = exp ((std::log (a * pp) + beta) / pp); - } - else - { - value = 1.0 - 2.0 / (t + 1.0); - } - } - } - - // Solve for X by a modified Newton-Raphson method, - // using the function BETAIN. - - r = 1.0 - pp; - t = 1.0 - qq; - yprev = 0.0; - sq = 1.0; - prev = 1.0; - - if (value < 0.0001) - { - value = 0.0001; - } - - if (0.9999 < value) - { - value = 0.9999; - } - - iex = std::max (- 5.0 / pp / pp - 1.0 / std::pow (a, 0.2) - 13.0, sae); - - acu = std::pow (10.0, iex); - - for ( ; ; ) - { - ycur = betain (value, pp, qq, beta, err); - - if (err) - { - return value; - } - - xin = value; - ycur = (ycur - a) * exp (beta + r * std::log (xin) - + t * std::log (1.0 - xin)); - - if (ycur * yprev <= 0.0) - { - prev = std::max (sq, fpu); - } - - g = 1.0; - - for ( ; ; ) - { - for ( ; ; ) - { - adj = g * ycur; - sq = adj * adj; - - if (sq < prev) - { - tx = value - adj; - - if (0.0 <= tx && tx <= 1.0) - { - break; - } - } - g /= 3.0; - } - - if (prev <= acu) - { - if (indx) - { - value = 1.0 - value; - } - return value; - } - - if (ycur * ycur <= acu) - { - if (indx) - { - value = 1.0 - value; - } - return value; - } - - if (tx != 0.0 && tx != 1.0) - { - break; - } - - g /= 3.0; - } - - if (tx == value) - { - break; - } - - value = tx; - yprev = ycur; - } - - if (indx) - { - value = 1.0 - value; - } - - return value; - } - - Array<double> - betaincinv (double x, double a, const Array<double>& b) - { - dim_vector dv = b.dims (); - octave_idx_type nel = dv.numel (); - - Array<double> retval (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betaincinv (x, a, b(i)); - - return retval; - } - - Array<double> - betaincinv (double x, const Array<double>& a, double b) - { - dim_vector dv = a.dims (); - octave_idx_type nel = dv.numel (); - - Array<double> retval (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betaincinv (x, a(i), b); - - return retval; - } - - Array<double> - betaincinv (double x, const Array<double>& a, const Array<double>& b) - { - Array<double> retval; - dim_vector dv = a.dims (); - - if (dv != b.dims ()) - err_betaincinv_nonconformant (dim_vector (0, 0), dv, b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betaincinv (x, a(i), b(i)); - - return retval; - } - - Array<double> - betaincinv (const Array<double>& x, double a, double b) - { - dim_vector dv = x.dims (); - octave_idx_type nel = dv.numel (); - - Array<double> retval (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betaincinv (x(i), a, b); - - return retval; - } - - Array<double> - betaincinv (const Array<double>& x, double a, const Array<double>& b) - { - Array<double> retval; - dim_vector dv = x.dims (); - - if (dv != b.dims ()) - err_betaincinv_nonconformant (dv, dim_vector (0, 0), b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betaincinv (x(i), a, b(i)); - - return retval; - } - - Array<double> - betaincinv (const Array<double>& x, const Array<double>& a, double b) - { - Array<double> retval; - dim_vector dv = x.dims (); - - if (dv != a.dims ()) - err_betaincinv_nonconformant (dv, a.dims (), dim_vector (0, 0)); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betaincinv (x(i), a(i), b); - - return retval; - } - - Array<double> - betaincinv (const Array<double>& x, const Array<double>& a, - const Array<double>& b) - { - Array<double> retval; - dim_vector dv = x.dims (); - - if (dv != a.dims () && dv != b.dims ()) - err_betaincinv_nonconformant (dv, a.dims (), b.dims ()); - - octave_idx_type nel = dv.numel (); - - retval.resize (dv); - - double *pretval = retval.fortran_vec (); - - for (octave_idx_type i = 0; i < nel; i++) - *pretval++ = betaincinv (x(i), a(i), b(i)); - - return retval; - } - Complex biry (const Complex& z, bool deriv, bool scaled, octave_idx_type& ierr) { @@ -2829,334 +1955,6 @@ return result; } - // FIXME: there is still room for improvement here... - - double - gammainc (double x, double a, bool& err) - { - if (a < 0.0 || x < 0.0) - (*current_liboctave_error_handler) - ("gammainc: A and X must be non-negative"); - - err = false; - - double retval; - - F77_XFCN (xgammainc, XGAMMAINC, (a, x, retval)); - - return retval; - } - - Matrix - gammainc (double x, const Matrix& a) - { - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.cols (); - - Matrix retval (nr, nc); - - bool err; - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - { - retval(i,j) = gammainc (x, a(i,j), err); - - if (err) - return Matrix (); - } - - return retval; - } - - Matrix - gammainc (const Matrix& x, double a) - { - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - Matrix retval (nr, nc); - - bool err; - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - { - retval(i,j) = gammainc (x(i,j), a, err); - - if (err) - return Matrix (); - } - - return retval; - } - - Matrix - gammainc (const Matrix& x, const Matrix& a) - { - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - - if (nr != a_nr || nc != a_nc) - (*current_liboctave_error_handler) - ("gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)", - nr, nc, a_nr, a_nc); - - Matrix retval (nr, nc); - - bool err; - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - { - retval(i,j) = gammainc (x(i,j), a(i,j), err); - - if (err) - return Matrix (); - } - - return retval; - } - - NDArray - gammainc (double x, const NDArray& a) - { - dim_vector dv = a.dims (); - octave_idx_type nel = dv.numel (); - - NDArray retval (dv); - - bool err; - - for (octave_idx_type i = 0; i < nel; i++) - { - retval(i) = gammainc (x, a(i), err); - - if (err) - return NDArray (); - } - - return retval; - } - - NDArray - gammainc (const NDArray& x, double a) - { - dim_vector dv = x.dims (); - octave_idx_type nel = dv.numel (); - - NDArray retval (dv); - - bool err; - - for (octave_idx_type i = 0; i < nel; i++) - { - retval(i) = gammainc (x(i), a, err); - - if (err) - return NDArray (); - } - - return retval; - } - - NDArray - gammainc (const NDArray& x, const NDArray& a) - { - dim_vector dv = x.dims (); - octave_idx_type nel = dv.numel (); - - if (dv != a.dims ()) - { - std::string x_str = dv.str (); - std::string a_str = a.dims ().str (); - - (*current_liboctave_error_handler) - ("gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)", - x_str.c_str (), a_str.c_str ()); - } - - NDArray retval (dv); - - bool err; - - for (octave_idx_type i = 0; i < nel; i++) - { - retval(i) = gammainc (x(i), a(i), err); - - if (err) - return NDArray (); - } - - return retval; - } - - float - gammainc (float x, float a, bool& err) - { - if (a < 0.0 || x < 0.0) - (*current_liboctave_error_handler) - ("gammainc: A and X must be non-negative"); - - err = false; - - float retval; - - F77_XFCN (xsgammainc, XSGAMMAINC, (a, x, retval)); - - return retval; - } - - FloatMatrix - gammainc (float x, const FloatMatrix& a) - { - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.cols (); - - FloatMatrix retval (nr, nc); - - bool err; - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - { - retval(i,j) = gammainc (x, a(i,j), err); - - if (err) - return FloatMatrix (); - } - - return retval; - } - - FloatMatrix - gammainc (const FloatMatrix& x, float a) - { - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - FloatMatrix retval (nr, nc); - - bool err; - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - { - retval(i,j) = gammainc (x(i,j), a, err); - - if (err) - return FloatMatrix (); - } - - return retval; - } - - FloatMatrix - gammainc (const FloatMatrix& x, const FloatMatrix& a) - { - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - - if (nr != a_nr || nc != a_nc) - (*current_liboctave_error_handler) - ("gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)", - nr, nc, a_nr, a_nc); - - FloatMatrix retval (nr, nc); - - bool err; - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - { - retval(i,j) = gammainc (x(i,j), a(i,j), err); - - if (err) - return FloatMatrix (); - } - - return retval; - } - - FloatNDArray - gammainc (float x, const FloatNDArray& a) - { - dim_vector dv = a.dims (); - octave_idx_type nel = dv.numel (); - - FloatNDArray retval (dv); - - bool err; - - for (octave_idx_type i = 0; i < nel; i++) - { - retval(i) = gammainc (x, a(i), err); - - if (err) - return FloatNDArray (); - } - - return retval; - } - - FloatNDArray - gammainc (const FloatNDArray& x, float a) - { - dim_vector dv = x.dims (); - octave_idx_type nel = dv.numel (); - - FloatNDArray retval (dv); - - bool err; - - for (octave_idx_type i = 0; i < nel; i++) - { - retval(i) = gammainc (x(i), a, err); - - if (err) - return FloatNDArray (); - } - - return retval; - } - - FloatNDArray - gammainc (const FloatNDArray& x, const FloatNDArray& a) - { - dim_vector dv = x.dims (); - octave_idx_type nel = dv.numel (); - - if (dv != a.dims ()) - { - std::string x_str = dv.str (); - std::string a_str = a.dims ().str (); - - (*current_liboctave_error_handler) - ("gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)", - x_str.c_str (), a_str.c_str ()); - } - - FloatNDArray retval (dv); - - bool err; - - for (octave_idx_type i = 0; i < nel; i++) - { - retval(i) = gammainc (x(i), a(i), err); - - if (err) - return FloatNDArray (); - } - - return retval; - } - Complex log1p (const Complex& x) { @@ -3575,6 +2373,7 @@ FloatComplexNDArray airy (const FloatComplexNDArray& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::airy (z, deriv, scaled, ierr); } FloatComplexNDArray biry (const FloatComplexNDArray& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::biry (z, deriv, scaled, ierr); } +<<<<<<< dest Array<double> betainc (double x, double a, const Array<double>& b) { return octave::math::betainc (x, a, b); } Array<double> betainc (double x, const Array<double>& a, double b) { return octave::math::betainc (x, a, b); } Array<double> betainc (double x, const Array<double>& a, const Array<double>& b) { return octave::math::betainc (x, a, b); } @@ -3593,22 +2392,6 @@ Array<float> betainc (const Array<float>& x, const Array<float>& a, float b) { return octave::math::betainc (x, a, b); } Array<float> betainc (const Array<float>& x, const Array<float>& a, const Array<float>& b) { return octave::math::betainc (x, a, b); } -Matrix gammainc (double x, const Matrix& a) { return octave::math::gammainc (x, a); } -Matrix gammainc (const Matrix& x, double a) { return octave::math::gammainc (x, a); } -Matrix gammainc (const Matrix& x, const Matrix& a) { return octave::math::gammainc (x, a); } - -NDArray gammainc (double x, const NDArray& a) { return octave::math::gammainc (x, a); } -NDArray gammainc (const NDArray& x, double a) { return octave::math::gammainc (x, a); } -NDArray gammainc (const NDArray& x, const NDArray& a) { return octave::math::gammainc (x, a); } - -FloatMatrix gammainc (float x, const FloatMatrix& a) { return octave::math::gammainc (x, a); } -FloatMatrix gammainc (const FloatMatrix& x, float a) { return octave::math::gammainc (x, a); } -FloatMatrix gammainc (const FloatMatrix& x, const FloatMatrix& a) { return octave::math::gammainc (x, a); } - -FloatNDArray gammainc (float x, const FloatNDArray& a) { return octave::math::gammainc (x, a); } -FloatNDArray gammainc (const FloatNDArray& x, float a) { return octave::math::gammainc (x, a); } -FloatNDArray gammainc (const FloatNDArray& x, const FloatNDArray& a) { return octave::math::gammainc (x, a); } - Array<double> betaincinv (double x, double a, const Array<double>& b) { return octave::math::betaincinv (x, a, b); } Array<double> betaincinv (double x, const Array<double>& a, double b) { return octave::math::betaincinv (x, a, b); } Array<double> betaincinv (double x, const Array<double>& a, const Array<double>& b) { return octave::math::betaincinv (x, a, b); }
--- a/liboctave/numeric/lo-specfun.h Sun Mar 18 18:43:04 2018 -0700 +++ b/liboctave/numeric/lo-specfun.h Mon Mar 19 13:50:10 2018 -0700 @@ -282,55 +282,6 @@ extern OCTAVE_API FloatComplexMatrix besselh2 (const FloatRowVector& alpha, const FloatComplexColumnVector& x, bool scaled, Array<octave_idx_type>& ierr); - extern OCTAVE_API double betainc (double x, double a, double b); - extern OCTAVE_API Array<double> betainc (double x, double a, - const Array<double>& b); - extern OCTAVE_API Array<double> betainc (double x, const Array<double>& a, - double b); - extern OCTAVE_API Array<double> betainc (double x, const Array<double>& a, - const Array<double>& b); - extern OCTAVE_API Array<double> betainc (const Array<double>& x, double a, - double b); - extern OCTAVE_API Array<double> betainc (const Array<double>& x, double a, - const Array<double>& b); - extern OCTAVE_API Array<double> betainc (const Array<double>& x, - const Array<double>& a, double b); - extern OCTAVE_API Array<double> betainc (const Array<double>& x, - const Array<double>& a, const Array<double>& b); - - extern OCTAVE_API float betainc (float x, float a, float b); - extern OCTAVE_API Array<float> betainc (float x, float a, - const Array<float>& b); - extern OCTAVE_API Array<float> betainc (float x, const Array<float>& a, - float b); - extern OCTAVE_API Array<float> betainc (float x, const Array<float>& a, - const Array<float>& b); - extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, - float b); - extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, - const Array<float>& b); - extern OCTAVE_API Array<float> betainc (const Array<float>& x, - const Array<float>& a, float b); - extern OCTAVE_API Array<float> betainc (const Array<float>& x, - const Array<float>& a, const Array<float>& b); - - extern OCTAVE_API double betaincinv (double x, double a, double b); - extern OCTAVE_API Array<double> betaincinv (double x, double a, - const Array<double>& b); - extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, - double b); - extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, - const Array<double>& b); - - extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, - double b); - extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, - const Array<double>& b); - extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, - const Array<double>& a, double b); - extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, - const Array<double>& a, const Array<double>& b); - extern OCTAVE_API Complex biry (const Complex& z, bool deriv, bool scaled, octave_idx_type& ierr); extern OCTAVE_API ComplexMatrix biry (const ComplexMatrix& z, bool deriv, @@ -771,42 +722,6 @@ OCTAVE_DEPRECATED (4.2, "use 'octave::math::biry' instead") extern OCTAVE_API FloatComplexNDArray biry (const FloatComplexNDArray& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -inline double betainc (double x, double a, double b) { return octave::math::betainc (x, a, b); } -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<double> betainc (double x, double a, const Array<double>& b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<double> betainc (double x, const Array<double>& a, double b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<double> betainc (double x, const Array<double>& a, const Array<double>& b); - -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -inline float betainc (float x, float a, float b) { return octave::math::betainc (x, a, b); } -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<double> betainc (const Array<double>& x, double a, const Array<double>& b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<double> betainc (const Array<double>& x, const Array<double>& a, double b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<double> betainc (const Array<double>& x, const Array<double>& a, const Array<double>& b); - -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API float betainc (float x, float a, float b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<float> betainc (float x, float a, const Array<float>& b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<float> betainc (float x, const Array<float>& a, float b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<float> betainc (float x, const Array<float>& a, const Array<float>& b); - -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, float b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, const Array<float>& b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<float> betainc (const Array<float>& x, const Array<float>& a, float b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead") -extern OCTAVE_API Array<float> betainc (const Array<float>& x, const Array<float>& a, const Array<float>& b); - OCTAVE_DEPRECATED (4.2, "use 'octave::math::gammainc' instead") inline double gammainc (double x, double a, bool& err) { return octave::math::gammainc (x, a, err); } OCTAVE_DEPRECATED (4.2, "use 'octave::math::gammainc' instead") @@ -887,27 +802,6 @@ OCTAVE_DEPRECATED (4.2, "use 'octave::math::dawson' instead") inline FloatComplex dawson (const FloatComplex& x) { return octave::math::dawson (x); } -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -inline double betaincinv (double x, double a, double b) { return octave::math::betaincinv (x, a, b); } - -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API double betaincinv (double x, double a, double b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API Array<double> betaincinv (double x, double a, const Array<double>& b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, double b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, const Array<double>& b); - -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, double b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, const Array<double>& b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, const Array<double>& a, double b); -OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead") -extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, const Array<double>& a, const Array<double>& b); - OCTAVE_DEPRECATED (4.2, "use 'octave::math::ellipj' instead") inline void ellipj (double u, double m, double& sn, double& cn, double& dn, double& err) { octave::math::ellipj (u, m, sn, cn, dn, err); } OCTAVE_DEPRECATED (4.2, "use 'octave::math::ellipj' instead")
--- a/scripts/help/bessel.m Sun Mar 18 18:43:04 2018 -0700 +++ b/scripts/help/bessel.m Mon Mar 19 13:50:10 2018 -0700 @@ -76,7 +76,7 @@ ## machine accuracy. ## ## @item -## Complete loss of significance by argument reduction, return @code{NaN}. +## Loss of significance by argument reduction, output may be inaccurate. ## ## @item ## Error---no computation, algorithm termination condition not met, return
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/specfun/betainc.m Mon Mar 19 13:50:10 2018 -0700 @@ -0,0 +1,240 @@ +## Copyright (C) 2018 Stefan Schlögl +## Copyright (C) 2018 Michele Ginesi +## +## This file is part of Octave. +## +## Octave is free software: you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <https://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {} {} betainc (@var{x}, @var{a}, @var{b}) +## @deftypefnx {} {} betainc (@var{x}, @var{a}, @var{b}, @var{tail}) +## Compute the incomplete beta function. +## +## This is defined as +## @tex +## $$ +## I_x (a, b) = {1 \over {B(a,b)}} \displaystyle{\int_0^x t^{a-1} (1-t)^{b-1} dt} +## $$ +## @end tex +## @ifnottex +## +## @example +## @group +## x +## / +## | +## I_x (a, b) = | t^(a-1) (1-t)^(b-1) dt +## | +## / +## 0 +## @end group +## @end example +## +## @end ifnottex +## +## with real @var{x} in the range [0,1]. The inputs @var{a} and @var{b} must +## be real and strictly positive (> 0). If one of the inputs is not a scalar +## then the other inputs must be scalar or of compatible dimensions. +## +## By default, @var{tail} is @qcode{"lower"} and the incomplete beta function +## integrated from 0 to @var{x} is computed. If @var{tail} is @qcode{"upper"} +## then the complementary function integrated from @var{x} to 1 is calculated. +## The two choices are related by +## +## betainc (@var{x}, @var{a}, @var{b}, @qcode{"upper"}) = +## 1 - betainc (@var{x}, @var{a}, @var{b}, @qcode{"lower"}). +## +## @code{betainc} uses a more sophisticated algorithm than subtraction to +## get numerically accurate results when the @qcode{"lower"} value is small. +## +## Reference: @nospell{A. Cuyt, V. Brevik Petersen, B. Verdonk, H. Waadeland, +## W.B. Jones}, @cite{Handbook of Continued Fractions for Special Functions}, +## ch.@: 18. +## +## @seealso{beta, betaincinv, betaln} +## @end deftypefn + +function y = betainc (x, a, b, tail = "lower") + + if (nargin < 3 || nargin > 4) + print_usage (); + endif + + [err, x, a, b] = common_size (x, a, b); + if (err > 0) + error ("betainc: X, A, and B must be of common size or scalars"); + endif + + if (iscomplex (x) || iscomplex (a) || iscomplex (b)) + error ("betainc: all inputs must be real"); + endif + + ## Remember original shape of data, but convert to column vector for calcs. + orig_sz = size (x); + x = x(:); + a = a(:); + b = b(:); + + if (any ((x < 0) | (x > 1))) + error ("betainc: X must be in the range [0, 1]"); + endif + + if (any (a <= 0)) + error ("betainc: A must be strictly positive"); + endif + + if (any (b <= 0)) + error ("betainc: B must be strictly positive"); + endif + + ## If any of the arguments is single then the output should be as well. + if (strcmp (class (x), "single") || strcmp (class (a), "single") + || strcmp (class (b), "single")) + a = single (a); + b = single (b); + x = single (x); + endif + + ## Convert to floating point if necessary + if (isinteger (x)) + y = double (x); + endif + if (isinteger (a)) + a = double (a); + endif + if (isinteger (b)) + b = double (b); + endif + + ## Initialize output array + y = zeros (size (x), class (x)); + + ## In the following, we use the fact that the continued fraction Octave uses + ## is more efficient when x <= a / (a + b). Moreover, to compute the upper + ## version, which is defined as I_x(a,b,"upper") = 1 - I_x(a,b) we use the + ## property I_x(a,b) + I_(1-x) (b,a) = 1. + + if (strcmpi (tail, "lower")) + fflag = (x > a./(a+b)); + x(fflag) = 1 - x(fflag); + [a(fflag), b(fflag)] = deal (b(fflag), a(fflag)); + elseif (strcmpi (tail, "upper")) + fflag = (x < (a ./ (a + b))); + x(! fflag) = 1 - x(! fflag); + [a(! fflag), b(! fflag)] = deal (b(! fflag), a(! fflag)); + else + error ("betainc: invalid value for TAIL"); + endif + + f = zeros (size (x), class (x)); + + ## Continued fractions: CPVWJ, formula 18.5.20, modified Lentz algorithm + ## implemented in a separate .cc file. This particular continued fraction + ## gives (B(a,b) * I_x(a,b)) / (x^a * (1-x)^b). + + f = __betainc__ (x, a, b); + + ## Divide continued fraction by B(a,b) / (x^a * (1-x)^b) to obtain I_x(a,b). + y = a .* log (x) + b .* log1p (-x) + gammaln (a + b) - ... + gammaln (a) - gammaln (b) + log (f); + y = real (exp (y)); + y(fflag) = 1 - y(fflag); + + ## Restore original shape + y = reshape (y, orig_sz); + +endfunction + + +## Double precision +%!test +%! a = [1, 1.5, 2, 3]; +%! b = [4, 3, 2, 1]; +%! v1 = betainc (1, a, b); +%! v2 = [1,1,1,1]; +%! x = [.2, .4, .6, .8]; +%! v3 = betainc (x, a, b); +%! v4 = 1 - betainc (1-x, b, a); +%! assert (v1, v2, sqrt (eps)); +%! assert (v3, v4, sqrt (eps)); + +## Single precision +%!test +%! a = single ([1, 1.5, 2, 3]); +%! b = single ([4, 3, 2, 1]); +%! v1 = betainc (1, a, b); +%! v2 = single ([1,1,1,1]); +%! x = single ([.2, .4, .6, .8]); +%! v3 = betainc (x, a, b); +%! v4 = 1 - betainc (1-x, b, a); +%! assert (v1, v2, sqrt (eps ("single"))); +%! assert (v3, v4, sqrt (eps ("single"))); + +## Mixed double/single precision +%!test +%! a = single ([1, 1.5, 2, 3]); +%! b = [4, 3, 2, 1]; +%! v1 = betainc (1,a,b); +%! v2 = single ([1,1,1,1]); +%! x = [.2, .4, .6, .8]; +%! v3 = betainc (x, a, b); +%! v4 = 1-betainc (1.-x, b, a); +%! assert (v1, v2, sqrt (eps ("single"))); +%! assert (v3, v4, sqrt (eps ("single"))); + +%!test <51157> +%! y = betainc ([0.00780;0.00782;0.00784],250.005,49750.995); +%! y_ex = [0.999999999999989; 0.999999999999992; 0.999999999999995]; +%! assert (y, y_ex, -1e-14); + +%!assert (betainc (0.001, 20, 30), 2.750687665855991e-47, -3e-14); +%!assert (betainc (0.0001, 20, 30), 2.819953178893307e-67, -3e-14); +%!assert (betainc (0.99, 20, 30, "upper"), 1.5671643161872703e-47, -3e-14); +%!assert (betainc (0.999, 20, 30, "upper"), 1.850806276141535e-77, -3e-14); +%!assert (betainc (0.5, 200, 300), 0.9999964565197356, -1e-15); +%!assert (betainc (0.5, 200, 300, "upper"), 3.54348026439253e-06, -1e-13); + +## Test trivial values +%!test +%! [a,b] = ndgrid (linspace (1e-4, 100, 20), linspace (1e-4, 100, 20)); +%! assert (betainc (0, a, b), zeros (20)); +%! assert (betainc (1, a, b), ones (20)); + +## Test input validation +%!error betainc () +%!error betainc (1) +%!error betainc (1,2) +%!error betainc (1,2,3,4,5) +%!error <must be of common size or scalars> betainc (ones (2,2), ones (1,2), 1) +%!error <all inputs must be real> betainc (0.5i, 1, 2) +%!error <all inputs must be real> betainc (0, 1i, 1) +%!error <all inputs must be real> betainc (0, 1, 1i) +%!error <X must be in the range \[0, 1\]> betainc (-0.1,1,1) +%!error <X must be in the range \[0, 1\]> betainc (1.1,1,1) +%!error <X must be in the range \[0, 1\]> +%! x = ones (1, 1, 2); +%! x(1,1,2) = -1; +%! betainc (x,1,1); +%!error <A must be strictly positive> betainc (0.5,0,1) +%!error <A must be strictly positive> +%! a = ones (1, 1, 2); +%! a(1,1,2) = 0; +%! betainc (1,a,1); +%!error <B must be strictly positive> betainc (0.5,1,0) +%!error <B must be strictly positive> +%! b = ones (1, 1, 2); +%! b(1,1,2) = 0; +%! betainc (1,1,b); +%!error <invalid value for TAIL> betainc (1,2,3, "foobar")
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/specfun/betaincinv.m Mon Mar 19 13:50:10 2018 -0700 @@ -0,0 +1,302 @@ +## Copyright (C) 2017 Michele Ginesi +## +## This file is part of Octave. +## +## Octave is free software: you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation, either version 3 of the License, or +## (at your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <https://www.gnu.org/licenses/>. + +## Author: Michele Ginesi <michele.ginesi@gmail.com> + +## -*- texinfo -*- +## @deftypefn {} {} betaincinv (@var{y}, @var{a}, @var{b}) +## @deftypefnx {} {} betaincinv (@var{y}, @var{a}, @var{b}, "lower") +## @deftypefnx {} {} betaincinv (@var{y}, @var{a}, @var{b}, "upper") +## Compute the inverse of the normalized incomplete beta function. +## +## The normalized incomplete beta function is defined as +## @tex +## $$ +## I_x (a, b) = {1 \over {B(a,b)}} \displaystyle{\int_0^x t^{a-1} (1-t)^{b-1} dt} +## $$ +## @end tex +## @ifnottex +## +## @example +## @group +## x +## / +## | +## I_x (a, b) = | t^(a-1) (1-t)^(b-1) dt +## | +## / +## 0 +## @end group +## @end example +## +## @end ifnottex +## +## If two inputs are scalar, then @code{betaincinv (@var{y}, @var{a}, @var{b})} +## is returned for each of the other inputs. +## +## If two or more inputs are not scalar, the sizes of them must agree, and +## @code{betaincinv} is applied element-by-element. +## +## The variable @var{y} must be in the interval [0,1], while @var{a} and +## @var{b} must be real and strictly positive. +## +## By default, @var{tail} is @qcode{"lower")} and the inverse of the incomplete +## beta function integrated from 0 to @var{x} is computed. If @var{tail} is +## @qcode{"upper"} then the complementary function integrated from @var{x} to 1 +## is inverted. +## +## The function is computed by standard Newton's method, by solving +## @tex +## $$ +## y - I_x (a, b) = 0 +## $$ +## @end tex +## @ifnottex +## +## @example +## @var{y} - betainc (@var{x}, @var{a}, @var{b}) = 0 +## @end example +## +## @end ifnottex +## +## @seealso{betainc, beta, betaln} +## @end deftypefn + +function x = betaincinv (y, a, b, tail = "lower") + + if (nargin < 3 || nargin > 4) + print_usage (); + endif + + [err, y, a, b] = common_size (y, a, b); + if (err > 0) + error ("betaincinv: Y, A, and B must be of common size or scalars"); + endif + + if (iscomplex (y) || iscomplex (a) || iscomplex (b)) + error ("betaincinv: all inputs must be real"); + endif + + ## FIXME: Should there be isnumeric checking? Right now it accepts char + ## arrays, but then produces a weird error later on. + + ## Remember original shape of data, but convert to column vector for calcs. + orig_sz = size (y); + y = y(:); + a = a(:); + b = b(:); + + if (any ((y < 0) | (y > 1))) + error ("betaincinv: Y must be in the range [0, 1]"); + endif + + if (any (a <= 0)) + error ("betaincinv: A must be strictly positive"); + endif + + if (any (b <= 0)) + error ("betaincinv: B must be strictly positive"); + endif + + ## If any of the arguments is single then the output should be as well. + if (strcmp (class (y), "single") || strcmp (class (a), "single") + || strcmp (class (b), "single")) + y = single (y); + a = single (a); + b = single (b); + endif + + ## Convert to floating point if necessary + if (isinteger (y)) + y = double (y); + endif + if (isinteger (a)) + a = double (a); + endif + if (isinteger (b)) + b = double (b); + endif + + ## Initialize output array + x = zeros (size (y), class (y)); + + ## Parameters for the Newton method + maxit = 20; + tol = eps (class (y)); + + if (strcmpi (tail, "lower")) + p = y; + q = 1 - y; + x(y == 0) = 0; + x(y == 1) = 1; + elseif (strcmpi (tail, "upper")) + p = 1 - y; + q = y; + x(y == 0) = 1; + x(y == 1) = 0; + else + error ("betaincinv: invalid value for TAIL") + endif + + ## Special values have been already computed. + todo = (y != 0) & (y != 1); + + ## We will invert the lower version for p < 0.5 and the upper otherwise. + i_low = (p < 0.5); + i_upp = (! i_low); + + idx = todo & i_low; + if (any (idx)); + n = nnz (idx); + ## Function and derivative of the lower version. + F = @(x, a, b, y) y - betainc (x, a, b); + JF = @(x, a, b) - real (exp ((a-1) .* log (x) + (b-1) .* log1p (-x) + ... + gammaln (a+b) - gammaln (a) - gammaln (b))); + + ## Compute the initial guess with a bisection method of 10 steps. + x0 = bisection_method (F, zeros (n,1), ones (n,1), ... + a(i_low), b(i_low), p(i_low), 10); + + ## Use Newton's method to iteratively find solution. + x(i_low) = newton_method (F, JF, x0, a(i_low), b(i_low), p(i_low), ... + tol, maxit); + endif + + idx = todo & i_upp; + if (any (idx)); + n = nnz (idx); + ## Function and derivative of the upper version. + F = @(x, a, b, y) y - betainc (x, a, b, "upper"); + JF = @(x, a, b) real (exp ((a-1) .* log (x) + (b-1) .* log1p (-x) + ... + gammaln (a+b) - gammaln (a) - gammaln (b))); + + ## Compute the initial guess with a bisection method of 10 steps. + x0 = bisection_method (F, zeros (n,1), ones (n,1), ... + a(i_upp), b(i_upp), q(i_upp), 10); + + ## Use Newton's method to iteratively find solution. + x(i_upp) = newton_method (F, JF, x0, a(i_upp), b(i_upp), q(i_upp), ... + tol, maxit); + endif + + ## Restore original shape + x = reshape (x, orig_sz); + +endfunction + + +## Subfunctions: Bisection and Newton Methods +function xc = bisection_method (F, xl, xr, a, b, y, maxit) + F_l = F (xl, a, b, y); + F_r = F (xr, a, b, y); + for it = 1:maxit + xc = (xl + xr) / 2; + F_c = F (xc, a, b, y); + flag_l = ((F_c .* F_r) < 0); + flag_r = ((F_c .* F_l) < 0); + flag_c = (F_c == 0); + xl(flag_l) = xc(flag_l); + xr(flag_r) = xc(flag_r); + xl(flag_c) = xr(flag_c) = xc(flag_c); + F_l(flag_l) = F_c(flag_l); + F_r(flag_r) = F_c(flag_r); + F_l(flag_c) = F_r(flag_c) = 0; + endfor +endfunction + +function x = newton_method (F, JF, x0, a, b, y, tol, maxit); + l = numel (y); + res = -F (x0, a, b, y) ./ JF (x0, a, b); + todo = (abs(res) >= tol * abs (x0)); + x = x0; + it = 0; + while (any (todo) && (it < maxit)) + it++; + x(todo) += res(todo); + res(todo) = -F(x(todo), a(todo), b(todo), y(todo)) ... + ./ JF (x(todo), a(todo), b(todo)); + todo = (abs(res) >= tol * abs (x)); + endwhile + x += res; +endfunction + + +%!test +%! x = linspace (0.1, 0.9, 11); +%! a = [2, 3, 4]; +%! [x,a,b] = ndgrid (x,a,a); +%! xx = betaincinv (betainc (x, a, b), a, b); +%! assert (xx, x, 3e-15); + +%!test +%! x = linspace (0.1, 0.9, 11); +%! a = [2, 3, 4]; +%! [x,a,b] = ndgrid (x,a,a); +%! xx = betaincinv (betainc (x, a, b, "upper"), a, b, "upper"); +%! assert (xx, x, 3e-15); + +%!test +%! x = linspace (0.1, 0.9, 11); +%! a = [0.1:0.1:1]; +%! [x,a,b] = ndgrid (x,a,a); +%! xx = betaincinv (betainc (x, a, b), a, b); +%! assert (xx, x, 3e-15); + +%!test +%! x = linspace (0.1, 0.9, 11); +%! a = [0.1:0.1:1]; +%! [x,a,b] = ndgrid (x,a,a); +%! xx = betaincinv (betainc (x, a, b, "upper"), a, b, "upper"); +%! assert (xx, x, 3e-15); + +## Test the conservation of the input class +%!assert (class (betaincinv (0.5, 1, 1)), "double") +%!assert (class (betaincinv (single (0.5), 1, 1)), "single") +%!assert (class (betaincinv (0.5, single (1), 1)), "single") +%!assert (class (betaincinv (int8 (0), 1, 1)), "double") +%!assert (class (betaincinv (0.5, int8 (1), 1)), "double") +%!assert (class (betaincinv (int8 (0), single (1), 1)), "single") +%!assert (class (betaincinv (single (0.5), int8 (1), 1)), "single") + +## Test input validation +%!error betaincinv () +%!error betaincinv (1) +%!error betaincinv (1,2) +%!error betaincinv (1,2,3,4,5) +%!error <must be of common size or scalars> +%! betaincinv (ones (2,2), ones (1,2), 1); +%!error <all inputs must be real> betaincinv (0.5i, 1, 2) +%!error <all inputs must be real> betaincinv (0, 1i, 1) +%!error <all inputs must be real> betaincinv (0, 1, 1i) +%!error <Y must be in the range \[0, 1\]> betaincinv (-0.1,1,1) +%!error <Y must be in the range \[0, 1\]> betaincinv (1.1,1,1) +%!error <Y must be in the range \[0, 1\]> +%! y = ones (1, 1, 2); +%! y(1,1,2) = -1; +%! betaincinv (y,1,1); +%!error <A must be strictly positive> betaincinv (0.5,0,1) +%!error <A must be strictly positive> +%! a = ones (1, 1, 2); +%! a(1,1,2) = 0; +%! betaincinv (1,a,1); +%!error <B must be strictly positive> betaincinv (0.5,1,0) +%!error <B must be strictly positive> +%! b = ones (1, 1, 2); +%! b(1,1,2) = 0; +%! betaincinv (1,1,b); +%!error <invalid value for TAIL> betaincinv (1,2,3, "foobar")
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/specfun/cosint.m Mon Mar 19 13:50:10 2018 -0700 @@ -0,0 +1,268 @@ +## Copyright (C) 2017 Michele Ginesi +## +## This file is part of Octave. +## +## Octave is free software: you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation, either version 3 of the License, or +## (at your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <https://www.gnu.org/licenses/>. + +## Author: Michele Ginesi <michele.ginesi@gmail.com> + +## -*- texinfo -*- +## @deftypefn {} {} cosint (@var{x}) +## Compute the cosine integral function: +## @tex +## $$ +## {\rm Ci} (x) = - \int_x^\infty {{\cos (t)} \over t} dt +## $$ +## @end tex +## @ifnottex +## +## @example +## @group +## +oo +## / +## Ci (x) = - | (cos (t)) / t dt +## / +## x +## @end group +## @end example +## +## @end ifnottex +## An equivalent definition is +## @tex +## $$ +## {\rm Ci} (x) = \gamma + \log (x) + \int_0^x {{\cos (t) - 1} \over t} dt +## $$ +## @end tex +## @ifnottex +## +## @example +## @group +## x +## / +## | cos (t) - 1 +## Ci (x) = gamma + log (x) + | ------------- dt +## | t +## / +## 0 +## @end group +## @end example +## +## @end ifnottex +## Reference: +## +## @nospell{M. Abramowitz and I.A. Stegun}, +## @cite{Handbook of Mathematical Functions} +## 1964. +## +## @seealso{sinint, expint, cos} +## +## @end deftypefn + +function y = cosint (x) + + if (nargin != 1) + print_usage (); + endif + + if (! isnumeric (x)) + error ("cosint: X must be numeric"); + endif + + ## Convert to floating point if necessary + if (isinteger (x)) + x = double (x); + endif + + ## Convert to column vector + orig_sz = size (x); + if (iscomplex (x)) + ## Work around reshape which narrows to real (bug #52953) + x = complex (real (x)(:), imag (x)(:)); + else + x = x(:); + end + + ## Initialize the result + y = zeros (size (x), class (x)); + tol = eps (class (x)); + + todo = true (size (x)); + + ## Special values + y(x == Inf) = 0; + y((x == -Inf) & !signbit (imag (x))) = 1i * pi; + y((x == -Inf) & signbit (imag (x))) = -1i * pi; + + todo(isinf (x)) = false; + + ## For values large in modulus, but not in the range (-oo,0), we use the + ## relation with expint. + + flag_large = (abs (x) > 2); + xx = x(flag_large); + + ## Abramowitz, relation 5.2.20 + ii_sw = (real (xx) <= 0 & imag (xx) < 0); + xx(ii_sw) = conj (xx(ii_sw)); + ii_nw = (real (xx) < 0); + xx(ii_nw) *= -1; + yy = -0.5 * (expint (1i * xx) + expint (-1i * xx)); + yy(ii_nw) += 1i * pi; + yy(ii_sw) = conj (yy(ii_sw)); + y(todo & flag_large) = yy; + + todo(flag_large) = false; + + ## For values small in modulus, use the series expansion (also near (-oo, 0]) + if (iscomplex (x)) + ## indexing can lose imag part: if it was -0, we could end up on the + ## wrong right side of the branch cut along the negative real axis. + xx = complex (real (x)(todo), imag (x)(todo)); + else + xx = x(todo); + end + ssum = - xx .^ 2 / 4; # First term of the series expansion + ## FIXME: This is way more precision than a double value can hold. + gma = 0.57721566490153286060651209008; # Euler gamma constant + yy = gma + log (complex (xx)) + ssum; # log(complex(...) handles signed zero + flag_sum = true (nnz (todo), 1); + it = 0; + maxit = 300; + while (any (flag_sum) && (++it < maxit)); + ssum .*= - xx .^ 2 * (2 * it) / ((2 * it + 2) ^ 2 * (2 * it + 1)); + yy(flag_sum) += ssum (flag_sum); + flag_sum = (abs (ssum) >= tol); + endwhile + y(todo) = yy; + + ## Clean up values which are purely real + flag_neg_zero_imag = (real (x) < 0) & (imag (x) == 0) & signbit (imag (x)); + y(flag_neg_zero_imag) = complex (real (y(flag_neg_zero_imag)), -pi); + + ## Restore original shape + y = reshape (y, orig_sz); + +endfunction + + +%!assert (cosint (1.1), 0.38487337742465081550, 2 * eps); + +%!test +%! x = [2, 3, pi; exp(1), 5, 6]; +%! A = cosint (x); +%! B = [0.422980828774864996, 0.119629786008000328, 0.0736679120464254860; ... +%! 0.213958001340379779, -0.190029749656643879, -0.0680572438932471262]; +%! assert (A, B, -5e-15); + +%!assert (cosint (0), - Inf) +%!assert (cosint (-0), -inf + 1i*pi) +%!assert (cosint (complex (-0, 0)), -inf + 1i*pi) +%!assert (cosint (complex (-0, -0)), -inf - 1i*pi) +%!assert (cosint (inf), 0) +%!assert (cosint (-inf), 1i * pi) +%!assert (cosint (complex (-inf, -0)), -1i * pi) +%!assert (isnan (cosint (nan))) + +%!assert (class (cosint (single (1))), "single") + +## tests against maple +%!assert (cosint (1), 0.337403922900968135, -2*eps) +%!assert (cosint (-1), 0.337403922900968135 + 3.14159265358979324*I, -2*eps) +%!assert (cosint (pi), 0.0736679120464254860, -2*eps) +%!assert (cosint (-pi), 0.0736679120464254860 + 3.14159265358979324*I, -2*eps) +%!assert (cosint (300), -0.00333219991859211178, -2*eps) +%!assert (cosint (1e4), -0.0000305519167244852127, -2*eps) +%!assert (cosint (20i), 1.28078263320282944e7 + 1.57079632679489662*I, -2*eps) + +%!test +%! x = (0:4).'; +%! y_ex = [-Inf +%! 0.337403922900968135 +%! 0.422980828774864996 +%! 0.119629786008000328 +%! -0.140981697886930412]; +%! assert (cosint (x), y_ex, -3e-15); + +%!test +%! x = -(1:4).'; +%! y_ex = [0.337403922900968135 + pi*1i +%! 0.422980828774864996 + pi*1i +%! 0.119629786008000328 + pi*1i +%! -0.140981697886930412 + pi*1i]; +%! assert (cosint (x), y_ex, -4*eps); + +%!test +%! x = complex (-(1:4).', 0); +%! y_ex = [0.337403922900968135 + pi*1i +%! 0.422980828774864996 + pi*1i +%! 0.119629786008000328 + pi*1i +%! -0.140981697886930412 + pi*1i]; +%! assert (cosint (x), y_ex, -4*eps); + +%!test +%! x = complex (-(1:4).', -0); +%! y_ex = [0.337403922900968135 - pi*1i +%! 0.422980828774864996 - pi*1i +%! 0.119629786008000328 - pi*1i +%! -0.140981697886930412 - pi*1i]; +%! assert (cosint (x), y_ex, -4*eps); + +%!test +%! x = 1i * (0:4).'; +%! y_ex = [-Inf +%! 0.837866940980208241 + 1.57079632679489662*I +%! 2.45266692264691452 + 1.57079632679489662*I +%! 4.96039209476560976 + 1.57079632679489662*I +%! 9.81354755882318556 + 1.57079632679489662*I]; +%! assert (cosint (x), y_ex, -4*eps); + +%!test +%! x = -1i * (1:4).'; +%! y_ex = [0.837866940980208241 - 1.57079632679489662*I +%! 2.45266692264691452 - 1.57079632679489662*I +%! 4.96039209476560976 - 1.57079632679489662*I +%! 9.81354755882318556 - 1.57079632679489662*I]; +%! assert (cosint (x), y_ex, -4*eps); + +%!test +%! x = [1+2i; -2+5i; 2-5i; 100; 10i; -1e-4 + 1e-6*1i; -20-1i]; +%! A = [ 2.03029639329172164 - 0.151907155175856884*I +%! 1.61538963829107749 + 19.7257540553382650*I +%! 1.61538963829107749 + 16.5841614017484717*I +%! -0.00514882514261049214 +%! 1246.11448604245441 + 1.57079632679489662*I +%! -8.63307471207423322 + 3.13159298695312800*I +%! 0.0698222284673061493 - 3.11847446254772946*I ]; +%! B = cosint (x); +%! assert (A, B, -3*eps); +%! B = cosint (single (x)); +%! assert (A, B, -3*eps ("single")); + +## Fails along negative real axis +%!test +%! x = [-25; -100; -1000]; +%! yex = [-0.0068485971797025909189 + pi*1i +%! -0.0051488251426104921444 + pi*1i +%! 0.000826315511090682282 + pi*1i]; +%! y = cosint (x); +%! assert (y, yex, -5*eps); + +## FIXME: Need a test for bug #52953 +%#!test <*52953> + +## Test input validation +%!error cosint () +%!error cosint (1,2) +%!error <X must be numeric> cosint ("1")
--- a/scripts/specfun/expint.m Sun Mar 18 18:43:04 2018 -0700 +++ b/scripts/specfun/expint.m Mon Mar 19 13:50:10 2018 -0700 @@ -1,5 +1,4 @@ -## Copyright (C) 2006-2017 Sylvain Pelissier -## Copyright (C) 2017 Michele Ginesi +## Copyright (C) 2018 Michele Ginesi ## ## This file is part of Octave. ## @@ -17,12 +16,10 @@ ## along with Octave; see the file COPYING. If not, see ## <https://www.gnu.org/licenses/>. -## Authors: Sylvain Pelissier <sylvain.pelissier@gmail.com> -## Michele Ginesi <michele.ginesi@gmail.com> - ## -*- texinfo -*- ## @deftypefn {} {} expint (@var{x}) ## Compute the exponential integral: +## ## @tex ## $$ ## {\rm E_1} (x) = \int_x^\infty {e^{-t} \over t} dt @@ -43,7 +40,8 @@ ## @end example ## ## @end ifnottex -## Note: For compatibility, this functions uses the @sc{matlab} definition +## +## Note: For compatibility, this function uses the @sc{matlab} definition ## of the exponential integral. Most other sources refer to this particular ## value as @math{E_1 (x)}, and the exponential integral as ## @tex @@ -75,6 +73,16 @@ ## @ifnottex ## @w{@code{E_1 (-x) = -Ei (x) - i*pi}}. ## @end ifnottex +## +## References: +## +## @nospell{M. Abramowitz and I.A. Stegun}, +## @cite{Handbook of Mathematical Functions}, 1964. +## +## @nospell{N. Bleistein and R.A. Handelsman}, +## @cite{Asymptotic expansions of integrals}, 1986. +## +## @seealso{cosint, sinint, exp} ## @end deftypefn function E1 = expint (x) @@ -85,11 +93,14 @@ if (! isnumeric (x)) error ("expint: X must be numeric"); - elseif (isinteger (x)) + endif + + ## Convert to floating point if necessary + if (isinteger (x)) x = double (x); endif - sparse_x = issparse (x); + orig_sparse = issparse (x); orig_sz = size (x); x = x(:); # convert to column vector @@ -102,19 +113,18 @@ tol = eps (class (x)); ## Divide the input into 3 regions and apply a different algorithm for each. - s_idx = (((real (x) + 19.5).^2 ./ (20.5^2) + imag (x).^2 ./ (10^2)) <= 1) ... + ## s = series expansion, cf = continued fraction, a = asymptotic series + s_idx = (((real (x) + 19.5).^ 2 ./ (20.5^2) + ... + imag (x).^2 ./ (10^2)) <= 1) ... | (real (x) < 0 & abs (imag (x)) <= 1e-8); - cf_idx = ((((real (x) + 1).^2 ./ (38^2) + imag (x).^2 ./ (40^2)) <= 1) ... + cf_idx = ((((real (x) + 1).^2 ./ (38^2) + ... + imag (x).^2 ./ (40^2)) <= 1) ... & (! s_idx)) & (real (x) <= 35); - a_idx = ! s_idx & ! cf_idx; + a_idx = (! s_idx) & (! cf_idx); x_s = x(s_idx); x_cf = x(cf_idx); x_a = x(a_idx); - ## FIXME: The performance of these routines need improvement. - ## There are numerous temporary variables created, some of which could - ## probably be eliminated. - ## Series expansion ## Abramowitz, Stegun, "Handbook of Mathematical Functions", ## formula 5.1.11, p 229. @@ -124,72 +134,40 @@ res = -x_s; ssum = res; k = 1; - fflag = true (size (res)); - while (k < 1e3 && any (fflag)) - res_tmp = res(fflag); - x_s_tmp = x_s(fflag); - ssum_tmp = ssum(fflag); - res_tmp .*= (k * -x_s_tmp/((k + 1)^2)); - ssum_tmp += res_tmp; + todo = true (size (res)); + while (k < 1e3 && any (todo)) + res(todo) .*= (k * (- x_s(todo)) / ((k + 1) ^ 2)); + ssum(todo) += res(todo); k += 1; - res(fflag) = res_tmp; - ssum(fflag) = ssum_tmp; - x_s(fflag) = x_s_tmp; - fflag = abs (res) > tol*abs (ssum); + todo = (abs (res) > (tol * abs (ssum))); endwhile e1_s -= ssum; - ## Continued fraction, + ## Continued fraction expansion, ## Abramowitz, Stegun, "Handbook of Mathematical Functions", ## formula 5.1.22, p 229. - ## modified Lentz's algorithm, from "Numerical recipes in Fortran 77" p.165. - f_new = 2^-100 * ones (size (x_cf), class (x_cf)); - C_new = f_new; - D_new = zeros (size (x_cf), class (x_cf)); + ## Modified Lentz's algorithm, from "Numerical recipes in Fortran 77" p.165. + + e1_cf = exp (-x_cf) .* __expint__ (x_cf); + + ## Remove spurious imaginary part if needed (__expint__ works automatically + ## with complex values) + + if (isreal (x_cf) && x_cf >= 0) + e1_cf = real (e1_cf); + endif + + ## Asymptotic series, from N. Bleistein and R.A. Handelsman + ## "Asymptotic expansion of integrals", pages 1-4. + e1_a = exp (-x_a) ./ x_a; + ssum = res = ones (size (x_a), class (x_a)); k = 0; - fflag = true (size (x_cf)); - Delta = C_old = D_old = f_old = ones (size (x_cf), class (x_cf)); - while (k < 1e3 && any (fflag)) - x_cf_tmp = x_cf(fflag); - C_new_tmp = C_new(fflag); - D_new_tmp = D_new(fflag); - f_old = f_new(fflag); - C_old = C_new_tmp; - D_old = D_new_tmp; - b = x_cf_tmp*(mod (k,2) == 0) + (mod (k,2) == 1); - a = exp (-x_cf_tmp)*(k == 0) + ceil (k/2)*(k >= 1); - D_new_tmp = b + a.*D_old; - D_new_tmp = D_new_tmp.*(D_new_tmp != 0) + 2^-100*(D_new_tmp == 0); - C_new_tmp = b + a./C_old; - C_new_tmp = C_new_tmp.*(C_new_tmp != 0) + 2^-100*(C_new_tmp == 0); - D_new_tmp = 1 ./ D_new_tmp; - Delta(fflag) = C_new_tmp.*D_new_tmp; - x_cf(fflag) = x_cf_tmp; - f_new(fflag) = f_old.*Delta(fflag); - C_new(fflag) = C_new_tmp; - D_new(fflag) = D_new_tmp; - fflag = abs (Delta-1) > tol; + todo = true (size (x_a)); + while (k < 1e3 && any (todo)) + res(todo) ./= (- x_a(todo) / (k + 1)); + ssum(todo) += res(todo); k += 1; - endwhile - - e1_cf = f_new; - ## Asymptotic series, from Fikioris, Tastsoglou, Bakas, - ## "Selected Asymptotic Methods with Application to Magnetics and Antennas" - ## formula A.10 p 161. - e1_a = exp (-x_a) ./ x_a; - oldres = ssum = res = ones (size (x_a), class (x_a)); - k = 0; - fflag = true (size (x_a)); - while (k < 1e3 && any (fflag)) - res_tmp = res(fflag); - oldres_tmp = res_tmp; - x_a_tmp = x_a(fflag); - res_tmp ./= (-x_a_tmp / (k+1)); - ssum(fflag) += res_tmp; - k += 1; - res(fflag) = res_tmp; - oldres(fflag) = oldres_tmp; - fflag = abs (oldres) > abs (res); + todo = abs (x_a) > k; endwhile e1_a .*= ssum; @@ -197,9 +175,11 @@ E1(s_idx) = e1_s; E1(cf_idx) = e1_cf; E1(a_idx) = e1_a; + + ## Restore shape and sparsity of input E1 = reshape (E1, orig_sz); - if (sparse_x) - E1 = sparse (E1); # if input was sparse format, output should be too. + if (orig_sparse) + E1 = sparse (E1); endif endfunction @@ -249,7 +229,7 @@ %! 9.40470496160143e-05 - 2.44265223110736e-04i, ... %! 6.64487526601190e-09 - 7.87242868014498e-09i, ... %! 3.10273337426175e-13 - 2.28030229776792e-13i]; -%! assert (expint (X), y_exp, -1e-12); +%! assert (expint (X), y_exp, -1e-14); ## Exceptional values (-Inf, Inf, NaN, 0, 0.37250741078) %!test @@ -270,6 +250,10 @@ %!assert (class (expint (uint64 (1))), "double") %!assert (issparse (expint (sparse (1)))) +## Test on the correct Image set +%!assert (isreal (expint (linspace (0, 100)))) +%!assert (! isreal (expint (-1))) + ## Test input validation %!error expint () %!error expint (1,2)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/specfun/gammainc.m Mon Mar 19 13:50:10 2018 -0700 @@ -0,0 +1,564 @@ +## Copyright (C) 2016 Marco Caliari +## Copyright (C) 2016 Nir Krakauer +## Copyright (C) 2018 Stefan Schlögl +## Copyright (C) 2018 Michele Ginesi +## +## This file is part of Octave. +## +## Octave is free software: you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <https://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {} {} gammainc (@var{x}, @var{a}) +## @deftypefnx {} {} gammainc (@var{x}, @var{a}, @var{tail}) +## Compute the normalized incomplete gamma function. +## +## This is defined as +## @tex +## $$ +## \gamma (x, a) = {1 \over {\Gamma (a)}}\displaystyle{\int_0^x t^{a-1} e^{-t} dt} +## $$ +## @end tex +## @ifnottex +## +## @example +## @group +## x +## 1 / +## gammainc (x, a) = --------- | exp (-t) t^(a-1) dt +## gamma (a) / +## t=0 +## @end group +## @end example +## +## @end ifnottex +## with the limiting value of 1 as @var{x} approaches infinity. +## The standard notation is @math{P(a,x)}, e.g., @nospell{Abramowitz} and +## @nospell{Stegun} (6.5.1). +## +## If @var{a} is scalar, then @code{gammainc (@var{x}, @var{a})} is returned +## for each element of @var{x} and vice versa. +## +## If neither @var{x} nor @var{a} is scalar then the sizes of @var{x} and +## @var{a} must agree, and @code{gammainc} is applied element-by-element. +## The elements of @var{a} must be non-negative. +## +## By default, @var{tail} is @qcode{"lower"} and the incomplete gamma function +## integrated from 0 to @var{x} is computed. If @var{tail} is @qcode{"upper"} +## then the complementary function integrated from @var{x} to infinity is +## calculated. +## +## If @var{tail} is @qcode{"scaledlower"}, then the lower incomplete gamma +## function is multiplied by +## @tex +## $\Gamma(a+1)\exp(x)x^{-a}$. +## @end tex +## @ifnottex +## @math{gamma(a+1)*exp(x)/(x^a)}. +## @end ifnottex +## If @var{tail} is @qcode{"scaledupper"}, then the upper incomplete gamma +## function is divided by the same quantity. +## +## References: +## +## @nospell{M. Abramowitz and I. Stegun}, +## @cite{Handbook of mathematical functions}, +## @nospell{Dover publications, Inc.}, 1972. +## +## @nospell{W. Gautschi}, +## @cite{A computational procedure for incomplete gamma functions}, +## @nospell{ACM Trans. Math Software}, pp. 466--481, Vol 5, No. 4, 2012. +## +## @nospell{W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery}, +## @cite{Numerical Recipes in Fortran 77}, ch.@: 6.2, Vol 1, 1992. +## +## @seealso{gamma, gammaincinv, gammaln} +## @end deftypefn + +## P(a,x) = gamma(a,x)/Gamma(a), upper +## 1-P(a,x)=Q(a,x)=Gamma(a,x)/Gamma(a), lower + +function y = gammainc (x, a, tail = "lower") + + if (nargin < 2 || nargin > 3) + print_usage (); + endif + + [err, x, a] = common_size (x, a); + if (err > 0) + error ("gammainc: X and A must be of common size or scalars"); + endif + + if (iscomplex (x) || iscomplex (a)) + error ("gammainc: all inputs must be real"); + endif + + ## Remember original shape of data, but convert to column vector for calcs. + x_sz = size (x); + x = x(:); + a = a(:); + + if (any (a < 0)) + error ("gammainc: A must be non-negative"); + endif + + if (nargin == 3 + && ! any (strcmpi (tail, {"lower","upper","scaledlower","scaledupper"}))) + error ("gammainc: invalid value for TAIL"); + endif + tail = tolower (tail); + + ## If any of the arguments is single then the output should be as well. + if (strcmp (class (x), "single") || strcmp (class (a), "single")) + x = single (x); + a = single (a); + endif + + ## Convert to floating point if necessary + if (isinteger (x)) + x = double (x); + endif + if (isinteger (a)) + a = double (a); + endif + + ## Initialize output array + y = zeros (x_sz, class (x)); + + ## Different x, a combinations are handled by different subfunctions. + todo = true (size (x)); # Track which elements need to be calculated. + + ## Case 0: x == Inf, a == Inf + idx = (x == Inf) & (a == Inf); + if (any (idx)) + y(idx) = NaN; + todo(idx) = false; + endif + + ## Case 1: x == 0, a == 0. + idx = (x == 0) & (a == 0); + if (any (idx)) + y(idx) = gammainc_00 (tail); + todo(idx) = false; + endif + + ## Case 2: x == 0. + idx = todo & (x == 0); + if (any (idx)) + y(idx) = gammainc_x0 (tail); + todo(idx) = false; + endif + + ## Case 3: x = Inf + idx = todo & (x == Inf); + if (any (idx)) + y(idx) = gammainc_x_inf (tail); + todo(idx) = false; + endif + + ## Case 4: a = Inf + idx = todo & (a == Inf); + if (any (idx)) + y(idx) = gammainc_a_inf (tail); + todo(idx) = false; + endif + + ## Case 5: a == 0. + idx = todo & (a == 0); + if (any (idx)) + y(idx) = gammainc_a0 (x(idx), tail); + todo(idx) = false; + endif + + ## Case 6: a == 1. + idx = todo & (a == 1); + if (any (idx)) + y(idx) = gammainc_a1 (x(idx), tail); + todo(idx) = false; + endif + + ## Case 7: positive integer a; exp (x) and a! both under 1/eps. + idx = (todo + & (a == fix (a)) & (a > 1) & (a <= 18) & (x <= 36) & (abs (x) >= .1)); + if (any (idx)) + y(idx) = gammainc_an (x(idx), a(idx), tail); + todo(idx) = false; + endif + + ## For a < 2, x < 0, we increment a by 2 and use a recurrence formula after + ## the computations. + + flag_a_small = todo & (abs (a) > 0) & (abs (a) < 2) & (x < 0); + a(flag_a_small) += 2; + + flag_s = (((x + 0.25 < a) | (x < 0) | (a < 5)) & (x > -20)) | (abs (x) < 1); + + ## Case 8: x, a relatively small. + idx = todo & flag_s; + if (any (idx)) + y(idx) = gammainc_s (x(idx), a(idx), tail); + todo(idx) = false; + endif + + ## Case 9: x positive and large relative to a. + idx = todo; + if (any (idx)) + y(idx) = gammainc_l (x(idx), a(idx), tail); + todo(idx) = false; + endif + + if (any (flag_a_small)) + if (strcmp (tail, "lower")) + y(flag_a_small) += D (x(flag_a_small), a(flag_a_small) - 1) + ... + D (x(flag_a_small), a(flag_a_small) - 2); + elseif (strcmp (tail, "upper")) + y(flag_a_small) -= D (x(flag_a_small), a(flag_a_small) - 1) + ... + D (x(flag_a_small), a(flag_a_small) - 2); + elseif (strcmp (tail, "scaledlower")) + y(flag_a_small) = y(flag_a_small) .* (x(flag_a_small) .^ 2) ./ ... + (a(flag_a_small) .* (a(flag_a_small) - 1)) + (x(flag_a_small) ./ ... + (a(flag_a_small) - 1)) + 1; + elseif (strcmp (tail, "scaledupper")) + y(flag_a_small) = y(flag_a_small) .* (x(flag_a_small) .^ 2) ./ ... + (a(flag_a_small) .* (a(flag_a_small) - 1)) - (x(flag_a_small) ./ ... + (a(flag_a_small) - 1)) - 1; + endif + endif + +endfunction + +## Subfunctions to handle each case: + +## x == 0, a == 0. +function y = gammainc_00 (tail) + if ((strcmp (tail, "upper")) || (strcmp (tail, "scaledupper"))) + y = 0; + else + y = 1; + endif +endfunction + +## x == 0. +function y = gammainc_x0 (tail) + if (strcmp (tail, "lower")) + y = 0; + elseif (strcmp (tail, "upper") || strcmp (tail, "scaledlower")) + y = 1; + else + y = Inf; + endif +endfunction + +## x == Inf. +function y = gammainc_x_inf (tail) + if (strcmp (tail, "lower")) + y = 1; + elseif (strcmp (tail, "upper") || strcmp (tail, "scaledupper")) + y = 0; + else + y = Inf; + endif +endfunction + +## a == Inf. +function y = gammainc_a_inf (tail) + if (strcmp (tail, "lower")) + y = 0; + elseif (strcmp (tail, "upper") || strcmp (tail, "scaledlower")) + y = 1; + else + y = Inf; + endif +endfunction + +## a == 0. +function y = gammainc_a0 (x, tail) + if (strcmp (tail, "lower")) + y = 1; + elseif (strcmp (tail, "scaledlower")) + y = exp (x); + else + y = 0; + endif +endfunction + +## a == 1. +function y = gammainc_a1 (x, tail) + if (strcmp (tail, "lower")) + y = 1 - exp (-x); + elseif (strcmp (tail, "upper")) + y = exp (-x); + elseif (strcmp (tail, "scaledlower")) + if (abs (x) < 1/2) + y = expm1 (x) ./ x; + else + y = (exp (x) - 1) ./ x; + endif + else + y = 1 ./ x; + endif +endfunction + +## positive integer a; exp (x) and a! both under 1/eps +## uses closed-form expressions for nonnegative integer a +## -- http://mathworld.wolfram.com/IncompleteGammaFunction.html. +function y = gammainc_an (x, a, tail) + y = t = ones (size (x), class (x)); + i = 1; + while (any (a(:) > i)) + jj = (a > i); + t(jj) .*= (x(jj) / i); + y(jj) += t(jj); + i++; + endwhile + if (strcmp (tail, "lower")) + y = 1 - exp (-x) .* y; + elseif (strcmp (tail, "upper")) + y .*= exp (-x); + elseif (strcmp (tail, "scaledlower")) + y = (1 - exp (-x) .* y) ./ D(x, a); + elseif (strcmp (tail, "scaledupper")) + y .*= exp (-x) ./ D(x, a); + endif +endfunction + +## x + 0.25 < a | x < 0 | x not real | abs(x) < 1 | a < 5. +## Numerical Recipes in Fortran 77 (6.2.5) +## series +function y = gammainc_s (x, a, tail) + if (strcmp (tail, "scaledlower") || strcmp (tail, "scaledupper")) + y = ones (size (x), class (x)); + term = x ./ (a + 1); + else + ## Of course it is possible to scale at the end, but some tests fail. + ## And try gammainc (1,1000), it take 0 iterations if you scale now. + y = D (x,a); + term = y .* x ./ (a + 1); + endif + n = 1; + while (any (abs (term(:)) > (abs (y(:)) * eps))) + ## y can be zero from the beginning (gammainc (1,1000)) + jj = abs (term) > abs (y) * eps; + n += 1; + y(jj) += term(jj); + term(jj) .*= x(jj) ./ (a(jj) + n); + endwhile + if (strcmp (tail, "upper")) + y = 1 - y; + elseif (strcmp (tail, "scaledupper")) + y = 1 ./ D (x,a) - y; + endif +endfunction + +## x positive and large relative to a +## NRF77 (6.2.7) +## Gamma (a,x)/Gamma (a) +## Lentz's algorithm +## __gammainc__ in libinterp/corefcn/__gammainc__.cc +function y = gammainc_l (x, a, tail) + y = __gammainc__ (x, a); + if (strcmp (tail, "lower")) + y = 1 - y .* D (x, a); + elseif (strcmp (tail, "upper")) + y .*= D (x, a); + elseif (strcmp (tail, "scaledlower")) + y = 1 ./ D (x, a) - y; + endif +endfunction + +## Compute exp(-x)*x^a/Gamma(a+1) in a stable way for x and a large. +## +## L. Knusel, Computation of the Chi-square and Poisson distribution, +## SIAM J. Sci. Stat. Comput., 7(3), 1986 +## which quotes Section 5, Abramowitz&Stegun 6.1.40, 6.1.41. +function y = D (x, a) + athresh = 10; # FIXME: can this be better tuned? + y = zeros (size (x), class (x)); + + todo = true (size (x)); + todo(x == 0) = false; + + ii = todo & (x > 0) & (a > athresh) & (a >= x); + if (any (ii)) + lnGa = log (2 * pi * a(ii)) / 2 + 1 ./ (12 * a(ii)) - ... + 1 ./ (360 * a(ii) .^ 3) + 1 ./ (1260 * a(ii) .^ 5) - ... + 1 ./ (1680 * a(ii) .^ 7) + 1 ./ (1188 * a(ii) .^ 9)- ... + 691 ./ (87360 * a(ii) .^ 11) + 1 ./ (156 * a(ii) .^ 13) - ... + 3617 ./ (122400 * a(ii) .^ 15) + ... + 43867 ./ (244188 * a(ii) .^ 17) - 174611 ./ (125400 * a(ii) .^ 19); + lns = log1p ((a(ii) - x(ii)) ./ x(ii)); + y(ii) = exp ((a(ii) - x(ii)) - a(ii) .* lns - lnGa); + todo(ii) = false; + endif + + ii = todo & (x > 0) & (a > athresh) & (a < x); + if (any (ii)) + lnGa = log (2 * pi * a(ii)) / 2 + 1 ./ (12 * a(ii)) - ... + 1 ./ (360 * a(ii) .^ 3) + 1 ./ (1260 * a(ii) .^ 5) - ... + 1 ./ (1680 * a(ii) .^ 7) + 1 ./ (1188 * a(ii) .^ 9)- ... + 691 ./ (87360 * a(ii) .^ 11) + 1 ./ (156 * a(ii) .^ 13) - ... + 3617 ./ (122400 * a(ii) .^ 15) + ... + 43867 ./ (244188 * a(ii) .^ 17) - 174611 ./ (125400 * a(ii) .^ 19); + lns = -log1p ((x(ii) - a(ii)) ./ a(ii)); + y(ii) = exp ((a(ii) - x(ii)) - a(ii) .* lns - lnGa); + todo(ii) = false; + endif + + ii = todo & ((x <= 0) | (a <= athresh)); + if (any (ii)) # standard formula for a not so large. + y(ii) = exp (a(ii) .* log (x(ii)) - x(ii) - gammaln (a(ii) + 1)); + todo(ii) = false; + endif + + ii = (x < 0) & (a == fix (a)); + if (any(ii)) # remove spurious imaginary part. + y(ii) = real (y(ii)); + endif + +endfunction + + +## Test: case 1,2,5 +%!assert (gammainc ([0, 0, 1], [0, 1, 0]), [1, 0, 1]) +%!assert (gammainc ([0, 0, 1], [0, 1, 0], "upper"), [0, 1, 0]) +%!assert (gammainc ([0, 0, 1], [0, 1, 0], "scaledlower"), [1, 1, exp(1)]) +%!assert (gammainc ([0, 0, 1], [0, 1, 0], "scaledupper"), [0, Inf, 0]) + +## Test: case 3,4 +%!assert (gammainc ([2, Inf], [Inf, 2]), [0, 1]) +%!assert (gammainc ([2, Inf], [Inf, 2], "upper"), [1, 0]) +%!assert (gammainc ([2, Inf], [Inf, 2], "scaledlower"), [1, Inf]) +%!assert (gammainc ([2, Inf], [Inf, 2], "scaledupper"), [Inf, 0]) + +## Test: case 5 +## Matlab fails for this test +%!assert (gammainc (-100,1,"upper"), exp (100), -eps) + +## Test: case 6 +%!assert (gammainc ([1, 2, 3], 1), 1 - exp (-[1, 2, 3])) +%!assert (gammainc ([1, 2, 3], 1, "upper"), exp (- [1, 2, 3])) +%!assert (gammainc ([1, 2, 3], 1, "scaledlower"), ... +%! (exp ([1, 2, 3]) - 1) ./ [1, 2, 3]) +%!assert (gammainc ([1, 2, 3], 1, "scaledupper"), 1 ./ [1, 2, 3]) + +## Test: case 7 +%!assert (gammainc (2, 2, "lower"), 0.593994150290162, -2e-15) +%!assert (gammainc (2, 2, "upper"), 0.406005849709838, -2e-15) +%!assert (gammainc (2, 2, "scaledlower"), 2.194528049465325, -2e-15) +%!assert (gammainc (2, 2, "scaledupper"), 1.500000000000000, -2e-15) +%!assert (gammainc ([3 2 36],[2 3 18], "upper"), ... +%! [4/exp(3) 5*exp(-2) (4369755579265807723 / 2977975)/exp(36)]) +%!assert (gammainc (10, 10), 1 - (5719087 / 567) * exp (-10), -eps) +%!assert (gammainc (10, 10, "upper"), (5719087 / 567) * exp (-10), -eps) + +## Test: case 8 +%!assert (gammainc (-10, 10), 3.112658265341493126871617e7, -2*eps) +## Matlab fails this next one%! %! +%!assert (isreal (gammainc (-10, 10)), true) +%!assert (gammainc (-10, 10.1, "upper"), ... +%! -2.9582761911890713293e7-1i * 9.612022339061679758e6, -30*eps) +%!assert (gammainc (-10, 10, "upper"), -3.112658165341493126871616e7, ... +%! -2*eps) +%!assert (gammainc (-10, 10, "scaledlower"), 0.5128019364747265, -1e-14); +%!assert (gammainc (-10, 10, "scaledupper"), -0.5128019200000000, -1e-14); +%!assert (gammainc (200, 201, "upper"), 0.518794309678684497, -2 * eps); +%!assert (gammainc (200, 201, "scaledupper"), +%! 18.4904360746560462660798514, -eps) +## Here we are very good (no D (x,a)) involved +%!assert (gammainc(1000, 1000.5, "scaledlower"), 39.48467539583672271, -2*eps) +%!assert (gammainc (709, 1000, "upper"), 0.99999999999999999999999954358, -eps) + +## Test: case 9 +%!test <47800> +%! assert (gammainc (60, 6, "upper"), 6.18022358081160257327264261e-20, +%! -10*eps); +## Matlab is better here than Octave +%!assert (gammainc (751, 750, "upper"), 0.4805914320558831327179457887, -12*eps) +%!assert (gammainc (200, 200, "upper"), 0.49059658199276367497217454, -4*eps) +%!assert (gammainc (200, 200), 0.509403418007236325027825459574527043, -3*eps) +%!assert (gammainc (200, 200, "scaledupper"), 17.3984438553791505135122900, +%! -eps) +%!assert (gammainc (200, 200, "scaledlower"), 18.065406676779221643065, -6*eps) +%!assert (gammainc (201, 200, "upper"), 0.46249244908276709524913736667, -7*eps) + +## Test small argument +%!assert (gammainc ([1e-05, 1e-07,1e-10,1e-14], 0.1), ... +%! [0.33239840504050, 0.20972940370977, 0.10511370061022, ... +%! 0.041846517936723], 1e-13); + +%!assert (gammainc ([1e-05, 1e-07,1e-10,1e-14], 0.2), ... +%! [0.10891226058559, 0.043358823442178, 0.010891244210402, ... +%! 0.0017261458806785], 1e-13); + +%!test +%!assert (gammainc ([1e-02, 1e-03, 1e-5, 1e-9, 1e-14], 0.9), ... +%! [0.016401189184068, 0.0020735998660840, 0.000032879756964708, ... +%! 8.2590606569241e-9, 2.6117443021738e-13], -1e-12); + +%!test +%!assert (gammainc ([1e-02, 1e-03, 1e-5, 1e-9, 1e-14], 2), ... +%! [0.0000496679133402659, 4.99666791633340e-7, 4.99996666679167e-11, ... +%! 4.99999999666667e-19, 4.99999999999997e-29], -1e-12); + +## FIXME: should this be tagged with a bug report number? +%!xtest +%! assert (gammainc (-20, 1.1, "upper"), ... +%! 6.50986687074979e8 + 2.11518396291149e8*i, -1e-13); + +## Test conservation of the class (five tests for each subroutine). +%!assert (class (gammainc (0, 1)) == "double") +%!assert (class (gammainc (single (0), 1)) == "single") +%!assert (class (gammainc (int8 (0), 1)) == "double") +%!assert (class (gammainc (0, single (1))) == "single") +%!assert (class (gammainc (0, int8 (1))) == "double") +%!assert (class (gammainc (1, 0)) == "double") +%!assert (class (gammainc (single (1), 0)) == "single") +%!assert (class (gammainc (int8 (1), 0)) == "double") +%!assert (class (gammainc (1, single (0))) == "single") +%!assert (class (gammainc (1, int8 (0))) == "double") +%!assert (class (gammainc (1, 1)) == "double") +%!assert (class (gammainc (single (1), 1)) == "single") +%!assert (class (gammainc (int8 (1), 1)) == "double") +%!assert (class (gammainc (1, single (1))) == "single") +%!assert (class (gammainc (1, int8(1))) == "double") +%!assert (class (gammainc (1, 2)) == "double") +%!assert (class (gammainc (single (1), 2)) == "single") +%!assert (class (gammainc (int8 (1), 2)) == "double") +%!assert (class (gammainc (1, single (2))) == "single") +%!assert (class (gammainc (1, int8 (2))) == "double") +%!assert (class (gammainc (-1, 0.5)) == "double") +%!assert (class (gammainc (single (-1), 0.5)) == "single") +%!assert (class (gammainc (int8 (-1), 0.5)) == "double") +%!assert (class (gammainc (-1, single (0.5))) == "single") +%!assert (class (gammainc (-1, int8 (0.5))) == "double") +%!assert (class (gammainc (1, 0.5)) == "double") +%!assert (class (gammainc (single (1), 0.5)) == "single") +%!assert (class (gammainc (int8 (1), 0.5)) == "double") +%!assert (class (gammainc (1, single (0.5))) == "single") +%!assert (class (gammainc (1, int8 (0.5))) == "double") + +## Test input validation +%!error gammainc () +%!error gammainc (1) +%!error gammainc (1,2,3,4) +%!error <must be of common size or scalars> gammainc ([0, 0],[0; 0]) +%!error <must be of common size or scalars> gammainc ([1 2 3], [1 2]) +%!error <all inputs must be real> gammainc (2+i, 1) +%!error <all inputs must be real> gammainc (1, 2+i) +%!error <A must be non-negative> gammainc (1, [0, -1, 1]) +%!error <A must be non-negative> +%! a = ones (2,2,2); +%! a(1,1,2) = -1; +%! gammainc (1, a); +%!error <invalid value for TAIL> gammainc (1,2, "foobar")
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/specfun/gammaincinv.m Mon Mar 19 13:50:10 2018 -0700 @@ -0,0 +1,315 @@ +## Copyright (C) 2017 Michele Ginesi +## +## This file is part of Octave. +## +## Octave is free software: you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <https://www.gnu.org/licenses/>. + +## Author: Michele Ginesi <michele.ginesi@gmail.com> + +## -*- texinfo -*- +## @deftypefn {} {} gammaincinv (@var{y}, @var{a}) +## @deftypefnx {} {} gammaincinv (@var{y}, @var{a}, @var{tail}) +## Compute the inverse of the normalized incomplete gamma function. +## +## The normalized incomplete gamma function is defined as +## @tex +## $$ +## \gamma (x, a) = {1 \over {\Gamma (a)}}\displaystyle{\int_0^x t^{a-1} e^{-t} dt} +## $$ +## @end tex +## @ifnottex +## +## @example +## @group +## x +## 1 / +## gammainc (x, a) = --------- | exp (-t) t^(a-1) dt +## gamma (a) / +## t=0 +## @end group +## @end example +## +## @end ifnottex +## +## and @code{gammaincinv (gammainc (@var{x}, @var{a}), @var{a}) = @var{x}} +## for each non-negative value of @var{x}. If @var{a} is scalar then +## @code{gammaincinv (@var{y}, @var{a})} is returned for each element of +## @var{y} and vice versa. +## +## If neither @var{y} nor @var{a} is scalar then the sizes of @var{y} and +## @var{a} must agree, and @code{gammaincinv} is applied element-by-element. +## The variable @var{y} must be in the interval @math{[0,1]} while @var{a} must +## be real and positive. +## +## By default, @var{tail} is @qcode{"lower"} and the inverse of the incomplete +## gamma function integrated from 0 to @var{x} is computed. If @var{tail} is +## @qcode{"upper"}, then the complementary function integrated from @var{x} to +## infinity is inverted. +## +## The function is computed with Newton's method by solving +## @tex +## $$ +## y - \gamma (x, a) = 0 +## $$ +## @end tex +## @ifnottex +## +## @example +## @var{y} - gammainc (@var{x}, @var{a}) = 0 +## @end example +## +## @end ifnottex +## +## Reference: @nospell{A. Gil, J. Segura, and N. M. Temme}, @cite{Efficient and +## accurate algorithms for the computation and inversion of the incomplete +## gamma function ratios}, @nospell{SIAM J. Sci. Computing}, pp. A2965--A2981, +## Vol 34, 2012. +## +## @seealso{gammainc, gamma, gammaln} +## @end deftypefn + +function x = gammaincinv (y, a, tail = "lower") + + if (nargin < 2 || nargin > 3) + print_usage (); + endif + + [err, y, a] = common_size (y, a); + if (err > 0) + error ("gammaincinv: Y and A must be of common size or scalars"); + endif + + if (iscomplex (y) || iscomplex (a)) + error ("gammaincinv: all inputs must be real"); + endif + + ## Remember original shape of data, but convert to column vector for calcs. + orig_sz = size (y); + y = y(:); + a = a(:); + + if (any ((y < 0) | (y > 1))) + error ("gammaincinv: Y must be in the range [0, 1]"); + endif + + if (any (a <= 0)) + error ("gammaincinv: A must be strictly positive"); + endif + + ## If any of the arguments is single then the output should be as well. + if (strcmp (class (y), "single") || strcmp (class (a), "single")) + y = single (y); + a = single (a); + endif + + ## Convert to floating point if necessary + if (isinteger (y)) + y = double (y); + endif + if (isinteger (a)) + a = double (a); + endif + + ## Initialize output array + x = zeros (size (y), class (y)); + + maxit = 20; + tol = eps (class (y)); + + ## Special cases, a = 1 or y = 0, 1. + + if (strcmpi (tail, "lower")) + x(a == 1) = - log (1 - y(a == 1)); + x(y == 0) = 0; + x(y == 1) = Inf; + p = y; + q = 1 - p; + elseif (strcmpi (tail, "upper")) + x(a == 1) = - log (y(a == 1)); + x(y == 0) = Inf; + x(y == 1) = 0; + q = y; + p = 1 - q; + else + error ("gammaincinv: invalid value for TAIL") + endif + + todo = (a != 1) & (y != 0) & (y != 1); + + ## Case 1: p small. + + i_flag_1 = todo & (p < ((0.2 * (1 + a)) .^ a) ./ gamma (1 + a)); + + aa = a(i_flag_1); + pp = p(i_flag_1); + + ## Initial guess. + + r = (pp .* gamma (1 + aa)) .^ (1 ./ aa); + + c2 = 1 ./ (aa + 1); + c3 = (3 * aa + 5) ./ (2 * (aa + 1) .^2 .* (aa + 2)); + c4 = (8 * aa .^ 2 + 33 * aa + 31) ./ (3 * (aa + 1) .^ 3 .* (aa + 2) .* ... + (aa + 3)); + c5 = (125 * aa .^ 4 + 1179 * aa .^ 3 + 3971 * aa.^2 + 5661 * aa + 2888) ... + ./ (24 * (1 + aa) .^4 .* (aa + 2) .^ 2 .* (aa + 3) .* (aa + 4)); + + ## FIXME: Would polyval() be better here for more accuracy? + x0 = r + c2 .* r .^ 2 + c3 .* r .^ 3 + c4 .* r .^4 + c5 .* r .^ 5; + + ## For this case we invert the lower version. + + F = @(p, a, x) p - gammainc (x, a, "lower"); + JF = @(a, x) - exp (-gammaln (a) - x + (a - 1) .* log (x)); + x(i_flag_1) = newton_method (F, JF, pp, aa, x0, tol, maxit); + + todo(i_flag_1) = false; + + ## Case 2: q small. + + i_flag_2 = (q < exp (- 0.5 * a) ./ gamma (1 + a)) & (a > 0) & (a < 10); + i_flag_2 &= todo; + + aa = a(i_flag_2); + qq = q(i_flag_2); + + ## Initial guess. + + x0 = (-log (qq) - gammaln (aa)); + + ## For this case, we invert the upper version. + + F = @(q, a, x) q - gammainc (x, a, "upper"); + JF = @(a, x) exp (- gammaln (a) - x) .* x .^ (a - 1); + x(i_flag_2) = newton_method (F, JF, qq, aa, x0, tol, maxit); + + todo(i_flag_2) = false; + + ## Case 3: a small. + + i_flag_3 = todo & ((a > 0) & (a < 1)); + + aa = a(i_flag_3); + pp = p(i_flag_3); + + ## Initial guess + + xl = (pp .* gamma (aa + 1)) .^ (1 ./ aa); + x0 = xl; + + ## For this case, we invert the lower version. + + F = @(p, a, x) p - gammainc (x, a, "lower"); + JF = @(a, x) - exp (-gammaln (a) - x) .* x .^ (a - 1); + x(i_flag_3) = newton_method (F, JF, pp, aa, x0, tol, maxit); + + todo(i_flag_3) = false; + + ## Case 4: a large. + + i_flag_4 = todo; + aa = a(i_flag_4); + qq = q(i_flag_4); + + ## Initial guess + + d = 1 ./ (9 * aa); + t = 1 - d + sqrt (2) * erfcinv (2 * qq) .* sqrt (d); + x0 = aa .* (t .^ 3); + + ## For this case, we invert the upper version. + + F = @(q, a, x) q - gammainc (x, a, "upper"); + JF = @(a, x) exp (- gammaln (a) - x + (a - 1) .* log (x)); + x(i_flag_4) = newton_method (F, JF, qq, aa, x0, tol, maxit); + + ## Restore original shape + x = reshape (x, orig_sz); + +endfunction + +## Subfunction: Newton's Method +function x = newton_method (F, JF, y, a, x0, tol, maxit); + l = numel (y); + res = -F (y, a, x0) ./ JF (a, x0); + todo = (abs (res) >= tol * abs (x0)); + x = x0; + it = 0; + while (any (todo) && (it++ < maxit)) + x(todo) += res(todo); + res(todo) = -F (y(todo), a(todo), x(todo)) ./ JF (a(todo), x(todo)); + todo = (abs (res) >= tol * abs (x)); + endwhile + x += res; +endfunction + + +%!test +%! x = [1e-10, 1e-09, 1e-08, 1e-07]; +%! a = [2, 3, 4]; +%! [x, a] = ndgrid (x, a); +%! xx = gammainc (gammaincinv (x, a), a); +%! assert (xx, x, -3e-14); + +%!test +%! x = [1e-10, 1e-09, 1e-08, 1e-07]; +%! a = [2, 3, 4]; +%! [x, a] = ndgrid (x, a); +%! xx = gammainc (gammaincinv (x, a, "upper"), a, "upper"); +%! assert (xx, x, -3e-14); + +%!test +%! x = linspace (0, 1)'; +%! a = [linspace(0.1, 1, 10), 2:5]; +%! [x, a] = ndgrid (x, a); +%! xx = gammainc (gammaincinv (x, a), a); +%! assert (xx, x, -1e-13); + +%!test +%! x = linspace (0, 1)'; +%! a = [linspace(0.1, 1, 10), 2:5]; +%! [x, a] = ndgrid (x, a); +%! xx = gammainc (gammaincinv (x, a, "upper"), a, "upper"); +%! assert (xx, x, -1e-13); + +## Test the conservation of the input class +%!assert (class (gammaincinv (0.5, 1)), "double") +%!assert (class (gammaincinv (single (0.5), 1)), "single") +%!assert (class (gammaincinv (0.5, single (1))), "single") +%!assert (class (gammaincinv (int8 (0), 1)), "double") +%!assert (class (gammaincinv (0.5, int8 (1))), "double") +%!assert (class (gammaincinv (int8 (0), single (1))), "single") +%!assert (class (gammaincinv (single (0.5), int8 (1))), "single") + +## Test input validation +%!error gammaincinv () +%!error gammaincinv (1) +%!error gammaincinv (1, 2, 3, 4) +%!error <must be of common size or scalars> +%! gammaincinv (ones (2,2), ones (1,2), 1); +%!error <all inputs must be real> gammaincinv (0.5i, 1) +%!error <all inputs must be real> gammaincinv (0, 1i) +%!error <Y must be in the range \[0, 1\]> gammaincinv (-0.1,1) +%!error <Y must be in the range \[0, 1\]> gammaincinv (1.1,1) +%!error <Y must be in the range \[0, 1\]> +%! y = ones (1, 1, 2); +%! y(1,1,2) = -1; +%! gammaincinv (y,1); +%!error <A must be strictly positive> gammaincinv (0.5, 0) +%!error <A must be strictly positive> +%! a = ones (1, 1, 2); +%! a(1,1,2) = 0; +%! gammaincinv (1,a,1); +%!error <invalid value for TAIL> gammaincinv (1,2, "foobar")
--- a/scripts/specfun/module.mk Sun Mar 18 18:43:04 2018 -0700 +++ b/scripts/specfun/module.mk Mon Mar 19 13:50:10 2018 -0700 @@ -2,11 +2,16 @@ %canon_reldir%_FCN_FILES = \ %reldir%/beta.m \ + %reldir%/betainc.m \ + %reldir%/betaincinv.m \ %reldir%/betaln.m \ + %reldir%/cosint.m \ %reldir%/ellipke.m \ %reldir%/expint.m \ %reldir%/factor.m \ %reldir%/factorial.m \ + %reldir%/gammainc.m \ + %reldir%/gammaincinv.m \ %reldir%/isprime.m \ %reldir%/lcm.m \ %reldir%/legendre.m \ @@ -17,7 +22,8 @@ %reldir%/primes.m \ %reldir%/reallog.m \ %reldir%/realpow.m \ - %reldir%/realsqrt.m + %reldir%/realsqrt.m \ + %reldir%/sinint.m %canon_reldir%dir = $(fcnfiledir)/specfun
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/specfun/sinint.m Mon Mar 19 13:50:10 2018 -0700 @@ -0,0 +1,207 @@ +## Copyright (C) 2017 Michele Ginesi +## +## This file is part of Octave. +## +## Octave is free software: you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation, either version 3 of the License, or +## (at your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <https://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {} {} sinint (@var{x}) +## Compute the sine integral function: +## @tex +## $$ +## {\rm Si} (x) = \int_0^x {\sin (t) \over t} dt +## $$ +## @end tex +## @ifnottex +## +## @example +## @group +## x +## / +## Si (x) = | sin (t) / t dt +## / +## 0 +## @end group +## @end example +## +## @end ifnottex +## +## Reference: @nospell{M. Abramowitz and I.A. Stegun}, +## @cite{Handbook of Mathematical Functions}, 1964. +## +## @seealso{cosint, expint, sin} +## @end deftypefn + +function y = sinint (x) + + if (nargin != 1) + print_usage (); + endif + + if (! isnumeric (x)) + error ("sinint: X must be numeric"); + endif + + ## Convert to floating point if necessary + if (isinteger (x)) + x = double (x); + endif + + ## Convert to column vector + orig_sz = size (x); + x = x(:); + if (iscomplex (x)) + ## Work around reshape which narrows to real (bug #52953) + x = complex (real (x)(:), imag (x)(:)); + else + x = x(:); + end + + ## Initialize the result + y = zeros (size (x), class (x)); + tol = eps (class (x)); + + todo = true (size (x)); + + ## Special values + y(x == 0) = x(x == 0); # correctly signed zero + y(x == Inf) = pi / 2; + y(x == - Inf) = - pi / 2; + + todo = ((todo) & (x != 0) & (x != Inf) & (x != - Inf)); + + ## For values large in modulus we use the relation with expint + + flag_large = abs (x) > 2; + xx = x(flag_large & todo); + ii_neg = (real (xx) < 0); + xx(ii_neg) *= -1; + ii_conj = (real (xx) == 0) & (imag (xx) < 0); + xx(ii_conj) = conj (xx(ii_conj)); + yy = -0.5i * (expint (1i * xx) - expint (-1i * xx)) + pi / 2; + yy(ii_neg) *= -1; + yy(ii_conj) = conj (yy(ii_conj)); + y(todo & flag_large) = yy; + + ## For values small in modulus we use the series expansion + + todo = (todo) & (! flag_large); + xx = x(todo); + ssum = xx; # First term of the series expansion + yy = ssum; + flag_sum = true (nnz (todo), 1); + it = 0; + maxit = 300; + while (any (flag_sum) && (it < maxit)) + ssum .*= - xx .^ 2 * (2 * it + 1) / ((2 * it + 3) ^ 2 * (2 * it + 2)); + yy(flag_sum) += ssum (flag_sum); + flag_sum = (abs (ssum) >= tol); + it++; + endwhile + + y(todo) = yy; + + y = reshape (y, orig_sz); + +endfunction + + +%!assert (sinint (1.1), 1.02868521867373, -5e-15) + +%!test +%! x = [2, 3, pi; exp(1), 5, 6]; +%! A = sinint (x); +%! B = [1.60541297680269, 1.84865252799947, 1.85193705198247e+00; ... +%! 1.82104026914757, 1.54993124494467, 1.42468755128051e+00]; +%! assert (A, B, -5e-15); + +## Test exceptional values +%!assert (sinint (0), 0) +%!assert (signbit (sinint (-0))) +%!assert (sinint (Inf), pi/2) +%!assert (sinint (-Inf), -pi/2) +%!assert (isnan (sinint (NaN))) + +## Check single data type is preserved +%!assert (class (sinint (single (1))), "single") + +## Tests against Maple +%!assert (sinint (1) , 0.9460830703671830149414, -2*eps) +%!assert (sinint (-1) , -0.9460830703671830149414, -2*eps) +%!assert (sinint (pi) , 1.851937051982466170361, -2*eps) +%!assert (sinint (-pi), -1.851937051982466170361, -2*eps) +%!assert (sinint (300), 1.5708810882137495193, -2*eps) +%!assert (sinint (1e4), 1.5708915453859619157, -2*eps) +%!assert (sinint (20i), 1.2807826332028294459e7*1i, -2*eps) + +%!test +%! x = (0:4)'; +%! y_ex = [0 +%! 0.946083070367183015 +%! 1.60541297680269485 +%! 1.84865252799946826 +%! 1.75820313894905306]; +%! assert (sinint (x), y_ex, -4*eps); + +%!test +%! x = -(0:4)'; +%! y_ex = - [0 +%! 0.946083070367183015 +%! 1.60541297680269485 +%! 1.84865252799946826 +%! 1.75820313894905306]; +%! assert (sinint (x), y_ex, -4*eps); + +%!test +%! x = 1i * (0:4).'; +%! y_ex = [0 +%! 1.05725087537572851*I +%! 2.50156743335497564*I +%! 4.97344047585980680*I +%! 9.81732691123303446*I]; +%! assert (sinint (x), y_ex, -4*eps); + +%!test +%! x = - 1i * (0:4).'; +%! y_ex = - [0 +%! 1.05725087537572851*I +%! 2.50156743335497564*I +%! 4.97344047585980680*I +%! 9.81732691123303446*I]; +%! assert (sinint (x), y_ex, -4*eps); + +%!test +%! % maple: +%! % > A := [1+2*I, -2 + 5*I, 100, 10*I, -1e-4 + 1e-6*I, -20 + I]; +%! % > for a in A do evalf(Si(a)) end do; +%! x = [1+2i; -2+5i; 100; 10i; -1e-4 + 1e-6*1i; -20-1i]; +%! A = [ 1.6782404878293681180 + 2.0396845546022061045*1i +%! -18.154174221650281533 + 1.6146414539230479060*1i +%! 1.5622254668890562934 +%! 1246.1144901994233444*1i +%! -0.000099999999944461111128 + 0.99999999833338888972e-6*1i +%! -1.5386156269726011209 - 0.053969388020443786229*1i ]; +%! B = sinint (x); +%! assert (A, B, -3*eps) +%! B = sinint (single (x)); +%! assert (A, B, -3*eps ("single")) + +## FIXME: Need a test for bug #52953 +%#!test <*52953> + +## Test input validation +%!error sinint () +%!error sinint (1,2) +%!error <X must be numeric> sinint ("1")
--- a/scripts/statistics/corrcoef.m Sun Mar 18 18:43:04 2018 -0700 +++ b/scripts/statistics/corrcoef.m Mon Mar 19 13:50:10 2018 -0700 @@ -292,3 +292,4 @@ %!error <"alpha" must be a number between 0 and 1> corrcoef (1,2, "alpha", 2) %!error <"rows" must be "all"...> corrcoef (1,2, "rows", "foobar") %!error <Unknown option "foobar"> corrcoef (1,2, "foobar", 1) +%!error <Unknown option "foobar"> corrcoef (1,2, "foobar", 1)