Google
 

Trailing-Edge - PDP-10 Archives - BB-H580E-SB_1985 - expon.mac
There are 7 other files named expon.mac in the archive. Click here to see a list.
; UPD ID= 2970 on 6/25/80 at 11:48 AM by NIXON                          
TITLE	EXPON FOR LIBOL V12C
SUBTTL	EXPONENTIATION		KERRY BENSMAN/ALB/CAM/DMN



	SEARCH	COPYRT
	SALL

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1974, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE.

;EDIT HISTORY
;V12*******************
;NAME	DATE		COMMENTS
;DMN	14-APR-79	      ADD DOUBLE PRECISION FUNCTIONS
;EHM	28-APR-78	[527] FIX ON SIZE ERROR FOR EXPONENTIATION

EXTERNAL OVFLO.

SEARCH	LBLPRM
	HISEG
	.COPYRIGHT		;Put COPYRIGHT statement in .REL file.

	SALL

;ACCUMULATOR DEFINITIONS
	A=	0
	B=	1

	C=	4
	D=	5
	E=	6
	F=	7
	G=	10
	H=	11

	K=	13
	L=	14

	PA=	16
	P=	17

DEFINE DOUBLE (A,B)<EXP A,B>
;THESE ROUTINES CALCULATE A**B
;THE GENERAL CALLING SEQUENCE IS:
;	MOVE	0,A
;	MOVEI	16,B
;	PUSHJ	17,<ROUTINE>

;FOR THE SINGLE PRECISION CASES, THE BASE IS IN ACCUMULATOR 0
;AND THE EXPONENT IS REFERENCED BY  ACCUMULATOR 16.
;FOR DOUBLE PRECISION BASES,THE HIGH ORDER WORD IS IN
;ACCUMULATOR 0 AND THE LOW ORDER WORD IS IN ACCUMULATOR 1.
;FOR DOUBLE PRECISION EXPONENTS, THE HIGH ORDER WORD IS
;REFERENCED BY 0(16) AND THE LOW ORDER WORD IS REFERENCED BY 1(16).
;NOTE ACCS 2 AND 3 CAN NOT BE USED AS SCRATCH SINCE ACC 16 MAY POINT TO THEM.

;ALL ANSWERS ARE RETURNED IN FLOATING POINT FORMAT.
SUBTTL	COMP-1 RAISED TO 1-WORD COMP POWER

;THIS ROUTINE CALCULATES A COMP-1 FLOATING POINT NUMBER
;TO A FIXED POINT (1-WORD) POWER.


;THE CALCULATION IS A**B, WHERE B IS OF THE FORM

;	B=Q(0) + Q(1)*2 + Q(2)*4 + ...WHERE Q(I)=0 OR 1

;THERE ARE NO RESTRICTIONS ON THE BASE OR EXPONENT


	ENTRY	E.C3C1

E.C3C1:	MOVE	C,0(PA)		;MOVE EXP TO SAFE PLACE

EC3C1:	JUMPE	A,ATEST		;SPECIAL CASE IF BASE = 0.

	MOVSI	D,(1.0)		;GET 1.0 IN AC D.
	JUMPGE	C, FEXP2	;IS EXPONENT POSITIVE?
	MOVMS	C		;NO, MAKE IT POSITIVE
	PUSHJ	P, FEXP2	;CALL MAIN PART OF PROGRAM.
	MOVSI	C, (1.0)	;GET 1.0 IN "B2.
	FDVM	C, A		;FORM 1/(A**B) FOR NEG. EXPONENT.
	POPJ	P,		;RETURN.

FEXP1:	FMP	A, A		;FORM A**N, FLOATING POINT.
	JFCL	1,AINF		;IF OVER/UNDERFLOW, GO TO AINF.
	LSH	C, -1		;SHIFT EXPONENT FOR NEXT BIT.
FEXP2:	JFCL	1,.+1		;CLEAR OVERFLOW FLAG
	TRZE	C, 1		;IS THE BIT ON?
	FMP	D, A		;YES, MULTIPLY ANSWER BY A**N.
	JFCL	1,AINF		;IF OVER/UNDERFLOW, GO TO AINF.
	JUMPN	C, FEXP1	;UPDATE A**N UNLESS ALL THROUGH.
	MOVE	A, D		;PICK UP RESULT FROM D.
	POPJ P,

ATEST:	JUMPGE	C,AZERO		;EXIT IF BASE =0,EXP>=0.
AINF:	HRLOI	A,377777	;ANS.=+INFINITY
	SETOM	OVFLO.		;[527] SET OVERFLOW ON
	POPJ	17,		;AND EXIT.

AZERO:	SETZ	A,		;ANSWER IS ZERO
	POPJ P,
SUBTTL	COMP-1 RAISED TO COMP-1 POWER

;THIS ROUTINE CALCULATES A COMP-1 FLOATING POINT NUMBER
;RAISED TO A COMP-1 FLOATING POINT POWER.

