function w = bilin(u,v)
% BILIN computes biliner form w = a(u,v) as convolution of fft's
% w_hat(jj) = \sum_{kk} lambda(-jj,kk,jj-kk) u_hat(kk)*v_hat(jj-kk)
n = min(length(u),length(v));
u_hat = fft(u,n);
v_hat = fft(v,n);
u_shift = fftshift(u_hat);
v_shift = fftshift(v_hat);
u_pad = [zeros(1,n/2) u_shift(2:n) zeros(1,n/2)];
v_pad = [zeros(1,n/2) v_shift(2:n) zeros(1,n/2)];
w_hat = zeros(1,n);
lcf = zeros(1,(2*n-1));
for j=2:n
    jj = j-n/2-1;
    for k=1:(2*n-1)
        kk = k-n;
        lcf(k) = lambda(-jj,kk,jj-kk);
    end
    v_pad_shift = fliplr(circshift(v_pad,[0 -jj]));
    w_hat(j) = sum(lcf.*u_pad.*v_pad_shift)/n;
end
w_hat = ifftshift(w_hat);
w = ifft(w_hat);
w = real(w);
end

function lcf = lambda(k1,k2,k3)
% LAMBDA computes three-wave interaction coefficient
mod = abs(k1) + abs(k2) + abs(k3);
if (mod==0)
    lcf=0;
else
    % coefficient for 0.5*[h^2]_{xx} - H[hu_{xx}] where h=H[u]
    lcf = 2*abs(k1*k2*k3)/mod;
end
end