A large commit.

[pdp8.git] / sw / os8 / v3d / sources / fortran / all / dexp.ra

/
/
/    SUBROUTINE DEXP
/
/ VERSION 5A 4-26-77 MH
/
/E^X=2^(X*LOG2(E))
/E^X=2^(M+F)
/M=INTEGER; F=FRACTION
/
/2^(M+F)=2^(M+N+R)
/WHERE 0<R<1/8
/AND M+N+R=M+F=X*LOG2(E)
/
/(2^M)*(2^N)*(2^R)=E^X
/
/2^M IS CALCULATED BY SUCCESSIVE MULTIPLIES
/2^N IS CALCULATED BY LOOK UP
/2^R=1+<A4/((B4/R)-C4+(D4*R)+(H4/(R+(B4/R))))>
/
/RESTRICTIONS:
/X=0 IMPLIES E^X=1
/
/X>88.028 IMPLIES E^X=3377/3377/3777/7777/777/7777
/
/X<-88.028 IMPLIES E^X=0
/
/
/
	SECT	DEXP
	JA	#DEXP
	DPCHK
	TEXT	+DEXP  +
/
DEXPXR,	SETX	XRDEXP
	SETB	BPDEXP
/
/BEGINNING OF BASE PAGE
/
BPDEXP,	F 0.0
XRDEXP,	F 0.0
X,	F 0.0
	F 0.0
/
	ORG	10*3+BPDEXP
	FNOP
	JA	DEXPXR
	0
DEXRTN,	JA	.
/
TOPLIM,	3377
	3377
	3777
	7777
	7777
	7777
M,	F 0.0
	F 0.0
N,	F 0.0
	F 0.0
R,	F 0.0
	F 0.0
LOG2E,	0001		/1.4426950408889634
	2705
	2435
	4512
	7013
	7603
DFP125,	7775		/.125
	3777
	7777
	7777
	7777
	7776
DEXFP1,	F 1.0
	F 0.0
/
DFR1S8,	0001		/2^1/8
	2134
	5340
	7437
	2505
	7302
DFP2S8,	0001		/2^2/8
	2301
	5770
	1214
	3334
	2524
DFP3S8,	0001		/2^3/8
	2457
	7553
	2515
	4250
	4720
DFP4S8,	0001		/2^4/8
	2650
	1171
	4637
	6357
	1425
DFP5S8,	0001		/2^5/8
	3053
	1625
	0212
	5174
	3070
DFP6S8,	0001		/2^6/8
	3272
	1176
	3126
	5516
	5532
DFP7S8,	0001		/2^7/8
	3526
	0143
	3476
	7222
	0722
/
/
DEXA4,	0006		/60.593191717336463
	3622
	7666
	6462
	2157
	5534
DEXB4,	0007		/87.417497202235527
	2566
	5341
	0613
	6705
	7214
DEXC4,	0005		/30.296595858668232
	3622
	7666
	6462
	2157
	5546
DEXD4,	0001		/1.0500
	2063
	1463
	1463
	1463
	1462
DEXH4,	0010		/214.17286814547704
	3261
	3040
	4261
	5654
	0240
DTEMP1,	F 0.0
	F 0.0
DFP2,	F 2.0
	F 0.0
/
	BASE	0
#DEXP,	STARTD
	FLDA	10*3
	FSTA	DEXRTN
	FLDA	0
	SETX	XRDEXP
	SETB	BPDEXP
	BASE	BPDEXP
	LDX	1,1
	LDX	73,2		/FOR ALIGNING
	FSTA	BPDEXP
	FLDA%	BPDEXP,1	/ADDRESS OF X
	FSTA	BPDEXP
	STARTE
	FLDA%	BPDEXP	/GET X
	LDX	0,0
	JGT	DEX1	/CHECK SIGN
	FNEG
	LDX	-1,0	/SET FLAG
DEX1,	JNE	DEX2	/X=0
	FLDA	DEXFP1		/E^0=1
	JA	DEXRTN
DEX2,	FSTA	X
	JA	DEX4
DEX3,	FCLA
	JA	DEXRTN		/RETURN 0 FOR TOO SMALL
/
/SET UP M+N+R=X*LOG2(E)
DEX4,	FLDA	LOG2E
	FMULM	X
	FLDA	X
	ALN	2	/FIX
	FNORM		/FLOAT
	FSTA	M	/INTEGER PART
	FLDA	X
	FSUB	M
	FSTA	N	/FRACTION
	JNE	DEX50	/0 IS SPECIAL CASE
	FLDA	DEXFP1	/1.0
	FSTA	N	/N
	FSTA	R	/R
	JA	DEX20	/SKIP 
/
/CALCULATE N+R
DEX50,	LDX	0,1
	FLDA	N
	FSTA	R	/IF < .125 ALREADY
DEX5,	FSUB	DFP125	/-.125
	JLT	DEX6	/DONE IF .LT.
	FSTA	R	/STORE REMAINDER
	ADDX	1,1	/NEXT POWER OF 2
	JA	DEX5	/AND AGAIN
/
/GET N FROM TABLE
DEX6,	FLDA	DEXFP1,1
	FSTA	N
/
/NOW CALCULATE R
	FLDA	R	/IF R=0
	JNE	DEX7
	FLDA	DEXFP1	/2^R=1
	FSTA	R
	JA 	DEX20	/NO CALCULATION
/
/
DEX7,	FLDA	DEXB4
	FDIV	R	/(B4/R)
	FSTA	X
	FLDA	DEXD4	/D4*R
	FMUL	R
	FADD	X	/+(B4/R)
	FSUB	DEXC4	/-C4
	FSTA	DTEMP1
	FLDA	R
	FADD	X	/R+(B4/R)
	FSTA	R
	FLDA	DEXH4
	FDIV	R	/H4/(R+B4/R)
	FADD	DTEMP1
	FSTA	DTEMP1
	FLDA	DEXA4
	FDIV	DTEMP1
	FADD	DEXFP1
	FSTA	R
/
/CALCULATE 2^M
/
DEX20,	FLDA	M
	JNE	DEX21
	FLDA	DEXFP1
	FSTA	M
	JA	DEX30
DEX21,	FNEG
	ATX	1
	FLDA	DEXFP1
	FSTA	M
	FLDA	DFP2
DEX22,	FMULM	M	/M*2
	JXN	DEX22,1+
/CALCULATE M*N*R
DEX30,	FLDA	M
	FMUL	N
	FMUL	R
	FSTA	X
	JXN	DEX31,0	/WAS X MINUS
	JA	DEXRTN
DEX31,	FLDA	DEXFP1	/.1/X IF -X
	FDIV	X
	JA	DEXRTN
\f\1e