function [ y ] = harmonicpoly( n,x )
%HARMONICPOLY harmonicpoly(n,x) computes the interpolating polynomial with
%   nodes 1,2,..,n and values 1/1, 1/2,...,1/n.
%   Uses the nested form of the polynomial
%   1 - (1/2)(x-1)[ 1 - (1/3)(x-2) [ ... 1 - (1/n-1)(x-(n-1))]...]
%
    y = ones(size(x));
    for i=(n-1):-1:1
        y = 1-(x-i).*y/(i+1);
    end
end