# HG changeset patch # User Michele Ginesi # Date 1484073310 -3600 # Tue Jan 10 19:35:10 2017 +0100 # Node ID a89ed6e4828432b4014a2cb9007c525cb73717e4 # Parent 8ba865cd470354a9afafe69144bc1e6c782ca93b expint.m: new strategies to compute the exponential integral (bug #47738) * expint now is computed using different strategies depending on the coordinates of the input in the complex plane. The strategies are now: ** series expansion ** continued fraction ** asymptotic series expansion diff -r 8ba865cd4703 -r a89ed6e48284 scripts/specfun/expint.m --- a/scripts/specfun/expint.m Mon Jan 09 18:05:13 2017 -0800 +++ b/scripts/specfun/expint.m Tue Jan 10 19:35:10 2017 +0100 @@ -1,267 +1,238 @@ -## Copyright (C) 2006-2016 Sylvain Pelissier -## -## 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 -## . - -## Author: Sylvain Pelissier - -## -*- texinfo -*- -## @deftypefn {} {} expint (@var{x}) -## Compute the exponential integral: -## @tex -## $$ -## {\rm E_1} (x) = \int_x^\infty {e^{-t} \over t} dt -## $$ -## @end tex -## @ifnottex -## -## @example -## @group -## infinity -## / -## E_1 (x) = | exp (-t)/t dt -## / -## x -## @end group -## @end example -## -## @end ifnottex -## Note: For compatibility, this functions 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 -## infinity -## / -## Ei (x) = - | exp (-t)/t dt -## / -## -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 -## @end deftypefn - -function y = 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; - endif - -endfunction - -## -*- 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 - -function y = expint_Ei (x) - - y = zeros (size (x)); - F = @(x) exp (-x)./x; - - 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 - -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); -%!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); - -## 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 = expint (x); -%! assert (y, y_exp, 5*eps); - -## Test input validation -%!error expint () -%!error expint (1,2) - - +## Copyright (C) 2006, 2013 Sylvain Pelissier +## Copyright (C) 2017 Michele Ginesi +## +## 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 +## . + +## Authors: Sylvain Pelissier +## Michele Ginesi + +## -*- texinfo -*- +## @deftypefn {Function File} {} expint (@var{z}) +## Compute the exponential integral: +## @tex +## $$ +## {\rm E_1} (z) = \int_z^\infty {e^{-t} \over t} dt +## $$ +## @end tex +## @ifnottex +## +## @example +## @group +## +oo +## / +## | exp(-t) +## E_1(z) = | -------- dt +## | t +## / +## z +## @end group +## @end example +## +## @end ifnottex +## Note: For compatibility, this functions uses the @sc{matlab} definition +## of the exponential integral. Most other sources refer to this particular +## value as @math{E_1 (z)}, and the exponential integral as +## @tex +## $$ +## {\rm Ei} (z) = - \int_{-z}^\infty {e^{-t} \over t} dt. +## $$ +## @end tex +## @ifnottex +## +## @example +## @group +## +oo +## / +## | exp(-t) +## Ei(z) = - | -------- dt +## | t +## / +## -z +## @end group +## @end example +## +## @end ifnottex +## The two definitions are related, for positive real values of @var{z}, by +## @tex +## $ +## E_1 (-z) = -{\rm Ei} (z) - i\pi. +## $ +## @end tex +## @ifnottex +## @w{@code{E_1 (-z) = -Ei (z) - i*pi}}. +## @end ifnottex +## @end deftypefn + +function [E1] = expint(z) + + if (nargin != 1) + print_usage (); + endif + + flag_sparse = issparse(z); # check if the input is sparse + sz = size (z); # Extract the size + z = z(:); # Write the input as column + # Control (and eventually conversion) of the class + if (ischar(z)) + error('expecting numeric argument') + elseif (isinteger (z) || islogical (z)) + z = double(z); + endif + # Initialize the result + if (isreal(z) && z >=0) + E1 = zeros (size (z), class (z)); + else + E1 = complex (zeros (size (z), class (z))); + endif + tol = eps (class (z)); + ## In the following, we extract the three zones in which we divide the input + flag_s = (((real (z) + 19.5).^2./(20.5^2) + imag (z).^2./(10^2)) <= 1) | ... + (real (z) < 0 & abs (imag (z)) <= 1e-08); + flag_cf = ((((real (z) + 1).^2./(38^2) + imag (z).^2./(40^2)) <= 1) & ... + (~flag_s)) & (real (z) <= 35); + flag_a = ~flag_s & ~flag_cf; + z_s = z(flag_s); + z_cf = z(flag_cf); + z_a = z(flag_a); + ## Series expansion + ## Abramowitz, Stegun 'Handbook of Mathematical function', + ## formula 5.1.11 p 229 + gm = 0.577215664901532860606512090082402431042159335; + e1_s = -gm - log (z_s); + res = -z_s; + ssum = res; + k = 1; + fflag = true (size (res)); + while (k < 1e3 && any(fflag)) + res_tmp = res(fflag); + z_s_tmp = z_s(fflag); + ssum_tmp = ssum(fflag); + res_tmp .*= k*(-z_s_tmp)/((k + 1)^2); + ssum_tmp += res_tmp; + k++; + res(fflag) = res_tmp; + ssum(fflag) = ssum_tmp; + z_s(fflag) = z_s_tmp; + fflag = abs (res) > tol*abs (ssum); + endwhile + e1_s -= ssum; + ## Continued fraction, + ## Abramowitz, Stegun 'Handbook of Mathematical function', + ## formula 5.1.22 p 229 + ## modified Lentz's algorithm, + ## from 'Numerical recipes in Fortrann 77' p.165 + f_new = 2^(-100)*ones (size (z_cf), class (z_cf)); + C_new = f_new; + D_new = zeros (size (z_cf), class (z_cf)); + k = 0; + fflag = true (size (z_cf)); + ## Variables initialization + Delta = C_old = D_old = f_old = ones (size (z_cf), class (z_cf)); + while (k < 1e3 && any(fflag)) + z_cf_tmp = z_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 = z_cf_tmp*(mod(k,2) == 0) + (mod(k,2) == 1); + a = exp(-z_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; + z_cf(fflag) = z_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++; + 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(-z_a)./z_a; + oldres = ssum = res = ones (size (z_a), class (z_a)); + k = 0; + fflag = true (size (z_a)); + while (k < 1e3 && any(fflag)) + res_tmp = res(fflag); + oldres_tmp = res_tmp; + z_a_tmp = z_a(fflag); + res_tmp ./= (-z_a_tmp/(k+1)); + ssum(fflag) += res_tmp; + k++; + res(fflag) = res_tmp; + oldres(fflag) = oldres_tmp; + fflag = abs (oldres) > abs (res); + endwhile + e1_a .*= ssum; + # Re-arranging + E1(flag_s) = e1_s; + E1(flag_cf) = e1_cf; + E1(flag_a) = e1_a; + E1 = reshape (E1, sz); + # if the input was in sparse format, also the output will be + if flag_sparse + E1 = sparse (E1); + endif +endfunction + +## Test +## 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 (y_exp,expint(X),-1e-12) + +## 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 input validation +%!error expint () +%!error expint (1,2) +%!error expint ('1') + +## Test on preservation or conversion of the class +%!test +%! 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 (class (expint (true)), 'double')