;THE CALCULATION IS


;	A**B= EXP(B*LOG(A))

;IF THE EXPONENT IS AN INTEGER < 2**35 IN MAGNITUDE, THE
;RESULT WILL BE COMPUTED USING  "E.C3C1"  AND THE ANSWER 
;WILL HAVE THE CORRECT SIGN. (REMEMBER THAT THE "INTEGER"
;HAS ONLY 27 SIGNIFCANT BITS.)


;SINCE NEGATIVE NUMBERS RAISED TO NON-INTEGER POWERS YIELD
;COMPLEX ANSWERS, THE MAIN ALGORITHM CALCULATES
;	EXP(B*LOG(ABSF(A)))


	ENTRY	E.C3C3

E.C3C3:	MOVE	C,0(PA)		;GET EXPONENT

EC3C3A:	JUMPE	A,ATEST		;SPECIAL IF BASE = 0

	MOVM	E,C		;SET EXP. POSITIVE.
	MOVEI	D,0		;CLEAR AC D TO ZERO
	LSHC	D,11		;SHIFT 9 PLACES LEFT
	SUBI	D,200		;TO OBTAIN SHIFTING FACTOR
	JUMPLE	D,EXP3GO	;IS D > 0
	HRR	F,D		;SET UP F AS AN INDEX REG.
	SETZ	D,		;CLEAR OUT AC D
	LSH	E,-1		;RIGHT ADJUST EXP TO BIT 1.
	ASHC	D,(F)		;SHIFT LFT BY CONTENTS OF F
	JFCL	1,EXP3GO	;IF OVERFLOW, GO TO EXP3GO.
	JUMPN	E,EXP3GO	;IS EXPONENT AN INTEGER ?
	SKIPGE	C		;YES, WAS IT NEGATIVE?
	MOVNS	D		;YES
	MOVE	C,D		;NO
	JRST	EC3C1		;USE INTEGER ROUTINE, THEN RETURN

EXP3GO:	MOVM	E,A		;GET ABS(BASE) IF NE 0 OR 1.
	MOVE	D,C		;SAVE AC C.
	PUSHJ	P,ALOG.		;CALCULATE LOG OF "B"
	FMPRM	E,D		;CALCULATE A*LOG(B)
	JRST	EXPON.		;CALCULATE EXP(A*LOG(B)) AND RETURN
SUBTTL ALOG.

;FLOATING POINT SINGLE PRECISION LOGARITHM FUNCTION
;LOG(ABSF(X)) IS CALCULATED BY THE SUBROUTINE, AND AN
;ARGUMENT OF ZERO IS RETURNED AS MINUS INFINITY.
;ALOG IS THE ENTRY POINT FOR LOGE(X).

;FOR LOGE(X), THE ALGORITHM IS:
;	LOGE(X) = (I + LOG2(F))*LOGE(2)
;	WHERE X = (F/2)*2^(I+1), AND LOG2(F) IS GIVEN BY
;	LOG2(F) = C1*Z + C3*Z^3 + C5*Z^5 - 1/2
;	AND Z = (F-SQRT(2))/(F+SQRT(2))

;THE CALLING SEQUENCE FOR THE ROUTINE IS AS FOLLOWS:
;	MOVE	E,ARGUMENT
;	PUSHJ	P,ALOG.
;THE ANSWER IS RETURNED IN ACCUMULATOR E
ALOG.:	JUMPE	E, LZERO	;CHECK FOR ZERO ARGUMENT
	CAMN	E, ONE		;CHECK FOR 1.0 ARGUMENT
	JRST	ZERANS		;IT IS 1.0 RETURN ZERO ANS.
	ASHC	E, -33		;SEPARATE FRACTION FROM EXPONENT
	ADDI	E, 211000	;FLOAT THE EXPONENT AND MULT. BY 2
	MOVSM	E, G		;NUMBER NOW IN CORRECT FL. FORMAT
	MOVSI	E, 567377	;SET UP -401.0 IN E
	FADM	E, G		;SUBTRACT 401 FROM EXP.*2
	ASH	F, -8		;SHIFT FRACTION FOR FLOATING
	TLC	F, 200000	;FLOAT THE FRACTION PART
	FAD	F, L1		;F = F-SQRT(2.0)/2.0
	MOVE	E, F		;PUT RESULTS IN E
	FAD	E, L2		;E = E+SQRT(2.0)
	FDV	F, E		;F = F/E
	MOVE	H, F		;STORE NEW VARIABLE IN H
	FMP	F, F		;CALCULATE Z^2
	MOVE	E, L3		;PICK UP FIRST CONSTANT
	FMP	E, F		;MULTIPLY BY Z^2
	FAD	E, L4		;ADD IN NEXT CONSTANT
	FMP	E, F		;MULTIPLY BY Z^2
	FAD	E, L5		;ADD IN NEXT CONSTANT
	FMP	E, H		;MULTIPLY BY Z
	FAD	E, G		;ADD IN EXPONENT TO FORM LOG2(X)
	FMP	E, L7		;MULTIPLY TO FORM LOGE(X)
	POPJ	P,

