GNU Octave  8.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
__expint__.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 2018-2023 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 "CNDArray.h"
31 #include "defun.h"
32 #include "fCNDArray.h"
33 
35 
36 DEFUN (__expint__, args, ,
37  doc: /* -*- texinfo -*-
38 @deftypefn {} {@var{y} =} __expint__ (@var{x})
39 Continued fraction expansion for the exponential integral.
40 @end deftypefn */)
41 {
42  if (args.length () != 1)
43  print_usage ();
44 
45  octave_value_list retval;
46 
47  bool is_single = args(0).is_single_type ();
48 
49  int numel_x = args(0).numel ();
50 
51  // Initialize output dimension vector
52  dim_vector output_dv (numel_x, 1);
53 
54  // Lentz's algorithm in two cases: single and double precision
55  if (is_single)
56  {
57  // Initialize output and inputs
58  FloatComplexColumnVector output (output_dv);
60 
61  if (numel_x == 1)
62  x = FloatComplexNDArray (output_dv, args(0).float_complex_value ());
63  else
64  x = args(0).float_complex_array_value ();
65 
66  // Initialize variables used in algorithm
67  static const FloatComplex tiny = math::exp2 (-50.0f);
68  static const float eps = std::numeric_limits<float>::epsilon ();
69  const FloatComplex cone (1.0, 0.0);
70  const FloatComplex czero (0.0, 0.0);
71  const int maxit = 100;
72 
73  // Loop over all elements
74  for (octave_idx_type i = 0; i < numel_x; ++i)
75  {
76  // Catch Ctrl+C
78 
79  // Variable initialization for the current element
80  FloatComplex xj = x(i);
81  FloatComplex y = tiny;
82  FloatComplex Cj = y;
83  FloatComplex Dj = czero;
84  FloatComplex alpha_j = cone;
85  FloatComplex beta_j = xj;
86  FloatComplex Deltaj = czero;
87  int j = 1;
88 
89  // Lentz's algorithm
90  while ((std::abs (Deltaj - cone) > eps) && (j < maxit))
91  {
92  Dj = beta_j + alpha_j * Dj;
93  if (Dj == czero)
94  Dj = tiny;
95  Cj = beta_j + alpha_j / Cj;
96  if (Cj == czero)
97  Cj = tiny;
98  Dj = cone / Dj;
99  Deltaj = Cj * Dj;
100  y *= Deltaj;
101  alpha_j = (j + 1) / 2;
102  if ((j % 2) == 0)
103  beta_j = xj;
104  else
105  beta_j = cone;
106  j++;
107  }
108 
109  output(i) = y;
110  }
111  retval(0) = output;
112  }
113  else
114  {
115  // Initialize output and inputs
116  ComplexColumnVector output (output_dv);
118 
119  if (numel_x == 1)
120  x = ComplexNDArray (output_dv, args(0).complex_value ());
121  else
122  x = args(0).complex_array_value ();
123 
124  // Initialize variables used in algorithm
125  static const Complex tiny = math::exp2 (-100.0);
126  static const double eps = std::numeric_limits<double>::epsilon ();
127  const Complex cone (1.0, 0.0);
128  const Complex czero (0.0, 0.0);
129  const int maxit = 200;
130 
131  // Loop over all scenarios
132  for (octave_idx_type i = 0; i < numel_x; ++i)
133  {
134  // Catch Ctrl+C
135  OCTAVE_QUIT;
136 
137  // Variable initialization for the current element
138  Complex xj = x(i);
139  Complex y = tiny;
140  Complex Cj = y;
141  Complex Dj = czero;
142  Complex alpha_j = cone;
143  Complex beta_j = xj;
144  Complex Deltaj = czero;
145  int j = 1;
146 
147  // Lentz's algorithm
148  while ((std::abs (Deltaj - cone) > eps) && (j < maxit))
149  {
150  Dj = beta_j + alpha_j * Dj;
151  if (Dj == czero)
152  Dj = tiny;
153  Cj = beta_j + alpha_j / Cj;
154  if (Cj == czero)
155  Cj = tiny;
156  Dj = cone / Dj;
157  Deltaj = Cj * Dj;
158  y *= Deltaj;
159  alpha_j = (j + 1) / 2;
160  if ((j % 2) == 0)
161  beta_j = xj;
162  else
163  beta_j = cone;
164  j++;
165  }
166 
167  output(i) = y;
168  }
169 
170  retval(0) = output;
171  }
172 
173  return retval;
174 }
175 
OCTAVE_END_NAMESPACE(octave)
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
T eps(const T &x)
Definition: data.cc:4942
OCTINTERP_API void print_usage(void)
Definition: defun-int.h:72
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition: defun.h:56
double exp2(double x)
Definition: lo-mappers.h:98
F77_RET_T const F77_DBLE * x
class OCTAVE_API ComplexNDArray
Definition: mx-fwd.h:39
class OCTAVE_API FloatComplexNDArray
Definition: mx-fwd.h:41
std::complex< double > Complex
Definition: oct-cmplx.h:33
std::complex< float > FloatComplex
Definition: oct-cmplx.h:34
static T abs(T x)
Definition: pr-output.cc:1678
#define OCTAVE_QUIT
Definition: quit.h:247