Mercurial > octave
changeset 23080:c56d30ea6cd4
expint.m: new strategies to compute the exponential integral (bug #47738).
* expint.m: output is now computed using three different strategies depending
on the coordinates of the input in the complex plane. The strategies are:
1) series expansion, 2) continued fraction, and 3) asymptotic series expansion.
New BIST tests added.
author | Michele Ginesi <michele.ginesi@gmail.com> |
---|---|
date | Tue, 10 Jan 2017 19:35:10 +0100 |
parents | 17dc6c7ef427 |
children | 7485462a6924 |
files | scripts/specfun/expint.m |
diffstat | 1 files changed, 184 insertions(+), 174 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/specfun/expint.m Sat Jan 21 16:19:38 2017 -0800 +++ b/scripts/specfun/expint.m Tue Jan 10 19:35:10 2017 +0100 @@ -1,4 +1,5 @@ -## Copyright (C) 2006-2016 Sylvain Pelissier +## Copyright (C) 2006, 2013 Sylvain Pelissier +## Copyright (C) 2017 Michele Ginesi ## ## This file is part of Octave. ## @@ -16,7 +17,8 @@ ## along with Octave; see the file COPYING. If not, see ## <http://www.gnu.org/licenses/>. -## Author: Sylvain Pelissier <sylvain.pelissier@gmail.com> +## Authors: Sylvain Pelissier <sylvain.pelissier@gmail.com> +## Michele Ginesi <michele.ginesi@gmail.com> ## -*- texinfo -*- ## @deftypefn {} {} expint (@var{x}) @@ -30,9 +32,11 @@ ## ## @example ## @group -## infinity +## +oo ## / -## E_1 (x) = | exp (-t)/t dt +## | exp (-t) +## E_1 (x) = | -------- dt +## | t ## / ## x ## @end group @@ -51,9 +55,11 @@ ## ## @example ## @group -## infinity +## +oo ## / -## Ei (x) = - | exp (-t)/t dt +## | exp (-t) +## Ei (x) = - | -------- dt +## | t ## / ## -x ## @end group @@ -71,197 +77,201 @@ ## @end ifnottex ## @end deftypefn -function y = expint (x) +function E1 = expint (x) if (nargin != 1) print_usage (); endif - y = x; # Copy over all values, including NaNs - - if (isreal (x)) - idx = (x >= 0); - y(idx) = -expint_Ei (-x(idx)); - - idx = (x < 0); - y(idx) = -expint_Ei (-x(idx)) - i*pi; - else - idx = (imag (x) > 0); - y(idx) = -expint_Ei (-x(idx)) - i*pi; - - idx = (imag (x) < 0); - y(idx) = -expint_Ei (-x(idx)) + i*pi; - - isreal_idx = (imag (x) == 0); - idx = (isreal_idx & real (x) >= 0); - y(idx) = -expint_Ei (-x(idx)); - - idx = (isreal_idx & real (x) < 0); - y(idx) = -expint_Ei (-x(idx)) - i*pi; + if (! isnumeric (x)) + error ("expint: X must be numeric"); + elseif (isinteger (x)) + x = double (x); endif -endfunction + sparse_x = 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_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); -## -*- texinfo -*- -## @deftypefn {} {@var{y} =} expint_Ei (@var{x}) -## Compute the exponential integral: -## @verbatim -## infinity -## / -## expint_Ei (x) = - | exp(-t)/t dt -## / -## -x -## @end verbatim -## @end deftypefn + ## FIXME: The performance of these routines need improvement. + ## There are numerous temporary variables created, some of which could + ## probably be eliminated. -function y = expint_Ei (x) - - y = zeros (size (x)); - F = @(x) exp (-x)./x; + ## 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; + fflag = true (size (res)); + while (k < 1e3 && any (fflag)) + res_tmp = res(fflag); + x_s_tmp = x_s(fflag); + ssum_tmp = ssum(fflag); + res_tmp .*= (k * -x_s_tmp/((k + 1)^2)); + ssum_tmp += res_tmp; + k += 1; + res(fflag) = res_tmp; + ssum(fflag) = ssum_tmp; + x_s(fflag) = x_s_tmp; + fflag = abs (res) > tol*abs (ssum); + endwhile + e1_s -= ssum; - for t = 1:numel (x) - xt = x(t); - if (xt < 0 && imag (xt) == 0) - ## Direct integration for most real inputs - y(t) = -quad (F, -xt, Inf, [0, 1e-10]); - elseif (xt > 2 && imag (xt) == 0) - persistent Ei_2 = 4.954234356001890; - y(t) = Ei_2 - quad (F, -xt, -2); - elseif (abs (xt) < 10) - ## Series Expansion for real (range [0,2]) or complex inputs (r < 10) - k = 1; - do - term = xt^k / (k*factorial (k)); - y(t) += term; - until (abs (term) < eps (abs (y(t))) / 2 || k++ >= 100) - y(t) = 0.57721566490153286 + log (xt) + y(t); - else - ## FIXME: This expansion is accurate to only 1e-13 at the beginning - ## near 10+i, although it becomes more accurate as the magnitude - ## of xt grows. - if (imag (xt) <= 0) - persistent a1 = 4.03640; - persistent a2 = 1.15198; - persistent b1 = 5.03637; - persistent b2 = 4.19160; - y(t) = -(xt^2 - a1*xt + a2) ... - / ((xt^2 - b1*xt + b2) * (-xt) * exp (-xt)) ... - - i*pi; - else - y(t) = conj (expint_Ei (conj (xt))); - endif; - endif - endfor + ## Continued fraction, + ## Abramowitz, Stegun, "Handbook of Mathematical Functions", + ## formula 5.1.22, p 229. + ## modified Lentz's algorithm, from "Numerical recipes in Fortran 77" p.165. + f_new = 2^-100 * ones (size (x_cf), class (x_cf)); + C_new = f_new; + D_new = zeros (size (x_cf), class (x_cf)); + k = 0; + fflag = true (size (x_cf)); + Delta = C_old = D_old = f_old = ones (size (x_cf), class (x_cf)); + while (k < 1e3 && any (fflag)) + x_cf_tmp = x_cf(fflag); + C_new_tmp = C_new(fflag); + D_new_tmp = D_new(fflag); + f_old = f_new(fflag); + C_old = C_new_tmp; + D_old = D_new_tmp; + b = x_cf_tmp*(mod (k,2) == 0) + (mod (k,2) == 1); + a = exp (-x_cf_tmp)*(k == 0) + ceil (k/2)*(k >= 1); + D_new_tmp = b + a.*D_old; + D_new_tmp = D_new_tmp.*(D_new_tmp != 0) + 2^-100*(D_new_tmp == 0); + C_new_tmp = b + a./C_old; + C_new_tmp = C_new_tmp.*(C_new_tmp != 0) + 2^-100*(C_new_tmp == 0); + D_new_tmp = 1 ./ D_new_tmp; + Delta(fflag) = C_new_tmp.*D_new_tmp; + x_cf(fflag) = x_cf_tmp; + f_new(fflag) = f_old.*Delta(fflag); + C_new(fflag) = C_new_tmp; + D_new(fflag) = D_new_tmp; + fflag = abs (Delta-1) > tol; + k += 1; + endwhile + + e1_cf = f_new; + ## Asymptotic series, from Fikioris, Tastsoglou, Bakas, + ## "Selected Asymptotic Methods with Application to Magnetics and Antennas" + ## formula A.10 p 161. + e1_a = exp (-x_a) ./ x_a; + oldres = ssum = res = ones (size (x_a), class (x_a)); + k = 0; + fflag = true (size (x_a)); + while (k < 1e3 && any (fflag)) + res_tmp = res(fflag); + oldres_tmp = res_tmp; + x_a_tmp = x_a(fflag); + res_tmp ./= (-x_a_tmp / (k+1)); + ssum(fflag) += res_tmp; + k += 1; + res(fflag) = res_tmp; + oldres(fflag) = oldres_tmp; + fflag = abs (oldres) > abs (res) + 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; + E1 = reshape (E1, orig_sz); + if (sparse_x) + E1 = sparse (E1); # if input was sparse format, output should be too. + endif endfunction -## Test against A&S Table 5.1 -%!test -%! x = [5:5:50]'/100; -%! gamma = 0.5772156649; -%! y_exp = [0.9876375971; -%! 0.9755453033; -%! 0.9637156702; -%! 0.9521414833; -%! 0.9408157528; -%! 0.9297317075; -%! 0.9188827858; -%! 0.9082626297; -%! 0.8978650778; -%! 0.8876841584 ]; -%! y = (expint (x) + log(x) + gamma) ./ x; -%! assert (y, y_exp, 1e-9); -%!test -%! x = [50:5:95]'/100; -%! y_exp = [0.559773595; -%! 0.503364081; -%! 0.454379503; -%! 0.411516976; -%! 0.373768843; -%! 0.340340813; -%! 0.310596579; -%! 0.284019269; -%! 0.260183939; -%! 0.238737524 ]; -%! y = expint (x); -%! assert (y, y_exp, 1e-9); -%!test -%! x = [100:5:145]'/100; -%! y_exp = [0.219383934; -%! 0.201872813; -%! 0.185990905; -%! 0.171555354; -%! 0.158408437; -%! 0.146413373; -%! 0.135450958; -%! 0.125416844; -%! 0.116219313; -%! 0.107777440 ]; -%! y = expint (x); -%! assert (y, y_exp, 1e-9); +## The following values were computed with the Octave symbolic package %!test -%! x = [150:5:200]'/100; -%! y_exp = [0.100019582; -%! 0.092882108; -%! 0.086308334; -%! 0.080247627; -%! 0.074654644; -%! 0.069488685; -%! 0.064713129; -%! 0.060294967; -%! 0.056204378; -%! 0.052414380; -%! 0.048900511 ]; -%! y = expint (x); -%! assert (y, y_exp, 1e-9); - -## Series expansion (-2 < x < 0) -## Expected values from Mathematica -%!test -%! x = [-0.1; -0.5; -1; -1.5; -2]; -%! y_exp = [ 1.6228128139692767 - i*pi; -%! -0.45421990486317358 - i*pi; -%! -1.8951178163559368 - i*pi; -%! -3.3012854491297978 - i*pi; -%! -4.9542343560018902 - i*pi]; -%! y = expint (x); -%! assert (y, y_exp, eps (real (y_exp))); - -## (x < -2, x real) -%!test -%! x = [-2.5; -3; -10;-15; -25]; -%! y_exp = [-7.0737658945786007 - i*pi; -%! -9.9338325706254165 - i*pi; -%! -2492.2289762418777 - i*pi; -%! -234955.85249076830 - i*pi; -%! -3.0059509065255486e9 - i*pi]; -%! y = expint (x); -%! assert (y, y_exp, 8*eps (real (y_exp))); - -## Complex values -%!test -%! x = [i; -1-i; 10-i; 10+i]; -%! y_exp = [-0.33740392290096813 - i*0.62471325642771360; -%! -1.7646259855638540 + i*0.75382280207927082; -%! 1.90746381979783120e-6 + i*3.67354374003294739e-6; -%! 1.90746381979783120e-6 - i*3.67354374003294739e-6]; -%! y = expint (x); -%! assert (y, y_exp, 1e-12); +%! 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-12); ## Exceptional values (-Inf, Inf, NaN, 0, 0.37250741078) %!test %! x = [-Inf; Inf; NaN; 0; -0.3725074107813668]; -%! y_exp = [-Inf - i*pi; -%! -Inf; # should be 0; -%! NaN; -%! Inf; -%! 0 - i*pi]; +%! y_exp = [-Inf - i*pi; 0; NaN; Inf; 0 - i*pi]; %! y = expint (x); %! assert (y, y_exp, 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 input validation %!error expint () %!error expint (1,2) +%!error <X must be numeric> expint ("1") -