Mercurial > octave
annotate scripts/specfun/expint.m @ 24927:c280560d9c96 stable
Overhaul special functions modified by GSOC2018 project.
* NEWS: Add note about new functions added. Add note explaining the changes
done to existing functions.
* __betainc__.cc: Renamed from __betainc_lentz__.cc. Use standard GPL v3
copyright block. Add missing #include "dNDArray.h". Add one-line texinfo
documentation for internal function. Remove fourth input argument to function.
Use is_single_type() to decide whether to operate with FloatNDArray or NDArray.
Delete temporary variables x_arg_s, a_arg_s, b_arg_s used in input validation.
Rename len_x, len_a, len_b to numel_[xab] for clarity. Remove input validation
that numel match (internal function, no need). Rename function output to
retval according to Octave coding conventions. Use more spacing and newlines
for readability of code.
* __expint__.cc: Renamed from __expint_lentz__.cc. Use standard GPL v3
copyright block. Cull list of #includes to just the necessary ones. Use
Complex and FloatComplex typedefs defined by Octave. Add one-line texinfo
documentation for internal function. Remove second input argument to function.
Use is_single_type() to decide whether to operate with FloatComplexNDArray or
ComplexNDArray. Delete temporary variables x_arg_s used in input validation.
Rename len_x to numel_x for clarity. Use constructor with dim_vector and
scalar value rather than fill() after creating array. Rename function output
to retval according to Octave coding conventions. Use more spacing and
newlines for readability of code.
* __gammainc__.cc: Renamed from __gammainc_lentz__.cc. Use standard GPL v3
copyright block. Add one-line texinfo documentation for internal function.
Remove third input argument to function. Remove input validation that numel
match (internal function, no need). Use is_single_type() to decide whether to
operate with FloatNDArray or NDArray. Delete temporary variables x_arg, a_arg
used in input validation. Use constructor with dim_vector and scalar value
rather than fill() after creating array. Rename function output to retval
according to Octave coding conventions. Use more spacing and newlines for
readability of code.
* __betainc_lentz__.cc, __expint_lentz__.cc, __gammainc_lentz__.cc: Removed.
* libinterp/corefcn/module.mk: Add renamed functions to build system.
* betainc.m: Use Octave standard GPL block. Rewrite parts of docstring.
Don't use array brackets around single output of function. Remove isscalar
checks on inputs because common_size() function will already handle it.
Use capital variable names in error messages to match documentation as displayed
in terminal. Reshape all inputs in to column vectors quickly so that input
validation tests that depend on all/any will pass with N-D arrays. Add
comments to code. Check for specific error messages in input validation BIST
tests.
* betaincinv.m: Use Octave standard GPL block. Rewrite parts of docstring.
Don't use array brackets around single output of function. Remove isscalar
checks on inputs because common_size() function will already handle it.
Use capital variable names in error messages to match documentation as displayed
in terminal. Reshape all inputs in to column vectors quickly so that input
validation tests that depend on all/any will pass with N-D arrays. Put most
common case of tail ("lower") first in if/elseif trees. Call functions directly
with function handle rather than using unnecessary feval() call. Use numel
in preference to length. Rename variable i_miss to todo for clarity. Add
comments to code. Check for specific error messages in input validation BIST
tests.
* cosint.m: Use Octave standard GPL block. Rewrite parts of docstring.
Don't use array brackets around single output of function. Remove isscalar
checks on inputs because common_size() function will already handle it.
Add input validation check for isnumeric value. Convert integer classes to
double before proceeding. Rename i_miss to todo for clarity. Use isinf to
detect both -Inf and +Inf rather than separate tests. Use ++it in while
loop conditional to shorten loop blocks. Add Input validation BIST tests.
* expint.m: Remove Sylvain who was not actually an author on this file.
Rewrite parts of docstring. Rename variable sparse_x to orig_sparse.
Eliminate temporary variables res_tmp, x_s_tmp, ssum_tmp. Rename i_miss to
todo. Use Octave coding conventions throughout. Add comments to code.
* gammainc.m: Use Octave standard GPL block. Rewrite parts of docstring.
Remove isscalar checks on inputs because common_size() function will already
handle it. Use capital variable names in error messages to match documentation
as displayed in terminal. Reshape all inputs in to column vectors quickly so
that input validation tests that depend on all/any will pass with N-D arrays.
Put most common case of tail ("lower") first in if/elseif trees. Add input
validation of tail. Rename variable ii to idx for clarity. Rename variable
i_done to todo and switch polarity so that the '!' operator is not required
every time the variable is updated. Use indexing and direct assignment to
update todo rather than logical operator '&' which is slower. Use tolower on
tail variable and then switch strcmpi calls to strcmp. Reformat %!test blocks
in to %!assert blocks to be more compact. Check for specific error messages in
input validation BIST tests.
* gammaincinv.m: Use Octave standard GPL block. Rewrite parts of docstring.
Don't use array brackets around single output of function. Remove isscalar
checks on inputs because common_size() function will already handle it.
Add input validation check for iscomplex value. Use capital variable names in
error messages to match documentation as displayed in terminal. Reshape all
inputs in to column vectors quickly so that input validation tests that depend
on all/any will pass with N-D arrays. Rename i_miss to todo. Use numel
in preference to length. Call functions directly with function handle rather
than using unnecessary feval() call.Rename variable i_miss to todo for clarity.
Use it++ in while loop conditional to shorten loop blocks. Add comments to
code. Check for specific error messages in input validation BIST tests.
Add input validation BIST tests for all error messages.
* sinint.m: Use Octave standard GPL block. Rewrite parts of docstring.
Add input validation for isnumeric. Convert integers to double for
calculation. Reshape input to column vector. Rename variable sz to orig_sz
for clarity. rename i_miss to todo. Reformat BIST tests to mak them more
compact. Add input validation BIST tests.
author | Rik <rik@octave.org> |
---|---|
date | Mon, 19 Mar 2018 10:01:48 -0700 |
parents | 40ab8be7d7ec |
children | baa7e37453b1 |
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 | 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 | 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 | 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 | 37 ## / |
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 | 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 | 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 | 47 ## @tex |
48 ## $$ | |
49 ## {\rm Ei} (x) = - \int_{-x}^\infty {e^{-t} \over t} dt. | |
50 ## $$ | |
51 ## @end tex | |
52 ## @ifnottex | |
53 ## | |
54 ## @example | |
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 | 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 | 61 ## / |
62 ## -x | |
63 ## @end group | |
64 ## @end example | |
17513
fedcd3717ebc
doc: grammarcheck of documentation before 3.8 release.
Rik <rik@octave.org>
parents:
17363
diff
changeset
|
65 ## |
17363 | 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 | 68 ## @tex |
69 ## $ | |
70 ## E_1 (-x) = -{\rm Ei} (x) - i\pi. | |
71 ## $ | |
72 ## @end tex | |
73 ## @ifnottex | |
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 | 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 | 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 | 233 |
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 | 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 | 238 %! y = expint (x); |
239 %! assert (y, y_exp, 5*eps); | |
240 | |
23080
c56d30ea6cd4
expint.m: new strategies to compute the exponential integral (bug #47738).
Michele Ginesi <michele.ginesi@gmail.com>
parents:
22755
diff
changeset
|
241 ## 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
|
242 %!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
|
243 %!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
|
244 %!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
|
245 %!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
|
246 %!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
|
247 %!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
|
248 %!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
|
249 %!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
|
250 %!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
|
251 %!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
|
252 |
24908
6c082a43abd8
expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents:
24534
diff
changeset
|
253 ## 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
|
254 %!assert (isreal (expint (linspace (0, 100)))) |
24927
c280560d9c96
Overhaul special functions modified by GSOC2018 project.
Rik <rik@octave.org>
parents:
24923
diff
changeset
|
255 %!assert (! isreal (expint (-1))) |
24908
6c082a43abd8
expint: moved the Lentz algorithm to .cc function.
Michele Ginesi <michele.ginesi@gmail.com>
parents:
24534
diff
changeset
|
256 |
17363 | 257 ## Test input validation |
16585
1a3bfb14b5da
Add and fix tests for ellipj, ellipke, and expint
Mike Miller <mtmiller@ieee.org>
parents:
16584
diff
changeset
|
258 %!error expint () |
1a3bfb14b5da
Add and fix tests for ellipj, ellipke, and expint
Mike Miller <mtmiller@ieee.org>
parents:
16584
diff
changeset
|
259 %!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
|
260 %!error <X must be numeric> expint ("1") |