function [w,c] = SummationFormula(n,alpha)

% Calculates weights w and nodes c for a summation formula that
% approximates sum(f(k),k=1..infinity) by sum(w(k)*f(c(k)),k=1..n).
% Works well if f(k) is asymptotic to k^(-alpha) for large k.
% Put alpha = 'exp' to use an exponential formula.

switch alpha
    case 'exp'
        k=n:-1:2;
        u=(k-1)/n;
        c1 = exp(2./(1-u).^2-1./u.^2/2); c = k + c1;
        w = 1 + c1.*(4./(1-u).^3+1./u.^3)/n;
        kinf = find(c==Inf);
        c(kinf)=[]; w(kinf)=[];
        c = [c 1]; w = [w 1];
    otherwise
        n = ceil(n/2); a1 = (alpha-1)/6; a6 = 1+1/a1;
        k2 = n-1:-1:0;
        w2 = n^a6./(n-k2).^a6-a6*k2/n;
        c2 = n+a1*(-n+n^a6./(n-k2).^(1/a1))-a6*k2.^2/2/n; 
        k1 = n-1:-1:1;
        w = [w2 ones(size(k1))];
        c = [c2           k1  ];
end

return