C
C	This is a collection of subroutines that perform extended precision
C	floating point arithmetic.
C
C	By  Mark P. Esplin
C	    16 Shawsheen Rd.
C	    Billerica, Ma 01821	
C
C	Whole programs do not have to be written using extended
C	precision.  Only the part that needs extra precision need use it,
C	then convert back to regular REAL*8 numbers.
C	The routines for addition and multiplication use the same algorithms 
C	that are taught in grade school but using base 10000 numbers.
C	This is not the most  efficient way to do it, but
C	it makes it possible to use standard FORTRAN77 (almost).  
C	These routines should work with any computer and compiler that supports
C	INTEGER*2 and INTEGER*4 arithmetic.
C	Each number is stored in an INTEGER*2 array of length N.
C	The first word, IA(1), is the exponent in base 10000 digits.
C	The least significant digit is in IA(2) and the most significant
C	digit is in IA(N).  If IA(N) > 5000 the number is negative.
C	When IA(N) = 5000 the number is assumed invalid.
C	Base 10000 digits have worse roundoff characteristics than binary.
C	Since roundoff causes the last base 10000 digit to be inaccurate, 
C	that really means the last 4 decimal digits are inaccurate.  Be
C	sure to use more digits than is necessary to allow for roundoff.
C	N must be at least 6 or some routines won't work correctly.
C
C	Subroutines with their FORTRAN equivalent in this package are:
C
C	ADD(A,B,C,N)      C=A+B
C	NEG(A,N)          A=-A
C	MUL(A,B,C,N)      C=A*B
C	INV(A,B,N,W)      B=1/A      W is work area of size 3*N
C	EXSQRT(A,B,N,W)   B=SQRT(A)  W is work area of size 4*N
C	EXHALF(A,N)       A=1/2      much faster than using INV
C	PRTEXT(A,N)       WRITE(*,*) A
C	WRTEXT(IU,A,N)    WRITE(IU,*) A
C	CVI2EX(I2,A,N)    A=I2       I2 is INTEGER*2
C	CVR8EX(R,A,N)     A=R        R is REAL*8
C	CVEXR8(A,N,R)     R=A        R is REAL*8
C	COPY(A,B,N)       B=A
C	IEXCMP(A,B,N,W)    Compares two positive numbers A and B returns 1, 0, -1
C	CHGPRE(A,NA,B,NB)  Changes the precision of A form NA to NB
C
C	There are other internal routines, but they would probably be of limited
C	value for most users.
C
C
	SUBROUTINE INCEX(IA,N)
	INTEGER*2 IA(*),N
C	This routine increments the mantissa of extended number IA.
C	This routine is used by other routines to perform rounding.
	INTEGER*2 IT,ICAR,I
	ICAR=1
	DO 10 I=2,N
	  IT=IA(I)+ICAR
	  ICAR=IT/10000
	  IA(I)=IT-ICAR*10000
C	In almost all cases this routine will end after a few passes.
	IF(ICAR.EQ.0) GOTO 20
10	CONTINUE
C	Shift back if necessary
C	The only time that a shift is necessary is when IA(N) = 5000 and
C	the rest of the mantissa is zero.
20	CONTINUE
	IF(IA(N).EQ.5000) THEN
	  IA(1)=IA(1)+1
	  DO 30 I=2,N-1
	    IA(I)=IA(I+1)
30	  CONTINUE
	  IA(N)=0
C	  Don't need to round again since if a shift was necessary
C	  the least significant digit is zero.
	ENDIF
	RETURN
	END
C
C
C
C
C
	FUNCTION ICKCON(IA,N)
C	This routine is used to check if an extended number IA is converging
C	to an integer value. 
C	This routine would probably not be very useful for the average user.
C	The value returned is 1 if IA has converged.  The list significant digit
C	is not checked since it might be just round off.
C	This function return a 1 in IA has converged to an integer otherwise
C	it returns a 0.
	INTEGER*2 IA(*)
	INTEGER*2 N,I
	INTEGER*2 ITEST
	ICKCON=0
C	The number must have already converged past the N-2 th digit before
C	this routine is called.
	ITEST=IA(N-2)
	DO 10 I=3,N-3
	  IF(IA(I).NE.ITEST) RETURN
