C "@(#)pirats.F 1.2 92/05/29 SMI"

#define STDERR 0
#define STDOUT 6
#define STDIN 5
#ifdef SP
#define REAL real*4
#define EXPO e0
#endif
#ifdef DP
#define REAL real*8
#define EXPO d0
#endif
#ifdef QP
#define REAL real*16
#define EXPO q0
#endif

#ifdef TRIGP
#define SIN sinp
#define COS cosp
#define TAN tanp
#define ASIN asinp
#define ACOS acosp
#define ATAN atanp
#else
#define SIN sin
#define COS cos
#define TAN tan
#define ASIN asin
#define ACOS acos
#define ATAN atan
#endif

      PROGRAM PIRATS
C
C.... PIRATS                     Copyright (C) by W. Kahan, Mar. 8, 1989
C
C.... This program finds close rational approximations
C.... A/B to PI, and computes the differences C = PI-A/B to
C.... working precision within errors no bigger than E.  Working
C.... precision arithmetic is presumed to carry no more than about
C.... 110 sig. dec., and to be rounded reasonably like a DEC VAX
C.... or an HP calculator or in conformity with IEEE 754/854.
C.... Then the program tests argument reductions in functions SIN,
C.... COS and TAN, checking to see whether they are consistent
C.... with the hypothesis that their actual periods are some single
C.... close approximation 2*P to 2*PI if not 2*PI itself.
C
      INTEGER NPI
      PARAMETER (NPI = 210)
C
C.... The following array of NPI+1 divisors D(J) > 0 in the continued
C.... fraction PI = D(0)+1/(D(1)+1/(D(2)+1/(D(3)+...))) will incur an
C.... error in PI no worse than 3.6E-234. Only D(NPI) isn't an integer.
C.... This data is based upon DJ's computed by Stewart McDonald
C.... in 1983 from an algorithm due to W. Gosper.  W. K.
C
      INTEGER I,J,INC
      REAL D(0:NPI),EPS,RDX,FL,A,B,C,E,R
      DATA (D(I),I = 0,99)/
     1   3. EXPO ,7. EXPO ,15. EXPO ,1. EXPO ,292. EXPO ,1. EXPO ,
     1	1. EXPO ,1. EXPO ,2. EXPO ,1. EXPO ,
     2   3. EXPO ,1. EXPO ,14. EXPO ,2. EXPO ,1. EXPO ,1. EXPO ,
     1	2. EXPO ,2. EXPO ,2. EXPO ,2. EXPO ,
     3   1. EXPO ,84. EXPO ,2. EXPO ,1. EXPO ,1. EXPO ,15. EXPO ,
     1	3. EXPO ,13. EXPO ,1. EXPO ,4. EXPO ,
     4   2. EXPO ,6. EXPO ,6. EXPO ,99. EXPO ,1. EXPO ,2. EXPO ,
     1	2. EXPO ,6. EXPO ,3. EXPO ,5. EXPO ,
     5   1. EXPO ,1. EXPO ,6. EXPO ,8. EXPO ,1. EXPO ,7. EXPO ,1. EXPO,
     1	2. EXPO ,3. EXPO ,7. EXPO ,
     6   1. EXPO ,2. EXPO ,1. EXPO ,1. EXPO ,12. EXPO ,1. EXPO ,
     1	1. EXPO ,1. EXPO ,3. EXPO ,1. EXPO ,
     7   1. EXPO ,8. EXPO ,1. EXPO ,1. EXPO ,2. EXPO ,1. EXPO ,6. EXPO,
     1	1. EXPO ,1. EXPO ,5. EXPO ,
     8   2. EXPO ,2. EXPO ,3. EXPO ,1. EXPO ,2. EXPO ,4. EXPO ,4. EXPO,
     1 	16. EXPO ,1. EXPO ,161. EXPO , 45. EXPO ,1. EXPO ,22. EXPO ,
     1	1. EXPO ,2. EXPO ,2. EXPO,1. EXPO,4. EXPO ,1. EXPO ,2. EXPO ,
     2  24. EXPO ,1. EXPO ,2. EXPO ,1. EXPO ,3. EXPO ,1. EXPO ,2. EXPO, 
     1  1. EXPO ,1. EXPO ,10. EXPO /
      DATA (D(I),I = 100,199)/
     1   2. EXPO ,5. EXPO ,4. EXPO ,1. EXPO ,2. EXPO ,2. EXPO ,8. EXPO ,
     1	1. EXPO ,5. EXPO ,2. EXPO ,
     2   2. EXPO ,26. EXPO ,1. EXPO ,4. EXPO ,1. EXPO ,1. EXPO ,8. EXPO 
     1	,2. EXPO ,42. EXPO ,2. EXPO ,
     3   1. EXPO ,7. EXPO ,3. EXPO ,3. EXPO ,1. EXPO ,1. EXPO ,7. EXPO ,
     1	2. EXPO ,4. EXPO ,9. EXPO ,
     4   7. EXPO ,2. EXPO ,3. EXPO ,1. EXPO ,57. EXPO ,1. EXPO ,18. EXPO
     1	 ,1. EXPO ,9. EXPO ,19. EXPO ,
     5   1. EXPO ,2. EXPO ,18. EXPO ,1. EXPO ,3. EXPO ,7. EXPO ,30. EXPO
     1	 ,1. EXPO ,1. EXPO ,1. EXPO ,
     6   3. EXPO ,3. EXPO ,3. EXPO ,1. EXPO ,2. EXPO ,8. EXPO ,1. EXPO ,
     1	1. EXPO ,2. EXPO ,1. EXPO ,
     7   15. EXPO ,1. EXPO ,2. EXPO ,13. EXPO ,1. EXPO ,2. EXPO ,1. EXPO
     1	 ,4. EXPO ,1. EXPO ,12. EXPO ,
     8   1. EXPO ,1. EXPO ,3. EXPO ,3. EXPO ,28. EXPO ,1. EXPO ,10. EXPO
     1	 ,3. EXPO ,2. EXPO ,20. EXPO ,
     9   1. EXPO ,1. EXPO ,1. EXPO ,1. EXPO ,4. EXPO ,1. EXPO ,1. EXPO ,
     1	1. EXPO ,5. EXPO ,3. EXPO , 2. EXPO ,1. EXPO ,6. EXPO ,1. EXPO ,
     1	4. EXPO ,1. EXPO ,120. EXPO ,2. EXPO ,1. EXPO ,1. EXPO /
      DATA (D(I),I = 200,NPI)/
     1   3. EXPO ,1. EXPO ,23. EXPO ,1. EXPO ,15. EXPO ,1. EXPO ,3. EXPO
     1	 ,7. EXPO ,1. EXPO ,16. EXPO ,
     2   1.338371961073448 EXPO 
     3   /
