	PROGRAM CALPI
C	This program calculates pi using a method described in 
C	Scientific American February 1988, page 115.
C
C	By  Mark P. Esplin
C	    16 Shawsheen Rd.
C	    Billerica, Ma 01821	
C
	INTEGER*2 Y(3000),ALP(3000)
	INTEGER*2 T1(3000),T2(3000),T3(3000),T4(3000),T5(3000)
	INTEGER*2 T6(3000)
	INTEGER*2 NT,N
	INTEGER*2 ITER,NITER,ITEMP,IFIL
	CHARACTER ANS*1,FILNAM*20
	DATA NT/3000/
10	WRITE(*,*) 'Enter number of desired digits'
	READ(*,*) NDIG
	N=NDIG/4+2
	IF(N.GT.NT) THEN
	  WRITE(*,*) 'Not enough memory available'
	  GOTO 10
	ENDIF
C	For now number of iterations is 
	WRITE(*,*) 'Enter number of iterations'
	READ(*,*) NITER
	IFIL=0
	WRITE(*,*) 'Log results to disk?'
	READ(*,'(A)') ANS
	IF(ANS.EQ.'y'.OR.ANS.EQ.'Y') THEN
	  IFIL=1
	  WRITE(*,*) 'File name?'
	  READ(*,'(A)') FILNAM
	  OPEN(1,FILE=FILNAM,STATUS='NEW')
	ENDIF
C	Set beginning values for Y and ALP.
C	Y=SQRT(2)-1
	ITEMP=2
	CALL CVI2EX(ITEMP,T1,N)
C	The square root use 4 temporaries, so it is actually using T3 - T6.
	CALL EXSQRT(T1,T2,N,T3)
	ITEMP=-1
	CALL CVI2EX(ITEMP,T1,N)
	CALL ADD(T1,T2,Y,N)
C	ALP=6-4*SQRT(2)
	ITEMP=4
	CALL CVI2EX(ITEMP,T1,N)
	CALL MUL(T1,T2,T3,N)
	CALL NEG(T3,N)
	ITEMP=6
	CALL CVI2EX(ITEMP,T1,N)
	CALL ADD(T1,T3,ALP,N)
C	Main loop
	DO 20 ITER=1,NITER
C	     T1=SQRT(SQRT(1-Y**4))
	  CALL MUL(Y,Y,T1,N)
	  CALL MUL(T1,T1,T2,N)
	  CALL NEG(T2,N)
	  ITEMP=1
	  CALL CVI2EX(ITEMP,T3,N)
	  CALL ADD(T2,T3,T1,N)
	  CALL EXSQRT(T1,T2,N,T3)
	  CALL EXSQRT(T2,T1,N,T3)	
C	     Y=(1-T1)/(1+T1)
	  CALL COPY(T1,T2,N)
	  CALL NEG(T2,N)
	  ITEMP=1
	  CALL CVI2EX(ITEMP,T3,N)
	  CALL ADD(T2,T3,T4,N)
	  CALL ADD(T1,T3,T2,N)
	  CALL COPY(T4,T1,N)
C	        INV using T4-T6
	  CALL INV(T2,T3,N,T4)
	  CALL MUL(T1,T3,Y,N)
C	     T1=1+Y
	  ITEMP=1
	  CALL CVI2EX(ITEMP,T2,N)
	  CALL ADD(T2,Y,T1,N)
C	    T2=ALP*(1+Y)**4
	  CALL MUL(T1,T1,T2,N)
	  CALL MUL(T2,T2,T3,N)
	  CALL MUL(ALP,T3,T2,N)
C	    T1=-(2**(2*(ITER-1)+3))*Y*(1+Y+Y**2)
	  CALL MUL(Y,Y,T3,N)
	  CALL ADD(T1,T3,T4,N)
	  CALL MUL(T4,Y,T5,N)
	  ITEMP=2**(2*(ITER-1)+3)
	  CALL CVI2EX(ITEMP,T6,N)
	  CALL MUL(T6,T5,T1,N)
	  CALL NEG(T1,N)
C	     ALP=T2+T1
	  CALL ADD(T2,T1,ALP,N)
C	     PI=1/ALP
	  CALL INV(ALP,T1,N,T2)
	  WRITE(*,*)
	  WRITE(*,900) ITER
900	  FORMAT(' Iteration number ',I3)
	  IF(IFIL.EQ.1) THEN
	    WRITE(1,*)
	    WRITE(1,900) ITER
	    ITEMP=1
	    CALL WRTEXT(ITEMP,T1,N)
	  ELSE
	    CALL PRTEXT(T1,N)
	  ENDIF
20	CONTINUE
	STOP
	END