10	CONTINUE
	ICKCON=1
	RETURN
	END
	SUBROUTINE NEG(IA,N)
C	This routine changes the sign of extended precision number IA.
	INTEGER*2 IA(*),N
	INTEGER*2 I,J
	J=2
	DO 10 I=2,N
	  IF(IA(J).NE.0) GOTO 20
	  J=J+1
10	CONTINUE
C	Number was zero
	RETURN
C	The first non-zero digit has been reached
20	IA(J)=10000-IA(J)
	J=J+1
	DO 30 I=J,N
	  IA(I)=9999-IA(I)
30	CONTINUE
	RETURN
	END
C
C
C
C
	SUBROUTINE PRTEXT(IA,N)
C	This routine writes extended precision numbers to output.
C	This routine could use some work to make a more readable output.
	INTEGER*2 IA(*),N
	INTEGER*2 I,INEG
	INEG=0
	WRITE(*,*) ' '
	IF(IA(N).GT.5000) THEN
	  WRITE(*,*) 'Minus'
	  CALL NEG(IA,N)
	  INEG=1
	ENDIF
	WRITE(*,900) IA(1)
900	FORMAT(' Exponent  ',I6)
	WRITE(*,910) (IA(I),I=N,2,-1)
910	FORMAT(10I5)
C	Restore IA to previous value.
	IF(INEG.NE.0) CALL NEG(IA,N)
	RETURN
	END
C
C
C
C
	SUBROUTINE WRTEXT(IUNIT,IA,N)
C	This routine writes extended precision numbers to file IUNIT.
C	This routine could use some work to make a more readable output.
	INTEGER*2 IUNIT,IA(*),N
	INTEGER*2 I,INEG
	INEG=0
	WRITE(IUNIT,*) ' '
	IF(IA(N).GT.5000) THEN
	  WRITE(IUNIT,*) 'Minus'
	  CALL NEG(IA,N)
	  INEG=1
	ENDIF
	WRITE(IUNIT,900) IA(1)
900	FORMAT(' Exponent  ',I6)
	WRITE(IUNIT,910) (IA(I),I=N,2,-1)
910	FORMAT(10I5)
C	Restore IA to previous value.
	IF(INEG.NE.0) CALL NEG(IA,N)
	RETURN
	END
C
C
C
C
C
	SUBROUTINE NORMIZ(IA,N)
C	This routine is used by other routines to normalize
C	floating point number after addition and other such processes.
	INTEGER*2 IA(*),N
	INTEGER*2 I,ISHFT
	IF(IA(N).GT.5000) THEN
C	  This is a negative number.
	  DO 10 I=N,2,-1
	    IF(IA(I).NE.9999) GOTO 20
10	  CONTINUE
	  I=I+1
20	  ISHFT=N-I
	  IF(IA(I).LE.5000) ISHFT=ISHFT-1
	ELSE
C	  This is a positive number.
	  DO 30 I=N,2,-1
	    IF(IA(I).NE.0) GOTO 40
30	  CONTINUE
C	  This number is zero
	  RETURN
40	  ISHFT=N-I
	  IF(IA(I).GE.5000) ISHFT=ISHFT-1
	ENDIF
	IF(ISHFT.EQ.0) RETURN
	IA(1)=IA(1)-ISHFT
C	Shift left the appropriate amount.
	DO 50 I=N,ISHFT+2,-1
	  IA(I)=IA(I-ISHFT)
50	CONTINUE
C	Zero unknown digits.
	DO 60 I=2,ISHFT+1
	  IA(I)=0
60	CONTINUE
	RETURN
	END
C
C
C
C
C
	SUBROUTINE ADD(IA,IB,IC,N)
C	This routine performs extended precision addition.
C	IC = IA + IB. IA and IB are not changed.
	INTEGER*2 IA(*),IB(*),IC(*)
	INTEGER*2 N,IPA,IPB,NML,I
	INTEGER*2 ISGNA,ISGNB,ISGNC,ISGEXT
	INTEGER*2 IT,ICAR
	ISGNA=1
	IF(IA(N).GT.5000) ISGNA=-1
	ISGNB=1
	IF(IB(N).GT.5000) ISGNB=-1
