Mercurial > octave
changeset 24927:c280560d9c96 stable
Overhaul special functions modified by GSOC2018 project.
* NEWS: Add note about new functions added. Add note explaining the changes
done to existing functions.
* __betainc__.cc: Renamed from __betainc_lentz__.cc. Use standard GPL v3
copyright block. Add missing #include "dNDArray.h". Add one-line texinfo
documentation for internal function. Remove fourth input argument to function.
Use is_single_type() to decide whether to operate with FloatNDArray or NDArray.
Delete temporary variables x_arg_s, a_arg_s, b_arg_s used in input validation.
Rename len_x, len_a, len_b to numel_[xab] for clarity. Remove input validation
that numel match (internal function, no need). Rename function output to
retval according to Octave coding conventions. Use more spacing and newlines
for readability of code.
* __expint__.cc: Renamed from __expint_lentz__.cc. Use standard GPL v3
copyright block. Cull list of #includes to just the necessary ones. Use
Complex and FloatComplex typedefs defined by Octave. Add one-line texinfo
documentation for internal function. Remove second input argument to function.
Use is_single_type() to decide whether to operate with FloatComplexNDArray or
ComplexNDArray. Delete temporary variables x_arg_s used in input validation.
Rename len_x to numel_x for clarity. Use constructor with dim_vector and
scalar value rather than fill() after creating array. Rename function output
to retval according to Octave coding conventions. Use more spacing and
newlines for readability of code.
* __gammainc__.cc: Renamed from __gammainc_lentz__.cc. Use standard GPL v3
copyright block. Add one-line texinfo documentation for internal function.
Remove third input argument to function. Remove input validation that numel
match (internal function, no need). Use is_single_type() to decide whether to
operate with FloatNDArray or NDArray. Delete temporary variables x_arg, a_arg
used in input validation. Use constructor with dim_vector and scalar value
rather than fill() after creating array. Rename function output to retval
according to Octave coding conventions. Use more spacing and newlines for
readability of code.
* __betainc_lentz__.cc, __expint_lentz__.cc, __gammainc_lentz__.cc: Removed.
* libinterp/corefcn/module.mk: Add renamed functions to build system.
* betainc.m: Use Octave standard GPL block. Rewrite parts of docstring.
Don't use array brackets around single output of function. Remove isscalar
checks on inputs because common_size() function will already handle it.
Use capital variable names in error messages to match documentation as displayed
in terminal. Reshape all inputs in to column vectors quickly so that input
validation tests that depend on all/any will pass with N-D arrays. Add
comments to code. Check for specific error messages in input validation BIST
tests.
* betaincinv.m: Use Octave standard GPL block. Rewrite parts of docstring.
Don't use array brackets around single output of function. Remove isscalar
checks on inputs because common_size() function will already handle it.
Use capital variable names in error messages to match documentation as displayed
in terminal. Reshape all inputs in to column vectors quickly so that input
validation tests that depend on all/any will pass with N-D arrays. Put most
common case of tail ("lower") first in if/elseif trees. Call functions directly
with function handle rather than using unnecessary feval() call. Use numel
in preference to length. Rename variable i_miss to todo for clarity. Add
comments to code. Check for specific error messages in input validation BIST
tests.
* cosint.m: Use Octave standard GPL block. Rewrite parts of docstring.
Don't use array brackets around single output of function. Remove isscalar
checks on inputs because common_size() function will already handle it.
Add input validation check for isnumeric value. Convert integer classes to
double before proceeding. Rename i_miss to todo for clarity. Use isinf to
detect both -Inf and +Inf rather than separate tests. Use ++it in while
loop conditional to shorten loop blocks. Add Input validation BIST tests.
* expint.m: Remove Sylvain who was not actually an author on this file.
Rewrite parts of docstring. Rename variable sparse_x to orig_sparse.
Eliminate temporary variables res_tmp, x_s_tmp, ssum_tmp. Rename i_miss to
todo. Use Octave coding conventions throughout. Add comments to code.
* gammainc.m: Use Octave standard GPL block. Rewrite parts of docstring.
Remove isscalar checks on inputs because common_size() function will already
handle it. Use capital variable names in error messages to match documentation
as displayed in terminal. Reshape all inputs in to column vectors quickly so
that input validation tests that depend on all/any will pass with N-D arrays.
Put most common case of tail ("lower") first in if/elseif trees. Add input
validation of tail. Rename variable ii to idx for clarity. Rename variable
i_done to todo and switch polarity so that the '!' operator is not required
every time the variable is updated. Use indexing and direct assignment to
update todo rather than logical operator '&' which is slower. Use tolower on
tail variable and then switch strcmpi calls to strcmp. Reformat %!test blocks
in to %!assert blocks to be more compact. Check for specific error messages in
input validation BIST tests.
* gammaincinv.m: Use Octave standard GPL block. Rewrite parts of docstring.
Don't use array brackets around single output of function. Remove isscalar
checks on inputs because common_size() function will already handle it.
Add input validation check for iscomplex value. Use capital variable names in
error messages to match documentation as displayed in terminal. Reshape all
inputs in to column vectors quickly so that input validation tests that depend
on all/any will pass with N-D arrays. Rename i_miss to todo. Use numel
in preference to length. Call functions directly with function handle rather
than using unnecessary feval() call.Rename variable i_miss to todo for clarity.
Use it++ in while loop conditional to shorten loop blocks. Add comments to
code. Check for specific error messages in input validation BIST tests.
Add input validation BIST tests for all error messages.
* sinint.m: Use Octave standard GPL block. Rewrite parts of docstring.
Add input validation for isnumeric. Convert integers to double for
calculation. Reshape input to column vector. Rename variable sz to orig_sz
for clarity. rename i_miss to todo. Reformat BIST tests to mak them more
compact. Add input validation BIST tests.
| author | Rik <rik@octave.org> |
|---|---|
| date | Mon, 19 Mar 2018 10:01:48 -0700 |
| parents | ff80c319e664 |
| children | f38ac278ff7d 996c30831e1a |
| files | NEWS configure.ac libinterp/corefcn/__betainc__.cc libinterp/corefcn/__betainc_lentz__.cc libinterp/corefcn/__expint__.cc libinterp/corefcn/__expint_lentz__.cc libinterp/corefcn/__gammainc__.cc libinterp/corefcn/__gammainc_lentz__.cc libinterp/corefcn/module.mk scripts/specfun/betainc.m scripts/specfun/betaincinv.m scripts/specfun/cosint.m scripts/specfun/expint.m scripts/specfun/gammainc.m scripts/specfun/gammaincinv.m scripts/specfun/sinint.m |
| diffstat | 16 files changed, 1294 insertions(+), 1256 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Tue Mar 13 10:29:52 2018 +0100 +++ b/NEWS Mon Mar 19 10:01:48 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.
--- a/configure.ac Tue Mar 13 10:29:52 2018 +0100 +++ b/configure.ac Mon Mar 19 10:01:48 2018 -0700 @@ -298,8 +298,10 @@ AC_MSG_NOTICE([configuring Octave to use system fonts in $SYSTEM_FREEFONT_DIR]) AC_DEFINE_UNQUOTED([SYSTEM_FREEFONT_DIR], ["$SYSTEM_FREEFONT_DIR"], [Define this to be the system directory containing the GNU FreeFont fonts.]) +fi AM_CONDITIONAL([AMCOND_INSTALL_INTERNAL_FONT_FILES], [test -z "$SYSTEM_FREEFONT_DIR"]) + ### Determine which C++ compiler to use (we expect to find g++). AC_PROG_CXX @@ -2924,26 +2926,6 @@ ## Add warning flags -ENABLE_HG_ID=yes -AC_ARG_ENABLE([hg-id], - [AS_HELP_STRING([--disable-hg-id], - [disable embedding of hg id in libraries])], - [if test "$enableval" = no; then ENABLE_HG_ID=no; fi], []) -AM_CONDITIONAL([AMCOND_ENABLE_HG_ID], [test $ENABLE_HG_ID = yes]) - -### Determine whether to install build logs with Octave. - -install_build_logs=no -AC_ARG_ENABLE([install-build-logs], - [AS_HELP_STRING([--enable-install-build-logs], - [install build logs (i.e. config.log) with Octave])], - [if test "$enableval" = yes; then install_build_logs=yes; fi]) -AM_CONDITIONAL([AMCOND_INSTALL_BUILD_LOGS], [test $install_build_logs = yes]) - -### Add extra compiler flags now that feature testing is complete. - -## Add warning flags - dnl Don't add -Wshadow for GCC 4.x. case "$GCC_VERSION" in *4*) ;;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/corefcn/__betainc__.cc Mon Mar 19 10:01:48 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; +}
--- a/libinterp/corefcn/__betainc_lentz__.cc Tue Mar 13 10:29:52 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,191 +0,0 @@ -// 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 -// <http://www.gnu.org/licenses/>. - -#if defined (HAVE_CONFIG_H) -# include "config.h" -#endif -#include "defun.h" -#include "fNDArray.h" - -DEFUN (__betainc_lentz__, args, , "Continued fraction for incomplete beta function") -{ - int nargin = args.length (); - octave_value_list outargs; - if (nargin != 4) - print_usage (); - else - { - // Values initialized in single precision - FloatNDArray x_arg_s = args(0).array_value (); - FloatNDArray a_arg_s = args(1).array_value (); - FloatNDArray b_arg_s = args(2).array_value (); - bool is_single = args(3).bool_value (); - - // total number of scenarios: get maximum of length of all vectors - int len_x = x_arg_s.rows (); - int len_a = a_arg_s.rows (); - int len_b = b_arg_s.rows (); - int len = std::max(len_x, len_a); - len = std::max(len, len_b); - - // input checks - if ((len_x != len) || (len_a != len) || (len_b != len)) - error("__betainc_lentz__: expecting arguments of same dimension"); - - // initialize scenario dependent output: - dim_vector dim_scen (len, 1); - ColumnVector f (dim_scen); - - // Lentz's algorithm in two cases: double and single precision - if (! is_single) - { - NDArray x_arg = args(0).array_value (); - NDArray a_arg = args(1).array_value (); - NDArray b_arg = args(2).array_value (); - NDArray x (dim_scen); - NDArray a (dim_scen); - NDArray b (dim_scen); - - // initialize scenario dependent input values (idx either 0 or ii) - if (len_x == 1) - x.fill (x_arg(0)); - else - x = x_arg; - // - if (len_a == 1) - a.fill (a_arg(0)); - else - a = a_arg; - // - if (len_b == 1) - b.fill (b_arg(0)); - else - b = b_arg; - // Variables initialization - 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 via all scenarios - for (octave_idx_type ii = 0; ii < len; ++ii) - { - // catch ctrl + c - OCTAVE_QUIT; - // Variable initialization for the current element - xj = x(ii); - y = tiny; - Cj = y; - Dj = 0; - aj = a(ii); - bj = b(ii); - 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++; - } - if (! error_state) - f(ii) = y; - } - outargs(0) = f; - } - else - { - FloatNDArray x_s (dim_scen); - FloatNDArray a_s (dim_scen); - FloatNDArray b_s (dim_scen); - - // initialize scenario dependent input values (idx either 0 or ii) - if (len_x == 1) - x_s.fill (x_arg_s(0)); - else - x_s = x_arg_s; - // - if (len_a == 1) - a_s.fill (a_arg_s(0)); - else - a_s = a_arg_s; - // - if (len_b == 1) - b_s.fill (b_arg_s(0)); - else - b_s = b_arg_s; - // Variables initialization - 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 via all scenarios - for (octave_idx_type ii = 0; ii < len; ++ii) - { - // catch ctrl + c - OCTAVE_QUIT; - // Variable initialization for the current element - xj = x_s(ii); - y = tiny; - Cj = y; - Dj = 0; - aj = a_s(ii); - bj = b_s(ii); - 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++; - } - if (! error_state) - f(ii) = y; - } - outargs(0) = f; - } - } - return octave_value (outargs); -}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/corefcn/__expint__.cc Mon Mar 19 10:01:48 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; +}
--- a/libinterp/corefcn/__expint_lentz__.cc Tue Mar 13 10:29:52 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,175 +0,0 @@ -// 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 "defun.h" -#include "fNDArray.h" -#include "CNDArray.h" -#include "fCNDArray.h" -#include <complex> - -DEFUN (__expint_lentz__, args, , "Continued fraction for the exponential integral") -{ - int nargin = args.length (); - octave_value_list outargs; - if (nargin != 2) - print_usage (); - else - { - // Value initialized in single precision - FloatComplexNDArray x_arg_s = args(0).complex_array_value (); - bool is_single = args(1).bool_value (); - - int len_x = x_arg_s.rows (); - - // initialize scenario dependent output: - dim_vector dim_scen (len_x, 1); - ComplexColumnVector f (dim_scen); - - // Lentz's algorithm in two cases: double and single precision - - if (! is_single) - { - ComplexNDArray x_arg = args(0).complex_array_value (); - ComplexNDArray x (dim_scen); - - // initialize scenario dependent input values (idx either 0 or ii) - if (len_x == 1) - x.fill (x_arg(0)); - else - x = x_arg; - // Variables initialization - static const std::complex<double> tiny = pow (2, -100); - static const double eps = std::numeric_limits<double>::epsilon(); - std::complex<double> cone(1.0, 0.0); - std::complex<double> czero(0.0, 0.0); - std::complex<double> xj = x(0); - std::complex<double> y = tiny; - std::complex<double> Cj = y; - std::complex<double> Dj = czero; - std::complex<double> alpha_j = cone; - std::complex<double> beta_j = xj; - std::complex<double> Deltaj = czero; - int j = 1; - int maxit = 200; - // loop via all scenarios - for (octave_idx_type ii = 0; ii < len_x; ++ii) - { - // catch ctrl + c - OCTAVE_QUIT; - // Variable initialization for the current element - xj = x(ii); - 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++; - } - if (! error_state) - f(ii) = y; - } - outargs(0) = f; - } - else - { - FloatComplexNDArray x_s (dim_scen); - - FloatComplexColumnVector f (dim_scen); - - // initialize scenario dependent input values (idx either 0 or ii) - if (len_x == 1) - x_s.fill (x_arg_s(0)); - else - x_s = x_arg_s; - // Variables initialization - static const std::complex<float> tiny = pow (2, -50); - static const float eps = std::numeric_limits<float>::epsilon(); - std::complex<float> cone(1.0, 0.0); - std::complex<float> czero(0.0, 0.0); - std::complex<float> xj = x_s(0); - std::complex<float> y = tiny; - std::complex<float> Cj = y; - std::complex<float> Dj = czero; - std::complex<float> alpha_j = cone; - std::complex<float> beta_j = czero; - std::complex<float> Deltaj = czero; - int j = 1; - int maxit = 100; - // loop via all scenarios - for (octave_idx_type ii = 0; ii < len_x; ++ii) - { - // catch ctrl + c - OCTAVE_QUIT; - // Variable initialization for the current element - xj = x_s(ii); - 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++; - } - if (! error_state) - f(ii) = y; - } - outargs(0) = f; - } - } - return octave_value (outargs); -}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/corefcn/__gammainc__.cc Mon Mar 19 10:01:48 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/__gammainc_lentz__.cc Tue Mar 13 10:29:52 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,161 +0,0 @@ -// Copyright (C) 2017 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 -// <http://www.gnu.org/licenses/>. - -#if defined (HAVE_CONFIG_H) -# include "config.h" -#endif -#include "defun.h" -#include "fNDArray.h" - -DEFUN (__gammainc_lentz__, args, , "Continued fraction for incomplete gamma function") -{ - int nargin = args.length (); - octave_value_list outargs; - if (nargin != 3) - print_usage (); - else - { - // Values initialized in single precision - FloatNDArray x_arg_s = args(0).array_value (); - FloatNDArray a_arg_s = args(1).array_value (); - bool is_single = args(2).bool_value (); - - // total number of scenarios: get maximum of length of all vectors - int len_x = x_arg_s.rows (); - int len_a = a_arg_s.rows (); - int len = std::max(len_x,len_a); - - // input checks - if (len_x > 1 && len_x != len) - error("gammainc_lentz_vec: expecting x to be of length 1 or %d", len); - if (len_a > 1 && len_a != len) - error("gammainc_lentz_vec: expecting a to be of length 1 or %d", len); - - // initialize scenario dependent output: - dim_vector dim_scen (len, 1); - ColumnVector f (dim_scen); - - // Lentz's algorithm in two cases: double and single precision - - if (! is_single) - { - NDArray x_arg = args(0).array_value (); - NDArray a_arg = args(1).array_value (); - NDArray x (dim_scen); - NDArray a (dim_scen); - - // initialize scenario dependent input values (idx either 0 or ii) - if (len_x == 1) - x.fill (x_arg(0)); - else - x = x_arg; - // - if (len_a == 1) - a.fill (a_arg(0)); - else - a = a_arg; - // Variables initialization - 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 via all scenarios - for (octave_idx_type ii = 0; ii < len; ++ii) - { - // catch ctrl + c - OCTAVE_QUIT; - // Variable initialization for the current element - y = tiny; - Cj = y; - Dj = 0; - bj = x(ii) - a(ii) + 1; - aj = a(ii); - 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(ii) - j); - j++; - } - if (! error_state) - f(ii) = y; - } - outargs(0) = f; - } - else - { - FloatNDArray x_s (dim_scen); - FloatNDArray a_s (dim_scen); - - // initialize scenario dependent input values (idx either 0 or ii) - if (len_x == 1) - x_s.fill (x_arg_s(0)); - else - x_s = x_arg_s; - // - if (len_a == 1) - a_s.fill (a_arg_s(0)); - else - a_s = a_arg_s; - // Variables initialization - 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 via all scenarios - for (octave_idx_type ii = 0; ii < len; ++ii) - { - // catch ctrl + c - OCTAVE_QUIT; - // Variable initialization for the current element - y = tiny; - Cj = y; - Dj = 0; - bj = x_s(ii) - a_s(ii) + 1; - aj = a_s(ii); - 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_s(ii) - j); - j++; - } - if (! error_state) - f(ii) = y; - } - outargs(0) = f; - } - } - return octave_value (outargs); -}
--- a/libinterp/corefcn/module.mk Tue Mar 13 10:29:52 2018 +0100 +++ b/libinterp/corefcn/module.mk Mon Mar 19 10:01:48 2018 -0700 @@ -107,11 +107,11 @@ COREFCN_SRC = \ %reldir%/Cell.cc \ - %reldir%/__betainc_lentz__.cc \ + %reldir%/__betainc__.cc \ %reldir%/__contourc__.cc \ %reldir%/__dsearchn__.cc \ - %reldir%/__expint_lentz__.cc \ - %reldir%/__gammainc_lentz__.cc \ + %reldir%/__expint__.cc \ + %reldir%/__gammainc__.cc \ %reldir%/__ichol__.cc \ %reldir%/__ilu__.cc \ %reldir%/__lin_interpn__.cc \
--- a/scripts/specfun/betainc.m Tue Mar 13 10:29:52 2018 +0100 +++ b/scripts/specfun/betainc.m Mon Mar 19 10:01:48 2018 -0700 @@ -3,7 +3,7 @@ ## ## This file is part of Octave. ## -## Octave is free software; you can redistribute it and/or modify it +## 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. @@ -15,14 +15,12 @@ ## ## 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/>. - -## Authors: Michele Ginesi <michele.ginesi@gmail.com> +## <https://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {} {} betainc (@var{x}, @var{a}, @var{b}) +## @deftypefn {} {} betainc (@var{x}, @var{a}, @var{b}) ## @deftypefnx {} {} betainc (@var{x}, @var{a}, @var{b}, @var{tail}) -## Compute the incomplete beta function ratio. +## Compute the incomplete beta function. ## ## This is defined as ## @tex @@ -46,147 +44,145 @@ ## ## @end ifnottex ## -## with @var{x} real in [0,1], @var{a} and @var{b} real and strictly positive. -## If one of the input has more than one components, then the others must be -## scalar or of compatible dimensions. +## 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 or if @var{tail} is @qcode{"lower"} the incomplete beta function -## ratio 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 as +## 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{"lower"}) = -## 1 - betainc (@var{x}, @var{a}, @var{b}, @qcode{"upper"}). +## betainc (@var{x}, @var{a}, @var{b}, @qcode{"upper"}) = +## 1 - betainc (@var{x}, @var{a}, @var{b}, @qcode{"lower"}). ## -## Reference +## @code{betainc} uses a more sophisticated algorithm than subtraction to +## get numerically accurate results when the @qcode{"lower"} value is small. ## -## @nospell{A. Cuyt, V. Brevik Petersen, B. Verdonk, H. Waadeland, W.B. Jones} -## @cite{Handbook of Continued Fractions for Special Functions}, -## ch. 18. +## 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") +function y = betainc (x, a, b, tail = "lower") - if ((nargin > 4) || (nargin < 3)) + if (nargin < 3 || nargin > 4) print_usage (); endif - if ((! isscalar (x)) || (! isscalar (a)) || (! isscalar (b))) - [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 + [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 - if ((iscomplex (x)) || (iscomplex (a)) || (iscomplex (b))) - error ("betainc: inputs must be real or integer"); + ## 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"); + error ("betainc: A must be strictly positive"); + endif + + if (any (b <= 0)) + error ("betainc: B must be strictly positive"); endif - if (any (b <= 0)) - error ("betainc: b must be strictly positive"); + ## 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 - if (any ((x > 1) | (x < 0))) - error ("betainc: x must be between 0 and 1"); - 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 - sz = size (x); - x = x(:); - a = a(:); - b = b(:); - y = zeros (size (x)); - - # If any of the argument is single, also the output should be + ## Initialize output array + y = zeros (size (x), class (x)); - if ((strcmpi (class (y), "single")) || (strcmpi (class (a), "single")) || (strcmpi (class (b), "single"))) - a = single (a); - b = single (b); - x = single (x); - y = single (y); - endif + ## 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. - # In the following, we use the fact that the continued fraction we will - # use 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 used the property - # I_x(a,b) + I_(1-x) (b,a) = 1. - - if (strcmpi (tail, "upper")) + 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)); - elseif (strcmpi (tail, "lower")) - 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 fflag") + 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 + ## 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_lentz__ (x, a, b, strcmpi (class (x), "single")); + f = __betainc__ (x, a, b); - # We divide the continued fraction by B(a,b) / (x^a * (1-x)^b) - # to obtain I_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); + gammaln (a) - gammaln (b) + log (f); y = real (exp (y)); y(fflag) = 1 - y(fflag); - y = reshape (y, sz); + + ## Restore original shape + y = reshape (y, orig_sz); endfunction -## Tests from betainc.cc -# Double precision +## Double precision %!test %! a = [1, 1.5, 2, 3]; %! b = [4, 3, 2, 1]; -%! v1 = betainc (1,a,b); +%! 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); +%! v4 = 1 - betainc (1-x, b, a); %! assert (v1, v2, sqrt (eps)); %! assert (v3, v4, sqrt (eps)); -# Single precision +## Single precision %!test %! a = single ([1, 1.5, 2, 3]); %! b = single ([4, 3, 2, 1]); -%! v1 = betainc (1,a,b); +%! 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); +%! v4 = 1 - betainc (1-x, b, a); %! assert (v1, v2, sqrt (eps ("single"))); %! assert (v3, v4, sqrt (eps ("single"))); -# Mixed double/single precision +## Mixed double/single precision %!test %! a = single ([1, 1.5, 2, 3]); %! b = [4, 3, 2, 1]; @@ -198,10 +194,8 @@ %! assert (v1, v2, sqrt (eps ("single"))); %! assert (v3, v4, sqrt (eps ("single"))); -## New test - -%!test #<51157> -%! y = betainc([0.00780;0.00782;0.00784],250.005,49750.995); +%!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); @@ -212,25 +206,35 @@ %!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 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)); +%! [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.1,1,1) -%!error betainc (-0.1,1,1) -%!error betainc (0.5,0,1) -%!error betainc (0.5,1,0) -%!error betainc (1,2,3,4) -%!error betainc (1,2) %!error betainc (1,2,3,4,5) -%!error betainc (0.5i, 1, 2) -%!error betainc (0, 1i, 1) -%!error betainc (0, 1, 1i) +%!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")
--- a/scripts/specfun/betaincinv.m Tue Mar 13 10:29:52 2018 +0100 +++ b/scripts/specfun/betaincinv.m Mon Mar 19 10:01:48 2018 -0700 @@ -2,19 +2,19 @@ ## ## This file is part of Octave. ## -## Octave is free software; you can redistribute it and/or modify it +## 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. +## 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. +## 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/>. +## <https://www.gnu.org/licenses/>. ## Author: Michele Ginesi <michele.ginesi@gmail.com> @@ -52,12 +52,13 @@ ## 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. +## 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 the inverse of the incomplete beta function integrated from 0 -## to @var{x} is computed. If @qcode{"upper"} is given then the complementary -## function integrated from @var{x} to 1 is inverted. +## 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 @@ -68,143 +69,144 @@ ## @ifnottex ## ## @example -## @group -## ## @var{y} - betainc (@var{x}, @var{a}, @var{b}) = 0 -## -## @end group ## @end example ## ## @end ifnottex ## -## -## @seealso{beta, betainc, betaln} +## @seealso{betainc, beta, betaln} ## @end deftypefn -function [x] = betaincinv (y, a, b, tail = "lower") +function x = betaincinv (y, a, b, tail = "lower") - if (nargin > 4 || nargin < 3) + if (nargin < 3 || nargin > 4) print_usage (); endif - if (! isscalar (y) || ! isscalar (a) || ! isscalar (b)) - [err, y, a, b] = common_size (y, a, b); - if (err > 0) - error ("betaincinv: y, a must be of common size or scalars"); - 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: inputs must be real or integer"); + 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) | any (b <= 0)) - error ("betaincinv: a must be strictly positive"); + 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 (y > 1 | y < 0)) - error ("betaincinv: y must be between 0 and 1"); + ## 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 - - # Extract the size. - sz = size (y); - # Write the inputs as two column vectors. - y = y(:); - a = a(:); - b = b(:); - l = length (y); - # Initialise the output. - x = zeros (l, 1); + ## Initialize output array + x = zeros (size (y), class (y)); - # If one of the inputs is of single type, also the output should be - - if (strcmpi (class (y), "single") || strcmpi (class (a), "single") || strcmpi (class (b), "single")) - a = single (a); - b = single (b); - y = single (y); - x = single (x); - endif - - # Parameters of the Newton method + ## Parameters for the Newton method maxit = 20; tol = eps (class (y)); - if (strcmpi (tail, "upper")) - p = 1 - y; - q = y; - x(y == 0) = 1; - x(y == 1) = 0; - elseif (strcmpi (tail, "lower")) + 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") + error ("betaincinv: invalid value for TAIL") endif - # Trivial values have been already computed. - i_miss = ((y != 0) & (y != 1)); + ## 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. + ## We will invert the lower version for p < 0.5 and the upper otherwise. i_low = (p < 0.5); - i_upp = (!i_low); + i_upp = (! i_low); - len_low = nnz (i_miss & i_low); - len_upp = nnz (i_miss & i_upp); - - if (any (i_miss & i_low)); - # Function and derivative of the lower version. + 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))); - # We compute the initial guess with a bisection method of 10 steps. - x0 = bisection_method (F, zeros (len_low,1), ones (len_low,1),... - a(i_low), b(i_low), p(i_low), 10); + 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); + tol, maxit); endif - if (any (i_miss & i_upp)); - # Function and derivative of the lower version. + 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))); - # We compute the initial guess with a bisection method of 10 steps. - x0 = bisection_method (F, zeros (len_upp,1), ones (len_upp,1),... - a(i_upp), b(i_upp), q(i_upp), 10); + 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); + tol, maxit); endif - ## Re-organize the output. - - x = reshape (x, sz); + ## 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 = feval (F, xl, a, b, y); - F_r = feval (F, xr, a, b, y); - for it = 1: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 = feval (F, xc, a, b, y); + 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); @@ -214,57 +216,55 @@ 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 + endfor endfunction function x = newton_method (F, JF, x0, a, b, y, tol, maxit); - l = length (y); - res = -feval (F, x0, a, b, y) ./ feval (JF, x0, a, b); - i_miss = (abs(res) >= tol * abs (x0)); + 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 (i_miss) && (it < maxit)) + while (any (todo) && (it < maxit)) it++; - x(i_miss) += res(i_miss); - res(i_miss) = - feval (F, x(i_miss), a(i_miss), b(i_miss), ... - y(i_miss)) ./ feval (JF, x(i_miss), a(i_miss), b(i_miss)); - i_miss = (abs(res) >= tol * abs (x)); + 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 - %!test -%! x = linspace(0.1, 0.9, 11); +%! 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); +%! 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); +%! 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); +%! 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 on the conservation of the class +## 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") @@ -276,8 +276,27 @@ ## Test input validation %!error betaincinv () %!error betaincinv (1) -%!error betaincinv (1,2,3,4) -%!error betaincinv (1, "2") -%!error betaincinv (0.5i, 1, 2) -%!error betaincinv (0, 1i, 1) -%!error betaincinv (0, 1, 1i) +%!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")
--- a/scripts/specfun/cosint.m Tue Mar 13 10:29:52 2018 +0100 +++ b/scripts/specfun/cosint.m Mon Mar 19 10:01:48 2018 -0700 @@ -2,19 +2,19 @@ ## ## This file is part of Octave. ## -## Octave is free software; you can redistribute it and/or modify it +## 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. +## 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. +## 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/>. +## <https://www.gnu.org/licenses/>. ## Author: Michele Ginesi <michele.ginesi@gmail.com> @@ -35,9 +35,9 @@ ## Ci (x) = - | (cos (t)) / t dt ## / ## x -## ## @end group ## @end example +## ## @end ifnottex ## An equivalent definition is ## @tex @@ -49,16 +49,16 @@ ## ## @example ## @group -## ## x ## / -## | cos (t) - 1 -## Ci (x) = gamma + log (x) + | --------------- dt -## | t +## | cos (t) - 1 +## Ci (x) = gamma + log (x) + | ------------- dt +## | t ## / ## 0 ## @end group ## @end example +## ## @end ifnottex ## Reference: ## @@ -66,39 +66,49 @@ ## @cite{Handbook of Mathematical Functions} ## 1964. ## -## @seealso{expint, cos, sinint} +## @seealso{sinint, expint, cos} ## ## @end deftypefn -function [y] = cosint (x) +function y = cosint (x) if (nargin != 1) print_usage (); endif - sz = size (x); - #x = reshape (x, numel (x), 1); + 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)) - ## workaround reshape narrowing to real (#52953) + ## 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)); - i_miss = true (length (x), 1); + todo = true (size (x)); - ## special values + ## Special values y(x == Inf) = 0; - y((x == -Inf) & !signbit (imag(x))) = 1i * pi; - y((x == -Inf) & signbit (imag(x))) = -1i * pi; + y((x == -Inf) & !signbit (imag (x))) = 1i * pi; + y((x == -Inf) & signbit (imag (x))) = -1i * pi; - i_miss = ((i_miss) & (x != Inf) & (x != -Inf)); + todo(isinf (x)) = false; - ## For values large in modulus and not in (-oo,0), we use the relation - ## with expint + ## 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); @@ -111,35 +121,38 @@ yy = -0.5 * (expint (1i * xx) + expint (-1i * xx)); yy(ii_nw) += 1i * pi; yy(ii_sw) = conj (yy(ii_sw)); - y(i_miss & flag_large) = yy; + y(todo & flag_large) = yy; + + todo(flag_large) = false; ## For values small in modulus, use the series expansion (also near (-oo, 0]) - i_miss = ((i_miss) & (! flag_large)); 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)(i_miss), imag (x)(i_miss)); + xx = complex (real (x)(todo), imag (x)(todo)); else - xx = x(i_miss); + 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 (i_miss), 1); - it = 1; + 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)); + 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); - it++; endwhile + y(todo) = yy; - y(i_miss) = yy; - 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); - y = reshape (y, sz); + ## 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 @@ -164,7 +177,7 @@ %!assert (class (cosint (single (1))), "single") -##tests against maple +## 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) @@ -180,7 +193,7 @@ %! 0.422980828774864996 %! 0.119629786008000328 %! -0.140981697886930412]; -%! assert (cosint(x), y_ex, -3e-15); +%! assert (cosint (x), y_ex, -3e-15); %!test %! x = -(1:4).'; @@ -191,7 +204,7 @@ %! assert (cosint (x), y_ex, -4*eps); %!test -%! x = complex(-(1:4).', 0); +%! x = complex (-(1:4).', 0); %! y_ex = [0.337403922900968135 + pi*1i %! 0.422980828774864996 + pi*1i %! 0.119629786008000328 + pi*1i @@ -233,15 +246,23 @@ %! -8.63307471207423322 + 3.13159298695312800*I %! 0.0698222284673061493 - 3.11847446254772946*I ]; %! B = cosint (x); -%! assert (A, B, -3*eps) +%! assert (A, B, -3*eps); %! B = cosint (single (x)); -%! assert (A, B, -3*eps ("single")) +%! assert (A, B, -3*eps ("single")); -## fails along negative real axis +## 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) +%! 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 Tue Mar 13 10:29:52 2018 +0100 +++ b/scripts/specfun/expint.m Mon Mar 19 10:01:48 2018 -0700 @@ -1,4 +1,3 @@ -## Copyright (C) 2006-2017 Sylvain Pelissier ## 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 @@ -79,15 +77,12 @@ ## References: ## ## @nospell{M. Abramowitz and I.A. Stegun}, -## @cite{Handbook of Mathematical Functions} -## 1964. +## @cite{Handbook of Mathematical Functions}, 1964. ## ## @nospell{N. Bleistein and R.A. Handelsman}, -## @cite{Asymptotic expansions of integrals} -## 1986. +## @cite{Asymptotic expansions of integrals}, 1986. ## -## @seealso{cosint, exp, sinint} -## +## @seealso{cosint, sinint, exp} ## @end deftypefn function E1 = expint (x) @@ -98,16 +93,19 @@ 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 ## Initialize the result - if ((isreal (x)) && (x >= 0)) + if (isreal (x) && x >= 0) E1 = zeros (size (x), class (x)); else E1 = complex (zeros (size (x), class (x))); @@ -115,77 +113,61 @@ 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) ... - | (real (x) < 0 & abs (imag (x)) <= 1e-08); - cf_idx = ((((real (x) + 1) .^ 2 ./ (38 ^ 2) + ... - imag (x) .^ 2 ./ (40 ^ 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) ... & (! 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. ## FIXME: Why so long? IEEE double doesn't have this much precision. gm = 0.577215664901532860606512090082402431042159335; - e1_s = - gm - log (x_s); - res = - x_s; + e1_s = -gm - log (x_s); + res = -x_s; ssum = res; k = 1; - fflag = true (size (res)); - while ((k < 1e03) && (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. + ## Modified Lentz's algorithm, from "Numerical recipes in Fortran 77" p.165. - e1_cf = exp(- x_cf); - e1_cf .*= __expint_lentz__ (x_cf, strcmpi (class (x_cf), "single")); + e1_cf = exp (-x_cf) .* __expint__ (x_cf); - # Remove spurious imaginary part if needed (__expint_lentz__ works - # automathically with complex values) + ## 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. + ## "Asymptotic expansion of integrals", pages 1-4. e1_a = exp (-x_a) ./ x_a; - oldres = ssum = res = ones (size (x_a), class (x_a)); + 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; + todo = true (size (x_a)); + while (k < 1e3 && any (todo)) + res(todo) ./= (- x_a(todo) / (k + 1)); + ssum(todo) += res(todo); k += 1; - res(fflag) = res_tmp; - oldres(fflag) = oldres_tmp; - fflag = abs (x_a) > k; + todo = abs (x_a) > k; endwhile e1_a .*= ssum; @@ -193,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 @@ -268,7 +252,7 @@ ## Test on the correct Image set %!assert (isreal (expint (linspace (0, 100)))) -%!assert (!isreal (expint (-1))) +%!assert (! isreal (expint (-1))) ## Test input validation %!error expint ()
--- a/scripts/specfun/gammainc.m Tue Mar 13 10:29:52 2018 +0100 +++ b/scripts/specfun/gammainc.m Mon Mar 19 10:01:48 2018 -0700 @@ -5,7 +5,7 @@ ## ## This file is part of Octave. ## -## Octave is free software; you can redistribute it and/or modify it +## 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. @@ -17,7 +17,7 @@ ## ## 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/>. +## <https://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {} {} gammainc (@var{x}, @var{a}) @@ -50,14 +50,14 @@ ## 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, the sizes of @var{x} and +## 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 nonnegative. +## The elements of @var{a} must be non-negative. ## -## By default or if @var{tail} is @qcode{"lower"} 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. +## 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 @@ -73,20 +73,17 @@ ## References: ## ## @nospell{M. Abramowitz and I. Stegun}, -## @cite{Handbook of mathematical functions} -## @nospell{Dover publications, INC.}, -## 1972. +## @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. +## @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. +## @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, gammainc, gammaln} +## @seealso{gamma, gammaincinv, gammaln} ## @end deftypefn ## P(a,x) = gamma(a,x)/Gamma(a), upper @@ -94,134 +91,145 @@ function y = gammainc (x, a, tail = "lower") - if ((nargin >= 4) || (nargin <= 1)) + if (nargin < 2 || nargin > 3) print_usage (); endif - if ((! isscalar (x)) || (! isscalar (a))) - [err, x, a] = common_size (x, a); - if (err > 0) - error ("gammainc: x, a must be of common size or scalars"); - 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 ((any (a < 0)) || (any (imag (a) != 0))) - error ("gammainc: a must be real and non negative"); + if (iscomplex (x) || iscomplex (a)) + error ("gammainc: all inputs must be real"); endif - if (any (imag (x) != 0)) - error ("gammainc: x must be real"); + ## 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 - ## Initialize output array. - if (isinteger (x)) - x = double (x); + if (nargin == 3 + && ! any (strcmpi (tail, {"lower","upper","scaledlower","scaledupper"}))) + error ("gammainc: invalid value for TAIL"); endif + tail = tolower (tail); - if ((strcmpi (class (a), "single")) || (strcmpi (class (x), "single"))) + ## 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 - y = zeros (size (x), class (x)); + ## 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. - i_done = false (size (x)); # Track which elements have been calculated. + todo = true (size (x)); # Track which elements need to be calculated. ## Case 0: x == Inf, a == Inf - - ii = ((x == Inf) & (a == Inf)); - if (any (ii(:))) - y(ii) = NaN; - i_done(ii) = true; + idx = (x == Inf) & (a == Inf); + if (any (idx)) + y(idx) = NaN; + todo(idx) = false; endif ## Case 1: x == 0, a == 0. - ii = ((x == 0) & (a == 0)); - if (any (ii(:))) - y(ii) = gammainc_00 (tail); - i_done(ii) = true; + idx = (x == 0) & (a == 0); + if (any (idx)) + y(idx) = gammainc_00 (tail); + todo(idx) = false; endif ## Case 2: x == 0. - ii = ((! i_done) & (x == 0)); - if (any (ii(:))) - y(ii) = gammainc_x0 (tail); - i_done(ii) = true; + idx = todo & (x == 0); + if (any (idx)) + y(idx) = gammainc_x0 (tail); + todo(idx) = false; endif ## Case 3: x = Inf - ii = ((! i_done) & (x == Inf)); - if (any (ii(:))) - y(ii) = gammainc_x_inf (tail); - i_done(ii) = true; + idx = todo & (x == Inf); + if (any (idx)) + y(idx) = gammainc_x_inf (tail); + todo(idx) = false; endif ## Case 4: a = Inf - ii = ((! i_done) & (a == Inf)); - if (any (ii(:))) - y(ii) = gammainc_a_inf (tail); - i_done(ii) = true; + idx = todo & (a == Inf); + if (any (idx)) + y(idx) = gammainc_a_inf (tail); + todo(idx) = false; endif ## Case 5: a == 0. - ii = ((! i_done) & (a == 0)); - if (any (ii(:))) - y(ii) = gammainc_a0 (x(ii), tail); - i_done(ii) = true; + idx = todo & (a == 0); + if (any (idx)) + y(idx) = gammainc_a0 (x(idx), tail); + todo(idx) = false; endif ## Case 6: a == 1. - ii = ((! i_done) & (a == 1)); - if (any (ii(:))) - y(ii) = gammainc_a1 (x(ii), tail); - i_done(ii) = true; + idx = todo & (a == 1); + if (any (idx)) + y(idx) = gammainc_a1 (x(idx), tail); + todo(idx) = false; endif - flag_an = ((a > 1) & (a == fix (a)) & (x <= 36) & (abs (x) >= 1e-01) & ... - (a <= 18)); - - ## Case 7: positive integer a; exp (x) and a! both under 1/eps. - ii = ((! i_done) & flag_an); - if (any (ii(:))) - y(ii) = gammainc_an (x(ii), a(ii), tail); - i_done(ii) = true; + ## 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 = ((abs (a) < 2) & (abs(a) > 0) & (! i_done) & (x < 0)); + 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)); + flag_s = (((x + 0.25 < a) | (x < 0) | (a < 5)) & (x > -20)) | (abs (x) < 1); ## Case 8: x, a relatively small. - ii = ((! i_done) & flag_s); - if (any (ii(:))) - y(ii) = gammainc_s (x(ii), a(ii), tail); - i_done(ii) = true; + 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. - ii = (! i_done); - if (any (ii(:))) - y(ii) = gammainc_l (x(ii), a(ii), tail); - i_done(ii) = true; + idx = todo; + if (any (idx)) + y(idx) = gammainc_l (x(idx), a(idx), tail); + todo(idx) = false; endif if (any (flag_a_small)) - if (strcmpi (tail, "lower")) + 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 (strcmpi (tail, "upper")) + 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 (strcmpi (tail, "scaledlower")) + 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 (strcmpi (tail, "scaledupper")) + 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; @@ -234,7 +242,7 @@ ## x == 0, a == 0. function y = gammainc_00 (tail) - if ((strcmpi (tail, "upper")) || (strcmpi (tail, "scaledupper"))) + if ((strcmp (tail, "upper")) || (strcmp (tail, "scaledupper"))) y = 0; else y = 1; @@ -243,10 +251,10 @@ ## x == 0. function y = gammainc_x0 (tail) - if ((strcmpi (tail, "upper")) || (strcmpi (tail, "scaledlower"))) + if (strcmp (tail, "lower")) + y = 0; + elseif (strcmp (tail, "upper") || strcmp (tail, "scaledlower")) y = 1; - elseif (strcmpi (tail, "lower")) - y = 0; else y = Inf; endif @@ -254,9 +262,9 @@ ## x == Inf. function y = gammainc_x_inf (tail) - if (strcmpi (tail, "lower")) + if (strcmp (tail, "lower")) y = 1; - elseif ((strcmpi (tail, "upper")) || (strcmpi (tail, "scaledupper"))) + elseif (strcmp (tail, "upper") || strcmp (tail, "scaledupper")) y = 0; else y = Inf; @@ -265,9 +273,9 @@ ## a == Inf. function y = gammainc_a_inf (tail) - if (strcmpi (tail, "lower")) + if (strcmp (tail, "lower")) y = 0; - elseif ((strcmpi (tail, "upper")) || (strcmpi (tail, "scaledlower"))) + elseif (strcmp (tail, "upper") || strcmp (tail, "scaledlower")) y = 1; else y = Inf; @@ -276,9 +284,9 @@ ## a == 0. function y = gammainc_a0 (x, tail) - if (strcmpi (tail, "lower")) + if (strcmp (tail, "lower")) y = 1; - elseif (strcmpi (tail, "scaledlower")) + elseif (strcmp (tail, "scaledlower")) y = exp (x); else y = 0; @@ -287,16 +295,16 @@ ## a == 1. function y = gammainc_a1 (x, tail) - if (strcmpi (tail, "lower")) + if (strcmp (tail, "lower")) y = 1 - exp (-x); - elseif (strcmpi (tail, "scaledlower")) + 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 - elseif (strcmpi (tail, "upper")) - y = exp (-x); else y = 1 ./ x; endif @@ -314,14 +322,14 @@ y(jj) += t(jj); i++; endwhile - if (strcmpi (tail, "upper")) - y .*= exp (-x); - elseif (strcmpi (tail, "lower")) + if (strcmp (tail, "lower")) y = 1 - exp (-x) .* y; - elseif (strcmpi (tail, "scaledupper")) + 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); - elseif (strcmpi (tail, "scaledlower")) - y = (1 - exp (-x) .* y) ./ D(x, a); endif endfunction @@ -329,7 +337,7 @@ ## Numerical Recipes in Fortran 77 (6.2.5) ## series function y = gammainc_s (x, a, tail) - if ((strcmpi (tail, "scaledlower")) || (strcmpi (tail, "scaledupper"))) + if (strcmp (tail, "scaledlower") || strcmp (tail, "scaledupper")) y = ones (size (x), class (x)); term = x ./ (a + 1); else @@ -346,9 +354,9 @@ y(jj) += term(jj); term(jj) .*= x(jj) ./ (a(jj) + n); endwhile - if (strcmpi (tail, "upper")) + if (strcmp (tail, "upper")) y = 1 - y; - elseif (strcmpi (tail, "scaledupper")) + elseif (strcmp (tail, "scaledupper")) y = 1 ./ D (x,a) - y; endif endfunction @@ -357,30 +365,32 @@ ## NRF77 (6.2.7) ## Gamma (a,x)/Gamma (a) ## Lentz's algorithm -## __gammainc_lentz__ in libinterp/corefcn/__gammainc_lentz__.cc +## __gammainc__ in libinterp/corefcn/__gammainc__.cc function y = gammainc_l (x, a, tail) - y = __gammainc_lentz__ (x, a, strcmpi (class (x), "single")); - if (strcmpi (tail, "upper")) - y .*= D (x, a); - elseif (strcmpi (tail, "lower")) - y = 1 - y .* D (x, a); - elseif (strcmpi (tail, "scaledlower")) - y = 1 ./ D (x, a) - y; - endif - endfunction + 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) - ## 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. - athresh = 10; ## FIXME: can be better tuned? + athresh = 10; # FIXME: can this be better tuned? y = zeros (size (x), class (x)); - i_done = false (size (x)); - i_done(x == 0) = true; - ii = ((! i_done) & (x > 0) & (a > athresh) & (a >= x)); - if (any (ii(:))) + + 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)- ... @@ -389,10 +399,11 @@ 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); - i_done(ii) = true; + todo(ii) = false; endif - ii = ((! i_done) & (x > 0) & (a > athresh) & (a < x)); - if (any (ii(:))) + + 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)- ... @@ -401,153 +412,111 @@ 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); - i_done(ii) = true; + todo(ii) = false; endif - ii = ((! i_done) & ((x <= 0) | (a <= athresh))); - if (any (ii(:))) ## standard formula for a not so large. + + 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)); - i_done(ii) = true; + todo(ii) = false; endif - ii = ((x < 0) & (a == fix (a))); - if any(ii(:)) ## remove spurious imaginary part. + + 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 -%!test -%! assert (gammainc ([0, 0, 1], [0, 1, 0]), [1, 0, 1]); -%!test -%! assert (gammainc ([0, 0, 1], [0, 1, 0], "upper"), [0, 1, 0]); -%!test -%! assert (gammainc ([0, 0, 1], [0, 1, 0], "scaledlower"), [1, 1, exp(1)]); -%!test -%! assert (gammainc ([0, 0, 1], [0, 1, 0], "scaledupper"), [0, Inf, 0]); +%!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 -%!test -%! assert (gammainc ([2, Inf], [Inf, 2]), [0, 1]); -%!test -%! assert (gammainc ([2, Inf], [Inf, 2], "upper"), [1, 0]); -%!test -%! assert (gammainc ([2, Inf], [Inf, 2], "scaledlower"), [1, Inf]); -%!test -%! assert (gammainc ([2, Inf], [Inf, 2], "scaledupper"), [Inf, 0]); +%!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 -%!test -%! % Here Matlab fails -%! assert (gammainc (-100,1,"upper"),exp (100),-eps); +## Matlab fails for this test +%!assert (gammainc (-100,1,"upper"), exp (100), -eps) ## Test: case 6 -%!test -%! assert (gammainc ([1, 2, 3], 1), 1 - exp (-[1, 2, 3])); -%!test -%! assert (gammainc ([1, 2, 3], 1, "upper"), exp (- [1, 2, 3])); -%!test -%! 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]); +%!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 -%!test -%! assert (gammainc (2, 2, "lower"), 0.593994150290162, -2e-15); -%!test -%! assert (gammainc (2, 2, "upper"), 0.406005849709838, -2e-15); -%!test -%! assert (gammainc (2, 2, "scaledlower"), 2.194528049465325, -2e-15); -%!test -%! assert (gammainc (2, 2, "scaledupper"), 1.500000000000000, -2e-15); -%!test -%! assert (gammainc ([3 2 36],[2 3 18], "upper"),... -%! [4/exp(3) 5*exp(-2) (4369755579265807723 / 2977975)/exp(36)]); -%!test -%! assert (gammainc (10, 10), 1 - (5719087 / 567) * exp (-10), -eps); -%!test -%! assert (gammainc (10, 10, "upper"), (5719087 / 567) * exp (-10), -eps); +%!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 -%!test -%! assert (gammainc (-10, 10), 3.112658265341493126871617e7, -2 * eps); -%!test -%! % Here Matlab fails -%! assert (isreal (gammainc (-10, 10)), true); -%!test -%! assert (gammainc (-10, 10.1, "upper"),... -%! -2.9582761911890713293e7-1i * 9.612022339061679758e6,-30 * eps); -%!test -%! assert (gammainc (-10, 10, "upper"), -3.112658165341493126871616e7, ... -%! -2 * eps); -%!test -%! assert (gammainc (-10, 10, "scaledlower"), ... -%! 0.5128019364747265,-1e-14); -%!test -%! assert (gammainc (-10, 10, "scaledupper"), ... -%! -0.5128019200000000, -1e-14); -%!test -%! assert (gammainc (200, 201, "upper"),... -%! 0.518794309678684497, -2 * eps); -%!test -%! assert (gammainc (200, 201, "scaledupper"), ... -%! 18.4904360746560462660798514, -eps); -%!test -%! % here we are very good (no D (x,a)) involved -%! assert (gammainc(1000, 1000.5, "scaledlower"), 39.48467539583672271, ... -%! -2* eps); -%!test -%! assert (gammainc (709, 1000, "upper"), ... -%! 0.99999999999999999999999954358, -eps); +%!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); -%!test -%! % Here Matlab is better -%! assert (gammainc (751, 750, "upper"),... -%! 0.4805914320558831327179457887, -12 * eps); -%!test -%! assert (gammainc (200, 200, "upper"), 0.49059658199276367497217454, ... -%! -4 * eps); -%!test -%! assert (gammainc (200, 200), 0.509403418007236325027825459574527043, ... -%! -3 * eps); -%!test -%! assert (gammainc (200, 200, "scaledupper"), 17.3984438553791505135122900, ... -%! -eps); -%!test -%! assert (gammainc (200, 200, "scaledlower"), 18.065406676779221643065, ... -%! -6 * eps); -%!test -%! assert (gammainc (201, 200, "upper"), 0.46249244908276709524913736667,... -%! -7 * eps); +%! 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 -%!test -%! 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.1), ... +%! [0.33239840504050, 0.20972940370977, 0.10511370061022, ... +%! 0.041846517936723], 1e-13); -%!test -%! assert (gammainc ([1e-05, 1e-07,1e-10,1e-14], 0.2), ... -%! [0.10891226058559, 0.043358823442178, 0.010891244210402, ... -%! 0.0017261458806785], 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); +%!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); +%!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); +%! assert (gammainc (-20, 1.1, "upper"), ... +%! 6.50986687074979e8 + 2.11518396291149e8*i, -1e-13); -## Test on the conservation of the class, five tests for each subroutine +## 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") @@ -583,9 +552,13 @@ %!error gammainc () %!error gammainc (1) %!error gammainc (1,2,3,4) -%!error gammainc (1, [0, 1i,1]) -%!error gammainc (1, [0, -1, 1]) -%!error gammainc ([0 0],[0; 0]) -%!error gammainc ([1 2 3], [1 2]) -%!error gammainc (1i,1) -%!error gammainc (1,1i) +%!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")
--- a/scripts/specfun/gammaincinv.m Tue Mar 13 10:29:52 2018 +0100 +++ b/scripts/specfun/gammaincinv.m Mon Mar 19 10:01:48 2018 -0700 @@ -2,7 +2,7 @@ ## ## This file is part of Octave. ## -## Octave is free software; you can redistribute it and/or modify it +## 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. @@ -14,7 +14,7 @@ ## ## 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/>. +## <https://www.gnu.org/licenses/>. ## Author: Michele Ginesi <michele.ginesi@gmail.com> @@ -44,21 +44,21 @@ ## @end ifnottex ## ## and @code{gammaincinv (gammainc (@var{x}, @var{a}), @var{a}) = @var{x}} -## for each nonnegative 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. +## 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, the sizes of @var{y} and +## 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 elements of @var{y} must be in @math{[0,1]} and those of @var{a} -## must be positive. +## The variable @var{y} must be in the interval @math{[0,1]} while @var{a} must +## be real and positive. ## -## By default or if @var{tail} is @qcode{"lower"} 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. +## 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 by standard Newton's method, by solving +## The function is computed with Newton's method by solving ## @tex ## $$ ## y - \gamma (x, a) = 0 @@ -67,207 +67,195 @@ ## @ifnottex ## ## @example -## @group -## ## @var{y} - gammainc (@var{x}, @var{a}) = 0 -## -## @end group ## @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. +## 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{gamma, gammainc, gammaln} +## @seealso{gammainc, gamma, gammaln} ## @end deftypefn -function [x] = gammaincinv (y, a, tail = "lower") +function x = gammaincinv (y, a, tail = "lower") - if (nargin >= 4 || nargin <= 1) + if (nargin < 2 || nargin > 3) print_usage (); endif - if (! isscalar (y) || ! isscalar (a)) - [err, y, a] = common_size (y, a); - if (err > 0) - error ("gammaincinv: y, a must be of common size or scalars"); - 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"); + error ("gammaincinv: A must be strictly positive"); endif - if (any (y > 1 | y < 0)) - error ("gammaincinv: y must be between 0 and 1"); + ## 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; - # Extract the size. - sz = size (y); - # Write the inputs as two column vectors. - y = y(:); - a = a(:); - l = length (y); - # Initialise the output. - x = zeros (l, 1); - - if (strcmpi (class (y), "single") || strcmpi (class (a), "single")) - a = single (a); - y = single (y); - x = single (x); - endif - tol = eps (class (y)); - # special cases, a = 1 or y = 0, 1. + ## Special cases, a = 1 or y = 0, 1. - if strcmpi (tail, "lower") + 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") + 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") + error ("gammaincinv: invalid value for TAIL") endif - i_miss = ((y != 0) & (y != 1) & (a != 1)); + todo = (a != 1) & (y != 0) & (y != 1); ## Case 1: p small. - i_flag_1 = p < ((0.2 * (1 + a)) .^ a) ./ gamma (1 + a); - i_flag_1 = ((i_flag_1) & (i_miss)); + 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. + ## Initial guess. - r = ((pp .* gamma (1 + aa)) .^ (1 ./ aa)); + 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)); + (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)); + ./ (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. + ## 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); - i_miss = ((i_miss) & (! i_flag_1)); + todo(i_flag_1) = false; ## Case 2: q small. - i_flag_2 = ((q < exp (- 0.5 * a) ./ gamma (1 + a)) & (a < 10) & (a > 0)); - i_flag_2 = ((i_flag_2) & (i_miss)); + 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. + ## Initial guess. x0 = (-log (qq) - gammaln (aa)); - # For this case, we invert the upper version. + ## 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); - i_miss = ((i_miss) & (! i_flag_2)); + todo(i_flag_2) = false; ## Case 3: a small. - i_flag_3 = ((a > 0) & (a < 1)); - i_flag_3 = ((i_flag_3) & (i_miss)); + i_flag_3 = todo & ((a > 0) & (a < 1)); aa = a(i_flag_3); pp = p(i_flag_3); - # Initial guess + ## Initial guess - xl = ((pp .* gamma (aa + 1)) .^ (1 ./ aa)); + xl = (pp .* gamma (aa + 1)) .^ (1 ./ aa); x0 = xl; - # For this case, we invert the lower version. + ## 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); - i_miss = ((i_miss) & (! i_flag_3)); + todo(i_flag_3) = false; ## Case 4: a large. - i_flag_4 = i_miss; - + i_flag_4 = todo; aa = a(i_flag_4); qq = q(i_flag_4); - # Initial guess + ## Initial guess d = 1 ./ (9 * aa); - t = 1 - d - norminv (qq) .* sqrt(d); + t = 1 - d + sqrt (2) * erfcinv (2 * qq) .* sqrt (d); x0 = aa .* (t .^ 3); - # For this case, we invert the upper version. + ## 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); - ## Reshape the output. + ## Restore original shape + x = reshape (x, orig_sz); - x = reshape (x, sz); endfunction -## Subfunction: Newton Method - +## Subfunction: Newton's Method function x = newton_method (F, JF, y, a, x0, tol, maxit); - l = length (y); - res = -feval (F, y, a, x0) ./ feval (JF, a, x0); - i_miss = (abs(res) >= tol * abs (x0)); + l = numel (y); + res = -F (y, a, x0) ./ JF (a, x0); + todo = (abs (res) >= tol * abs (x0)); x = x0; it = 0; - while (any (i_miss) && (it < maxit)) - it++; - x(i_miss) += res(i_miss); - res(i_miss) = - feval (F, y(i_miss), a(i_miss), x(i_miss)) ./ ... - feval (JF, a(i_miss), x(i_miss)); - i_miss = (abs(res) >= tol * abs (x)); + 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 - %!test %! x = [1e-10, 1e-09, 1e-08, 1e-07]; %! a = [2, 3, 4]; @@ -296,8 +284,7 @@ %! xx = gammainc (gammaincinv (x, a, "upper"), a, "upper"); %! assert (xx, x, -1e-13); - -## Test on the conservation of the class +## 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") @@ -310,4 +297,19 @@ %!error gammaincinv () %!error gammaincinv (1) %!error gammaincinv (1, 2, 3, 4) -%!error gammaincinv (1, "2") +%!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/sinint.m Tue Mar 13 10:29:52 2018 +0100 +++ b/scripts/specfun/sinint.m Mon Mar 19 10:01:48 2018 -0700 @@ -2,21 +2,19 @@ ## ## This file is part of Octave. ## -## Octave is free software; you can redistribute it and/or modify it +## 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. +## 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. +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see -## <http://www.gnu.org/licenses/>. - -## Author: Michele Ginesi <michele.ginesi@gmail.com> +## <https://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {} {} sinint (@var{x}) @@ -37,107 +35,116 @@ ## 0 ## @end group ## @end example +## ## @end ifnottex ## -## Reference: -## -## @nospell{M. Abramowitz and I.A. Stegun}, -## @cite{Handbook of Mathematical Functions} -## 1964. +## Reference: @nospell{M. Abramowitz and I.A. Stegun}, +## @cite{Handbook of Mathematical Functions}, 1964. ## ## @seealso{cosint, expint, sin} -## ## @end deftypefn -function [y] = sinint (x) +function y = sinint (x) if (nargin != 1) print_usage (); endif - sz = size (x); + 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)); - i_miss = true (length (x), 1); + todo = true (size (x)); - ## Trivial values - y(x == 0) = x(x == 0); # correctly signed zero + ## Special values + y(x == 0) = x(x == 0); # correctly signed zero y(x == Inf) = pi / 2; y(x == - Inf) = - pi / 2; - i_miss = ((i_miss) & (x != 0) & (x != Inf) & (x != - Inf)); + 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 & i_miss); + xx = x(flag_large & todo); ii_neg = (real (xx) < 0); xx(ii_neg) *= -1; - ii_conj = ((real (xx) == 0) & (imag (xx) < 0)); + 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(i_miss & flag_large) = yy; + y(todo & flag_large) = yy; ## For values small in modulus we use the series expansion - i_miss = ((i_miss) & (! flag_large)); - xx = x(i_miss); - ssum = xx; # First term of the series expansion + todo = (todo) & (! flag_large); + xx = x(todo); + ssum = xx; # First term of the series expansion yy = ssum; - flag_sum = true (nnz (i_miss), 1); + flag_sum = true (nnz (todo), 1); it = 0; maxit = 300; - while ((any (flag_sum)) && (it < maxit)); + 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(i_miss) = yy; + y(todo) = yy; - y = reshape (y, sz); + y = reshape (y, orig_sz); endfunction -%!test -%! x = 1.1; -%! %y = sym(11)/10; -%! A = sinint (x); -%! %B = double (sinint (y)); -%! B = 1.02868521867373; -%! assert (A, B, -5e-15); +%!assert (sinint (1.1), 1.02868521867373, -5e-15) %!test -%! %y = [2 3 sym(pi); exp(sym(1)) 5 6]; %! x = [2, 3, pi; exp(1), 5, 6]; %! A = sinint (x); -%! %B = double (sinint (y)); %! 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))) +%!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) +## 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) +%!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)'; @@ -190,3 +197,11 @@ %! 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")