% function [z, e2] = fastpoly2rank1pt(m, s_max, omega, gamma, beta)
%
% Polynomial lattice rule over Z_2/f with deg(f) = m.
%
% Note the beta is only in the product up front for the squared worst-case
% error formula, not in the matrix-vector multiplication.
% -> changed as for the scalar case
%
% inputs
%   m           degree of the irreducible polynomial
%   s_max       number of dimensions, scalar
%   omega       a function handle for the varying kernel part of the
%               shift-invariant kernel function (assumed to be symmetric)
%   gamma       gamma parameters for AC weighting per dimension, vector [s_max x 1]
%   beta        beta parameters for DC weighting per dimension, vector [s_max x 1]
% outputs
%   z           generating vector of the lattice rule, vector [s_max x 1]
%   e2          optimal square error per dimension (= one for each iteration),
%               vector [s_max x 1]
%
% e.g.
%   m = 10; s_max = 100;
%   usob = @(x) 1/6 - pow2(floor(log2(x))-1);
%   gamma = ones(s_max, 1) / s_max; beta = ones(s_max, 1);
%   [z, e2] = fastrank1ppt(m, s_max, usob, gamma, beta);
%   % plot the square error in function of the dimension
%   semilogy(e2);
%
% (C) 2007, <dirk.nuyens@cs.kuleuven.ac.be>
function [z, e2] = fastpoly2rank1pt(m, s_max, omega, gamma, beta)

z = zeros(s_max, 1);
e2 = zeros(s_max, 1);
n = pow2(m);

g = gf(2, m); % generator $g = 2 = (10)_2 = x$ since $f$ is primitive 
perm = gf(zeros(n-1, 1), m);
perm(1) = 1; for j=1:n-2, perm(j+1) = perm(j)*g; end;
psi = omega(vm(perm.x, g.prim_poly, m)');
psi0 = omega(0);
fft_psi = fft(psi);

E2 = zeros(n, 1);
cumbeta = cumprod(beta);
q = ones(n-1, 1);
q0 = 1;

for s = 1:s_max
    E2 = real(ifft(fft_psi .* fft(q)));
    [min_E2, w] = min(E2);
    if s == 1, w = 1; noise = abs(E2(1) - min_E2), end;
    z(s) = getfield(perm(w), 'x');
    e2(s) = -cumbeta(s) + ( beta(s) * (q0 + sum(q)) + ...
               gamma(s) * (psi0*q0 + min_E2) ) / n;
    q = (beta(s) + gamma(s) * psi([w:-1:1 n-1:-1:w+1])) .* q;
    q0 = (beta(s) + gamma(s) * psi0) * q0;
    fprintf('s=%4d, z=%6d, w=%6d, e2=%.4e, e=%.4e\n', ...
             s, z(s), w, e2(s), sqrt(e2(s)));
end;