function [M, t, h, evals] = maxmod(p, epsilon, verbose,nint0)
%MAXMOD
%
%   The maximum modulus of a complex polynomial p on the unit
%   disc is estimated using a method of repeated subdivision
%   of [0,2pi] and rejection of non-maximising subintervals.
%   The rejection criterion is based on a lemma of S.B. Steckin.
%
%   [M, t, h, evals] = MAXMOD(p, epsilon, verbose, nint0)
%
%   p        complex vector of polynomial coefficients
%   epsilon  relative precision required [0.5e-8]
%   verbose  set non-zero for run-time information
%   nint0    initial number of intervals [5*N]
%
%   M        estimate of the maximum modulus of p
%   t        vector of midpoints of non-rejected subintervals
%   h        half-width of non-rejected intervals
%   evals    number of evaluations of the polynomial
%
%   The algorithm is described in
%
%   J.J. Green. Calculating the maximimum modulus of a polynomial
%   with Steckin's lemma. SIAM J. Numer. Anal., 83:211-222, 2000.
%
%   J.J.Green 1999, 2008, 2012, 2015

%   This code is contributed to the public domain.
%
%   The author requests that derived code is named in such a way as
%   to avoid confusion with the original version.
%
%   Changelog:
%   14 Mar 2013 - replaced octave rows() by matlab size(,1)
%   23 Mar 2012 - fixed incorrect calulation of beta
%   14 May 2012 - replaced poly evaluation by built-in polyval
%                 and some cosmetic tidying of the code

% get options

    if nargin < 1
        error('too few args')
    end

    p = p(:);
    N = size(p, 1);

    if nargin < 2
        epsilon = 0.5e-8;
    end

    if nargin < 3
        verbose = 0;
    end

    if nargin < 4
        nint0 = 5*N;
    end

    if nargin > 4
        error('too many args');
    end

    % sanity bounds on input

    if nint0 < 2 * N
        error('nint0 must be at least twice the order of p');
    end

    beta  = conv(p(N:-1:1), conj(p));
    beta0 = real(beta(N));
    normq = sum(abs(beta));

    evals = 0;

    if normq - beta0 < 4 * epsilon * beta0

        M2min = beta0;
        M2max = normq;

        M = sqrt((M2min + M2max) / 2);
        t = [];
        h = 2 * pi;
        evals = 0;

        return;
    end

    M2min = 0;
    M2max = inf;

    A = max([2 * beta0 - normq, 0]);

    % initial interval-set - V is a matrix
    %
    % [z0 q0
    %  z1 q1
    %    :
    %  zn qn]
    %
    % where qi = |p(zi)|^2

    h = pi/nint0;

    V = zeros(nint0, 2);

    V(:, 1) = (2 * h) * [0:(nint0 - 1)]';
    V(:, 2) = abs2(fft(p, nint0));

    eval0 = nint0;
    evals = nint0;

    if (verbose)
        fprintf('lower bound\tupper bound\tevals\trejected\n');
    end

    while 1

        % evaluate bounds

        M2min = max(V(:, 2));
        SC = cos(N * h);
        M2max = (M2min - A) ./ SC + A;

        if M2max - M2min < 4 * epsilon * M2min
            if verbose
                fprintf('%f\t%f\t%i\t-\n',...
                        sqrt(M2min),...
                        sqrt(M2max),...
                        eval0);
            end
            break;
        end

        % rejection

        rT = (M2min - A) * SC + A;
        ind = find(V(:, 2) > rT);
        nrej = size(V, 1) - numel(ind);
        V = V(ind, :);

        if verbose
            fprintf('%f\t%f\t%i\t%i\n',...
                    sqrt(M2min),...
                    sqrt(M2max),...
                    eval0, nrej);
        end

        % subdivision

        h = h / 3;
        offset = 2 * h;
        nt = size(V, 1);
        t = V(:, 1);

        % expand V and assign indices of the bottom end

        V = [V ; zeros(2 * nt, 2)];
        ind = (nt + 1):(3 * nt);

        % fill bottom end with new t valuee

        V(ind, 1) = [ t + offset ; t - offset ];

        % new z-values

        z  = exp(V(ind, 1).*1i);

        % assign new q-values

        V(ind, 2) = abs2(polyval(p, z));

        % update statistics

        eval0 = numel(ind);
        evals = evals + eval0;

    end

    % we can only get here from the break above

    M = sqrt((M2min + M2max) / 2);
    t = sort(V(:, 1));

end

% element-wise caculation of the magnitude
% squared of the complex vector z using only
% (vectorised) multiply and add

function a2z = abs2(z)
    rz = real(z);
    iz = imag(z);
    a2z = rz .* rz + iz .* iz;
end