C	Figure out pointers for main loop
	IF(IA(1).GE.IB(1)) THEN
	  IC(1)=IA(1)
	  IPA=2
	  IPB=IA(1)-IB(1)+2
	  NML=N-IPB+2
	  IF(NML.LT.2) THEN
C	    The difference between the two numbers is too large.
	    CALL COPY(IA,IC,N)
	    RETURN
	  ENDIF
	  IF(IPB.EQ.2) THEN
	    ICAR=0
	  ELSE
	    ICAR=(IB(IPB-1)+5000)/10000
	  ENDIF
	ELSE
	  IC(1)=IB(1)
	  IPB=2
	  IPA=IB(1)-IA(1)+2
	  NML=N-IPA+2
	  IF(NML.LT.2) THEN
C	    The difference between the two numbers is too large.
	    CALL COPY(IB,IC,N)
	    RETURN
	  ENDIF
	  ICAR=(IA(IPA-1)+5000)/10000
	ENDIF
C	Main addition loop.
	DO 10 I=2,NML
	  IT=IA(IPA)
	  IT=IT+IB(IPB)+ICAR
	  ICAR=IT/10000
	  IC(I)=IT-ICAR*10000
	  IPA=IPA+1
	  IPB=IPB+1
10	CONTINUE	
C	Do sign extended part.
	IF(IA(1).GT.IB(1)) THEN
	  IF(ISGNB.GE.1) THEN
	    ISGEXT=0
	  ELSE
	    ISGEXT=9999
	  ENDIF
	  DO 20 I=NML+1,N
	    IT=IA(I)
	    IT=IT+ISGEXT+ICAR
	    ICAR=IT/10000
	    IC(I)=IT-ICAR*10000
20	  CONTINUE
	ELSE
	  IF(ISGNA.GE.1) THEN
	    ISGEXT=0
	  ELSE
	    ISGEXT=9999
	  ENDIF
	  DO 30 I=NML+1,N
	    IT=IB(I)
	    IT=IT+ISGEXT+ICAR
	    ICAR=IT/10000
	    IC(I)=IT-ICAR*10000
30	  CONTINUE
	ENDIF  
C	Shift right one place if it has overflowed during addition.
	IF(ISGNA*ISGNB.EQ.1) THEN
	  ISGNC=1
	  IF(IC(N).GE.5000) ISGNC=-1
	  IF(ISGNA*ISGNC.EQ.-1) THEN
C	    An overflow has occurred.
	    IC(1)=IC(1)+1
	    ICAR=IC(2)
	    DO 40 I=2,N-1
	      IC(I)=IC(I+1)
40	    CONTINUE
	    IC(N)=ISGEXT
C	    Do rounding and return.
	    IF(ICAR.GE.5000) CALL INCEX(IC,N)
	    RETURN
	  ENDIF
	ENDIF
C	Normalize answer.
	CALL NORMIZ(IC,N)
	RETURN
	END
C
C
C
C
	SUBROUTINE MUL(IA,IB,IC,N)
C	This routine multiplies two extended precision numbers.
C	N must be at least 4 or this routine won't work.
C	IC = IA*IB
	INTEGER*2 IA(*),IB(*),IC(*),N
	INTEGER*2 ISNA,ISNB,IXC1,IXC2,ISHFT
C	The computation is carried out with two extra placesIXC1 and IXC2
C	This makes it so precision is not lost if the answer has to be
C	shifted left.  IXC1 is the least significant.
	INTEGER*4 IT,ICAR
	INTEGER*2 I,J,K,IAS
	ISNA=1
	ISNB=1
	IF(IA(N).GT.5000) THEN
	  CALL NEG(IA,N)
	  ISNA=-1
	ENDIF
	IF(IB(N).GT.5000) THEN
	  CALL NEG(IB,N)
	  ISNB=-1
	ENDIF
	IC(1)=IA(1)+IB(1)
	IXC1=0
	IXC2=0
	IAS=N+1
C	Step through each IB
	DO 30 I=2,N
	  ICAR=0
	  IF(I.LE.(N-2)) THEN
