function J=jacobian_map_phi_even_int(beven)
global mu sigma
i2=intval(2);
i6=intval(6);
isigma=intval(sigma);

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

m=length(beven);
J=intval(zeros(m,m));

c2fft=quadratic_sumFFT_real_int([beven;intval(zeros(m,1))]);
c2fft=c2fft.';
 
 for n=1:m
     j=1:n;       
     J(n,j)=(ipi^i2*(i2*intval(n))^i2+isigma*imu)*(n==j)+real(isigma*i6*(c2fft(n-j+1)-c2fft(n+j+1)));       
     J(j,n)=J(n,j);     
 end
end

function [s]=quadratic_sumFFT_real_int(a1)

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

a=[intval(0);a1];
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







