function r=radii_polynomials_phi_odd(x,s)

global mu sigma
is=intval(s);
i0=intval(0);
i1=intval(1);
i2=intval(2);
i3=intval(3);
i4=intval(4);
i6=intval(6);

ipi=midrad(3.14159265358979323,1e-16);
imu=intval(mu);
isigma=intval(sigma);

ix=intval(x);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m=length(x)+1;

normx=max(sup(abs(x).*intval([1:1:m-1]).'.^is));
p=3;

M=p*(m-1)+1;
iM=intval(M);

Dfm=jacobian_map_phi_odd_int(ix);
muk=ipi^i2*(i2*[1:M-1]'-i1).^i2+isigma*imu;

Am=inv(Dfm.mid);

f=map_phi_odd_int(ix,M-1);
fm=f(1:m-1);
Ws=intval([1:1:M]).'.^(-is);

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% compute alpha2 alpha3 alphatilde

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

alpha3zero=alpha2zero+i2*sum(alpha2([1:1:M-1]).*(intval([1:1:M-1]).'.^(-i2*is)))+i2*alpha2(M)/((iM-i1)^(i2*is-i1)*(i2*is-i1));
alpha3=alpha2zero*ones(M,1);
for k=1:M-1
    k1=[1:M-k-1].';
    k2=[1:1:k-1].';
    k3=[1:1:M-1].';
    ik1=intval(k1);
    ik2=intval(k2);
    ik3=intval(k3);
    ik=intval(k);
   alpha3(k)=alpha3(k)+sum(ik^is*alpha2(k1+k)./(ik1.^is.*(ik1+ik).^is))+alpha2(M)*ik^is*(i1/(iM^is*(iM-ik)^is)+i1/(iM^is*(iM-ik)^is*(is-i1)))+...
        alpha2(k)+sum(ik^is*alpha2(k2)./(ik2.^is.*(ik-ik2).^is))+sum(ik^is*alpha2(k3)./(ik3.^is.*(ik3+ik).^is))+alpha2(M)/((iM-i1)^(is-i1)*(is-i1));    
end
k1=[1:1:M-1].';
ik1=intval(k1);
alpha3(M)=alpha3(M)+alpha2(M)*(i2+i2^(-is)+i3^(-is)+i1/(i3^(is-i1)*(is-i1))+i1/((iM-i1)^(is-i1)*(is-i1))+gammaM)+sum(alpha2(k1)./(ik1.^is).*(i1+iM^is*(iM-ik1).^(-is)));


alpha3tilde=max(sup([alpha3zero;alpha3]));
mutilde=ipi^i2*(i2*iM-i1)^i2+isigma*imu;

PtildeM0=i6*intval(normx)^i2*alpha3tilde/mutilde-i1;
if sup(PtildeM0)>0
    disp('the tail fails');
    r=Nan
    return
end

eps3=i2*alpha2(M)/((is-i1)*(iM-i1)^(is-i1))*(intval([1:1:M-1]+M)).^(-is);
for k=1:M-1
   eps3(k)=eps3(k)+sum(alpha2([M:1:M+k-1]-k)./(intval([M:1:M+k-1]).'.^is.*(intval([M:1:M+k-1]-k)).'.^is)); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%  Y   %%%%%%%%%%%
Y=[abs(intval(Am)*fm);abs(muk(m:M-1).^(-i1).*f(m:M-1))];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Z %%%%%%%%%%%%%%%%%

Z0=abs(intval(eye(m-1))-intval(Am)*Dfm)*Ws(1:m-1);

cub1=cubic_sumFFT_int([abs(ix);i0],[abs(ix);i0],[intval(zeros(m-1,1));Ws(m:M-1)]);
Z1(1:m-1,1)=cub1(2:m)+normx^i2*eps3(1:m-1)';

cub2=cubic_sumFFT_int([abs(ix);i0],[abs(ix);i0],Ws(1:M-1));
Z1(m:M-1,1)=cub2(m+1:M)+normx^i2*eps3(m:M-1)';

cub3=cubic_sumFFT_int(abs(ix),Ws(1:M-1),Ws(1:M-1));
Z2=cub3(2:M)+i2*normx*eps3(1:M-1)';

cub4=cubic_sumFFT_int(Ws(1:M-1),Ws(1:M-1),Ws(1:M-1));
Z3=cub4(2:M)+i3*eps3(1:M-1)';

%%%%%%%%%%%%%%%%  Tail


Z0=[Z0;zeros(M-m,1)];
rd=[Y,Z0+i6*([abs(Am)*Z1(1:m-1,1);Z1(m:M-1)./abs(muk(m:M-1))])-Ws(1:M-1),i6*i2*[abs(Am)*Z2(1:m-1,1);Z2(m:M-1,1)./abs(muk(m:M-1))],i6*[abs(Am)*Z3(1:m-1);Z3(m:M-1)./abs(muk(m:M-1))]];
   
rd=real([rd;[i6*normx^i2*alpha3tilde/mutilde-i1,i6*i2*normx*alpha3tilde/mutilde,i6*alpha3tilde/mutilde,0]]);
r=compute_intersection(rd.sup);



    function r=compute_intersection(rd)
        neg_int=[];
        for k=1:M-1
       Root=roots(rd(k,end:-1:1));
       Realroots=[-Inf;sort(Root(imag(Root)==0),'ascend')];
       neg_int=[neg_int;[Realroots(end-1),Realroots(end)]];
        end
       Root=roots(rd(end,end:-1:1));
       if isreal(Root)==1
           neg_int=[neg_int;sort(Root,'ascend')'];
       else
           disp('error:  positive tail polynomial');
           return
       end
       intersection=[max(neg_int(:,1)),min(neg_int(:,2))];
       if intersection(1)<intersection(2)
           disp('Enclosure successful')
                r=max(intersection(1),0);
       else
           r=NaN
       end      
       done=0;
       attemps=0;
       while done==0
       ird=intval(rd);
       ir=intval(r);
       Evaluation=ird(:,1)+ird(:,2)*ir+ird(:,3)*ir^i2+rd(:,4)*ir^i3;
       if sum(sup(Evaluation)>0,1)==0
           done=1;
       end
       if attemps==20
           return
       end
       r=r+10^-16;
       attemps=attemps+1;
       end
    end     
    
end

function [s]=cubic_sumFFT_int(a1,a2,a3)
% the convolution term using FFT

i2=intval('2');
i4=intval('4');

Ma1=size(a1,1);
Ma2=size(a2,1);
Ma3=size(a3,1);

M=max([Ma1,Ma2,Ma3]);
a1=[a1;intval(zeros(M-Ma1,1))];
a2=[a2;intval(zeros(M-Ma2,1))];
a3=[a3;intval(zeros(M-Ma3,1))];

M=size(a1,1)+1;
M1=(2^ceil(log(4*M-1)/log(2))-4*M)/2; 

b1=flipud(a1(1:end));
b2=flipud(a2(1:end));
b3=flipud(a3(1:end));
ta1=[intval(zeros(M+M1+1,1));intval(0);b1;a1;intval(zeros(M+M1,1))];
ta2=[intval(zeros(M+M1+1,1));intval(0);b2;a2;intval(zeros(M+M1,1))];
ta3=[intval(zeros(M+M1+1,1));intval(0);b3;a3;intval(zeros(M+M1,1))];

tu1=verifyfft(ifftshift(ta1),-1);
tu2=verifyfft(ifftshift(ta2),-1);
tu3=verifyfft(ifftshift(ta3),-1);

F=fftshift(verifyfft(tu1.*tu2.*tu3,1));
s=(i4*intval(M)+i2*intval(M1))^i2*F(M1+2*M+1:M1+3*M,1);

end