function F=map_phi_even_int(beven,M)

global mu sigma
i0=intval(0);
i1=intval(1);
i2=intval(2);
%i3=intval(3);
i4=intval(4);


ipi=midrad(3.14159265358979323,1e-16);
imu=intval(mu);
%aodd=2*bodd;
m=length(beven);

%aodd=[aodd;zeros(M-m,1)];
beven=[beven;intval(zeros(M-m,1))];
              % all combinations with p,k,l in [-m+1,...,m]
 %PRa=perms_rep(m,3);

 %for n=1:M
%     F(n,1)=(ipi^i2*(i2*intval(n))^i2+sigma*imu)*beven(n)+i2*sigma*(-i3*convol(1,1,-1,n)+i3*convol(1,-1,-1,n)+convol(1,1,1,n));   
% end
cf=cubic_sumFFT_int([i0;beven],[i0;beven],[i0;beven]);
%cs=cubic_sum_loops(beven)
F=(ipi^i2*(i2*[1:1:M].').^i2+sigma*imu).*beven(1:M)+i2*sigma*cf(2:end);
    
function [s]=cubic_sumFFT_real(a1)

m=length(a1);

b1=-(flipdim(a1(2:m),1));
ta1=[intval(zeros(m,1));b1;a1;intval(zeros(m,1))];
tu1=ifft_int(ifftshift(ta1));


Ft=fftshift(fft_int(tu1.*tu1.*tu1));

s=(i4*im-i1)^i2*real(Ft(2*m:3*m-1));
end

function Conv=convol(a,b,c,n)
Sumpkl=a*PRa(:,1)+b*PRa(:,2)+c*PRa(:,3);
PKL=PRa(Sumpkl==n,:);
Conv=sum(beven(PKL(:,1)).*beven(PKL(:,2)).*beven(PKL(:,3)));
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))];
% 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));-a1(M:-1:2,1);a1;intval(zeros(M+M1,1))];
ta2=[intval(zeros(M+M1+1,1));-a2(M:-1:2,1);a2;intval(zeros(M+M1,1))];
ta3=[intval(zeros(M+M1+1,1));-a3(M:-1:2,1);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