ZERANS:	TDZA	E, E		;MAKE ANSWER ZERO
LZERO:	MOVE	E,MIFI		;PICK UP MINUS INFINITY
	POPJ	P,


;CONSTANTS

L1:	577225754146		;-0.707106781187
L2:	201552023632		;1.414213562374
L3:	200462532521		;0.5989786496
L4:	200754213604		;0.9614706323
L5:	202561251002		;2.8853912903
L7:	200542710300		;0.69314718056
MIFI:	400000000001		;LARGEST NEGATIVE FLOATING NUMBER
SUBTTL SINGLE PRECISION EXPONENTIAL FUNCTION

;FLOATING POINT SINGLE PRECISION EXPONENTIAL FUNCTION
;CALCULATE EXP(C).

;IF X<=-89.415..., THE PROGRAM RETURNS ZERO AS THE ANSWER
;IF X>=88.029..., THE PROGRAM RETURNS 377777777777 AS THE ANSWER
;THE RANGE OF THE ARGUMENT IS REDUCED AS FOLLOWS:
;EXP(X) = 2**(X*LOG(E)BASE2) = 2**(M+F)
;WHERE M IS AN INTEGER AND F IS B FRACTION
;2**M IS CALCULATED BY ALGEBRAICALLY ADDING M TO THE EXPONENT
;OF THE RESULT OF 2**F. 2**F IS CALCULATED AS

;2**F = 2(0.5+F(B+A*F^2 - F-C(F^2 +D)**-1)**-1)

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
;	MOVE	PA,ARGUEMENT
;	PUSHJ	P,EXPON.
;THE ANSWER IS RETURNED IN ACCUMULATOR 0


EXPON.:	CAMGE	D,E77		;IS EXP. < -89.41...?
	JRST	AZERO		;YES, GO TO EXIT.
	CAMLE	D,E7		;IS EXP. > +88.029..?
	JRST	AINF		;YES, GET LARGEST FLOATING NUMBER
				;FALL INTO STANDARD ALGORITHM

EXP1:	SETZ	G,		;INITIALIZE G
	MULI	D, 400		;SEPARATE FRACTION AND EXPONENT
	TSC	D, D		;GET "B" POSITIVE EXPONENT
	MUL	E, E5		;FIXED POINT MULTIPLY BY LOG2(E)
	ASHC	E, -242(D)	;SEPARATE FRACTION AND INTEGER	
	AOSG	E		;ALGORITHM CALLS FOR MULT. BY 2
	ADDI	E,1		;ADJUST IF FRACTION WAS NEGATIVE
	HRRZM	E,H		;SAVE FOR FUTURE SCALING
	JUMPG	F,EXP2		;GO AHEAD IF ARG > 0.
	TRNE	F,377		;ARE ALL THESE BITS 0?
	ADDI	F,200		;NO. FIX UP.

EXP2:	ASH	F, -8		;MAKE ROOM FOR EXPONENT
	TLC	F, 200000	;PUT 200 IN EXPONENT BITS
	FADB	F, G		;NORMALIZE, RESULTS TO F AND G
	FMP 	F, F		;FORM X^2
	MOVE	A, E2		;GET FIRST CONSTANT
	FMP	A, F		;E2*X^2 IN A
	FAD	F, E4		;ADD E4 TO RESULTS IN F
	MOVE	D, E3		;PICK UP E3
	FDV	D, F		;CALCULATE E3/(F^2 + E4)
	FSB	A, D		;E2*F^2-E3(F^2 + E4)**-1
	MOVE	E, G		;GET G AGAIN
	FSB	A, E		;SUBTRACT FROM PARTIAL SUM
	FAD	A, E1		;ADD IN E1
	FDVM	E, A		;DIVIDE BY F
	FAD	A, E6		;ADD O.5
	FSC	A, (H)		;SCALE THE RESULTS
	POPJ	P,

E1:	204476430062		;9.95459578
E2:	174433723400		;0.03465735903
E3:	212464770715		;617.97226953
E4:	207535527022		;87.417497202
E5:	270524354513		;LOG(E),  BASE 2
E6:	0.5
E7:	207540074636		;88.029...
E77:	570232254037		;-89.415986

IFE BIS,<
ONE:	1.0

	ENTRY	E.F2D1,E.F2D2,E.F2FP,E.F2F2
E.F2D1:E.F2D2:E.F2FP:E.F2F2:
	OUTSTR	[ASCIZ	/?EXPON is not assembled for double precision floating point.
/]
	JRST	STOPR.##

	END
>
SUBTTL	COMP-2 RAISED TO COMP POWER

