Mercurial > octave

annotate scripts/specfun/expint.m @ 29789:c8585732ce08 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
tests: Relax tolerance for tests on macOS. * libinterp/corefcn/mappers.cc (asin): Use larger tolerance for tests on macOS. * scripts/specfun/expint.m: Slightly relax tolerance for the test to also pass on macOS.
author Markus Mützel <markus.muetzel@gmx.de>
date Sun, 20 Jun 2021 10:03:38 +0200
parents 0a5b15007766
children 2d17a87740dd
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
1 ########################################################################
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
2 ##
29358
0a5b15007766 update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 27923
diff changeset
3 ## Copyright (C) 2018-2021 The Octave Project Developers
27918
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 26946
diff changeset
4 ##
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
5 ## See the file COPYRIGHT.md in the top-level directory of this
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
6 ## distribution or <https://octave.org/copyright/>.
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
7 ##
16586
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
8 ## 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
9 ##
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23220
diff changeset
10 ## 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
11 ## 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
12 ## 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
13 ## (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
14 ##
16586
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
15 ## 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
16 ## WITHOUT ANY WARRANTY; without even the implied warranty of
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22323
diff changeset
17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22323
diff changeset
18 ## 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
19 ##
16586
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
20 ## 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
21 ## 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
22 ## <https://www.gnu.org/licenses/>.
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
23 ##
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
24 ########################################################################
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
25
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
26 ## -*- texinfo -*-
20852
516bb87ea72e 2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents: 20162
diff changeset
27 ## @deftypefn {} {} expint (@var{x})
25106
d7ad543255c5 doc: Shorten very long first sentences of docstrings (bug #53388).
Rik <rik@octave.org>
parents: 25015
diff changeset
28 ## Compute the exponential integral.
d7ad543255c5 doc: Shorten very long first sentences of docstrings (bug #53388).
Rik <rik@octave.org>
parents: 25015
diff changeset
29 ##
d7ad543255c5 doc: Shorten very long first sentences of docstrings (bug #53388).
Rik <rik@octave.org>
parents: 25015
diff changeset
30 ## The exponential integral is defined as:
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
31 ##
16587
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
32 ## @tex
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
33 ## $$
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
34 ## {\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
35 ## $$
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
36 ## @end tex
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
37 ## @ifnottex
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
38 ##
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
39 ## @example
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
40 ## @group
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
41 ## +oo
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
42 ## /
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
43 ## | exp (-t)
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
44 ## 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
45 ## | t
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
46 ## /
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
47 ## x
16587
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
48 ## @end group
a3fdd6041e64 Fix ellipj, ellipke, and expint docstrings
Mike Miller <mtmiller@ieee.org>
parents: 16586
diff changeset
49 ## @end example
17513
fedcd3717ebc doc: grammarcheck of documentation before 3.8 release.
Rik <rik@octave.org>
parents: 17363
diff changeset
50 ##
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
51 ## @end ifnottex
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
52 ##
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
53 ## Note: For compatibility, this function uses the @sc{matlab} definition
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
54 ## 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
55 ## 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
56 ## @tex
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
57 ## $$
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
58 ## {\rm Ei} (x) = - \int_{-x}^\infty {e^{-t} \over t} dt.
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
59 ## $$
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
60 ## @end tex
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
61 ## @ifnottex
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
62 ##
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
63 ## @example
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
64 ## @group
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
65 ## +oo
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
66 ## /
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
67 ## | exp (-t)
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
68 ## 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
69 ## | t
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
70 ## /
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
71 ## -x
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
72 ## @end group
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
73 ## @end example
17513
fedcd3717ebc doc: grammarcheck of documentation before 3.8 release.
Rik <rik@octave.org>
parents: 17363
diff changeset
74 ##
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
75 ## @end ifnottex
17514
5b916efea542 doc: spellcheck of documentation before 3.8 release.
Rik <rik@octave.org>
parents: 17513
diff changeset
76 ## 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
77 ## @tex
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
78 ## $
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
79 ## E_1 (-x) = -{\rm Ei} (x) - i\pi.
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
80 ## $
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
81 ## @end tex
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
82 ## @ifnottex
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
83 ## @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
84 ## @end ifnottex
24908
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
85 ##
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
86 ## References:
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
87 ##
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
88 ## @nospell{M. Abramowitz and I.A. Stegun},
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
89 ## @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
90 ##
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
91 ## @nospell{N. Bleistein and R.A. Handelsman},
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
92 ## @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
93 ##
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
94 ## @seealso{cosint, sinint, exp}
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
95 ## @end deftypefn
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
96
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
97 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
98
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
99 if (nargin != 1)
16586
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
100 print_usage ();
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
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
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
103 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
104 error ("expint: X must be numeric");
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
105 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
106
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
107 ## Convert to floating point if necessary
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
108 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
109 x = double (x);
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
110 endif
16586
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
111
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
112 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
113 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
114 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
115
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
116 ## Initialize the result
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
117 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
118 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
119 else
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
120 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
121 endif
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
122 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
123
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
124 ## 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
125 ## 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
126 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
127 imag (x).^2 ./ (10^2)) <= 1) ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
128 | (real (x) < 0 & abs (imag (x)) <= 1e-8);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
129 cf_idx = ((((real (x) + 1).^2 ./ (38^2) + ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
130 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
131 & (! s_idx)) & (real (x) <= 35);
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
132 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
133 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
134 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
135 x_a = x(a_idx);
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
136
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
137 ## Series expansion
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
138 ## 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
139 ## 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
140 ## 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
141 gm = 0.577215664901532860606512090082402431042159335;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
142 e1_s = -gm - log (x_s);
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
143 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
144 ssum = res;
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
145 k = 1;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
146 todo = true (size (res));
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
147 while (k < 1e3 && any (todo))
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
148 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
149 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
150 k += 1;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
151 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
152 endwhile
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
153 e1_s -= ssum;
16586
f423873d3275 Style fixes for ellipj.cc, ellipke.m, and expint.m
Mike Miller <mtmiller@ieee.org>
parents: 16585
diff changeset
154
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
155 ## 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
156 ## 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
157 ## formula 5.1.22, p 229.
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
158 ## 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
159
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
160 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
161
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
162 ## 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
163 ## with complex values)
24908
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
164
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
165 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
166 e1_cf = real (e1_cf);
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
167 endif
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
168
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
169 ## 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
170 ## "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
171 e1_a = exp (-x_a) ./ x_a;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
172 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
173 k = 0;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
174 todo = true (size (x_a));
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
175 while (k < 1e3 && any (todo))
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
176 res(todo) ./= (- x_a(todo) / (k + 1));
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
177 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
178 k += 1;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
179 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
180 endwhile
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
181 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
182
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
183 ## 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
184 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
185 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
186 E1(a_idx) = e1_a;
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
187
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
188 ## 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
189 E1 = reshape (E1, orig_sz);
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
190 if (orig_sparse)
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
191 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
192 endif
21758
ffad2baa90f7 maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents: 20852
diff changeset
193
16584
2f766ceeb03e Add ellipj, ellipke, and expint functions from Octave Forge
Mike Miller <mtmiller@ieee.org>
parents:
diff changeset
194 endfunction
16585
1a3bfb14b5da Add and fix tests for ellipj, ellipke, and expint
Mike Miller <mtmiller@ieee.org>
parents: 16584
diff changeset
195
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
196
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
197 ## 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
198 %!test
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
199 %! 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
200 %! -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
201 %! -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
202 %! -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
203 %! -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
204 %! -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
205 %! 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
206 %! -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
207 %! 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
208 %! 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
209 %! 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
210 %! 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
211 %! 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
212 %! 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
213 %! 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
214 %! 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
215 %! 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
216 %! 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
217 %! 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
218 %! 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
219 %! 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
220 %! 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
221 %! -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
222 %! -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
223 %! -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
224 %! -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
225 %! 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
226 %! 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
227 %! 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
228 %! 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
229 %! 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
230 %! 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
231 %! -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
232 %! -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
233 %! -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
234 %! -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
235 %! -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
236 %! -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
237 %! -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
238 %! 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
239 %! 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
240 %! 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
241 %! assert (expint (X), y_exp, -1e-14);
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
242
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
243 ## 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
244 %!test
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
245 %! 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
246 %! 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
247 %! y = expint (x);
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
248 %! assert (y, y_exp, 5*eps);
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
249
26946
04e5cb5e2cb3 update bug status in tests
John W. Eaton <jwe@octave.org>
parents: 26376
diff changeset
250 %!test <*53351>
25015
baa7e37453b1 Added more tests for betainc and expint.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24927
diff changeset
251 %! assert (expint (32.5 + 1i),
baa7e37453b1 Added more tests for betainc and expint.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24927
diff changeset
252 %! 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
253 %! assert (expint (44 + 1i),
baa7e37453b1 Added more tests for betainc and expint.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24927
diff changeset
254 %! 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
255
26946
04e5cb5e2cb3 update bug status in tests
John W. Eaton <jwe@octave.org>
parents: 26376
diff changeset
256 %!test <*47738>
29789
c8585732ce08 tests: Relax tolerance for tests on macOS.
Markus Mützel <markus.muetzel@gmx.de>
parents: 29358
diff changeset
257 %! assert (expint (10i), 0.0454564330044554 + 0.0875512674239774i, -5*eps);
25015
baa7e37453b1 Added more tests for betainc and expint.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24927
diff changeset
258
23080
c56d30ea6cd4 expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents: 22755
diff changeset
259 ## 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
260 %!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
261 %!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
262 %!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
263 %!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
264 %!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
265 %!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
266 %!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
267 %!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
268 %!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
269 %!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
270
24908
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
271 ## 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
272 %!assert (isreal (expint (linspace (0, 100))))
24927
c280560d9c96 Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents: 24923
diff changeset
273 %!assert (! isreal (expint (-1)))
24908
6c082a43abd8 expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents: 24534
diff changeset
274
17363
9aca7020c89f expint.m: Overhaul function.
Rik <rik@octave.org>
parents: 17338
diff changeset
275 ## Test input validation
16585
1a3bfb14b5da Add and fix tests for ellipj, ellipke, and expint
Mike Miller <mtmiller@ieee.org>
parents: 16584
diff changeset
276 %!error expint ()
1a3bfb14b5da Add and fix tests for ellipj, ellipke, and expint
Mike Miller <mtmiller@ieee.org>
parents: 16584
diff changeset
277 %!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
278 %!error <X must be numeric> expint ("1")