C	    Round digit that is past least significant.
	    IT=IA(IAS-3)
	    IT=IT*IB(I)+5000
	    ICAR=IT/10000
	  ENDIF
	  IF(I.LE.(N-1)) THEN
C	    Do extended precision 1
	    IT=IA(IAS-2)
	    IT=IT*IB(I)+ICAR+IXC1
	    ICAR=IT/10000
	    IXC1=IT-ICAR*10000
	  ENDIF
C	  Do extended precision 2
	  IT=IA(IAS-1)
	  IT=IT*IB(I)+ICAR+IXC2
	  ICAR=IT/10000
	  IXC2=IT-ICAR*10000
C	     The first time through the I loop the J loop is not executed
	  K=2
	  DO 20 J=IAS,N
	    IT=IA(J)
	    IT=IT*IB(I)+ICAR+IC(K)
	    ICAR=IT/10000
	    IC(K)=IT-ICAR*10000
	    K=K+1
20	  CONTINUE
	  IC(K)=ICAR
	  IAS=IAS-1
30	CONTINUE
C	Normalize output if needed rounding when possible.
	IF((IC(N).EQ.0).AND.(IC(N-1).LT.5000)) THEN
	  IF((IC(N-1).EQ.0).AND.(IC(N-2).LT.5000)) THEN
	    ISHFT=2
	  ELSE
	    ISHFT=1
	  ENDIF
	  IC(1)=IC(1)-ISHFT
	  DO 40 I=N,ISHFT+2,-1
	    IC(I)=IC(I-ISHFT)
40	  CONTINUE
	  IF(ISHFT.EQ.2) THEN
C	    No rounding is possible in this case.
	    IC(3)=IXC2
	    IC(2)=IXC1
	  ELSE
	    IC(2)=IXC2
	    IF(IXC1.GE.5000) CALL INCEX(IC,N)
	  ENDIF
	ELSE
	  IF(IXC2.GE.5000) CALL INCEX(IC,N)
	ENDIF
C	Set sign
	IF(ISNA*ISNB.LT.0) CALL NEG(IC,N)
	RETURN
	END
C
C
C
C
	SUBROUTINE EXSQRT(IA,IB,N,IW)
C	This routine does an extended precision square root.
C	IB = SQRT(IA)
C	It uses an iterative method.  See Arithmetic Operations in Digital
C	Computers by R. K. Richards  D. Van Nostrand Company, Inc. Copyright 1955.
C	B(k+1) = B(k)/2 * (3 - B(k)^2/A)  This series only requires one division.
C	IW is a work area of size 4*N
	INTEGER*2 IA(*),IB(*),IW
	INTEGER*2 N,NUSED,NUSEDL
	INTEGER*2 IALD,IMONE,ITHREE
	DIMENSION IW(N,4)
	REAL*8 R
	IALD=0
	IMONE=-1
	ITHREE=3
	IF(IA(N).GE.5000) THEN
	  WRITE(*,*) 'Square root of negative number?'
	  STOP
	ENDIF
C	Find first guess.
	CALL CVEXR8(IA,N,R)
	R=SQRT(R)
	NUSED=6
	CALL CVR8EX(R,IB,NUSED)
C	IW(1) is 1/IA at full precision.
	CALL INV(IA,IW(1,1),N,IW(1,2))
C	Beginning of iterative loop
10	NUSEDL=NUSED
	NUSED=(NUSED-1)*2
	IF(NUSED.GT.N) THEN
	  NUSED=N
	  IALD=1
	ENDIF
	CALL CHGPRE(IB,NUSEDL,IB,NUSED)
	CALL MUL(IB,IB,IW(1,2),NUSED)
C	Return if convergence is complete
	IF(IALD.EQ.1) THEN
	  IF(IEXCMP(IA,IW(1,2),N,IW(1,3)).EQ.0) RETURN
	ENDIF
	CALL CHGPRE(IW(1,1),N,IW(1,3),NUSED)
	CALL MUL(IW(1,2),IW(1,3),IW(1,4),NUSED)
	CALL NEG(IW(1,4),NUSED)
	CALL CVI2EX(ITHREE,IW(1,2),NUSED)
	CALL ADD(IW(1,2),IW(1,4),IW(1,3),NUSED)
