annotate libinterp/corefcn/__gammainc__.cc @ 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
children 1db0b81efafe
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
1 /*
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
2
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
3 Copyright (C) 2017 Nir Krakauer
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
4 Copyright (C) 2018 Michele Ginesi
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
5 Copyright (C) 2018 Stefan Schlögl
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
6
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
7 This file is part of Octave.
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
8
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
9 Octave is free software; you can redistribute it and/or modify it
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
10 under the terms of the GNU General Public License as published by
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
11 the Free Software Foundation; either version 3 of the License, or
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
12 (at your option) any later version.
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
13
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
14 Octave is distributed in the hope that it will be useful, but
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
15 WITHOUT ANY WARRANTY; without even the implied warranty of
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
17 GNU General Public License for more details.
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
18
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
19 You should have received a copy of the GNU General Public License
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
20 along with Octave; see the file COPYING. If not, see
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
21 <http://www.gnu.org/licenses/>.
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
22
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
23 */
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
24
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
25 #if defined (HAVE_CONFIG_H)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
26 # include "config.h"
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
27 #endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
28
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
29 #include "defun.h"
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
30 #include "fNDArray.h"
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
31
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
32 DEFUN (__gammainc__, args, ,
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
33 doc: /* -*- texinfo -*-
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
34 @deftypefn {} {@var{y} =} __gammainc__ (@var{x}, @var{a})
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
35 Continued fraction for incomplete gamma function.
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
36 @end deftypefn */)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
37 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
38 int nargin = args.length ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
39
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
40 if (nargin != 2)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
41 print_usage ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
42
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
43 bool is_single = args(0).is_single_type () || args(1).is_single_type ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
44
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
45 // Total number of scenarios: get maximum of length of all vectors
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
46 int numel_x = args(0).numel ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
47 int numel_a = args(1).numel ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
48 int len = std::max (numel_x, numel_a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
49
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
50 octave_value_list retval;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
51 // Initialize output dimension vector
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
52 dim_vector output_dv (len, 1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
53
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
54 // Lentz's algorithm in two cases: single and double precision
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
55 if (is_single)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
56 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
57 // Initialize output and inputs
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
58 FloatColumnVector output (output_dv);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
59 FloatNDArray x, a;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
60
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
61 if (numel_x == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
62 x = FloatNDArray (output_dv, args(0).float_scalar_value ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
63 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
64 x = args(0).float_array_value ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
65
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
66 if (numel_a == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
67 a = FloatNDArray (output_dv, args(1).float_scalar_value ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
68 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
69 a = args(1).float_array_value ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
70
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
71 // Initialize variables used in algorithm
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
72 static const float tiny = pow (2, -50);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
73 static const float eps = std::numeric_limits<float>::epsilon();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
74 float y, Cj, Dj, bj, aj, Deltaj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
75 int j, maxit;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
76 maxit = 200;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
77
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
78 // Loop over all elements
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
79 for (octave_idx_type i = 0; i < len; ++i)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
80 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
81 // Catch Ctrl+C
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
82 OCTAVE_QUIT;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
83
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
84 // Variable initialization for the current element
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
85 y = tiny;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
86 Cj = y;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
87 Dj = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
88 bj = x(i) - a(i) + 1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
89 aj = a(i);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
90 Deltaj = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
91 j = 1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
92
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
93 // Lentz's algorithm
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
94 while ((std::abs ((Deltaj - 1) / y) > eps) && (j < maxit))
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
95 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
96 Cj = bj + aj/Cj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
97 Dj = 1 / (bj + aj*Dj);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
98 Deltaj = Cj * Dj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
99 y *= Deltaj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
100 bj += 2;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
101 aj = j * (a(i) - j);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
102 j++;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
103 }
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
104
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
105 output(i) = y;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
106 }
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
107
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
108 retval(0) = output;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
109 }
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
110 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
111 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
112 // Initialize output and inputs
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
113 ColumnVector output (output_dv);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
114 NDArray x, a;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
115
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
116 if (numel_x == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
117 x = NDArray (output_dv, args(0).scalar_value ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
118 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
119 x = args(0).array_value ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
120
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
121 if (numel_a == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
122 a = NDArray (output_dv, args(1).scalar_value ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
123 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
124 a = args(1).array_value ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
125
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
126 // Initialize variables used in algorithm
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
127 static const double tiny = pow (2, -100);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
128 static const double eps = std::numeric_limits<double>::epsilon();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
129 double y, Cj, Dj, bj, aj, Deltaj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
130 int j, maxit;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
131 maxit = 200;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
132
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
133 // Loop over all elements
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
134 for (octave_idx_type i = 0; i < len; ++i)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
135 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
136 // Catch Ctrl+C
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
137 OCTAVE_QUIT;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
138
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
139 // Variable initialization for the current element
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
140 y = tiny;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
141 Cj = y;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
142 Dj = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
143 bj = x(i) - a(i) + 1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
144 aj = a(i);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
145 Deltaj = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
146 j = 1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
147
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
148 // Lentz's algorithm
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
149 while ((std::abs ((Deltaj - 1) / y) > eps) && (j < maxit))
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
150 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
151 Cj = bj + aj/Cj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
152 Dj = 1 / (bj + aj*Dj);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
153 Deltaj = Cj * Dj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
154 y *= Deltaj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
155 bj += 2;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
156 aj = j * (a(i) - j);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
157 j++;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
158 }
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
159
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
160 output(i) = y;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
161 }
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
162
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
163 retval(0) = output;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
164 }
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
165
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
166 return retval;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
167 }