From: RIKA fazi on
% sum[(-1)^i* factorial(M-1,i)* sum ( prod(k-m) )]
% i=0->M-1 j=0toM-1,j<>i k=0toM-1,k<>j
% sum [of sum of (prod)]
% for fixed M; m varies from 0 to M-1;

% example M=5
% i=0; k=0,1,2,3,4; m=1

% (-1)^i*bin(4,i)*sum( j=0 (0-m)(1-m)(2-m)(3-m)(4-m)/(0-m) is nul because i=j
% ( j=1 (0-m)(1-m)(2-m)(3-m)(4-m)/(1-m) is not nul i and j differe
% ( j=2 (0-m)(1-m)(2-m)(3-m)(4-m)/(2-m) is not nul i and j differe
% ( j=3 (0-m)(1-m)(2-m)(3-m)(4-m)/(3-m) is not nul i and j differe
% ( j=4 (0-m)(1-m)(2-m)(3-m)(4-m)/(4-m) is not nul i and j differe




% prod(k-m) for k varying from 0 to M-1 with k not equal to m
%we can use the function nonzeros for force k not equal to m
c21=1;
for k=0:M-1
c21=c21*prod(nonzeros(k-m));
end
c1c=c21;
C1=0; %%% intialize the external sum

for i=0:M-1

c212=0;c210=0;
for j=0:M-1
djm=eq(j,m);
dij=eq(i,j);
if dij>=1 %% if i=j then the product is nul
cc1=0; %% that is cc1=0
elseif (dij<=0) % if i and j are different
if djm>=1 % then if j and m are equal we take the same product
c210=c1c;
elseif djm<=0 % else that is j and m different we divide the product by j-m
c211=c1c./(j-m) ;
c212=c212+sum(c211); % sum all product
cc1=(-1).^i*bin(M-1,i)*(c210+c212); % multiply the sum by (-1)^i*bin(M-1,i)
C11=cc1;
end
end
end
C1=C1+C11;

end
C1=C1