Mercurial > octave-libtiff
view scripts/specfun/expint.m @ 30564:796f54d4ddbf stable
update Octave Project Developers copyright for the new year
In files that have the "Octave Project Developers" copyright notice,
update for 2021.
In all .txi and .texi files except gpl.txi and gpl.texi in the
doc/liboctave and doc/interpreter directories, change the copyright
to "Octave Project Developers", the same as used for other source
files. Update copyright notices for 2022 (not done since 2019). For
gpl.txi and gpl.texi, change the copyright notice to be "Free Software
Foundation, Inc." and leave the date at 2007 only because this file
only contains the text of the GPL, not anything created by the Octave
Project Developers.
Add Paul Thomas to contributors.in.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Tue, 28 Dec 2021 18:22:40 -0500 |
parents | 2d17a87740dd |
children | 5d3faba0342e |
line wrap: on
line source
######################################################################## ## ## Copyright (C) 2018-2022 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. ## ## This file is part of Octave. ## ## Octave is free software: you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation, either version 3 of the License, or ## (at your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## <https://www.gnu.org/licenses/>. ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {} expint (@var{x}) ## Compute the exponential integral. ## ## The exponential integral is defined as: ## ## @tex ## $$ ## {\rm E_1} (x) = \int_x^\infty {e^{-t} \over t} dt ## $$ ## @end tex ## @ifnottex ## ## @example ## @group ## +oo ## / ## | exp (-t) ## E_1 (x) = | -------- dt ## | t ## / ## x ## @end group ## @end example ## ## @end ifnottex ## ## Note: For compatibility, this function uses the @sc{matlab} definition ## of the exponential integral. Most other sources refer to this particular ## value as @math{E_1 (x)}, and the exponential integral as ## @tex ## $$ ## {\rm Ei} (x) = - \int_{-x}^\infty {e^{-t} \over t} dt. ## $$ ## @end tex ## @ifnottex ## ## @example ## @group ## +oo ## / ## | exp (-t) ## Ei (x) = - | -------- dt ## | t ## / ## -x ## @end group ## @end example ## ## @end ifnottex ## The two definitions are related, for positive real values of @var{x}, by ## @tex ## $ ## E_1 (-x) = -{\rm Ei} (x) - i\pi. ## $ ## @end tex ## @ifnottex ## @w{@code{E_1 (-x) = -Ei (x) - i*pi}}. ## @end ifnottex ## ## References: ## ## @nospell{M. Abramowitz and I.A. Stegun}, ## @cite{Handbook of Mathematical Functions}, 1964. ## ## @nospell{N. Bleistein and R.A. Handelsman}, ## @cite{Asymptotic expansions of integrals}, 1986. ## ## @seealso{cosint, sinint, exp} ## @end deftypefn function E1 = expint (x) if (nargin < 1) print_usage (); endif if (! isnumeric (x)) error ("expint: X must be numeric"); endif ## Convert to floating point if necessary if (isinteger (x)) x = double (x); endif orig_sparse = issparse (x); orig_sz = size (x); x = x(:); # convert to column vector ## Initialize the result if (isreal (x) && x >= 0) E1 = zeros (size (x), class (x)); else E1 = complex (zeros (size (x), class (x))); endif tol = eps (class (x)); ## Divide the input into 3 regions and apply a different algorithm for each. ## s = series expansion, cf = continued fraction, a = asymptotic series s_idx = (((real (x) + 19.5).^ 2 ./ (20.5^2) + ... imag (x).^2 ./ (10^2)) <= 1) ... | (real (x) < 0 & abs (imag (x)) <= 1e-8); cf_idx = ((((real (x) + 1).^2 ./ (38^2) + ... imag (x).^2 ./ (40^2)) <= 1) ... & (! s_idx)) & (real (x) <= 35); a_idx = (! s_idx) & (! cf_idx); x_s = x(s_idx); x_cf = x(cf_idx); x_a = x(a_idx); ## Series expansion ## Abramowitz, Stegun, "Handbook of Mathematical Functions", ## formula 5.1.11, p 229. ## FIXME: Why so long? IEEE double doesn't have this much precision. gm = 0.577215664901532860606512090082402431042159335; e1_s = -gm - log (x_s); res = -x_s; ssum = res; k = 1; todo = true (size (res)); while (k < 1e3 && any (todo)) res(todo) .*= (k * (- x_s(todo)) / ((k + 1) ^ 2)); ssum(todo) += res(todo); k += 1; todo = (abs (res) > (tol * abs (ssum))); endwhile e1_s -= ssum; ## Continued fraction expansion, ## Abramowitz, Stegun, "Handbook of Mathematical Functions", ## formula 5.1.22, p 229. ## Modified Lentz's algorithm, from "Numerical recipes in Fortran 77" p.165. e1_cf = exp (-x_cf) .* __expint__ (x_cf); ## Remove spurious imaginary part if needed (__expint__ works automatically ## with complex values) if (isreal (x_cf) && x_cf >= 0) e1_cf = real (e1_cf); endif ## Asymptotic series, from N. Bleistein and R.A. Handelsman ## "Asymptotic expansion of integrals", pages 1-4. e1_a = exp (-x_a) ./ x_a; ssum = res = ones (size (x_a), class (x_a)); k = 0; todo = true (size (x_a)); while (k < 1e3 && any (todo)) res(todo) ./= (- x_a(todo) / (k + 1)); ssum(todo) += res(todo); k += 1; todo = abs (x_a) > k; endwhile e1_a .*= ssum; ## Combine results from each region into final output E1(s_idx) = e1_s; E1(cf_idx) = e1_cf; E1(a_idx) = e1_a; ## Restore shape and sparsity of input E1 = reshape (E1, orig_sz); if (orig_sparse) E1 = sparse (E1); endif endfunction ## The following values were computed with the Octave symbolic package %!test %! X = [-50 - 50i -30 - 50i -10 - 50i 5 - 50i 15 - 50i 25 - 50i %! -50 - 30i -30 - 30i -10 - 30i 5 - 30i 15 - 30i 25 - 30i %! -50 - 10i -30 - 10i -10 - 10i 5 - 10i 15 - 10i 25 - 10i %! -50 + 5i -30 + 5i -10 + 5i 5 + 5i 15 + 5i 25 + 5i %! -50 + 15i -30 + 15i -10 + 15i 5 + 15i 15 + 15i 25 + 15i %! -50 + 25i -30 + 25i -10 + 25i 5 + 25i 15 + 25i 25 + 25i]; %! y_exp = [ -3.61285286166493e+19 + 6.46488018613387e+19i, ... %! -4.74939752018180e+10 + 1.78647798300364e+11i, ... %! 3.78788822381261e+01 + 4.31742823558278e+02i, ... %! 5.02062497548626e-05 + 1.23967883532795e-04i, ... %! 3.16785290137650e-09 + 4.88866651583182e-09i, ... %! 1.66999261039533e-13 + 1.81161508735941e-13i; %! 3.47121527628275e+19 + 8.33104448629260e+19i, ... %! 1.54596484273693e+11 + 2.04179357837414e+11i, ... %! 6.33946547999647e+02 + 3.02965459323125e+02i, ... %! 2.19834747595065e-04 - 9.25266900230165e-06i, ... %! 8.49515487435091e-09 - 2.95133588338825e-09i, ... %! 2.96635342439717e-13 - 1.85401806861382e-13i; %! 9.65535916388246e+19 + 3.78654062133933e+19i, ... %! 3.38477774418380e+11 + 8.37063899960569e+10i, ... %! 1.57615042657685e+03 - 4.33777639047543e+02i, ... %! 2.36176542789578e-05 - 5.75861972980636e-04i, ... %! -6.83624588479039e-09 - 1.47230889442175e-08i, ... %! -2.93020801760942e-13 - 4.03912221595793e-13i; %! -1.94572937469407e+19 - 1.03494929263031e+20i, ... %! -4.22385087573180e+10 - 3.61103191095041e+11i, ... %! 4.89771220858552e+02 - 2.09175729060712e+03i, ... %! 7.26650666035639e-04 + 4.71027801635222e-04i, ... %! 1.02146578536128e-08 + 1.51813977370467e-08i, ... %! 2.41628751621686e-13 + 4.66309048729523e-13i; %! 5.42351559144068e+19 + 8.54503231614651e+19i, ... %! 1.22886461074544e+11 + 3.03555953589323e+11i, ... %! -2.13050339387819e+02 + 1.23853666784218e+03i, ... %! -3.68087391884738e-04 + 1.94003994408861e-04i, ... %! -1.39355838231763e-08 + 6.57189276453356e-10i, ... %! -4.55133112151501e-13 - 8.46035902535333e-14i; %! -7.75482228205081e+19 - 5.36017490438329e+19i, ... %! -1.85284579257329e+11 - 2.08761110392897e+11i, ... %! -1.74210199269860e+02 - 8.09467914953486e+02i, ... %! 9.40470496160143e-05 - 2.44265223110736e-04i, ... %! 6.64487526601190e-09 - 7.87242868014498e-09i, ... %! 3.10273337426175e-13 - 2.28030229776792e-13i]; %! assert (expint (X), y_exp, -1e-14); ## Exceptional values (-Inf, Inf, NaN, 0, 0.37250741078) %!test %! x = [-Inf; Inf; NaN; 0; -0.3725074107813668]; %! y_exp = [-Inf - i*pi; 0; NaN; Inf; 0 - i*pi]; %! y = expint (x); %! assert (y, y_exp, 5*eps); %!test <*53351> %! assert (expint (32.5 + 1i), %! 1.181108930758065e-16 - 1.966348533426658e-16i, -4*eps); %! assert (expint (44 + 1i), %! 9.018757389858152e-22 - 1.475771020004195e-21i, -4*eps); %!test <*47738> %! assert (expint (10i), 0.0454564330044554 + 0.0875512674239774i, -5*eps); ## Test preservation or conversion of the class %!assert (class (expint (single (1))), "single") %!assert (class (expint (int8 (1))), "double") %!assert (class (expint (int16 (1))), "double") %!assert (class (expint (int32 (1))), "double") %!assert (class (expint (int64 (1))), "double") %!assert (class (expint (uint8 (1))), "double") %!assert (class (expint (uint16 (1))), "double") %!assert (class (expint (uint32 (1))), "double") %!assert (class (expint (uint64 (1))), "double") %!assert (issparse (expint (sparse (1)))) ## Test on the correct Image set %!assert (isreal (expint (linspace (0, 100)))) %!assert (! isreal (expint (-1))) ## Test input validation %!error <Invalid call> expint () %!error <X must be numeric> expint ("1")