function y = zeta(p)
% ZETA - compute the sum of 1/k^p from k=1 to infinity.
%
% y = zeta(p) returns the zeta function at p
%
% Inputs:
%  p - the power of the reciprocal sum (p>1)
% Outputs:
%  y - the answer
%
% Barnett 4/1/13

if p<=1, error('p is not greater than 1'); end       % catch bad input data
n = 1e9;  % this is lame: really you should figure out how many terms needed!
y = 0;
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 1 number y
end