;THIS ROUTINE CALCULATES A DOUBLE PRECISION NUMBER RAISED
;TO A FIXED POINT POWER. THE CALCULATION IS A**N, WHERE N
;IS AN INTEGER OF THE FORM
;	N= Q(0) + Q(1)*2 + Q(2)*4 + ... WHERE Q(I) = 0 OR 1
;THE ONLY RESTRICTION ON THE BASE OR EXPONENT IS THAT
;AN EXPONENT OF 400000000000 IS NOT HANDLED CORRECTLY.

	ENTRY	E.F2D1		;COMP-2 TO 1-WORD COMP
	ENTRY	E.F2D2		;COMP-2 TO 2-WORD COMP
	

E.F2D2:	SKIPN	G,(PA)		;GET HIGH WORD
	JRST	[SKIPN	G,1(PA)		;0, GET LOW WORD
		JRST	AB.ONE		;BOTH ZERO
		JRST	E.F2G]		;NO, CONTINUE WITH G SETUP
	AOJE	G,[SKIPGE G,1(PA)	;-1, GET LOW WORD
		JRST	E.F2G		;SMALL NEGATIVE NUMBER
		JRST	.+1]		;DON'T KNOW WHAT TO USE
	SETZB	A,B		;ERROR, RETURN 0 OR +INFINITY

E.F2D1:	SKIPN	G,(PA)		;IS EXPONENT 0?
	JRST	AB.ONE		;YES, A**0 GIVES 1
	JUMPE	A,ABTEST	;IF BASE = 0

E.F2G:	JUMPGE	G,[DMOVE D,A		;IS EXPONENT NEGATIVE?
		JRST	DEX2]		;NO, PUT ARG IN D,D+1 AND START MAIN LOOP
	MOVMS	G		;GET POSITIVE VALUE
	DMOVE	D,ONE		;GET DOUB. PRECISION 1.0
	DFDV	D,A		;CALCULATE (1/X)**N, SINCE N .L. 0
	JOV	OVR1		;OVER/UNDERFLOW OCCURED

DEX2:	DMOVE	A,ONE		;GET DOUB. PREC. 1.0
	JRST	DEX4		;START CALCULATING POWERS OF X (OR 1/X)

DEX3:	DFMP	D,D		;SQUARE X (OR 1/X) AGAIN
	JOV	OVR		;TEST FOR OVER/UNDERFLOW
	LSH	G,-1		;LOOK AT NEXT BIT IN EXPONENT
DEX4:	TRZN	G,1		;IS LO BIT ON? (I.E. NEXT POWER OF 2 WANTED)
	JRST	DEX5		;NO, DON'T MULTIPLY INTO ANSWER
	DFMP	A,D		;MULTIPLY POWER OF X INTO ANSWER
	JOV	OVR		;TEST FOR OVER/UNDERFLOW

DEX5:	JUMPN	G,DEX3		;IF G .N. 0, GET MORE POWERS OF X (OR 1/X)
DEX6:	POPJ	P,

OVR:	DMOVE	D,A		;PUT CURRENT RESULT IN SAFE PLACE
OVR1:	JSP	E,.+1		;GET THE FLAGS
				;NOTE IN EXTENDED SECTIONS USE SFM (JRST 14,)
	TLNE	E,(1B11)	;TEST FOR UNDERFLOW
	JRST	UNDERF		;YES, SET ANS TO ZERO
	PUSHJ	PP,ABTEST	;OVERFLOW, RETURN LARGE NUMBER
	JUMPGE	D,DEX6		;IF THE ARG IS NEGATIVE
	TRNN	G,1		;AND THE EXPONENT IS ODD
	DMOVN	A,A		;THEN RETURN -INFINITY
	POPJ	P,

AB.ONE:	DMOVE	A,ONE		;A**0 GIVES 1
	POPJ	P,

ABTEST:	SKIPL	(PA)		;IS BASE 0 WITH POSITIVE EXPONENT?
	JRST	ABZERO		;YES, RETURN 0
AB.INF:	HRLOI	A,377777	;BASE IS 0 WITH NEG. EXP- OVERFLOW
	HRLOI	B,377777	;RETURN LARGEST POSITIVE NUMBER
	SETOM	OVFLO.		;SET OVERFLOW ON
	POPJ P,

UNDERF:	SETOM	OVFLO.		;UNDERFLOW
ABZERO:	SETZB	A,B		;OUTPUT ANS = 0
	POPJ	P,		;RETURN
SUBTTL	COMP-2 RAISED TO FLOATING POINT POWER

;THIS PROGRAM CALCULATES A**B, WHERE A AND B ARE DOUBLE
;PRECISION FLOATING POINT NUMBERS. THE ALGORITHM USED IS
;A**B = DEXP(B*DLOG(/A/)).  THE ABSOLUTE VALUE OF A IS USED
;IN THIS CALCULATION BECAUSE A NEGATIVE NUMBER TO A
;NON-INTEGER POWER PRODUCES A COMPLEX ANSWER, AND B IS
;PRESUMED TO BE NON-INTEGER.
;CHECKS TO SEE IF B IS AN INTEGER AND TREATS IT
;ACCORDINGLY BY CALLING E.F2D1.
;GIVES ERROR FOR NEGATIVE NUMBER TO NON-INTEGER POWER.

	ENTRY	E.F2FP		;COMP-2 RAISED TO COMP-1 POWER
	ENTRY	E.F2F2		;COMP-2 RAISED TO COMP-2 POWER