C	IW(3) is now equal to   3 - B(k)^2/A
	CALL EXHALF(IW(1,2),NUSED)
	CALL MUL(IW(1,2),IB,IW(1,4),NUSED)
	CALL MUL(IW(1,4),IW(1,3),IB,NUSED)
	GOTO 10
	END
C
C
C
C
C
C
	SUBROUTINE CVI2EX(I2,IA,N)
C	This routine converts an INTEGER*2 number to extended precision.
	INTEGER*2 I2,IA(*)
	INTEGER*2 N,I
	INTEGER*2 ISGN,ITEMP
C	Zero array before start
	DO 10 I=1,N
	  IA(I)=0
10	CONTINUE
	ITEMP=I2
	ISGN=1
	IF(ITEMP.LT.0) THEN
	  ISGN=-1
	  ITEMP=-ITEMP
	ENDIF
	IF(ITEMP.GE.5000) THEN
	  IA(1)=2
	  IA(N)=ITEMP/10000
	  IA(N-1)=ITEMP-IA(N)*10000
	ELSE
	  IA(1)=1
	  IA(N)=ITEMP
	ENDIF
	  IF(ISGN.EQ.-1) CALL NEG(IA,N)
	RETURN
	END
C
C
	SUBROUTINE COPY(IA,IB,N)
C	This routine copies extended number IA to IB.
	INTEGER*2 IA(*),IB(*),N
	INTEGER*2 I
	DO 10 I=1,N
	  IB(I)=IA(I)
10	CONTINUE
	END
C
C
C
C
	SUBROUTINE EXHALF(IA,N)
C	This routine returns an extended precision 1/2.
	INTEGER*2 IA(*)
	INTEGER*2 N,I
	IA(1)=1
	IA(N)=0
	IA(N-1)=5000
	DO 10 I=2,N-2
	  IA(I)=0
10	CONTINUE
	RETURN
	END
C
C
C
C
C
	SUBROUTINE CVR8EX(R,IA,N)
C	Convert REAL*8 R to an extended precision number IA.
C	N must be greater than or equal to 6.
	REAL*8 R,RT
	INTEGER*2 IA(*)
	INTEGER*2 N,I
	INTEGER*2 ISIGN
	ISIGN=1
	RT=R
	IF(R.LT.0) THEN
	  ISIGN=-1
	  RT=-RT
	ENDIF
	IA(1)=LOG(RT)/LOG(10000.)
	RT=RT/(10000.**IA(1))
	IA(1)=IA(1)+1
	IF(RT.GE.5000) THEN
	  IA(1)=IA(1)+1
	  RT=RT/10000.
	ENDIF
	DO 10 I=N,N-4,-1
	  IA(I)=RT
	  RT=RT-IA(I)
	  RT=RT*10000
10	CONTINUE
	DO 20 I=N-5,2,-1
	  IA(I)=0
20	CONTINUE
	IF(ISIGN.LT.0) CALL NEG(IA,N)
	RETURN
	END
C
C
C
	SUBROUTINE CVEXR8(IA,N,R)
C	Convert an extended precision number IA to REAL*8 R.
C	N must be greater than or equal to 6.
	REAL*8 R,RT,EXP
	INTEGER*2 IA(*)
	INTEGER*2 N,I
	INTEGER*2 ISIGN
	ISIGN=1
	IF(IA(N).GE.5000) THEN
	  ISIGN=-1
	  CALL NEG(IA,N)
	ENDIF
	RT=0
	EXP=1
	EXP=(EXP/10000)**5
	DO 10 I=N-4,N
	  RT=RT+IA(I)*EXP
	  EXP=EXP*10000
10	CONTINUE
	EXP=10000
	R=RT*(EXP**IA(1))
	IF(ISIGN.LE.0) R=-R
	RETURN
	END
C
C
C
C
	SUBROUTINE CHGPRE(IA,NA,IB,NB)
C	This routine changes the precision of an extended precision number.
C	This routine is useful for many iterative solutions where the full
C	precision is not needed until the final iteration.
C	IA and IB can be the same array.
C	Input is IA of size NA, output is IB of size NB.
	INTEGER*2 IA(*),IB(*),ICAR
	INTEGER*2 NA,NB,I,J
	IB(1)=IA(1)
	IF(NA.GT.NB) THEN
