% Fairly loose checking:
tol = 1e-12;

rel_err = @(v1, v2) abs(v1-v2)/max(realmin, abs(mean([v1,v2])));

fprintf('Testing small values\n');

test_vals = log([2, 1, 0.1, 0.01, 1e-3, 1e-4, 1e-5, 1e-6, 1e-9, ...
    eps, eps/2, eps/3, eps/10, eps/1e3, ...
    realmin*10, realmin*2, realmin, 0]);

fprintf('Should see some failures from naive method\n');
fprintf('Otherwise there would not be any point in my code\n');
for x = test_vals
    mine = log1pexp(x);
    ml = log1p(exp(x));
    naive = log(1 + exp(x));
    if exp(x) < 1
        taylor = log1ptaylor(exp(x));
        assert(rel_err(mine, taylor) < tol);
    end
    assert(rel_err(mine, ml) < tol);
    if (ml~=naive) && ~(tol > rel_err(ml, naive))
        fprintf('    Naive method fails on x = %g\n', x);
    end
end
fprintf('OK\n');


fprintf('Testing large values\n');

test_vals = [1, 10, 100, 1000, 10000, 1e99, realmax, Inf];

contender = @(x) x + log1p(exp(-x));

fprintf('Should see some failures from naive method\n');
fprintf('Otherwise there would not be any point in my code\n');
for x = test_vals
    mine = log1pexp(x);
    compare = contender(x);
    naive = log(1 + exp(x));
    if mine ~= compare % Might be both Inf, which rel_err barfs on
        assert(rel_err(mine, compare) < tol);
    end
    if (compare~=naive) && ~(tol > rel_err(compare, naive))
        fprintf('    Naive method fails on x = %g\n', x);
    end
end
fprintf('OK\n');

assert(isnan(log1pexp(NaN)));
fprintf('OK\n');

test_exit