E.F2FP:	MOVE	F,(PA)		;GET EXPONENT
	TDZA	G,G		;ZERO LOW WORD

E.F2F2:	DMOVE	F,(PA)		;GET EXPONENT
	JUMPE	F,AB.ONE	;IF EXPONENT = 0 RETURN ANSWER OF 1.0
	JUMPE	A,ABTEST	;IF BASE = 0
	MOVM	D,F		;MAGNITUDE OF WORD 1 OF EXP
	SETZ	C,		;CLEAR AC C
	LSHC	C,11		;SHIFT 9 BITS
	CAILE	C,233		;INTEGER MUST FIT IN 27 BITS
	JRST	DEXPD		;NOT AN INTEGER
	SUBI	C,200		;200 COMPLEMENT ON EXP
	MOVM	E,G		;WORD 2 OF EXP
	JUMPN	E,DEXPD		;IS WORD 2 ZERO?
	LSH	D,-1		;COULD STILL BE INTEGER
	HRRZ	E,C		;SET E AS INDEX REG
	SETZ	C,		;CLEAR C
	ASHC	C,(E)		;SHIFT LEFT BY CONTENT S OF E
	JOV	DEXPD		;OVERFLOW
	JUMPN	D,DEXPD		;IS IT INTEGER??
	SKIPL	(PA)		;YES-NEED SIGN OF EX
	SKIPA	K,C		;GET EXPONENT IN SAFE PLACE
	MOVN	K,C		;NEGATIVE-NEGATE INTEGER
	MOVEI	PA,K		;POINT TO EXPONENT
	JRST	E.F2D1		;CALL INTEGER ROUTINE

DEXPD:	JUMPGE	A,DEXPD2	;NEGATIVE??
	DMOVN	A,A		;YES-GET ABS VALUE
	OUTSTR	[ASCIZ /%Attempt to raise negative double precision number to non-integer power
/]
DEXPD2:	PUSHJ	P,DLOG.		;CALC. LOG(/A/).
	DFMP	A,(PA)		;CALC. B*LOG(/A/)
	JOV	OVUNFL
	JRST	DEXP.		;CALC. EXP[B*LOG(/A/)] AND
				;LEAVE IT IN AC'S 0 AND 1.

OVUNFL:	JUMPE	A,AB.ONE	;IF EXP= 0, ANS = 1.0.
	JUMPL	A,UNDERF	;GO TO UNDERFLOW.
	JRST	AB.INF		;OVERFLOW, OUTPUT ANS = INFINITY
SUBTTL	DLOG. DOUBLE PRECISION LOGARITHM FUNCTION

;THIS PROGRAM CALCULATES THE LOGARITHM OF A DOUBLE PRECISION
;ARGUMENT. THE ALGORITHM USED IS DESCRIBED ON PAGES 29-30 OF
;RALSTON AND WILF, "MATHEMATICAL METHODS FOR DIGITAL COMPUTERS".
;THE ARGUMENT X IS WRITTEN AS

;	X = (2**N)*F	WHERE 1/2 < F < 1
;THEN LOG(X) = (N*LOGE(2)) + LOG(F)
;F IS REDUCED BY FIXED POINT MULTIPLICATION BY NOT MORE THAN
;THREE CONSTANTS. THIS YIELDS

;	0 < T = A1*A2*A3*F - 1.0 < (2**-7)/5
;NOTE THAT NOT ALL THE A1,A2,A3 NEED BE INCLUDED IN THE PRODUCT.
;FINALLY, 
;	LOG(F) = LOG(1+T) - LOG(A1) - LOG(A2) - LOG(A3)
;LOG(1+T) IS CALCULATED AS A TAYLOR SERIES IN T.

;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
;	DMOVE	0,BASE
;	PUSHJ	17,DLOG.

;THE RESULT IS LEFT IN ACCUMULATOR 0 AND 1.
;AVAILABLE TABLES FOR LOGARITHMS LIST 16 SIGNIFICANT DIGITS,
;ALL OF WHICH AGREE WITH VALUES PRODUCED BY THIS PROGRAM. AN
;ESTIMATE OF THE MAXIMUM THEORETICAL ERROR IS OBTAINED BY OB-
;SERVING THAT, AFTER REDUCTION, THE ARGUMENT IS LESS THAN
;(2**-7)/5 = .00625. HENCE, THE ERROR IN THE TAYLOR SERIES
;IS LESS THAN
;	(.00625**9)/9 = 1.4573*10**-20/9 = .1619*10**-20

;CONSTANTS AND TEMPORARY LOCATIONS AND STUFF

