C INTPI3.F OF N.TAJIMA'S FORTRAN BENCH-MARK TESTS (VER.2)
C CALCULATES INTEGRAL OVER 3-DIM SPHERE REGION R**2=X**2+Y**2+Z**2<N**2
C OF INTEGRAND (N**2 -R**2 +3/4) WITH MID-POINT RULE (STEP SIZE=1)
C USING ONLY 4-BYTE INTEGERS. "PARM1" IS DUMMY AND MUST BE ZERO.
C HISTORY: 98/6/25-7/1
      PROGRAM INTPI3
      IMPLICIT INTEGER (A-Z)
      PARM1=0
      CALL INTPI(PARM1)
      END
C
      SUBROUTINE INTPI(PARM1)
      IMPLICIT INTEGER (A-Z)
      PARAMETER(N=3000,NM1=N-1,M=1000000,MM=1000000000)
      N1=0
      N2=0
      N3=0
      S1=N**2
      DO 30 I=0,NM1
        S2=S1-I**2-I
        F=6+PARM1
        DO 20 J=0,I
          S3=S2-J**2-J
          IF(J .EQ. I) F=3
          DO 10 K=0,J-1
            S3=S3-2*K
            IF(S3 .LE. 0) GO TO 15
            N1=N1+S3*F
            IF(N1 .GT. MM) THEN
              NT=N1/M
              N2=N2+NT
              N1=N1-NT*M
            END IF
   10     CONTINUE
   15     S3=S3-2*J
          IF(S3 .GT. 0) N1=N1+S3*(F/2)
          NT=N1/M
          N2=N2+NT
          N1=N1-NT*M
   20   CONTINUE
        NT=N2/M
        N3=N3+NT
        N2=N2-NT*M
   30 CONTINUE
      WRITE(6,900) N3,N2,N1,N
  900 FORMAT(' PI=',I12,I6.6,I6.6,'*120/(8*N**5+15*N**3), N=',I7)
C*    WRITE(6,'(A,F20.16,A,1PE12.3)') ' PI='
C*   & ,((DBLE(N3)*M+N2)*M+N1)*120/(DBLE(N)**5*8+DBLE(N)**3*15)
C*   & ,'  ERROR='
C*   & ,((DBLE(N3)*M+N2)*M+N1)*120/(DBLE(N)**5*8+DBLE(N)**3*15)
C*   & -4*ATAN(1.0D0)
      END