function siam3(n,ndig)

% run g_siam3(200,50) for forty digits of precision

if nargin<2
   ndig = 64;
end
sol = mp('1.27422415282122818821234063972507809947246945663850391770381031749793967672676606068189521687585587774980110443371332140781337107034873203689684677183560860062297082957415624325822189608066128846250820485366710992484957006452161261548043926902807588055894387353313875269029');

L = 5; table = zeros(L-1,2);
alpha = {4,'exp'};


        [w,c]=SummationFormula2(n,ndig);
        m = length(w); w = sqrt(w);
        A = zeros(m,m);
        A = mp(A);
        for j=1:m
            for k=1:m
                A(j,k) = w(j)*w(k)/...
                    ((c(j)+c(k)-1)*(c(j)+c(k))/2-(c(k)-1));
            end
        end
r = mp(0);
r = normest(A,mp(10^(-ndig-1)));


fprintf(['%d | %.',num2str(ndig),'f | %.1e \n'],n,r,abs(sol-r));

