N = 20;
gamma = 3;
a = .1 + 3 * rand(N, 1);

mina = min(a);

x = linspace(0, gamma / mina, 1000);

f = @(beta) sum(min(gamma, a * beta(:)')) - N;

corners = gamma ./ a;
corners = [corners; 0; gamma / mina];
values = f(corners);
[values, reorder] = sort(values);
corners = corners(reorder);

low = 1;
upp = length(corners);

iters = 0;
while upp - low > 1
    iters = iters + 1;
    mid = floor((upp + low) / 2 + 0.1);
    % low < mid < upp
    if values(mid) > N
        upp = mid;
    else
        low = mid;
    end
end

fprintf('Binary search in %d iters, N at [%f, %f]n', ...
    iters, corners(low), corners(upp));

root = corners(low) + ...
    (corners(upp) - corners(low)) / (values(upp) - values(low)) * ...
    (N - values(low));

fprintf('The root is %fn', root);

plot(x, f(x), '-', corners, values, 'r.', ...
    x, N, 'k-', [root], [N], 'g*');