C
C.... Environmental constants EPS, RDX, FL
C
      CALL ENVCNS(EPS,RDX,FL)
C
C.... for adequate accuracy
C
      INC = INT(LOG(FL)/LOG(2.6 EXPO ))
C
      PRINT 1000
1000  FORMAT(1X,'PIRATS: computes  PI = A/B+C and R = PI-P.'/
     1   /1X,1X,'J',1X,18X,'A',18X,'/',18X,'B'
     2   /1X,2X,1X,27X,'+',5X,'C'
     3   /1X,2X,1X,41X,'R',5X,'+/-',5X,'E'
     4   )
C
C.... Initialize A, B
C
      A = D(0)
      B = 1. EXPO 
C
      DO 100 J = 1,NPI-INC+5
C
C.... Get next pair A, B
C
         CALL NEXTAB(D(J),A,B)
         IF (A.GT.FL.OR.B.GT.FL) THEN
	    GOTO 110
	 ENDIF
C
C....    Get C = PI-A/B+/-E
C
         CALL GETCE(NPI,J,INC,EPS,FL,D,C,E)
C
C....    Get R = PI-P+/-E
C
         CALL TSTRIG(EPS,A,B,C,E,R)
C
C....    Display these numbers
C
         CALL DISPLA(J,A,B,C,E,R)
C
C....    Test them for consistency
C
         CALL COMPAR(E,R)
100   CONTINUE
110   CONTINUE
C
C.... Normal termination.
C
      call ucbpass
      STOP
      END
C
	REAL function volatilizer(x)
	REAL x
	volatilizer=x
	end

      SUBROUTINE ENVCNS(EPS,RDX,FL)
      REAL EPS,RDX,FL
C
C.... Environmental constants.  This subroutine computes
C....    EPS = nextafter(1,2)-1 = 1.000...001-1
C....    RDX = RaDiX of floating-point arithmetic (2, 10, 16)
C....    FL = RDX/EPS = last of consecutive floating-point integers
C....                       among 1, 2, 3, ..., FL-2, FL-1, FL.
C
C.... local variables
C
	REAL volatilizer
      REAL H,X,Y,T,U,V
C
C.... First seek EPS
C
      T = 4. EXPO /3. EXPO 
      X = T-1. EXPO 
      Y = ABS((X+X-1. EXPO )+X)/64. EXPO 
      EPS = 0. EXPO 
      IF (Y.EQ.0. EXPO ) THEN
	 PRINT *,'Is 4/3 exact?'
         GOTO 299
      ENDIF
200   IF (EPS.NE.0. EXPO ) GOTO 210
         U = volatilizer(1. EXPO +Y)
         EPS = volatilizer(U-1. EXPO) 
         Y = Y+Y
      GOTO 200
