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