% function [z, e2, x, CBC_optimal] = fastrank1expt(p, m1, m2, ...
%   s_max, omega, gamma, beta, CBC_optimal, verbose)
function [z, e2, x, CBC_optimal] = fastrank1expt(p, m1, m2, ...
  s_max, omega, gamma, beta, CBC_optimal, verbose)

if (nargin < 9) || isempty(verbose), verbose = 1; end;
if ~isprime(p), error('p must be prime'); end;
if (m1 < 1) || (m1 > m2), error('must have 1 <= m1 <= m2'); end;

if (nargin < 8) || isempty(CBC_optimal) % if precalculations needed
  CBC_optimal = {}; CBC_optimal{m2} = ones(s_max, 1);
  if m1 ~= m2
    for m=m1:m2
      CBC_optimal{m} = CBC_optimal{m2}; % alias using Matlab's lazy copy
      [opt_z, CBC_optimal{m}] = fastrank1expt(p, m, m, s_max, ...
          omega, gamma, beta, CBC_optimal, verbose-1);
      end;
    end;
end;

n = p^m2;                                   % total points $n = p^{m_2}$
if p == 2, g = 5; else g = generator(p*p); end; % pseudo when $p=2$
phip0 = 1; phip = zeros(m2, 1);            % calculate $\eulerphi(p^m)$
phip(1) = p - 1; for m=2:m2, phip(m) = p*phip(m-1); end;
phip(phip >= 2) = phip(phip >= 2) / 2; % the size for $d=p^m$ is in \texttt{phip(m)}
perm = zeros(phip(m2), 1); perm(1) = 1;% $g^i \bmod{p^m}$ sequence
for i=2:phip(m2), perm(i) = mod(perm(i-1)*g, n); end;
if m2~=1, sidx = [1; cumsum(phip(1:m2-1))+1]; else sidx = [1]; end;
eidx = sidx + phip - 1;                     % start and ending indices

psi0 = omega(0);                            % build the $\vec{\psi}$ vector
psi = zeros(sum(phip), 1); fft_psi = zeros(sum(phip), 1);
for m=1:m2
  psi(sidx(m):eidx(m)) = omega(mod(perm(1:phip(m)), p^m) / p^m);
  fft_psi(sidx(m):eidx(m)) = fft(psi(sidx(m):eidx(m)));
end;

q0 = 1; q = ones(sum(phip), 1);
cumbeta = cumprod(beta);
z = zeros(s_max, 1); e2 = zeros(s_max, 1); x = zeros(s_max, 1);

if (verbose >= 1) && (m1 ~= m2)
  fprintf('Calculating rule for %d^%d to %d^%d\n', p, m1, p, m2);
elseif (verbose >= 1)
  fprintf('Calculating optimal CBC for %d^%d\n', p, m1);
end;
for s=1:s_max
  E2 = zeros(phip(m2), 1); X = zeros(phip(m2), m2-m1+1);
  C = beta(s) * q0; % accumulation of the constant parts goes in \texttt{C}
  E2(1) = gamma(s) * psi0 * q0;
  if p == 2, r = 1; else r = 2; end; % $p=2$ is special
  v = 1; % number of choices from the previous iteration for $m$ (next loop)
  for m=1:m2   % calculate error for each $m$
    C = C + beta(s) * r * sum(q(sidx(m):eidx(m)));
    E2(1:phip(m)) = repmat(E2(1:v), phip(m)/v, 1) ...
      + gamma(s) * r * real(ifft(fft_psi(sidx(m):eidx(m)) ...
        .* fft(q(sidx(m):eidx(m)))));
    if m >= m1 % store and calculate $X$ where needed
      X(:, m-m1+1) = repmat(sqrt((-cumbeta(s) + C/p^m ...
        + E2(1:phip(m))/p^m) / CBC_optimal{m}(s)), ...
          phip(m2)/phip(m), 1);
    end;
    r = 2; v = phip(m);
  end;
  [x(s), w] = min(max(X, [], 2)); % pick the good choice
  [dummy, mloc] = max(X(w,:)); mloc = mloc + m1 - 1;
  if s == 1, w = 1; end;
  z(s) = perm(w); z(s) = min(z(s), n-z(s)); % unpermute
  e2(s) = -cumbeta(s) + C/n + E2(w)/n;  % calculate real worst-case error
  q0 = (beta(s) + gamma(s) * psi0) * q0; % update $\vec{q}$...
  for m=1:m2                                 % in turn for each $m$
    ww = mod(w-1, phip(m)) + 1;
    q(sidx(m):eidx(m)) = (beta(s) + gamma(s) ...
      * psi([ww:-1:1 phip(m):-1:ww+1]+sidx(m)-1)) ...
        .* q(sidx(m):eidx(m));
  end;
  if verbose >= 2
    fprintf('s=%4d, z=%6d, e2=%.4e, x=%.4e, mloc=%4d\n', ...
             s, z(s), e2(s), x(s), mloc);
  end;
end;