C
C     Now seek EPS/RDX = 1-nextafter(1,0) = 1-0.999...999 :
C
210   H = 1. EXPO /2. EXPO 
      T = 2. EXPO /3. EXPO 
      X = T-H
      Y = ABS((X+X-H)+X)/64. EXPO 
      V = 0. EXPO 
      IF (Y.EQ.0. EXPO ) THEN
	 PRINT *,'Is 2/3 exact?'
         GOTO 299
      ENDIF
220   IF (V.NE.0. EXPO ) GOTO 230
         U = volatilizer((H-Y)+H)
         V = volatilizer((H-U)+H)
         Y = Y+Y
      GOTO 220
C
C.... in case Division is dirty
C
230   RDX = AINT(EPS/V+0.0001 EXPO )
      IF (RDX.LT.2. EXPO ) THEN
	 PRINT 5000,'Radix =',RDX
5000     FORMAT(1X,A,F4.0)
         GOTO 299
      ENDIF
C
C.... Confirm that RDX = Radix of Floating-point arithmetic.
C
      T = RDX
      X = 1. EXPO 
C
C.... until X.EQ.0 or X.EQ.RDX
C
240   IF (X.NE.1. EXPO ) GOTO 250
         T = T+T
         U = volatilizer(T+1. EXPO)
         X = volatilizer(U-T)
      GOTO 240
250   IF (X.EQ.0. EXPO ) THEN
         Y = 1. EXPO 
260      IF (X.NE.0. EXPO ) GOTO 270
            Y = Y+Y
            U = volatilizer(T+Y)
            X = volatilizer(U-T)
	 GOTO 260
      ENDIF
270   IF (X.NE.RDX) THEN
         PRINT 6000,'Is Radix ',X,' or ',RDX,'?'
6000     FORMAT(1X,A,F4.0,A,F4.0,A)
	 GOTO 299
      ENDIF
C
C.... Confirm that FL = RDX/EPS:
C
      FL = RDX
      X = 1. EXPO 
280   IF (X.NE.1. EXPO ) GOTO 290
         FL = volatilizer(FL*RDX)
	 U = volatilizer(FL+1. EXPO)
	 X = volatilizer(U-FL)
      GOTO 280
290   IF (FL*V.EQ.1. EXPO ) THEN
	 RETURN
      ENDIF
C
C.... ENVCNS cannot compute environmental constants correctly:
C
      PRINT 3000,'Is FL ',FL,' or ',1. EXPO /V,'?'
3000  FORMAT(1X,A6,1PE45.37E3/1X,2X,A4,1PE45.37E3,A1)
299   PRINT *,'Subroutine ENVCNS cannot compute correctly the'
      PRINT *,'Environmental constants'
      PRINT *,'EPS = 1.000...001 - 1 , and'
      PRINT *,'FL = Last consecutive Floating-point integer'
      PRINT *,'        among 1, 2, 3, ..., FL-2, FL-1, FL.'
      PRINT *,'Please substitute them for the subroutine.'
	call ucbfail
      STOP
      END
C
      SUBROUTINE NEXTAB(DJ,A,B)
      REAL DJ,A,B
C
C.... Get next pair A, B
C
C.... local variables
C
      REAL T,A0,B0
      SAVE A0,B0
      DATA A0,B0/1. EXPO ,0. EXPO /
C
      T = DJ*B+B0
      B0 = B
      B = T
      T = DJ*A+A0
      A0 = A
      A = T
C
C.... Now A/B = D0+1/(D1+1/(D2+...+1/DJ)).
C
      RETURN
      END
C
      SUBROUTINE GETCE(NPI,J,INC,EPS,FL,D,C,E)
      INTEGER NPI,J,INC
      REAL EPS,FL,D(0:NPI),C,E
C
C.... This subroutine computes the continued fraction's tail
C....    Z = D(J+1)+1/(D(J+2)+1/(D(J+3)+1/(D(J+4)+...)))
C.... to working accuracy by using INC terms of it, and then
C.... computes the effect C of cutting it off to get A/B .
C....
C.... Get  C = PI-A/B+/-E
C
C.... local variables
C
      INTEGER I,I9
      REAL X,Y,Z
C
      I9 = MIN(NPI,J+INC)
      Z = D(I9)
      DO 400 I = I9-1,J+1,-1
         Z = D(I)+1. EXPO /Z
400   CONTINUE
      X = FL*FL
C
C.... C = 1/Z-1/X always
C
      C = 1. EXPO /Z-1. EXPO /X
      DO 410 I = J,1,-1
	 Y = D(I)
         Z = Y+1. EXPO /Z
	 X = Y+1. EXPO /X
	 C = -C/(X*Z)
410   CONTINUE
C
C.... E > accumulated roundoff (mixed arithmetic)
C
      E = 4. EXPO *J*EPS*ABS(C)
      RETURN
      END
