function rad=checksimpleig(Beta,V,m,M,Xbar,mphi,F,Bint,rphi,s,sphi)

%m=finite dimensional parameter
%F= matrice dei coefficienti F(k,l)
%x= numerical solution of f^m=0
%bbar= bbar_k k=0...2mphi    (bbar(1)=bbar_0=0)
%Beta eigenvalue
global mu


radV=V(1).rad;
bbar=Bint(1:2*mphi+1).mid;
mphi=(length(bbar)-1)/2;


i1=intval(1);
i2=intval(2);
i3=intval(3);
i4=intval(4);
i6=intval(6);
ipi=midrad(3.14159265358979323,1e-16);
is=intval(s);
isphi=intval(sphi);
iM=intval(M);
irphi=intval(rphi);
im=intval(m);

imu=intval(mu);

format long

conv2bbar=real(quadratic_sumFFT_int([intval(bbar);zeros(4*mphi-size(bbar,1),1)]));
conv2bbar(2:2:end)=zeros(fix(size(conv2bbar,1)/2),1);
conv2bbar=[conv2bbar;zeros(4*M,1)];

%%%%%%% w_k k=0.....M

iW=[i1;intval([1:1:4*M].')];

%%%%%%%%%%%% compute alpha2 


alpha2zero=i4+i1/(i2^(i2*isphi-i1)*(i2*isphi-i1));
alpha2=i2*(i2+i2^(-isphi)+i3^(-isphi)+i1/(i3^(isphi-i1)*(isphi-i1)))*ones(2*M,1);
for k=1:2*M
    alpha2(k)=alpha2(k)+sum(intval(k)^isphi./(intval([1:1:k-1]).^isphi.*(intval(k-[1:1:k-1])).^isphi));
end
gammaM=i2*(iM/(iM-i1))^isphi+(i4*log(iM-i2)/iM+(ipi^i2-i6)/i3)*(i2/iM+i1/i2)^(isphi-i2);

%%%% \mathcal E(q)
bbarsym=[-flip(bbar(2:end),1);bbar];
Wsym=[flip(iW(2:end),1);iW];      %Wsym(M+1)=w_0
mEzero=i4^isphi*irphi^i2*alpha2zero+i2*irphi*i2^isphi*sum(abs(bbarsym).*Wsym(4*M+1+2*mphi:-1:4*M+1-2*mphi).^(-isphi));
mE=intval(zeros(1,2*M));
for q=1:2*M
    if even(q)==1
   mE(q)=i4^isphi*irphi^i2*alpha2(q)/(intval(q)^isphi)+i2*irphi*i2^isphi*sum(abs(bbarsym).*Wsym(4*M+1+q+2*mphi:-1:4*M+1+q-2*mphi).^(-isphi));
    end
end
mEtilde=[1:1:length(mE)].^(isphi).*mE;
% Computing C_Lambda
C_F=i2*(i2*norm(bbar)^i2+mEzero+mE(2*M));
C_Lambda=i1/(ipi^i2-(abs(Beta)+C_F)/(iM)^i2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x01=intval(Xbar(1));
Xk=Xbar(2:end);    % vector of ck,dk
iXk=intval(Xk);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Constant matrix, used to define f(x) and the Df(x)

Mat=[];
Constant=repmat([-i2,i1;i1,-i2],[1,m]);
for k=1:m
    lineFk=reshape([F(k,1:m);F(k,1:m)],1,[]);
    Matk=Constant.*repmat(lineFk,[2,1]);
    Mat=[Mat;Matk];
end
Mat=Mat+diag(reshape(repmat((ipi*intval([1:1:m])).^i2,[2,1]),[],1));
MatL=Mat+diag(repmat([imu-Beta,imu+Beta],[1,m]));


%%%%%%%%%
% g(xbar)  k=1....m

igm=x01*V+MatL*iXk;
igm=[sum(V(1:2*m).*iXk(1:2*m));igm];
%%%%%
% Dg(m)
iDgm=[intval('0'),V';V,MatL];

%%%%%%%%%%%%%%%%%
%%%% Am numerical approximation of Dfm
Dgm=iDgm.mid;
Am=inv(Dgm);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Buonds Y,Z
Y=abs(intval(Am)*igm);


%%%%%%%%%%%%
%% vector Wj 

iWm=[i1;reshape(intval([1:1:m;1:1:m]).^(-is),[],1)];

Z0=abs(intval(eye(2*m+1))-intval(Am)*iDgm)*iWm;

%%%%%%%%%%%%%%
H1=intval(zeros(M,1));
for k=1:m
    J=[m+1:1:4*mphi+k].';
   H1(k)=i3*( sum(i2*(abs((k+J<4*mphi+1).*(conv2bbar(1+k+J)))+abs(conv2bbar(1+J-k))).*intval(J).^(-is))+i2*(mE(m+1+k)+mE(m+1-k))/((is-i1)*im^(is-i1)));
end
for k=m+1:M
    J=[max(1,k-4*mphi):1:k-1,k+1:1:k+4*mphi].';
    H1(k)=i3*(sum(i2*(abs((k+J<4*mphi+1).*(conv2bbar(1+k+J)))+abs(conv2bbar(1+abs(J-k)))).*intval(J).^(-is))+i2*sum(mE(k-[1:1:k-1]).*[1:1:k-1].^(-s))+...
        i2*(mE(1)*i1/(intval(k)+i1)^is+mE(2)*i1/((is-i1)*(intval(k)+i1)^(is-i1))+mE(k+1)+mE(k+2)*i1/(is-i1)));
    H1(k)=H1(k)+radV*k^(-is);
    
end
C1=reshape([H1';H1'],[],1);

Z1=abs(intval(Am))*[intval(radV)*i1/((i2*is-i1)*intval(m)^(i2*is-i1));C1(1:2*m,1)];

for k=m+1:M

Fkl=i2*(infsup(conv2bbar(1+k+[1:1:m]).inf-(mE(k+[1:1:m]).sup)',conv2bbar(1+k+[1:1:m]).sup+(mE(k+[1:1:m]).sup)')-infsup(conv2bbar(1+k-[1:1:m]).inf-(mE(k-[1:1:m]).sup)',conv2bbar(1+k-[1:1:m]).sup+(mE(k-[1:1:m]).sup)'));
lineFk=reshape([Fkl.';Fkl.'],1,[]);
Matk=Constant.*repmat(lineFk,[2,1]);

Fkk=i2*(infsup(conv2bbar(1+k+k).inf-(mE(k+k).sup)',conv2bbar(1+k+k).sup+(mE(k+k).sup)')-infsup(conv2bbar(1).inf-(mEzero.sup)',conv2bbar(1).sup+(mEzero.sup)'));


inv_Lambdak=i1/((ipi^i2*intval(k)^i2+imu-Beta-i2*Fkk)*(ipi^i2*intval(k)^i2+imu+Beta-i2*Fkk)-Fkk*Fkk)*[ipi^i2*intval(k)^i2+imu+Beta-i2*Fkk,-Fkk;-Fkk,ipi^i2*intval(k)^i2+imu-Beta-i2*Fkk];
Y=[Y;abs(inv_Lambdak*Matk*iXk)];
Z1=[Z1;abs(inv_Lambdak)*[H1(k);H1(k)]];

end

%Tail bound

iXk=reshape(iXk,2,[]);


L=[1:1:m];
mH=i2*sum(repmat(mEtilde(M+L)+(i1-intval(L)/iM).^(-is).*mEtilde(M-L),[2,1]).*[abs(i2*iXk(1,1:m)-iXk(2,1:m));abs(iXk(1,1:m)-i2*iXk(2,1:m))],2);

YM=C_Lambda/(iM^(i2+isphi))*max(sup(mH));

P=[1:1:4*mphi].';
CM1=i6*(sum(conv2bbar(1+P).*(i1./((i1-intval(P)/iM).^is+i1)))+mEtilde(M+1)+mEtilde(M+2)/(is-i1)+mEtilde(1)*(gammaM+i1)+mEtilde(2)/(is-i1));
CM1=CM1+radV;
ZM=C_Lambda/(iM^(is+i2))*CM1;


% Tail bound



% Polynomials
iWm=[iWm;reshape(intval([m+1:1:M;m+1:1:M]).^(-is),[],1)];
Z_r=[Z0;zeros(2*(M-m),1)]+Z1;
rd=[Y,Z_r-iWm];
rd=[rd;[YM,ZM(1)-(M)^(-is),]];
rd=rd.sup;

if max(rd(:,2))>0
    disp('error: Simplicity of eigenvalue  cannot  be proved');
    return
else
    disp(['Eigenvalue ' num2str(Beta.mid) ' is simple']);
    rad=max(-rd(:,1)./rd(:,2));
end

end


function [s]=quadratic_sumFFT_int(a)

i2=intval('2');
i4=intval('4');
M=size(a,1);
% We make sure that the inputs are powers of 2 %
M1=(2^ceil(log(4*M-1)/log(2))-4*M)/2;   % Hence 4*M+2*M1 is a power of 2 %
ta1=[intval(zeros(M+M1+1,1));-a(M:-1:2,1);a;intval(zeros(M+M1,1))];
tu1=verifyfft(ifftshift(ta1),-1);
F=fftshift(verifyfft(tu1.*tu1,1));
s=(i4*intval(M)+i2*intval(M1))*F(M1+2*M+1:M1+3*M,1);
end
