function[e] = numEigContour(poly,L,C,R,N)
% This function approximates the number of eigenvalues of $P(\lambda)$ inside $\Gamma$.
% We used the paper: Stochastic estimation method of eigenvalue density for nonlinear problem on the comple plane
% by Y. Maeda, Y. Futamura, T, Sakurai.

% Input:  coeffs  -  cell array containing coefficients of matrix polynomial: P(x) = A_l x^l + ...+ A_2 x^2 + A_1 x + A_0
%              L  -  number of sample vectors
%              C  -  center of circle
%              R  -  radius of circle
%              N  -  number of nodes of the trapezoid rule to approximate the contour integral
% Output:      e  -  approximation to the number of eigenvalues inside contour Gamma

len = length(poly)-1; n = length(poly{1}); v = zeros(n,L);

for i=1:L
    %random vector v_i
    index = randi(0:1,[n 1]);
    %(rand(n,1) <= 0.5)
    for k=1:n
        if(index(k)==1)
            v(k,i) = 1;
        else
            v(k,i) = -1;
        end
    end
end

S = 0;
for j=0:N-1
    expon = exp(1i*j*2*pi/N); phi = (C+R*expon);
    % Computing F(phi) and DF(phi)
    F = zeros(n,n); DF = zeros(n,n);
    for m = 0:len
        F = F + phi^m*poly{m+1};        
        if(m>0)
            DF = DF + m*phi^(m-1)*poly{m+1};
        end
    end
    S1 = 0;
    for k=1:L
        S1 = S1 + v(:,k)'*(F\(DF*v(:,k)));
    end
    
    S = S + (R/N)*expon*S1;
end
e = round(real((1/L)*S));
end