% Relative error in the computation of the harmonic sum
npts=25;
Ns=zeros(1,npts); hsum=zeros(1,npts);
relerr=zeros(1,npts); relerrSP=zeros(1,npts);
nhsum=zeros(1,npts); nhsumSP=zeros(1,npts);
for i=1:npts
    Ns(i)=2^i;
    nhsum(i)=harmonic(Ns(i));
    nhsumSP(i)=harmonicSP(Ns(i));
    hsum(i)=(psi(Ns(i)+1)-psi(1)); % Theoretical result
    relerr(i)=abs(nhsum(i)-hsum(i))/hsum(i);
    relerrSP(i)=abs(nhsumSP(i)-hsum(i))/hsum(i);
end

figure(1);
loglog(Ns,relerr,'ro--', Ns,relerrSP,'bs-');
title('Error in harmonic sum');
xlabel('N'); ylabel('Relative error');
legend('double','single', 'Location','NorthWest');

figure(2);
semilogx(Ns,nhsum,'ro--', Ns,nhsumSP,'bs:', Ns,hsum,'g.-');
title('Harmonic sum');
xlabel('N'); ylabel('H(N)');
legend('double','single','"exact"', 'Location','NorthWest');
