function F=map_phi_odd_int(bodd,M)

global mu sigma

imu=intval(mu);
i1=intval(1);
i2=intval(2);
% i3=intval(3);
% i4=intval(4);
% i6=intval(6);
ipi=midrad(3.14159265358979323,1e-16);

%bodd=intval(bodd);

% aodd=i2*bodd;
m=length(bodd);

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

cf=cubic_sumFFT_int(bodd,bodd,bodd);


% for n=1:M
%     F(n,1)=(ipi^i2*(i2*intval(n)-i1)^i2+sigma*imu)*bodd(n)+(i2)*sigma*cf(n+1);
% end
F=(ipi^i2*(i2*intval([1:1:M].')-i1).^i2+sigma*imu).*bodd(1:M)+(i2)*sigma*cf(2:M+1);


% function Conv=convol(a,b,c,n)
% Sumpkl=a*PRa(:,1)+b*PRa(:,2)+c*PRa(:,3);
% PKL=PRa(Sumpkl==n,:);
% Conv=sum(bodd(PKL(:,1)).*bodd(PKL(:,2)).*bodd(PKL(:,3)));
% 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=-(flipdim(a1(1:end),1));
b2=-(flipdim(a2(1:end),1));
b3=-(flipdim(a3(1:end),1));
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+2:M1+3*M+1,1);

end
% function [s]=cubic_sumFFT_real(a1)
% 
% m=length(a1);
% 
% % 
%  b1=-(flipdim(a1(1:m),1));
%  ta1=[zeros(m,1);0;b1;a1;zeros(m,1)];
% 
% tu1=ifft(ifftshift(ta1));
% 
% Ft=fftshift(fft(tu1.*tu1.*tu1));
% 
% s=(4*m+1)^2*real(Ft(2*m+2:3*m+2));
% end