Mercurial > octave
changeset 24918:ea3edda05b66 stable
Lentz, gammainc: added single precision
Duplicated the code of the algorithm in __gammainc_lentz__.cc to
handle both single and double precision
--
changed libinterp/corefcn/__gammainc_lentz__.cc
changed scripts/specfun/gammainc.m
| author | Michele Ginesi <michele.ginesi@gmail.com> |
|---|---|
| date | Tue, 20 Feb 2018 20:37:25 +0100 |
| parents | 509e78c26e82 |
| children | ed6f6bbed604 |
| files | libinterp/corefcn/__gammainc_lentz__.cc scripts/specfun/gammainc.m |
| diffstat | 2 files changed, 106 insertions(+), 52 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/__gammainc_lentz__.cc Sat Feb 17 18:51:40 2018 +0100 +++ b/libinterp/corefcn/__gammainc_lentz__.cc Tue Feb 20 20:37:25 2018 +0100 @@ -22,78 +22,132 @@ # 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 != 2) + if (nargin != 3) print_usage (); else { - NDArray x_arg = args(0).array_value (); - NDArray a_arg = args(1).array_value (); + 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.rows (); - int len_a = a_arg.rows (); + 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); + 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); + 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); - NDArray x (dim_scen); - NDArray a (dim_scen); + 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; + // 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; + 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; + y = tiny; + Cj = y; + Dj = 0; + bj = x(ii) - a(ii) + 1; + aj = a(ii); + Deltaj = 0; + j = 1; - 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; - y = tiny; - Cj = y; - Dj = 0; - bj = x(ii) - a(ii) + 1; - aj = a(ii); - Deltaj = 0; - j = 1; + 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); - 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; + // 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; + 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; + y = tiny; + Cj = y; + Dj = 0; + bj = x_s(ii) - a_s(ii) + 1; + aj = a_s(ii); + Deltaj = 0; + j = 1; + + 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/scripts/specfun/gammainc.m Sat Feb 17 18:51:40 2018 +0100 +++ b/scripts/specfun/gammainc.m Tue Feb 20 20:37:25 2018 +0100 @@ -118,8 +118,9 @@ x = double (x); endif - if (strcmpi (class (a), "single")) + if (strcmpi (class (a), "single") || strcmpi (class (x), "single")) x = single (x); + a = single (a); endif y = zeros (size (x), class (x)); @@ -358,8 +359,7 @@ ## Lentz's algorithm ## __gammainc_lentz__ in libinterp/corefcn/__gammainc_lentz__.cc function y = gammainc_l (x, a, tail) - % calling vectorizied version of c++ code - y = __gammainc_lentz__ (x, a); + y = __gammainc_lentz__ (x, a, strcmpi (class (x), "single")); if (strcmpi (tail, "upper")) y .*= D (x, a); elseif (strcmpi (tail, "lower"))