Mercurial > octave

annotate libinterp/corefcn/__betainc__.cc @ 27919:1891570abac8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update Octave Project Developers copyright for the new year In files that have the "Octave Project Developers" copyright notice, update for 2020.
author John W. Eaton <jwe@octave.org>
date Mon, 06 Jan 2020 22:29:51 -0500
parents b442ec6dda5c
children bd51beb6205e
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
27919
1891570abac8 update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 27918
diff changeset
3 Copyright (C) 2018-2020 The Octave Project Developers
27918
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 26376
diff changeset
4
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 26376
diff changeset
5 See the file COPYRIGHT.md in the top-level directory of this distribution
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 26376
diff changeset
6 or <https://octave.org/COPYRIGHT.html/>.
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 26376
diff changeset
7
24927
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 This file is part of Octave.
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
10
25059
85867171af5f maint: update GPL license header
Mike Miller <mtmiller@octave.org>
parents: 24927
diff changeset
11 Octave is free software: you can redistribute it and/or modify it
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
12 under the terms of the GNU General Public License as published by
25059
85867171af5f maint: update GPL license header
Mike Miller <mtmiller@octave.org>
parents: 24927
diff changeset
13 the Free Software Foundation, either version 3 of the License, or
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
14 (at your option) any later version.
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
15
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
16 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
17 WITHOUT ANY WARRANTY; without even the implied warranty of
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
19 GNU General Public License for more details.
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
20
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
21 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
22 along with Octave; see the file COPYING. If not, see
25059
85867171af5f maint: update GPL license header
Mike Miller <mtmiller@octave.org>
parents: 24927
diff changeset
23 <https://www.gnu.org/licenses/>.
24927
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 */
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
26
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
27 #if defined (HAVE_CONFIG_H)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
28 # include "config.h"
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
29 #endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
30
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
31 #include "defun.h"
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
32 #include "dNDArray.h"
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
33 #include "fNDArray.h"
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
34
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
35 DEFUN (__betainc__, args, ,
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
36 doc: /* -*- texinfo -*-
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
37 @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
38 Continued fraction for incomplete beta function.
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
39 @end deftypefn */)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
40 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
41 int nargin = args.length ();
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 if (nargin != 3)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
44 print_usage ();
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 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
47 || args(2).is_single_type ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
48
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
49 // 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
50 int numel_x = args(0).numel ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
51 int numel_a = args(1).numel ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
52 int numel_b = args(2).numel ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
53 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
54
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
55 octave_value_list retval;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
56 // Initialize output dimension vector
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
57 dim_vector output_dv (len, 1);
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 // 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
60 if (is_single)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
61 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
62 // Initialize output and inputs
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
63 FloatColumnVector output (output_dv);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
64 FloatNDArray x, a, b;
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_x == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
67 x = FloatNDArray (output_dv, args(0).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 x = args(0).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
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
72 if (numel_a == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
73 a = FloatNDArray (output_dv, args(1).float_scalar_value ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
74 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
75 a = args(1).float_array_value ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
76
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
77 if (numel_b == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
78 b = FloatNDArray (output_dv, args(2).float_scalar_value ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
79 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
80 b = args(2).float_array_value ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
81
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
82 // Initialize variables used in algorithm
25542
e5208e98ab92 fix ambiguous overload build error on Solaris (bug #54217)
Mike Miller <mtmiller@octave.org>
parents: 25059
diff changeset
83 static const float tiny = octave::math::exp2 (-50.0f);
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
84 static const float eps = std::numeric_limits<float>::epsilon ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
85 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
86 int j, maxit;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
87 maxit = 200;
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 // Loop over all elements
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
90 for (octave_idx_type i = 0; i < len; ++i)
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 // Catch Ctrl+C
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
93 OCTAVE_QUIT;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
94
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
95 // Variable initialization for the current element
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
96 xj = x(i);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
97 y = tiny;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
98 Cj = y;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
99 Dj = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
100 aj = a(i);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
101 bj = b(i);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
102 Deltaj = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
103 alpha_j = 1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
104 beta_j = aj - (aj * (aj + bj)) / (aj + 1) * xj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
105 x2 = xj * xj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
106 j = 1;
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 // Lentz's algorithm
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
109 while ((std::abs ((Deltaj - 1)) > eps) && (j < maxit))
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
110 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
111 Dj = beta_j + alpha_j * Dj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
112 if (Dj == 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
113 Dj = tiny;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
114 Cj = beta_j + alpha_j / Cj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
115 if (Cj == 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
116 Cj = tiny;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
117 Dj = 1 / Dj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
118 Deltaj = Cj * Dj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
119 y *= Deltaj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
120 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
121 / ((aj + 2 * j - 1) * (aj + 2 * j - 1)) * x2;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
122 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
123 - ((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
124 j++;
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 output(i) = y;
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
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
130 retval(0) = output;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
131 }
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
132 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
133 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
134 // Initialize output and inputs
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
135 ColumnVector output (output_dv);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
136 NDArray x, a, b;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
137
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
138 if (numel_x == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
139 x = NDArray (output_dv, args(0).scalar_value ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
140 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
141 x = args(0).array_value ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
142
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
143 if (numel_a == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
144 a = NDArray (output_dv, args(1).scalar_value ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
145 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
146 a = args(1).array_value ();
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 if (numel_b == 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
149 b = NDArray (output_dv, args(2).scalar_value ());
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
150 else
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
151 b = args(2).array_value ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
152
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
153 // Initialize variables used in algorithm
25542
e5208e98ab92 fix ambiguous overload build error on Solaris (bug #54217)
Mike Miller <mtmiller@octave.org>
parents: 25059
diff changeset
154 static const double tiny = octave::math::exp2 (-100.0);
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
155 static const double eps = std::numeric_limits<double>::epsilon ();
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
156 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
157 int j, maxit;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
158 maxit = 200;
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 // Loop over all elements
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
161 for (octave_idx_type i = 0; i < len; ++i)
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 // Catch Ctrl+C
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
164 OCTAVE_QUIT;
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 // Variable initialization for the current element
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
167 xj = x(i);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
168 y = tiny;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
169 Cj = y;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
170 Dj = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
171 aj = a(i);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
172 bj = b(i);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
173 Deltaj = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
174 alpha_j = 1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
175 beta_j = aj - (aj * (aj + bj)) / (aj + 1) * xj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
176 x2 = xj * xj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
177 j = 1;
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 // Lentz's algorithm
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
180 while ((std::abs ((Deltaj - 1)) > eps) && (j < maxit))
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
181 {
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
182 Dj = beta_j + alpha_j * Dj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
183 if (Dj == 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
184 Dj = tiny;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
185 Cj = beta_j + alpha_j / Cj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
186 if (Cj == 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
187 Cj = tiny;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
188 Dj = 1 / Dj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
189 Deltaj = Cj * Dj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
190 y *= Deltaj;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
191 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
192 / ((aj + 2 * j - 1) * (aj + 2 * j - 1)) * x2;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
193 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
194 - ((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
195 j++;
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 output(i) = y;
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 retval(0) = output;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
202 }
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
203
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
204 return retval;
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
diff changeset
205 }