PROGRAM KYUSU
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (PI=3.14159265358979323846,ICT=4)
      DOUBLE PRECISION T(0:ICT-1,ICT),M
      N=10
      DO 400 K=1,4
          T(0,K)=0
          DO 300 I=0,N
              DO 200 J=0,N
                  T(0,K)=T(0,K)+1.0/(2*I+1)/
     +             (2*J+1)/((2*I+1)**2+(2*J+1)**2)*(-1)**(I+J)
  200         CONTINUE
  300     CONTINUE
          T(0,K)=T(0,K)*80/PI**4
          WRITE(*,'(A2,I6,A6,F18.15)')'N=',N,', SUM=',T(0,K)
          N=N*10
  400 CONTINUE
      WRITE(*,*)'Richardson extrapolation of 1st order:'
      M=10
      DO 3000 I=1,ICT-1
         DO 2000 I2=1,ICT-I
            T(I,I2)=(M*T(I-1,I2+1)-T(I-1,I2))/(M-1)
 2000    CONTINUE
         M=M*10
 3000 CONTINUE
      DO 5000 I=0,ICT-1
          WRITE(*,'(5F18.15)')(T(I-I2+1,I2),I2=1,I+1)
 5000 CONTINUE
      WRITE(*,*)'Error:'
      WRITE(*,'(4F19.15)')((T(0,I)-T(ICT-1,1)),I=1,ICT)
      WRITE(*,*)'Ratio:'
      WRITE(*,'(4F19.15)')((T(0,I+1)-T(ICT-1,1))/
     +            (T(0,I)-T(ICT-1,1)),I=1,ICT-1)
      M=1000
      DO 8000 I=1,ICT-1
         DO 7000 I2=1,ICT-I
            T(I,I2)=(M*T(I-1,I2+1)-T(I-1,I2))/(M-1)
 7000    CONTINUE
         M=M*1000
 8000 CONTINUE
      WRITE(*,*)'Richardson extrapolation of 3rd order:'
      DO 9000 I=0,ICT-1
          WRITE(*,'(5F18.15)')(T(I-I2+1,I2),I2=1,I+1)
 9000 CONTINUE
      WRITE(*,*)'Error:'
      WRITE(*,'(4F19.15)')((T(0,I)-T(ICT-1,1)),I=1,ICT)
      WRITE(*,*)'Ratio:'
      WRITE(*,'(4F19.15)')((T(0,I+1)-T(ICT-1,1))/
     +                 (T(0,I)-T(ICT-1,1)),I=1,ICT-1)
      END