LOGLST:	DOUBLE 200471174064,325425031470	;0.61180 15411 05992 8976
	DOUBLE 200402252251,151350376610	;0.50455 60107 52395 2859
	DOUBLE 177637144373,057714113734	;0.40546 51081 08164 3810
	DOUBLE 177506061360,207057302360	;0.31845 37311 18534 6147
	DOUBLE 176710776761,346515041520	;0.22314 35513 14209 7553
	DOUBLE 176537746034,051711723600	;0.17185 02569 26659 2214
	DOUBLE 175557032242,271265512760	;0.08961 21586 89687 12374
	DOUBLE 173770123303,236031377700	;0.03077 16586 66753 68689

LIST:	354000000000		;1.11011 BINARY
	324000000000		;1.10101 BINARY
	300000000000		;1.10000 BINARY
	260000000000		;1.01100 BINARY
	240000000000		;1.01000 BINARY
	230000000000		;1.00110 BINARY
	214000000000		;1.00011 BINARY
	204000000000		;1.00001 BINARY

DLOGE2:	DOUBLE 200542710277,276434757153	;0.69314 71805 59945 30941 72321

A9:	DOUBLE 175707070707,034343434344	;1/9
A8:	DOUBLE 601400000000,0			;-1/8
A7:	DOUBLE 176444444444,222222222222	;1/7
A6:	DOUBLE 601252525252,652525252526	;-1/6
A5:	DOUBLE 176631463146,146314631463	;1/5
A4:	DOUBLE 600400000000,0			;-1/4
A3:	DOUBLE 177525252525,125252525253	;1/3
A2:	DOUBLE 577400000000,0			;-1/2
ONE:
A1:	1.0
DZERO:	DOUBLE 0,0
DLOG.:	JUMPG	A,DLOGP		;ARGUMENT IS GREATER THAN 0
	JUMPE	A,LGZERO	;CHECK FOR ZERO ARGUMENT
	OUTSTR	[ASCIZ /%Attempt to take DLOG of Negative Arg
/]
	DMOVN	A,A		;GET ABSF(X)
DLOGP:	CAMN	A,ONE		;X PRECISELY 1.0?
	JUMPE	B,ABZERO	;YES, RETURN ZERO
DLOG1:	LDB	E,[POINT 8,A,8]	;NO,PICK UP EXPONENT FROM HI ORDER
	SUBI	E,200		;GET EXPONENT EXCESS 200
	FSC	E,233		;MAKE FLOATING POINT NUMBER
	SETZ	F,		;SET UP LO HALF
	DFMP	E,DLOGE2	;CALCULATE N*LOGE(2)
	DMOVE	K,E		;SAVE AWAY IN SAFE PLACE

	TLZ	A,777000	;MASK OUT EXPONENT
	ASHC	A,8		;NORMALIZE FRACTION TO BIT 1
	SETZB	E,F		;INITIALIZE REDUCTION CONSTANT TO 0
DLOG2:	LDB	H,[POINT 3,A,4]	;GET HI 3 BITS BELOW 1/2
	MOVE	C,B		;GET COPY OF LOW HALF
	MUL	C,LIST(H)	;FIXED POINT MULTIPLY FOR LO HALF
	MUL	A, LIST(H)	;MULTIPLY HI ORDER, TOO
	TLO	B,(1B0)		;SET BIT 0, TO AVOID OVERFLOW
	ADD	B,C		;COMBINE RESULTS OF MULTIPLY
	TLZN	B,(1B0)		;CLEAR BIT 0, WAS THERE OVERFLOW?
	ADDI	A,1		;YES, PROPOGATE CARRY
	ASH	H,1		;TURN BITS INTO D.P. POINTER
	DFAD	E,LOGLST(H)
	TLZE	A,200000	;IS THE PRODUCT .GE. 1.0?
	JRST	DLOG3		;YES
	ASHC	A,1		;NO, GET RID OF EXTRANEOUS ZERO
	JRST	DLOG2		;TRY ANOTHER MULTIPLICATION

DLOG3:	ASHC	A,-7		;MAKE ROOM FOR THE EXPONENT
	TLC	A,200000	;INSERT EXPONENT
	DFAD	A,DZERO		;NORMALIZE
	DMOVN	E,E		;NEGATE LOG OF MAGIC NUMBERS
	DFAD	E,K		;AND ADD INTO FINAL SUMMATION
	DMOVEM	E,K
	DMOVE	E,A9		;PICK UP CONSTANT TO START
	MOVEI	H,A8		;INIT TABLE POINTER AT A8
DLOG4:	DFMP	E,A		;MULTIPLY ACCUMULATED SUM BY X
	DFAD	E,0(H)		;ADD NEXT CONSTANT INTO PARTIAL SUM
	ADDI	H,2		;UPDATE POINTER TO NEXT CONSTANT
	CAIG	H,A1		;ARE WE DONE YET?
	JRST	DLOG4		;NO, LOOP BACK FOR MORE TAYLOR SERIES
	DFMP	A,E		;YES, ONE LAST MULTIPLICATION
	DFAD	A,K		;AND ADD SERIES SUM INTO FINAL ANSWER
	POPJ	P,		;RETURN