C
      SUBROUTINE TSTRIG(EPS,A,B,C,E,R)
      REAL EPS,A,B,C,E,R
C
C.... Get R = PI-P+/-E
C
C     This subroutine tests whether the programs that compute
C     TRIG(X) = trig(X*PI/P) for TRIG = SIN, COS or TAN
C     always use the same approximation P to PI during their
C     argument reduction.  If so, 3 or 4 values R = Q+C
C     derived from A = B*(PI-C+/-E) ought to agree.
C
C.... local variables
C
      REAL Q,Q0,S,W,W0,X,Y
      CHARACTER*11 TS
C
      REAL FNODD
      CHARACTER FNS
C
      CHARACTER*10 SI,CO,TA
      DATA SI,CO,TA/'arcsin(sin','arcsin(cos','arctan(tan'/
C
C.... FNODD(floating-point integer X) = (-1)**X
C
      FNODD(X) = (AINT(ABS(X)*0.5 EXPO +0.125 EXPO )*2. EXPO -ABS(X))*
     1	2. EXPO +1. EXPO 
C
C.... FNS(1) = '+', FNS(-1) = '-'
C
      FNS(X) = CHAR(44-INT(X))
C
      Q = ATAN(TAN(A))/B
      R = Q+C
      W = 3. EXPO *EPS*ABS(Q)
      E = E+W
      S = FNODD(B)
      X = A
      Y = B
      TS = FNS(S)//SI
      Q0 = ASIN(SIN(X))/Y*S
      W0 = W+6. EXPO *EPS*ABS(Q0)
      CALL WHET(Q,Q0,W0,X,Y,A,B,TS,TA)
      X = A*0.5 EXPO 
      Y = B*0.5 EXPO 
      IF (S.LT.0. EXPO ) THEN
C
C....    (B+1) is even
C
         S = FNODD((B+1. EXPO )*0.5 EXPO )
         TS = FNS(S)//CO
         Q0 = ASIN(COS(X))/Y*S
         W0 = W+6. EXPO *EPS*ABS(Q0)
      ELSE
C
C....    B = 2y is even
C
         TS = ' '//TA
         Q0 = ATAN(TAN(X))/Y
         W0 = 3. EXPO *EPS*ABS(Q0)
	 CALL WHET(Q,Q0,W0,X,Y,A,B,TS,TA)
         S = FNODD(Y)
         TS = FNS(S)//SI
	 Q0 = ASIN(SIN(X))/Y*S
         W0 = W+6. EXPO *EPS*ABS(Q0)
      ENDIF
      CALL WHET(Q,Q0,W0,X,Y,A,B,TS,TA)
      RETURN
      END
C
      SUBROUTINE WHET(Q,Q0,W0,X,Y,A,B,TS,TA)
      REAL Q,Q0,W0,X,Y,A,B
      CHARACTER*(*) TS,TA
C
C.... Test whether Q0.EQ.Q within +/- W0 (.GE.->.GT.)
C
      IF (ABS(Q-Q0).GT.W0) THEN
C
C....    Difference too big suggests P is not a constant after all.
C
         PRINT 4000,TS,X,')/',Y,' =',Q0,' differs from ',
     1      TA,A,')/',B,' =',Q,' too much.'
4000  FORMAT(/1X,4X,70('%')
     1   /1X,4X,A11,0PF38.1,A
     2   /1X,4X,11X,0PF38.1,A
     3   /1X,4X,11X,1PE45.37E3,A
     4   /1X,4X,1X,A10,0PF38.1,A
     5   /1X,4X,11X,0PF38.1,A
     6   /1X,4X,11X,1PE45.37E3,A
     7   /1X,4X,70('%')
     8   /)
	call ucbfail
      ENDIF
      RETURN
      END
C
      SUBROUTINE DISPLA(J,A,B,C,E,R)
      INTEGER J
      REAL A,B,C,E,R
C
C     Display Formatting
C
C.... display  J, A, B, C, R, E
C
      PRINT 2000,J,A,B,C,R,E
2000  FORMAT(1X,I2,1X,F37.0,'/',F37.0
     1   /1X,2X,1X,27X,'+',1X,1PE45.37E3
     2   /1X,2X,1X,1X,1PE45.37E3,1X,'+/-',1PE16.8E3
     3   )
      RETURN
      END
C
      SUBROUTINE COMPAR(E,R)
      REAL E,R
C
C.... test for consistency
C
C.... local variables
C
      REAL E0,R0
      SAVE E0,R0
      DATA E0,R0/0.1416 EXPO ,0. EXPO /
C
C.... .LT.->.LE.
C
      IF (ABS(R-R0).LE.E+E0) THEN
	 E0 = E
         R0 = R
      ELSE
         PRINT *,'Varying R implies defective argument reduction.'
	call ucbfail
         STOP
      ENDIF
      RETURN
      END
