GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
__betainc__.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 2018-2025 The Octave Project Developers
4//
5// See the file COPYRIGHT.md in the top-level directory of this
6// distribution or <https://octave.org/copyright/>.
7//
8// This file is part of Octave.
9//
10// Octave is free software: you can redistribute it and/or modify it
11// under the terms of the GNU General Public License as published by
12// the Free Software Foundation, either version 3 of the License, or
13// (at your option) any later version.
14//
15// Octave is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18// GNU General Public License for more details.
19//
20// You should have received a copy of the GNU General Public License
21// along with Octave; see the file COPYING. If not, see
22// <https://www.gnu.org/licenses/>.
23//
24////////////////////////////////////////////////////////////////////////
25
26#if defined (HAVE_CONFIG_H)
27# include "config.h"
28#endif
29
30#include "defun.h"
31#include "dNDArray.h"
32#include "fNDArray.h"
33
34OCTAVE_BEGIN_NAMESPACE(octave)
35DEFUN (__betainc__, args, ,
36 doc: /* -*- texinfo -*-
37@deftypefn {} {@var{y} =} __betainc__ (@var{x}, @var{a}, @var{b})
38Continued fraction for incomplete beta function.
39@end deftypefn */)
40{
41 if (args.length () != 3)
42 print_usage ();
43
44 bool is_single = (args(0).is_single_type () || args(1).is_single_type ()
45 || args(2).is_single_type ());
46
47 // Total number of scenarios: get maximum of length of all vectors
48 int numel_x = args(0).numel ();
49 int numel_a = args(1).numel ();
50 int numel_b = args(2).numel ();
51 int len = std::max (std::max (numel_x, numel_a), numel_b);
52
53 octave_value_list retval;
54 // Initialize output dimension vector
55 dim_vector output_dv (len, 1);
56
57 // Lentz's algorithm in two cases: single and double precision
58 if (is_single)
59 {
60 // Initialize output and inputs
61 FloatColumnVector output (output_dv);
62 FloatNDArray x, a, b;
63
64 if (numel_x == 1)
65 x = FloatNDArray (output_dv, args(0).float_scalar_value ());
66 else
67 x = args(0).float_array_value ();
68
69
70 if (numel_a == 1)
71 a = FloatNDArray (output_dv, args(1).float_scalar_value ());
72 else
73 a = args(1).float_array_value ();
74
75 if (numel_b == 1)
76 b = FloatNDArray (output_dv, args(2).float_scalar_value ());
77 else
78 b = args(2).float_array_value ();
79
80 // Initialize variables used in algorithm
81 static const float tiny = math::exp2 (-50.0f);
82 static constexpr float eps = std::numeric_limits<float>::epsilon ();
83 float xj, x2, y, Cj, Dj, aj, bj, Deltaj, alpha_j, beta_j;
84 int j, maxit;
85 maxit = 200;
86
87 // Loop over all elements
88 for (octave_idx_type i = 0; i < len; ++i)
89 {
90 // Catch Ctrl+C
91 OCTAVE_QUIT;
92
93 // Variable initialization for the current element
94 xj = x(i);
95 y = tiny;
96 Cj = y;
97 Dj = 0;
98 aj = a(i);
99 bj = b(i);
100 Deltaj = 0;
101 alpha_j = 1;
102 beta_j = aj - (aj * (aj + bj)) / (aj + 1) * xj;
103 x2 = xj * xj;
104 j = 1;
105
106 // Lentz's algorithm
107 while ((std::abs ((Deltaj - 1)) > eps) && (j < maxit))
108 {
109 Dj = beta_j + alpha_j * Dj;
110 if (Dj == 0)
111 Dj = tiny;
112 Cj = beta_j + alpha_j / Cj;
113 if (Cj == 0)
114 Cj = tiny;
115 Dj = 1 / Dj;
116 Deltaj = Cj * Dj;
117 y *= Deltaj;
118 alpha_j = ((aj + j - 1) * (aj + bj + j -1) * (bj - j) * j)
119 / ((aj + 2 * j - 1) * (aj + 2 * j - 1)) * x2;
120 beta_j = aj + 2 * j + ((j * (bj - j)) / (aj + 2 * j - 1)
121 - ((aj + j) * (aj + bj + j)) / (aj + 2 * j + 1)) * xj;
122 j++;
123 }
124
125 output(i) = y;
126 }
127
128 retval(0) = output;
129 }
130 else
131 {
132 // Initialize output and inputs
133 ColumnVector output (output_dv);
134 NDArray x, a, b;
135
136 if (numel_x == 1)
137 x = NDArray (output_dv, args(0).scalar_value ());
138 else
139 x = args(0).array_value ();
140
141 if (numel_a == 1)
142 a = NDArray (output_dv, args(1).scalar_value ());
143 else
144 a = args(1).array_value ();
145
146 if (numel_b == 1)
147 b = NDArray (output_dv, args(2).scalar_value ());
148 else
149 b = args(2).array_value ();
150
151 // Initialize variables used in algorithm
152 static const double tiny = math::exp2 (-100.0);
153 static constexpr double eps = std::numeric_limits<double>::epsilon ();
154 double xj, x2, y, Cj, Dj, aj, bj, Deltaj, alpha_j, beta_j;
155 int j, maxit;
156 maxit = 200;
157
158 // Loop over all elements
159 for (octave_idx_type i = 0; i < len; ++i)
160 {
161 // Catch Ctrl+C
162 OCTAVE_QUIT;
163
164 // Variable initialization for the current element
165 xj = x(i);
166 y = tiny;
167 Cj = y;
168 Dj = 0;
169 aj = a(i);
170 bj = b(i);
171 Deltaj = 0;
172 alpha_j = 1;
173 beta_j = aj - (aj * (aj + bj)) / (aj + 1) * xj;
174 x2 = xj * xj;
175 j = 1;
176
177 // Lentz's algorithm
178 while ((std::abs ((Deltaj - 1)) > eps) && (j < maxit))
179 {
180 Dj = beta_j + alpha_j * Dj;
181 if (Dj == 0)
182 Dj = tiny;
183 Cj = beta_j + alpha_j / Cj;
184 if (Cj == 0)
185 Cj = tiny;
186 Dj = 1 / Dj;
187 Deltaj = Cj * Dj;
188 y *= Deltaj;
189 alpha_j = ((aj + j - 1) * (aj + bj + j - 1) * (bj - j) * j)
190 / ((aj + 2 * j - 1) * (aj + 2 * j - 1)) * x2;
191 beta_j = aj + 2 * j + ((j * (bj - j)) / (aj + 2 * j - 1)
192 - ((aj + j) * (aj + bj + j)) / (aj + 2 * j + 1)) * xj;
193 j++;
194 }
195
196 output(i) = y;
197 }
198
199 retval(0) = output;
200 }
201
202 return retval;
203}
204
205OCTAVE_END_NAMESPACE(octave)
dim_vector
Vector representing the dimensions (size) of an Array.
Definition dim-vector.h:90
octave_idx_type
octave_value_list
Definition ovl.h:42
OCTAVE_BEGIN_NAMESPACE
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
eps
T eps(const T &x)
Definition data.cc:5008
print_usage
void print_usage()
Definition defun-int.h:72
DEFUN
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition defun.h:56
x
F77_RET_T const F77_DBLE * x
Definition lo-slatec-proto.h:38
OCTAVE_QUIT
#define OCTAVE_QUIT
Definition quit.h:250
len
F77_RET_T len
Definition xerbla.cc:61