function y = zeta(p)
% ZETA - estimate the sum of 1/k^p from k=1 to infinity.
%
% y = zeta(p) returns the zeta function for p>1. It approximates with direct
%  summation using the minimum number of terms needed, but limited to 1e9
%  terms, and hence is not very accurate for p around 2 or less.
%  
% Inputs:
%  p - the power of the reciprocal sum (p>1)
% Outputs:
%  y - the answer
%
% Barnett 4/1/13, auto choose n 4/6/14

nmax = 1e9;                       % max number of terms you're willing to do
tol = 1e-16;                      % desired accuracy (maxn may prevent it)
if p<=1, error('p is not greater than 1'); end       % catch bad input data

n = ceil(tol^(1/(-p+1))); % use n=O(eps^-1/(p-1)) integral estimate as in HW1
if n>nmax                % but make sure not too big; if so clip & warn user
  n=nmax;
  warning(sprintf('small p: choosing n=%d for error %.3g',n,n^(-p+1)))
end

y = 0;               % running sum
for k=n:-1:1         % notice summing from smallest to largest
  y = y + 1/k^p;     % notice we don't use a bunch of RAM, just one number y
end
