function INT=NLS_control_int(m,M,mphi,F,x,index_max,Bint,rphi,s,sphi)
%m=finite dimensional parameter
%F= matrix  with coefficients  F(k,l)
%x= numerical solution of f^m=0
%bbar= bbar_k k=0...2mphi    (bbar(1)=bbar_0=0)   center of Bint

%Then uses bbar(-n)=-bbar(n)

global mu sigma

beta=x(1);
bbar=Bint(1:2*mphi+1).mid;
mphi=(length(bbar)-1)/2;

i0=intval(0);
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);
ibeta=intval(beta);

isigma=intval(sigma);

imu=intval(mu);

format long



persistent conv2bbar alpha2zero alpha2 mEzero mE mEtilde Mat H1 C1 Constant

% The next if... end is calculated only the first time: the variables there
% do not depend on the solution x but only on the data bbar. hence once they
% are computed , they are the same for the next computation

if isempty(conv2bbar)==1

    
    % conv2bbar  convolution of bbar : sum_k1+k2=k bbar_k1 bbar_k2    for k=0....
    
 
    conv2bbar=quadratic_sumFFT_int(intval(bbar));
    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:m-1
           alpha2(k)=alpha2(k)+sum(intval(k)^isphi./(intval([1:1:k-1]).^isphi.*(intval(k-[1:1:k-1])).^isphi));
    end
    for k=m:2*M
    alpha2(k)=alpha2(k)+gamma_fun(k,isphi);
    end


%%%% \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:m-1
        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
    for q=m:2*M
        mE(q)=i4^isphi*irphi^i2*alpha2(q)/(intval(q)^isphi)+i2*irphi*i2^isphi*sum(abs(bbarsym).*Wsym(4*M+1+2*mphi+q:-1:4*M+1+q-2*mphi).^(-isphi));
    end
    mEtilde=[1:1:length(mE)].^(isphi).*mE;
    

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

    
    %%%%%% construct H1
    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)));
        H1(k)=H1(k)+i3*(i2*sum(mE((m+1:1:m+k)-k).*(m+1:1:m+k).^(-is))+i2*mE(m+1+k)/((is-i1)*im^(is-i1))+i2*mE(m+1)/((is-i1)*(im+k)^(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)));
        H1(k)=H1(k)+i3*(i2*sum(mE(k-[1:1:k-1]).*[1:1:k-1].^(-is))+i2*sum(mE([k+1:1:m+k]-k).*[k+1:1:m+k].^(-s))...
                        +i2*(mE(k+1)+mE(k+2)*i1/(is-i1)+mE(m+1)*i1/((is-i1)*(m+k)^(is-i1))));%
  
    end
    C1=reshape([H1';H1'],[],1);

end

% Computing C_Lambda
C_F=i2*(i2*norm(bbar)^i2+mEzero+mE(2*M));
C_Lambda=i1/((iM)^i2*(ipi^i2+(isigma*imu-abs(ibeta)-C_F)/(iM)^i2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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

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


MatL=Mat+diag(repmat([isigma*imu-ibeta,isigma*imu+ibeta],[1,m]));

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

ifm=MatL*iXk;

%%%%%
% Df(m)
iDfm=[((-1).^[1:1:2*m])'.*iXk,MatL(:,1:2*(index_max-1)),MatL(:,2*index_max:end)];

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


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



%%%%%%%%%%%%
%% vector Wj 
iWm=[i1;reshape(intval([1:1:index_max-1;1:1:index_max-1]).^(-is),[],1);intval(index_max)^(-is);reshape(intval([index_max+1:1:m;index_max+1:1:m]).^(-is),[],1)];


Z0=abs(intval(eye(2*m))-intval(Am)*iDfm)*iWm;

%%%%%%%%%%%%%%

C2=[i2*reshape(repmat([1:1:index_max-1].^(-is),[2,1]),[],1);index_max^(-is)*[intval(10^-15);2];i2*reshape(repmat([index_max+1:1:M].^(-is),[2,1]),[],1)];

Z1=abs(intval(Am))*C1(1:2*m,1);
Z2=abs(intval(Am))*C2(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-ibeta-i2*Fkk)*(ipi^i2*intval(k)^i2+imu+ibeta-i2*Fkk)-Fkk*Fkk)*[ipi^i2*intval(k)^i2+imu+ibeta-i2*Fkk,-Fkk;-Fkk,ipi^i2*intval(k)^i2+imu-ibeta-i2*Fkk];
Y=[Y;abs(inv_Lambdak*Matk*iXk)];
Z1=[Z1;abs(inv_Lambdak)*[H1(k);H1(k)]];
Z2=[Z2;abs(inv_Lambdak)*C2(2*(k-1)+1:2*(k-1)+2)];
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^isphi-is)*max(sup(mH));

P=[1:1:4*mphi].';
Z1M=i6*(sum(conv2bbar(1+P).*(i1./((i1-intval(P)/iM).^is+i1)))+mEtilde(M+1)+mEtilde(M+2)/(is-i1)+max(sup(mEtilde(1:m)))*(gamma_fun(M,is)+i1+i1/(is-i1)));
Z2M=i2;
%ZM=C_Lambda/(iM^(is+i2))*[CM1,CM2];
ZM=C_Lambda*[Z1M,Z2M];


% 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,Z2];
rd=real([rd;[YM,ZM(1)-1,ZM(2)]]);


rd=fliplr(rd); 
INT=evaluate_intersection_int(real(rd.sup));
end


function gam=gamma_fun(M,isphi)
i1=intval(1);
i2=intval(2);
i3=intval(3);
i4=intval(4);
i6=intval(6);
iM=intval(M);
ipi=midrad(3.14159265358979323,1e-16);
gam=i2*(iM/(iM-i1))^isphi+(i4*log(iM-i2)/iM+(ipi^i2-i6)/i3)*(i2/iM+i1/i2)^(isphi-i2);
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

function neg_inter=evaluate_intersection_int(rd)
i2=intval(2);
i4=intval(4);
rd=intval(rd);
neg_inter=NaN;
Rm=(-rd(:,2)-sqrt(rd(:,2).^i2-i4*rd(:,1).*rd(:,3)))./(i2*rd(:,1));  % root of less magnitude
Rp=-Rm-rd(:,2)./rd(:,1);                                            % root of larger magnitude
ISCOMP=find(isreal(Rp)+isreal(Rm)<2, 1);
if isempty(ISCOMP)==0
    disp('error: a polynomial is everywhere positive');
    return
end
intersection=[max(Rm.sup),min(Rp.inf)];
if intersection(1)<intersection(2)
    neg_inter=intersection;
    disp('Enclosure successful');
end

end