LGZERO:	MOVSI	A,400000	;SET UP LARGE NEG. NUM.
	DFAD	A,A		;CAUSE NEGATIVE OVERFLOW
				;TRAP ROUTINE FIXED IT UP AND
				;PRINTED MESSAGE
	POPJ	P,
SUBTTL	DEXP. DOUBLE PRECISION EXPONENTIAL FUNCTION

;ARGUMENT. AN ARGUMENT OF ZERO CAUSES AN IMMEDIATE EXIT WITH AN
;ANSWER OF 1.0 . AN ARGUMENT WHOSE MAGNITUDE EXCEEDS 88.028
;CAUSES THE ROUTINE TO EXIT WITH 0 IF THE ARGUMENT WAS NEGATIVE.
;AND 377777777777 IF THE ARGUMENT WAS POSITIVE. THIS IS
;BECAUSE DIRECT CALCULATION OF EXP(X) FOR ABSF(X)>88.028 WOULD
;CAUSE EXPONENT OVERFLOW OR UNDERFLOW.

;THE ROUTINE USES THE FOLLOWING ALGORITHM:
;EXP(X)	= 2**(X*LOG2(E))
;	= 2**(M+F) WHERE M IS AN INTEGER AND 0<F<1
;	= 2**(M+N+R) WHERE 0<R<1/8 AND M+N+R=X*LOG2(E)

;	= 2**(M+N) * EXP(R*LOG(2))
;2**M IS CALCULATED EASILY WITH THE FLOATING SCALE INSTRUCTION.
;2**N IS CALCULATED BY DETERMINING THE CORRECT INTERVAL OF N AND
;USING A TABLE OF POWERS OF TWO FROM 2**1/8 TO 2**7/8.

;FINALLY, EXP(R*LOG(2)) IS CALCULATED BY A CONTINUED FRACTION
;TAKEN FROM RALSTON AND WILF, "METHODS FOR DIGITAL COMPUTERS" :
;EXP(R*LOG(2)) = 1+AA4/((B4/R) -C4 + D4*R + H4/(R + B4/R))

;THE FOLLOWING ERRORS HAVE BEEN OBSERVED WITH DEXP:
;	1. WITH -10.0<X<10.0, ERRORS RANGED FROM 0 TO 48
;	   UNITS IN THE 19TH SIGNIFICANT DIGIT. THE MEAN ERROR
;	   FOR 20 READINGS WAS 5.4 UNITS IN THE 19TH DIGIT.
;	2. WITH 10.0<X<40.0, ERRORS RANGED FROM 0 TO 41 UNITS
;	   IN THE 19TH SIGNIFICANT DIGIT. THE MEAN ERROR FOR
;	   7 READINGS WAS 16 UNITS IN THE 19TH SIGNIFICANT DIGIT.
;	3. WITH 40.0<X<88.0, ERRORS RANGED FROM 5 TO 57 UNITS IN
;	   THE 19TH SIGNIFICANT DIGIT. THE MEAN ERROR FOR 12
;	   READINGS WAS 44.4 UNITS IN THE 19TH SIGNIFICANT DIGIT.
;THE ERRORS REFERRED TO ABOVE ARE ABSOLUTE ERRORS. IT SHOULD
;BE NOTED THAT ADDITIONAL ERRORS ARE INTRODUCED BY ERRORS IN
;THE DOUBLE PRECISION INPUT AND OUTPUT ROUTINES.

;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
;	DMOVE	0,BASE
;	PUSHJ	17,DEXP.

;CONSTANTS AND TEMPORARY LOCATIONS AND STUFF

DCON1:	88.028
DLOG2E:	DOUBLE	201561250731,112701376057	;LOG2(E) = 1.44269 50408 88963 40740
TABLE:	0			;0
	040000000000		;1/8
	100000000000		;2/8
	140000000000		;3/8
	200000000000		;4/8
	240000000000		;5/8
	300000000000		;6/8
	340000000000		;7/8

POWERS:	DOUBLE 1.0,0				;2**0 = 1.0
	DOUBLE 201427127017,037250572672	;2**1/8 = 1.09050 77326 65257 65919
	DOUBLE 201460337602,214333425134	;2**2/8 = 1.18920 71150 02721 06671
	DOUBLE 201513773265,115425047073	;2**3/8 = 1.29683 95546 51009 66590
	DOUBLE 201552023631,237635714441	;2**4/8 = 1.41421 35623 73095 04878
	DOUBLE 201612634520,212520333270	;2**5/8 = 1.54221 08254 07940 824
	DOUBLE 201656423746,126551655275	;2**6/8 = 1.68179 28305 07429 086
	DOUBLE 201725403067,076722207113	;2**7/8 = 1.83400 80864 09342 463

