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

% function [w,c]=SummationFormula(n,alpha)
%
% calculatesweights 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 asy,ptotic 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