Mercurial > octave
changeset 24924:10d2fc9edaff stable
Added comments in Lentz's algorithm files
--
changed libinterp/corefcn/__betainc_lentz__.cc
changed libinterp/corefcn/__expint_lentz__.cc
changed libinterp/corefcn/__gammainc_lentz__.cc
| author | Michele Ginesi <michele.ginesi@gmail.com> |
|---|---|
| date | Sun, 25 Feb 2018 00:03:24 +0100 |
| parents | 40ab8be7d7ec |
| children | 24ae3461fb85 |
| files | libinterp/corefcn/__betainc_lentz__.cc libinterp/corefcn/__expint_lentz__.cc libinterp/corefcn/__gammainc_lentz__.cc |
| diffstat | 3 files changed, 31 insertions(+), 6 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/__betainc_lentz__.cc Sun Feb 25 00:02:44 2018 +0100 +++ b/libinterp/corefcn/__betainc_lentz__.cc Sun Feb 25 00:03:24 2018 +0100 @@ -31,6 +31,7 @@ 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 (); @@ -51,6 +52,7 @@ 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 (); @@ -70,10 +72,12 @@ 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; @@ -84,6 +88,7 @@ { // catch ctrl + c OCTAVE_QUIT; + // Variable initialization for the current element xj = x(ii); y = tiny; Cj = y; @@ -95,7 +100,7 @@ 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; @@ -132,10 +137,12 @@ 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; @@ -146,6 +153,7 @@ { // catch ctrl + c OCTAVE_QUIT; + // Variable initialization for the current element xj = x_s(ii); y = tiny; Cj = y; @@ -157,7 +165,7 @@ 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;
--- a/libinterp/corefcn/__expint_lentz__.cc Sun Feb 25 00:02:44 2018 +0100 +++ b/libinterp/corefcn/__expint_lentz__.cc Sun Feb 25 00:03:24 2018 +0100 @@ -33,6 +33,7 @@ print_usage (); else { + // Value initialized in single precision FloatComplexNDArray x_arg_s = args(0).complex_array_value (); bool is_single = args(1).bool_value (); @@ -40,19 +41,21 @@ // 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); - ComplexColumnVector f (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); @@ -71,6 +74,7 @@ { // catch ctrl + c OCTAVE_QUIT; + // Variable initialization for the current element xj = x(ii); y = tiny; Cj = y; @@ -79,6 +83,7 @@ beta_j = xj; Deltaj = czero; j = 1; + //Lentz's algorithm while((std::abs (Deltaj - cone) > eps) & (j < maxit)) { Dj = beta_j + alpha_j * Dj; @@ -113,6 +118,7 @@ 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); @@ -131,6 +137,7 @@ { // catch ctrl + c OCTAVE_QUIT; + // Variable initialization for the current element xj = x_s(ii); y = tiny; Cj = y; @@ -139,6 +146,7 @@ beta_j = xj; Deltaj = czero; j = 1; + //Lentz's algorithm while((std::abs (Deltaj - cone) > eps) & (j < maxit)) { Dj = beta_j + alpha_j * Dj;
--- a/libinterp/corefcn/__gammainc_lentz__.cc Sun Feb 25 00:02:44 2018 +0100 +++ b/libinterp/corefcn/__gammainc_lentz__.cc Sun Feb 25 00:03:24 2018 +0100 @@ -32,6 +32,7 @@ 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 (); @@ -51,6 +52,8 @@ 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 (); @@ -63,10 +66,12 @@ 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; @@ -77,6 +82,7 @@ { // catch ctrl + c OCTAVE_QUIT; + // Variable initialization for the current element y = tiny; Cj = y; Dj = 0; @@ -84,7 +90,7 @@ aj = a(ii); Deltaj = 0; j = 1; - + //Lentz's algorithm while((std::abs ((Deltaj - 1) / y) > eps) & (j < maxit)) { Cj = bj + aj/Cj; @@ -110,10 +116,12 @@ 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; @@ -124,6 +132,7 @@ { // catch ctrl + c OCTAVE_QUIT; + // Variable initialization for the current element y = tiny; Cj = y; Dj = 0; @@ -131,7 +140,7 @@ aj = a_s(ii); Deltaj = 0; j = 1; - + //Lentz's algorithm while((std::abs ((Deltaj - 1) / y) > eps) & (j < maxit)) { Cj = bj + aj/Cj;