annotate libinterp/corefcn/__betainc__.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 85867171af5f
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) 2018 Michele Ginesi
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
4 Copyright (C) 2018 Stefan Schlögl
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
5
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
6 This file is part of Octave.
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
7
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
8 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
9 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
10 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
11 (at your option) any later version.
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
12
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
13 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
14 WITHOUT ANY WARRANTY; without even the implied warranty of
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
16 GNU General Public License for more details.
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
17
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
18 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
19 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
20 <http://www.gnu.org/licenses/>.
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
21
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 #if defined (HAVE_CONFIG_H)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
25 # include "config.h"
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
26 #endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
27
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
28 #include "defun.h"
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
29 #include "dNDArray.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 (__betainc__, 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} =} __betainc__ (@var{x}, @var{a}, @var{b})
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
35 Continued fraction for incomplete beta 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 != 3)
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 || args(2).is_single_type ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
45
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
46 // 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
47 int numel_x = args(0).numel ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
48 int numel_a = args(1).numel ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
49 int numel_b = args(2).numel ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
50 int len = std::max (std::max (numel_x, numel_a), numel_b);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
51
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
52 octave_value_list retval;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
53 // Initialize output dimension vector
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
54 dim_vector output_dv (len, 1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
55
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
56 // 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
57 if (is_single)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
58 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
59 // Initialize output and inputs
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
60 FloatColumnVector output (output_dv);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
61 FloatNDArray x, a, b;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
62
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
63 if (numel_x == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
64 x = FloatNDArray (output_dv, args(0).float_scalar_value ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
65 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
66 x = args(0).float_array_value ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
67
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
68
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
69 if (numel_a == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
70 a = FloatNDArray (output_dv, args(1).float_scalar_value ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
71 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
72 a = args(1).float_array_value ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
73
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
74 if (numel_b == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
75 b = FloatNDArray (output_dv, args(2).float_scalar_value ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
76 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
77 b = args(2).float_array_value ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
78
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
79 // Initialize variables used in algorithm
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
80 static const float tiny = pow (2, -50);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
81 static const float eps = std::numeric_limits<float>::epsilon ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
82 float xj, x2, y, Cj, Dj, aj, bj, Deltaj, alpha_j, beta_j;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
83 int j, maxit;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
84 maxit = 200;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
85
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
86 // Loop over all elements
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
87 for (octave_idx_type i = 0; i < len; ++i)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
88 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
89 // Catch Ctrl+C
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
90 OCTAVE_QUIT;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
91
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
92 // Variable initialization for the current element
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
93 xj = x(i);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
94 y = tiny;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
95 Cj = y;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
96 Dj = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
97 aj = a(i);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
98 bj = b(i);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
99 Deltaj = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
100 alpha_j = 1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
101 beta_j = aj - (aj * (aj + bj)) / (aj + 1) * xj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
102 x2 = xj * xj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
103 j = 1;
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 // Lentz's algorithm
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
106 while ((std::abs ((Deltaj - 1)) > eps) && (j < maxit))
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 Dj = beta_j + alpha_j * Dj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
109 if (Dj == 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
110 Dj = tiny;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
111 Cj = beta_j + alpha_j / Cj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
112 if (Cj == 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
113 Cj = tiny;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
114 Dj = 1 / Dj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
115 Deltaj = Cj * Dj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
116 y *= Deltaj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
117 alpha_j = ((aj + j - 1) * (aj + bj + j -1) * (bj - j) * j)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
118 / ((aj + 2 * j - 1) * (aj + 2 * j - 1)) * x2;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
119 beta_j = aj + 2 * j + ((j * (bj - j)) / (aj + 2 * j - 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
120 - ((aj + j) * (aj + bj + j)) / (aj + 2 * j + 1)) * xj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
121 j++;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
122 }
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
123
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
124 output(i) = y;
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
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
127 retval(0) = output;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
128 }
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
129 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
130 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
131 // Initialize output and inputs
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
132 ColumnVector output (output_dv);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
133 NDArray x, a, b;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
134
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
135 if (numel_x == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
136 x = NDArray (output_dv, args(0).scalar_value ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
137 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
138 x = args(0).array_value ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
139
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
140 if (numel_a == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
141 a = NDArray (output_dv, args(1).scalar_value ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
142 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
143 a = args(1).array_value ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
144
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
145 if (numel_b == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
146 b = NDArray (output_dv, args(2).scalar_value ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
147 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
148 b = args(2).array_value ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
149
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
150 // Initialize variables used in algorithm
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
151 static const double tiny = pow (2, -100);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
152 static const double eps = std::numeric_limits<double>::epsilon ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
153 double xj, x2, y, Cj, Dj, aj, bj, Deltaj, alpha_j, beta_j;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
154 int j, maxit;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
155 maxit = 200;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
156
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
157 // Loop over all elements
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
158 for (octave_idx_type i = 0; i < len; ++i)
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 // Catch Ctrl+C
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
161 OCTAVE_QUIT;
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 // Variable initialization for the current element
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
164 xj = x(i);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
165 y = tiny;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
166 Cj = y;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
167 Dj = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
168 aj = a(i);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
169 bj = b(i);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
170 Deltaj = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
171 alpha_j = 1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
172 beta_j = aj - (aj * (aj + bj)) / (aj + 1) * xj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
173 x2 = xj * xj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
174 j = 1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
175
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
176 // Lentz's algorithm
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
177 while ((std::abs ((Deltaj - 1)) > eps) && (j < maxit))
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
178 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
179 Dj = beta_j + alpha_j * Dj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
180 if (Dj == 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
181 Dj = tiny;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
182 Cj = beta_j + alpha_j / Cj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
183 if (Cj == 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
184 Cj = tiny;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
185 Dj = 1 / Dj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
186 Deltaj = Cj * Dj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
187 y *= Deltaj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
188 alpha_j = ((aj + j - 1) * (aj + bj + j - 1) * (bj - j) * j)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
189 / ((aj + 2 * j - 1) * (aj + 2 * j - 1)) * x2;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
190 beta_j = aj + 2 * j + ((j * (bj - j)) / (aj + 2 * j - 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
191 - ((aj + j) * (aj + bj + j)) / (aj + 2 * j + 1)) * xj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
192 j++;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
193 }
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
194
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
195 output(i) = y;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
196 }
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
197
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
198 retval(0) = output;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
199 }
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
200
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
201 return retval;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
202 }