Mercurial > octave

annotate scripts/specfun/expint.m @ 25015:baa7e37453b1 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added more tests for betainc and expint. * betainc.m: Added to show that bug #34405 is solved. * expint.m: Added tests to show that bugs #53351 and #47738 are solved.
author Michele Ginesi <michele.ginesi@gmail.com>
date Fri, 23 Mar 2018 15:43:59 +0100
parents c280560d9c96
children d7ad543255c5
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
24923
40ab8be7d7ec Fixed style in specfun scripts
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24922
diff changeset
1 ## Copyright (C) 2018 Michele Ginesi
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
2 ##
16586
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
3 ## This file is part of Octave.
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
4 ##
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23220
diff changeset
5 ## Octave is free software: you can redistribute it and/or modify it
16586
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
6 ## under the terms of the GNU General Public License as published by
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23220
diff changeset
7 ## the Free Software Foundation, either version 3 of the License, or
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22323
diff changeset
8 ## (at your option) any later version.
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
9 ##
16586
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
10 ## Octave is distributed in the hope that it will be useful, but
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22323
diff changeset
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22323
diff changeset
13 ## GNU General Public License for more details.
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
14 ##
16586
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
15 ## You should have received a copy of the GNU General Public License
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
16 ## along with Octave; see the file COPYING. If not, see
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23220
diff changeset
17 ## <https://www.gnu.org/licenses/>.
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
18
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
19 ## -*- texinfo -*-
20852
516bb87ea72e 2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents: 20162
diff changeset
20 ## @deftypefn {} {} expint (@var{x})
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
21 ## Compute the exponential integral:
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
22 ##
16587
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
23 ## @tex
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
24 ## $$
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
25 ## {\rm E_1} (x) = \int_x^\infty {e^{-t} \over t} dt
16587
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
26 ## $$
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
27 ## @end tex
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
28 ## @ifnottex
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
29 ##
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
30 ## @example
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
31 ## @group
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
32 ## +oo
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
33 ## /
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
34 ## | exp (-t)
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
35 ## E_1 (x) = | -------- dt
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
36 ## | t
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
37 ## /
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
38 ## x
16587
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
39 ## @end group
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
40 ## @end example
17513
fedcd3717ebc doc: grammarcheck of documentation before 3.8 release.
Rik <rik@octave.org>
parents: 17363
diff changeset
41 ##
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
42 ## @end ifnottex
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
43 ##
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
44 ## Note: For compatibility, this function uses the @sc{matlab} definition
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
45 ## of the exponential integral. Most other sources refer to this particular
20162
2645f9ef8c88 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 19593
diff changeset
46 ## value as @math{E_1 (x)}, and the exponential integral as
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
47 ## @tex
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
48 ## $$
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
49 ## {\rm Ei} (x) = - \int_{-x}^\infty {e^{-t} \over t} dt.
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
50 ## $$
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
51 ## @end tex
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
52 ## @ifnottex
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
53 ##
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
54 ## @example
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
55 ## @group
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
56 ## +oo
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
57 ## /
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
58 ## | exp (-t)
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
59 ## Ei (x) = - | -------- dt
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
60 ## | t
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
61 ## /
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
62 ## -x
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
63 ## @end group
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
64 ## @end example
17513
fedcd3717ebc doc: grammarcheck of documentation before 3.8 release.
Rik <rik@octave.org>
parents: 17363
diff changeset
65 ##
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
66 ## @end ifnottex
17514
5b916efea542 doc: spellcheck of documentation before 3.8 release.
Rik <rik@octave.org>
parents: 17513
diff changeset
67 ## The two definitions are related, for positive real values of @var{x}, by
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
68 ## @tex
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
69 ## $
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
70 ## E_1 (-x) = -{\rm Ei} (x) - i\pi.
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
71 ## $
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
72 ## @end tex
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
73 ## @ifnottex
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
74 ## @w{@code{E_1 (-x) = -Ei (x) - i*pi}}.
16587
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
75 ## @end ifnottex
24908
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
76 ##
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
77 ## References:
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
78 ##
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
79 ## @nospell{M. Abramowitz and I.A. Stegun},
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
80 ## @cite{Handbook of Mathematical Functions}, 1964.
24908
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
81 ##
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
82 ## @nospell{N. Bleistein and R.A. Handelsman},
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
83 ## @cite{Asymptotic expansions of integrals}, 1986.
24908
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
84 ##
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
85 ## @seealso{cosint, sinint, exp}
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
86 ## @end deftypefn
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
87
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
88 function E1 = expint (x)
16586
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
89
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
90 if (nargin != 1)
16586
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
91 print_usage ();
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
92 endif
16586
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
93
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
94 if (! isnumeric (x))
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
95 error ("expint: X must be numeric");
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
96 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
97
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
98 ## Convert to floating point if necessary
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
99 if (isinteger (x))
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
100 x = double (x);
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
101 endif
16586
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
102
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
103 orig_sparse = issparse (x);
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
104 orig_sz = size (x);
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
105 x = x(:); # convert to column vector
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
106
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
107 ## Initialize the result
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
108 if (isreal (x) && x >= 0)
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
109 E1 = zeros (size (x), class (x));
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
110 else
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
111 E1 = complex (zeros (size (x), class (x)));
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
112 endif
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
113 tol = eps (class (x));
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
114
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
115 ## Divide the input into 3 regions and apply a different algorithm for each.
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
116 ## s = series expansion, cf = continued fraction, a = asymptotic series
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
117 s_idx = (((real (x) + 19.5).^ 2 ./ (20.5^2) + ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
118 imag (x).^2 ./ (10^2)) <= 1) ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
119 | (real (x) < 0 & abs (imag (x)) <= 1e-8);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
120 cf_idx = ((((real (x) + 1).^2 ./ (38^2) + ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
121 imag (x).^2 ./ (40^2)) <= 1) ...
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
122 & (! s_idx)) & (real (x) <= 35);
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
123 a_idx = (! s_idx) & (! cf_idx);
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
124 x_s = x(s_idx);
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
125 x_cf = x(cf_idx);
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
126 x_a = x(a_idx);
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
127
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
128 ## Series expansion
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
129 ## Abramowitz, Stegun, "Handbook of Mathematical Functions",
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
130 ## formula 5.1.11, p 229.
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
131 ## FIXME: Why so long? IEEE double doesn't have this much precision.
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
132 gm = 0.577215664901532860606512090082402431042159335;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
133 e1_s = -gm - log (x_s);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
134 res = -x_s;
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
135 ssum = res;
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
136 k = 1;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
137 todo = true (size (res));
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
138 while (k < 1e3 && any (todo))
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
139 res(todo) .*= (k * (- x_s(todo)) / ((k + 1) ^ 2));
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
140 ssum(todo) += res(todo);
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
141 k += 1;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
142 todo = (abs (res) > (tol * abs (ssum)));
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
143 endwhile
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
144 e1_s -= ssum;
16586
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
145
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
146 ## Continued fraction expansion,
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
147 ## Abramowitz, Stegun, "Handbook of Mathematical Functions",
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
148 ## formula 5.1.22, p 229.
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
149 ## Modified Lentz's algorithm, from "Numerical recipes in Fortran 77" p.165.
24908
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
150
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
151 e1_cf = exp (-x_cf) .* __expint__ (x_cf);
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
152
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
153 ## Remove spurious imaginary part if needed (__expint__ works automatically
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
154 ## with complex values)
24908
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
155
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
156 if (isreal (x_cf) && x_cf >= 0)
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
157 e1_cf = real (e1_cf);
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
158 endif
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
159
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
160 ## Asymptotic series, from N. Bleistein and R.A. Handelsman
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
161 ## "Asymptotic expansion of integrals", pages 1-4.
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
162 e1_a = exp (-x_a) ./ x_a;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
163 ssum = res = ones (size (x_a), class (x_a));
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
164 k = 0;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
165 todo = true (size (x_a));
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
166 while (k < 1e3 && any (todo))
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
167 res(todo) ./= (- x_a(todo) / (k + 1));
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
168 ssum(todo) += res(todo);
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
169 k += 1;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
170 todo = abs (x_a) > k;
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
171 endwhile
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
172 e1_a .*= ssum;
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
173
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
174 ## Combine results from each region into final output
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
175 E1(s_idx) = e1_s;
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
176 E1(cf_idx) = e1_cf;
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
177 E1(a_idx) = e1_a;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
178
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
179 ## Restore shape and sparsity of input
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
180 E1 = reshape (E1, orig_sz);
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
181 if (orig_sparse)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
182 E1 = sparse (E1);
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
183 endif
21758
ffad2baa90f7 maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents: 20852
diff changeset
184
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
185 endfunction
16585
1a3bfb14b5da Add and fix tests for ellipj, ellipke, and expint
Mike Miller <mtmiller@ieee.org>
parents: 16584
diff changeset
186
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
187
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
188 ## The following values were computed with the Octave symbolic package
16585
1a3bfb14b5da Add and fix tests for ellipj, ellipke, and expint
Mike Miller <mtmiller@ieee.org>
parents: 16584
diff changeset
189 %!test
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
190 %! X = [-50 - 50i -30 - 50i -10 - 50i 5 - 50i 15 - 50i 25 - 50i
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
191 %! -50 - 30i -30 - 30i -10 - 30i 5 - 30i 15 - 30i 25 - 30i
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
192 %! -50 - 10i -30 - 10i -10 - 10i 5 - 10i 15 - 10i 25 - 10i
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
193 %! -50 + 5i -30 + 5i -10 + 5i 5 + 5i 15 + 5i 25 + 5i
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
194 %! -50 + 15i -30 + 15i -10 + 15i 5 + 15i 15 + 15i 25 + 15i
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
195 %! -50 + 25i -30 + 25i -10 + 25i 5 + 25i 15 + 25i 25 + 25i];
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
196 %! y_exp = [ -3.61285286166493e+19 + 6.46488018613387e+19i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
197 %! -4.74939752018180e+10 + 1.78647798300364e+11i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
198 %! 3.78788822381261e+01 + 4.31742823558278e+02i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
199 %! 5.02062497548626e-05 + 1.23967883532795e-04i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
200 %! 3.16785290137650e-09 + 4.88866651583182e-09i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
201 %! 1.66999261039533e-13 + 1.81161508735941e-13i;
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
202 %! 3.47121527628275e+19 + 8.33104448629260e+19i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
203 %! 1.54596484273693e+11 + 2.04179357837414e+11i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
204 %! 6.33946547999647e+02 + 3.02965459323125e+02i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
205 %! 2.19834747595065e-04 - 9.25266900230165e-06i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
206 %! 8.49515487435091e-09 - 2.95133588338825e-09i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
207 %! 2.96635342439717e-13 - 1.85401806861382e-13i;
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
208 %! 9.65535916388246e+19 + 3.78654062133933e+19i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
209 %! 3.38477774418380e+11 + 8.37063899960569e+10i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
210 %! 1.57615042657685e+03 - 4.33777639047543e+02i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
211 %! 2.36176542789578e-05 - 5.75861972980636e-04i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
212 %! -6.83624588479039e-09 - 1.47230889442175e-08i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
213 %! -2.93020801760942e-13 - 4.03912221595793e-13i;
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
214 %! -1.94572937469407e+19 - 1.03494929263031e+20i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
215 %! -4.22385087573180e+10 - 3.61103191095041e+11i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
216 %! 4.89771220858552e+02 - 2.09175729060712e+03i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
217 %! 7.26650666035639e-04 + 4.71027801635222e-04i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
218 %! 1.02146578536128e-08 + 1.51813977370467e-08i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
219 %! 2.41628751621686e-13 + 4.66309048729523e-13i;
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
220 %! 5.42351559144068e+19 + 8.54503231614651e+19i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
221 %! 1.22886461074544e+11 + 3.03555953589323e+11i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
222 %! -2.13050339387819e+02 + 1.23853666784218e+03i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
223 %! -3.68087391884738e-04 + 1.94003994408861e-04i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
224 %! -1.39355838231763e-08 + 6.57189276453356e-10i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
225 %! -4.55133112151501e-13 - 8.46035902535333e-14i;
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
226 %! -7.75482228205081e+19 - 5.36017490438329e+19i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
227 %! -1.85284579257329e+11 - 2.08761110392897e+11i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
228 %! -1.74210199269860e+02 - 8.09467914953486e+02i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
229 %! 9.40470496160143e-05 - 2.44265223110736e-04i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
230 %! 6.64487526601190e-09 - 7.87242868014498e-09i, ...
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
231 %! 3.10273337426175e-13 - 2.28030229776792e-13i];
24908
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
232 %! assert (expint (X), y_exp, -1e-14);
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
233
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
234 ## Exceptional values (-Inf, Inf, NaN, 0, 0.37250741078)
19593
446c46af4b42 strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents: 17744
diff changeset
235 %!test
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
236 %! x = [-Inf; Inf; NaN; 0; -0.3725074107813668];
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
237 %! y_exp = [-Inf - i*pi; 0; NaN; Inf; 0 - i*pi];
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
238 %! y = expint (x);
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
239 %! assert (y, y_exp, 5*eps);
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
240
25015
baa7e37453b1 Added more tests for betainc and expint.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24927
diff changeset
241 %!test <53351>
baa7e37453b1 Added more tests for betainc and expint.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24927
diff changeset
242 %! assert (expint (32.5 + 1i),
baa7e37453b1 Added more tests for betainc and expint.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24927
diff changeset
243 %! 1.181108930758065e-16 - 1.966348533426658e-16i, -4*eps);
baa7e37453b1 Added more tests for betainc and expint.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24927
diff changeset
244 %! assert (expint (44 + 1i),
baa7e37453b1 Added more tests for betainc and expint.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24927
diff changeset
245 %! 9.018757389858152e-22 - 1.475771020004195e-21i, -4*eps);
baa7e37453b1 Added more tests for betainc and expint.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24927
diff changeset
246
baa7e37453b1 Added more tests for betainc and expint.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24927
diff changeset
247 %!test <47738>
baa7e37453b1 Added more tests for betainc and expint.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24927
diff changeset
248 %! assert (expint (10i), 0.0454564330044554 + 0.0875512674239774i, -4*eps);
baa7e37453b1 Added more tests for betainc and expint.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24927
diff changeset
249
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
250 ## Test preservation or conversion of the class
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
251 %!assert (class (expint (single (1))), "single")
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
252 %!assert (class (expint (int8 (1))), "double")
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
253 %!assert (class (expint (int16 (1))), "double")
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
254 %!assert (class (expint (int32 (1))), "double")
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
255 %!assert (class (expint (int64 (1))), "double")
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
256 %!assert (class (expint (uint8 (1))), "double")
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
257 %!assert (class (expint (uint16 (1))), "double")
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
258 %!assert (class (expint (uint32 (1))), "double")
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
259 %!assert (class (expint (uint64 (1))), "double")
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
260 %!assert (issparse (expint (sparse (1))))
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
261
24908
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
262 ## Test on the correct Image set
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
263 %!assert (isreal (expint (linspace (0, 100))))
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
264 %!assert (! isreal (expint (-1)))
24908
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
265
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
266 ## Test input validation
16585
1a3bfb14b5da Add and fix tests for ellipj, ellipke, and expint
Mike Miller <mtmiller@ieee.org>
parents: 16584
diff changeset
267 %!error expint ()
1a3bfb14b5da Add and fix tests for ellipj, ellipke, and expint
Mike Miller <mtmiller@ieee.org>
parents: 16584
diff changeset
268 %!error expint (1,2)
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
269 %!error <X must be numeric> expint ("1")