C	  IA has the most precision.
	  ICAR=0
	  J=NA-NB+1
	  IF(IA(J).GT.5000) ICAR=1
	  DO 10 I=2,NB
	    J=J+1
	    IB(I)=IA(J)
10	  CONTINUE
	  IF(ICAR.EQ.1) CALL INCEX(IB,NB)
	ELSE
C	  IB has the most precision.
	  J=NB
	  DO 20 I=NA,2,-1
	    IB(J)=IA(I)
	    J=J-1
20	  CONTINUE
	  DO 30 I=2,J
	    IB(I)=0
30	  CONTINUE
	ENDIF
	RETURN
	END
C
C
C
C
C
	SUBROUTINE INV(IX,IY,N,IW)
C	This routine does an extended precision inverse.
C	IY=1/IX
C	It uses an iterative method.  A routine more like MUL would probably 
C	be more efficient and would certainly take less storage, but I haven't
C	gotten around to writing it yet.
C	Y(I+1)=Y(I)*(2-X*Y(I))
C	IW is a work area of size 3*N
	INTEGER*2 IX(*),IY(*),IW
	INTEGER*2 N,NUSED,NUSEDL
	INTEGER*2 IALD,ISGN,INEGT
	DIMENSION IW(N,3)
	REAL*8 R
	INEGT=-2
	IALD=0
	ISGN=1
	IF(IX(N).GE.5000) THEN
	  CALL NEG(IX,N)
	  ISGN=-1
	ENDIF
C	Find first guess.
	CALL CVEXR8(IX,N,R)
	R=1/R
	NUSED=6
	CALL CVR8EX(R,IY,NUSED)
10	NUSEDL=NUSED
	NUSED=(NUSED-1)*2
	IF(NUSED.GT.N) THEN
	  NUSED=N
	  IALD=1
	ENDIF
	CALL CHGPRE(IX,N,IW(1,1),NUSED)
	CALL CHGPRE(IY,NUSEDL,IY,NUSED)
	CALL MUL(IW(1,1),IY,IW(1,2),NUSED)
C	Return if convergence is complete
	IF(IALD.EQ.1) THEN
	  IF(ICKCON(IW(1,2),N).EQ.1) THEN
	    IF(ISGN.LT.0) CALL NEG(IY,N)
	    RETURN
	  ENDIF
	ENDIF
	CALL CVI2EX(INEGT,IW(1,1),NUSED)
	CALL ADD(IW(1,1),IW(1,2),IW(1,3),NUSED)
	CALL NEG(IW(1,3),NUSED)
	CALL MUL(IY,IW(1,3),IW(1,1),NUSED)
	CALL COPY(IW(1,1),IY,NUSED)
	GOTO 10
	END
C
C
C
C
C
C
C
	FUNCTION IEXCMP(IA,IB,N,IW)
C	This routine compares two extended precision numbers.
C	The returned value is 1 if IA > IB, 0 if IA is within roundoff of IB, etc.
C	IW is a work area of size N.
C	Only works for positive numbers
	INTEGER*2 IA(*),IB(*),IW(*)
	INTEGER*2 N,I
	REAL*8 A,B
	REAL DIF
C	First check if the numbers are close.
	CALL CVEXR8(IA,N,A)
	CALL CVEXR8(IB,N,B)
	DIF=(A-B)/(A+B)
	IF(ABS(DIF).GT.1E-12) THEN
	  IEXCMP=1
	  IF(DIF.LT.0) IEXCMP=-1
	  RETURN
	ENDIF
C	The values are close
	CALL NEG(IB,N)
	CALL ADD(IA,IB,IW,N)
	CALL NEG(IB,N)
C	Check for zero
	IF(IA(1)-IW(1).GE.N-2) THEN
C	  Within roundoff of zero
	  IEXCMP=0
	  RETURN
	ENDIF
	DO 10 I=2,N
	  IF(IW(I).NE.0) GOTO 20
10	CONTINUE
	IEXCMP=0
	RETURN
C	Results are not zero
20	IEXCMP=1
	IF(IW(N).GT.5000) IEXCMP=-1
	RETURN
	END

