.TITLE	$EXP		;EXP FORTRAN LIBRARY FUNCTION
	.IDENT	/V1.0/

; COPYRIGHT (C) 1975,1976,1977,
; DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.
;
;	L. WAN
;	MODIFIED FOR USE WITH DECUS BASIC
;
;+
;  FUNCTION:
;	EVALUATE FORTRAN EXP FUNCTION
;
;  METHODS:
;	EXP(X)=1	IF ABS(X) < 2**-60
;	EXP(X)=0	IF X < -88.7
;	EXP(X)=
;	2**Y*2**W*2**Z 	WHERE Y=INTEGER(X*LOG2(E))
;			LET V=16(FRACTION(X*LOG2(E)))
;			WHERE Z=INTEGER(V/16)
;			WHERE 0 < W < 1/16
;			2**W = (P(W**2) + WQ(W**2))/(P(W**2) - WQ(W**2))
;			WHERE P AND Q ARE FIRST DEGREE POLYNOMIALS IN W**2
;			THE COEFFICIENTS OF P AND Q ARE DRAWN FROM HART
;			#1121.  THE RELATIVE ERROR IS  <10**-16.4
;  NOTE:
;	SINGLE PRECISION IS COMPUTED AS DOUBLE PRECISION WITH
;	CONVERSION TO AND FROM SINGLE PRECISION.
;
;  CALL SEQUENCE:
;
;	R5 =>	ARG COUNT
;		ADDRESS OF ARGUMENT
;
;
;-

	F0=%0
	F1=%1
	F2=%2
	F3=%3

	.PSECT	OTS$I

EXP::	SETF
	LDF	@2(R5),F0
	JSR	PC,$$DEXP
	SETF
	STF	F0,-(SP)
	MOV	(SP)+,R0
	MOV	(SP)+,R1
	RTS	PC


$$DEXP:	SETD
	SETI

	TSTD	F0		;TEST THE ARGUMENT
	CFCC
	BGT	POS		;JUMP IF POSITIVE
	CMPD	MIN,F0		;CHECK LOWER BOUND
	CFCC	
	BGT	ZERO		;RETURN ZERO
	BR	SMTST		;
	
POS:	CMPD	MAX,F0		;CHECK UPPER BOUND
	CFCC	
	BLT	OVER		;ERROR OVERFLOW
	
SMTST:	LDD	F0,F2		;MOVE ARG TO F2 FOR LATER USE
	ABSD	F0		;GET THE ABSOLUTE VALUE OF THE ARGUMENT
	CMPD	MAG,F0		;CHECK THE MAGNITUDE
	CFCC
	BGT	ONE

	MOV	#FCONST,R0	;POINTER TO CONSTANTS
	MODD	(R0)+,F2	;F2=FRACT(X*LOG2(E))
	STCDI	F3,R4		;Z=INT(X*LOG2(E))

	TSTD	F2		;FRACT NEGATIVE?
	CFCC	
	BGE	CKONE
	ADDD	#1.0,F2		;MAKE FRACT POSITIVE
	DEC	R4		;AND ADJUST Z=Z-1

CKONE:	CMPD	#1.0,F2		;FRACT = 1.0 NOW?
	CFCC
	BGT	M16
	CLRD	F2		;RESET TO ZERO AND
	INC	R4		;BUMP Z=Z+1

M16:	MODD	#16.0,F2	;F2=FRACT(16*(X*LOG2(E)-FLOAT(Z)))
	STCDI	F3,R3		;D=INT   (16*(...
	DIVD	#16.0,F2	;E=F2/16
	LDD	F2,F3
	MULD	F3,F3		;E*E

	LDD	F3,F1
	ADDD	(R0)+,F1	;A=E*E+Q0
	MULD	(R0)+,F3
	ADDD	(R0)+,F3
	MULD	F2,F3		;B=(E*E*P1 + P0)*E
	LDD	F1,F0
	ADDD	F3,F0		;A+B
	SUBD	F3,F1		;A-B
	DIVD	F1,F0		;(A+B)/(A-B)

SCALE:	ASR	R3		;SHIFT D
	BCC	NOMULT
	MULD	(R0)+,F0	;MULTIPLY BY ROOT OF 2
	BR	SCALE
NOMULT:	BEQ	SCALE1
	ADD	#8.,R0		;POINT TO NEXT ROOT OF 2
	BR	SCALE

SCALE1:
	STEXP	F0,R0		;GET EXPONENT
	ADD	R4,R0		;ADD IN THE MODIFIER
	LDEXP	R0,F0		;RESTORE THE EXPONENT

OUT:	RTS	PC		;EXIT

ONE:	LDD	#1.0,F0		;ANSWER IS ONE
	RTS	PC

ZERO:	.WORD	FUFERR		;TOO SMALL
ZRO:	CLRD	F0		;ANSWER IS ZERO
	RTS	PC		;EXIT
OVER:	.WORD	FOFERR		;TOO LARGE
	LDD	BIG,F0		;SET BIG
	RTS	PC

;	ORDER-DEPENDENT CONSTANTS
;	R0 POINTS AT NEXT CONSTANT IN FPP VERSION
BIG:	.WORD	77777,177777,177777,177777
MIN:	.FLT4	-89.416
MAX:	.FLT4	88.029693
MAG:	.WORD	021400,000000,000000,000000 ;2**-60
FCONST:	.WORD	040270,125073,024534,013761 ;LOG2(E)

	.WORD	041246,101232,074433,171042 ;Q0
	.WORD	037154,113360,153011,153703 ;P1
	.WORD	040746,152405,015345,033343 ;P0
	.WORD	040205,125303,063714,044173 ;2**1/16
	.WORD	040213,112701,161752,105727 ;2**1/8
	.WORD	040230,033760,050615,134251 ;2**1/4
	.WORD	040265,002363,031771,157145 ;2**1/2

	.END