AA4:	DOUBLE 206744575555,062215755376	;AA4= 60.59319 17173 36463 11080
B4:	DOUBLE 207535527021,213670572221	;B4 = 87.41749 72022 35527 474
MC4:	DOUBLE 572033202222,715562022402	;MC4=-C4 = -30.29659 58586 68231 555
D4:	DOUBLE 201414631463,063146314632	;D4 = 1.05
H4:	DOUBLE 210654261010,261565402456	;H4 = 214.17286 81454 77042 3113
DEXP.:	JUMPE	A,AB.ONE	;RETURN 1.0 FOR ARG OF ZERO
	MOVM	C,A		;GET POS VALUE OF EXPONENT
	CAML	C,DCON1		;TOO BIG TO COMPUTE?
	JRST	[JUMPGE	A,AB.INF	;YES, RETURN INFINITY
		JRST	ABZERO]		;OR ZERO IF NEGATIVE
	DFMP	A,DLOG2E
	HLRE	E,A		;EXTRACT EXPONENT
	ASH	E,-9		;...
	TSC	E,E		;TAKE 1'S COMPLEMENT IF NUM .L. 0
	SKIPGE	A		;CHANGE  EXPONENT BITS TO SIGN BITS
	TLOA	A,377000	;NUMBER NEGATIVE, SET BITS
	TLZ	A,377000	;NUMBER POSITIVE, CLEAR BITS
	ASHC	A,8		;LEFT JUSTIFY ARG FRACTION BITS
	DMOVE	C,A		;GET ANOTHER COPY OF FRACTION
	ASHC	A,-243(E)	;SIMULATE THREE-WORD SHIFT WITH...
				;TWO LONG SHIFTS. THIS LEAVES INTEGER
				;IN A, FRAC IN B AND C.
	LSH	D,1		;SQUEEZE OUT SIGN BIT
IFN C-B-1,<
	EXCH	B,C-1		;CROCK BECAUSE B AND C ARE NOT ADJACENT
>
	LSHC	C,43-200(E)	;THEN DO 2ND LONG SHIFT, (THE LSHC HERE
				;PREVENTS OVERFLOW GOING LEFT)
IFN C-B-1,<
	EXCH	B,C-1		;PUT ACCS BACK THE WAY WE WANT
>
	TLZ	B,(1B0)		;CLEAR SIGN BIT. IF FRAC WAS <0,
				;THIS GIVES 1-FRAC AND PROPER EXP.
	HRRE	H,A		;SAVE EXPONENT FOR FUTURE SCALING
	MOVEI	G,7		;GET INDEX REGISTER POINTER TO TABLE
REDUCE:	CAME	B,TABLE(G)	;DOES ARGUMENT MATCH TABLE ENTRY?
	JRST	REDUC2		; NO - KEEP LOOKING
	JFFO	C,.+2		; HOW MANY LEADING 0'S
	  JRST	REDUC1		; NONE - C CONTAINS 0
	CAIL	D,^D28		; ARE ANY LEFT OF BIT 28 ON?
REDUC1:	JRST	[LSH	G,1		; NO, IFF LO HALF=0!<377(OCT). CHANGE INDEX TO POINTER
		DMOVE	A,POWERS(G)	;PICK UP DOUBLE WORD ANSWER
		JRST	SCALE]		;SCALE RESULTS AND EXIT
REDUC2:	CAMGE	B,TABLE(G)	;IS ARGUMENT GREATER THAN ENTRY?
	SOJA	G,REDUCE	;NO, TRY NEXT LOWER ENTRY
	SUB	B,TABLE(G)	;YES, ALL DONE -REDUCE ARGUMENT
	LSH	G,1		;SAVE INDEX AS A POINTER
IFN C-B-1,<
	EXCH	C,B+1
>
	ASHC	B,-8		;MAKE ROOM FOR EXPONENT
	TLO	B,200000	;INSERT EXPONENT
	DFAD	B,DZERO		;NORMALIZE RESULT
	DMOVE	A,B		;PUT RESULT AN A,B
IFN C-B-1,<
	EXCH	C,B+1
>
	DMOVE	C,B4		;GET B4/F
	DFDV	C,A
	DMOVEM	C,K		;SAVE B4/F
	DFAD	C,A		;GET F+B4/F
	DMOVEM	C,E		;GET H4/(F+B4/F)
	DMOVE	C,H4
	DFDV	C,E
	DFMP	A,D4		;GET D4*F
	DFAD	C,A		;GET (B4/F)-C4+D4*F+(H4/(F+B4/F))
	DFAD	C,MC4
	DFAD	C,K
	DMOVE	A,AA4		;GET 1.0+AA4/REST
	DFDV	A,C
	DFAD	A,ONE
	JUMPE	G,SCALE		;MULTIPLY BY POWER OF TWO?
	DFMP	A,POWERS(G)

SCALE:	FSC	A,(H)		;SCALE RESULTS
	POPJ	P,		;RETURN


	END