function y = expint(x)
%MATLAB Code Generation Library Function

%   Copyright 2002-2014 The MathWorks, Inc.
%#codegen

eml_invariant(nargin == 1,'MATLAB:minrhs');
eml_invariant(isfloat(x),'Coder:toolbox:expint_2');
ytype = zeros('like',complex(real(x)));
y = coder.internal.applyScalarFunction(mfilename,ytype,@scalar_expint,x);

%--------------------------------------------------------------------------

function y = scalar_expint(x)
ONE = complex(ones('like',x));
ZERO = complex(zeros('like',x));
INF = complex(coder.internal.inf('like',x));
% figure out which algorithm to use by evaluating interpolating polynomial
% at real(z)
egamma = cast(0.57721566490153286061,'like',real(x));
p = cast( ...
    [-3.602693626336023e-09 -4.819538452140960e-07 -2.569498322115933e-05 ...
    -6.973790859534190e-04 -1.019573529845792e-02 -7.811863559248197e-02 ...
    -3.012432892762715e-01 -7.773807325735529e-01  8.267661952366478e+00], ...
    'like',real(x));
% series expansion
pk = polyval(p,real(x));
absimagxk = abs(imag(x));
if x == 0
    y = INF;
elseif absimagxk <= pk
    y = -egamma - log(complex(x));
    j = 1;
    pterm = x;
    term = x;
    while abs(term) > eps(y)
        y = y + term;
        j = j + 1;
        pterm = coder.internal.div(-x*pterm,j);
        term = coder.internal.div(pterm,j);
    end
elseif absimagxk > pk
    n = 1;
    am2 = ZERO;
    bm2 = ONE;
    am1 = ONE;
    bm1 = complex(x);
    f = coder.internal.div(am1,bm1);
    oldf = INF;
    j = 2;
    while abs(f - oldf) > 100*eps(f)
        % calculate the coefficients of the recursion formulas for j
        % even
        alpha1 = n - 1 + eml_rdivide(j,2); % note: beta1 = 1
        % calculate A(j), B(j), and f(j)
        a = am1 + alpha1*am2;
        b = bm1 + alpha1*bm2;
        % save new normalized variables for next pass through the loop
        % note: normalization to avoid overflow or underflow
        am2 = coder.internal.div(am1,b);
        bm2 = coder.internal.div(bm1,b);
        am1 = coder.internal.div(a,b);
        bm1 = ONE;
        f = am1;
        j = j + 1;
        % calculate the coefficients for j odd
        alpha1 = eml_rdivide(j-1,2);
        beta1 = x;
        a = beta1*am1 + alpha1*am2;
        b = beta1*bm1 + alpha1*bm2;
        am2 = coder.internal.div(am1,b);
        bm2 = coder.internal.div(bm1,b);
        am1 = coder.internal.div(a,b);
        bm1 = ONE;
        oldf = f;
        f = am1;
        j = j + 1;
    end
    % y = exp(-x)*f - 1i*pi*((real(x)<0)&&(imag(x)==0));
    y = exp(-x)*f;
    if real(x) < 0 && imag(x) == 0
        y = y - 1i*pi;
    end
else
    y = ZERO;
end

%--------------------------------------------------------------------------
