%EXAMPLE1
%
f0='%10.0f %4.0f %6.0f %18.12f\n';
dig=48; digits(dig);
N=20;
ab0=sr_jacobi01(dig,N,-1/2,0);
disp('       beta   m      n           I') 
for beta=[50 100 200 400]
  for m=[0 2 4 6 8]
    if m==0
      ab=subs(ab0);
    else
      z(1)=vpa(pi*i/beta); z(2)=conj(z(1));
      z(3)=vpa(3*pi*i/beta); z(4)=conj(z(3));
      z(5)=vpa(5*pi*i/beta); z(6)=conj(z(5));
      z(7)=vpa(7*pi*i/beta); z(8)=conj(z(7));
      ab=mod_md(N,m,z,ab0,@rho0);
      ab=subs(real(ab));
    end
    for n=4:4:N
      xw=gauss(n,ab); x=xw(:,1); w=xw(:,2);
      I=sum(w.*f1(x,m,beta));
      fprintf(f0,beta,m,n,I)
    end
    fprintf('\n')
  end
end
