.TITLE DTRIGS
/COPYRIGHT 1969, DIGITAL EQUIPMENT CORP., MAYNARD, MASS.
/DOUBLE PRECISION TRIG FUNCTIONS
/
/	DSQRT		SQUARE ROOT
/	DATAN		ARCTANGENT
/	DEXP		NATURAL EXPONENTIAL
/	DLOG		NATURAL LOG
/	DLOG10		COMMON LOG
/	DSIN		SINE
/	DCOS		COSINE
/DSQRT	/COMPUTE THE SQUARE ROOT
/OF A D.P. ARG.
/
/
/--------------------
/
	.IFUND	%NDSQ
DSQRT	0
	JMS	%FG.
%DSQ10	0
	JMS	DLAC
	XCT	%DSQ10
	LAC	%FAC2	/GET SIGN WORD
	SPA!CLL		/IF NEGATIVE DISPLAY
	JMP	%DSQB	/ERROR 1
	ISZ	DSQRT
	SNA		/IF ZERO
	JMP*	DSQRT	/EXIT IMMEDIATELY
	LAC	%FAC1	/GET ARG EXPONENT
	SPA		/IF .LT.0 EXTEND EXPONENT
	STL
	RAR
	DAC	%DSQ1	/SET EXP OF CONSTANT (K)=EXP/2
	SNL		/WAS EXPONENT ODD?
	TAD	(-1	/NO...FAC1=EXP/2-1
	DAC	%FAC1	/EXP OF ARG = X/2 OR X/2-1
	JMS	DFAD	/DOUBLE PRECISION ADD
	.DSA	%DSQ1	/(ARG+K=P(0))
	LAW	-4	/SET ITERATION COUNTER
	DAC	%DSQ1
%DSQA	JMS	DDAC	/DOUBLE PRECISION STORE
	.DSA	%DSQ20	/(P(N))
	JMS	DDVR	/DOUBLE PRECISION REVERSE DIVIDE
	XCT	%DSQ10	/ARG/P(N))
	JMS	DFAD
	.DSA	%DSQ20	/(P(N) + (ARG/P(N)))
	LAW	-1	/EXP FROM ARITHMETIC
	TAD	%FAC1	/ABOVE=EXP-1
	DAC	%FAC1	/(DIVIDE ARG BY 2)
	ISZ	%DSQ1	/BUMP ITERATION COUNTER
	JMP	%DSQA	/LOOP BACK AGAIN
	JMP*	DSQRT	/DONE-EXIT
/
%DSQB	STL!GLK		/ERROR 1
	JMP*	DSQRT	/TRY AGAIN
			/STORAGE FOR INPUT
%DSQ20	0
	0
	0
%DSQ1	0
	200000
	0
	.ENDC
/
/--------------------
/
/COMPUTE LOG BASE 2 FOR DOUBLE ARGUMENT
/
/
	.IFUND	%NDLG
	.IFUND	%NDL10
%DLOGS	0
	LAC	%FAC2
	SNA!SPA
	JMP	%DLGA	/ERROR 2 IF.LT.0
	ISZ	%DLOGS	/OK POINT TO NORMAL EXIT
	LAC	%FAC1
	SNA		/CHECK FOR 1.0(SPECIAL CASE)
	JMP	%DLGAA
	AND	(777776
	SZA
	JMP	%DLGAA
	LAC	%FAC3
	SZA
	JMP	%DLGAA
	LAC	%FAC2
	SAD	(200000
	JMP	%DLOGZ
/
/
%DLGAA	LAW	-1	/GET EXP-1
	TAD	%FAC1
	DAC	%DLG2	/FOR INT DETERMINATION
	LAC	(1	/(1)
	DAC	%FAC1	/SET ARG EXP = 1 (1 .LE. F .LT. 2)
	DAC	%DLG4	/SET EXP OF SQRT (2) TO 1
	JMS	DFAD	/ADD DOUBLE
	.DSA	%DLG4	/ (F + SQRT(2)
	ISZ	%DLG4	/BUMP EXP FOR 2 + SQRT(2)
	JMS	DDAC	/STORE DOUBLE
	.DSA	%DLG5	/ (F+SQRT(2),
	JMS	DSUB	/SUBTRACT DOUBLE
	.DSA	%DLG4	/ (F+SQRT(2)-2*SQRT(2)=F-SQRT(2))
	JMS	DDVD	/DIVIDE DOUBLE
	.DSA	%DLG5	/ (F-SQRT(2) / F+SQRT(2))
	JMS	DPOLY	/POLYNOMIAL EVALUATE ABOVE
	.DSA	%DLG7	/ADDR OF CALLING SEQUENCE
	JMS	DDAC	/STORE DOUBLE
	.DSA	%DLG5	/ (RESULT TO TEMP)
	LAC	%DLG2	/GET INTEGER
	STL		/ADD 0.5
	RAL		/ (INTEGER * 2)
	JMS	%FLOT.	/FLOAT
	LAW	-1	/EXP = EXP-1 (INTEGER/2)
	TAD	%FAC1
	DAC	%FAC1
	JMS	DFAD	/ADD DOUBLE
	.DSA	%DLG5	/ (RESULT+INTEGER + 0.5)
	JMP*	%DLOGS	/EXIT
%DLOGZ	DZM	%FAC1
	DZM	%FAC2
	JMP*	%DLOGS
%DLGA	LAC	(2	/ERROR 2
	JMP*	%DLOGS	/ERROR-EXIT
%DLG2		0	/INTEGER STORAGE
%DLG5		0	/TEMP STORAGE	(1)
		0	/		(2)
		0	/		(3)
%DLG7		777775	/NO OF COEFF	(4)
	777777;	336256;	455134	/0.4342597513D0
	000000;	223466;	040146	/0.5765843421D0
	000000;	366161;	114432	/0.9618007623D0
	000002;	270524;	354400	/2.885390073D0
%DLG4		1	/ (X * SQRT(2))	(1)/ (X - 1 OR 2
		265011	/	   	(2)/ AND SET BY
		714640	/		(3)/ PROGRAM)
	.ENDC
	.ENDC
/
/--------------------
/
/DOUBLE PRECISION POLYNOMIAL EVALUATE
/
	.IFUND	%NDPOL
DPOLY	0
	LAC*	DPOLY		/GET PARAM ADDR.
	DAC	%DPOL1		/STORE AS POINTER
	LAC*	%DPOL1		/PICK UP (-N)
	DAC	%DPOL2
	ISZ	%DPOL1		/BUMP POINTER FOR C(N)
	ISZ	DPOLY		/BUMP EXIT ADDR.
	JMS	DDAC		/STORE DOUBLE
	.DSA	%DPOL4		/ (ARG)
	JMS	DMPY		/SQUARE THE INPUT ARG.
	.DSA	%DPOL4
	JMS	DDAC		/STORE DOUBLE
	.DSA	%DPOL6		/ (ARG2)
	JMS	DLAC		/LOAD DOUBLE
	XCT	%DPOL1		/ (C(N))
%DPOLA	JMS	DMPY		/MULTIPLY DOUBLE
	.DSA	%DPOL6		/ (ARG2*C(N))
	ISZ	%DPOL1		/POINTER = POINTER+3
	ISZ	%DPOL1
	ISZ	%DPOL1
	JMS	DFAD		/ADD DOUBLE
	XCT	%DPOL1		/(C(N-1)+(ARGXC(N)))
	ISZ	%DPOL2		/N=N+1
	JMP	%DPOLA		/MORE TERMS-CYCLE AGAIN
	JMS	DMPY		/DONE-MULTIPLY DOUBLE
	.DSA	%DPOL4		/ CUBE ARG.
	JMP*	DPOLY		/EXIT
%DPOL1		0		/POINTER
%DPOL2		0		/COUNTER
%DPOL4		0		/INPUT ARG (1)
		0		/   (2)
		0		/    (3)
%DPOL6		0		/BUILD RESULT (1)
		0		/     (2)
		0		/      (3)
	.ENDC
/
/--------------------
/
/
/DLOG	/COMPUTE THE NATURAL LOGARITHM
/OF A D P FLTING POINT ARG IN THE
/FLTING ACCUMULATOR.
/
/METHOD:	LOG(E)X=LOG(2)X*LOG(E)2
/
	.IFUND	%NDLG
DLOG	0
	JMS	%FG.
	0
	JMS	DLAC
	XCT	.-2
	JMS	%DLOGS	/GET LOG (2) ARG
	JMP*	DLOG	/ERROR EXIT
	JMS	DMPY	/LOG2(ARG)*LOGE(2)
	.DSA	%DLOG1
	ISZ	DLOG	/NORMAL EXIT
	JMP*	DLOG	/EXIT
%DLOG1		0	/LOG(E)2	(1) 0.6931471806
		261344	/	(2)
		137700	/	(3)
	.ENDC
/
/--------------------
/
/DLOG10	/COMPUTE COMMON LOGARITHM
/OF D P ARG IN THE FLTING ACC.
/
/METHOD:	 LOG(10) X = LOG(2)X * LOG(10)2
/
	.IFUND	%NDL10
DLOG10	0
	JMS	%FG.
	0
	JMS	DLAC
	XCT	.-2
	JMS	%DLOGS	/COMPUTE LOG(2) ARG
	JMP*	DLOG10
	JMS	DMPY
	.DSA	%DLG10
	ISZ	DLOG10
	JMP*	DLOG10	/EXIT
%DLG10	777777		/LOG(10)2 (1) 0.3010299957
	232101		/	 (2)
	152052		/	 (3)
	.ENDC
/
/--------------------
/
/DATAN	/COMPUTE THE ARCTANGENT
/OF A D.P. FLTING POINT VALUE.
/
	.IFUND	%NDTAN
DATAN	0
	JMS	%FG.
	0
	JMS	DLAC
	XCT	.-2
	LAC	%FAC2	/GET SIGN WORD
	DAC	%DTN2	/SAVE AS ANS SIGN
	AND	(377777	/STRIP SIGN
	DAC	%FAC2
	LAC	%FAC1	/GET EXPONENT
	DAC	%DTN5	/SAVE FOR LATER TEST
	SNA!SPA		/IF (ARG .GT. 1) TAKE 1/ARG
	JMP	%DATNA	/OR SKIP TO POLY EVALUTE
	JMS	DDVR	/REVERSE DIVIDE DOUBLE
	.DSA	%DTN8	/1/ARG
%DATNA	JMS	DDAC	/SAVE ARGUMENT
	.DSA	%DTN20
	JMS	DPOLY	/POLYNOMIAL EVALUATE
	.DSA	%DTN11
	JMS	DDAC	/SAVE INTERMEDIATE ANGLE
	.DSA	%DTN21
	JMS	%DSIN	/CALCULATE TENGENT OF INTERMEDIATE ANGLE
	NOP
	JMS	DDAC	/TAN (A)=SIN ((PI/2)-A)
	.DSA	%DTN22
	JMS	DLAC
	.DSA	%DTN14
	JMS	DSUB
	.DSA	%DTN21
	JMS	%DSIN
	NOP
	JMS	DDVR
	.DSA	%DTN22
	JMS	DDAC	/SAVE TANGENT OF INTERMEDIATE ANGLE
	.DSA	%DTN22
	JMS	DMPY	/CALCULATE CORRECTION ANGLE TO BRING
	.DSA	%DTN20	/INTERMEDIATE ANGLE TO DESIRED ACCURACY
	JMS	DFAD	/C.A.=(ARG-TAN(A))/(1+ARG*TAN(A))
	.DSA	%DTN8
	JMS	DDAC
	.DSA	%DTN24
	JMS	DLAC
	.DSA	%DTN20
	JMS	DSUB
	.DSA	%DTN22
	JMS	DDVD	/C.A. NOW IN FLOATING AC.
	.DSA	%DTN24
	~JMS	DFAD	/ADD INTERMEDIATE ANGLE TO CORRECTION
	.DSA	%DTN21	/ANGLE TO GET FINAL ANGLE
	LAC	%DTN5	/TEST ARG AGAIN IF .GT. 1 (EXP .GT. 0)
	SNA!SPA		/ (EXP POSITIVE AND NON-ZERO)
	JMP	%DATNB	/IS NOT .GT. 1, SKIP
	JMS	DSBR	/REVERSE SUBTRACT DOUBLE
	.DSA	%DTN14	/ (PI/2-ANS)
%DATNB	LAC	%DTN2	/GET ORIGINAL SIGN WORD
	AND	(400000	/KEEP BIT 0 ONLY
	XOR	%FAC2	/SET SIGN WORD WITH NEW SIGN
	DAC	%FAC2	/RESTORE SIGN WORD
	JMS	%FNOR.
	JMP*	DATAN
%DTN2	0	/ANS SIGN
%DTN5	0	/EXP STORAGE
%DTN8	1	/FLOATING 1 	(1)
	200000	/	    (2)
	0	/		(3)
%DTN14	1	/PI/2(1) (1.57079632679)
	311037	/      (2)
	552420	/      (3)
%DTN20	0	/TEMP STORAGE
	0
	0
%DTN21	0
	0
	0
%DTN22	0
	0
	0
%DTN24	0
	0
	0
%DTN11	777775		/NO OF COEFFICIENTS
	777774;	637556;	023444	/-0.389929D-1
	777776;	225623;	041646	/ 0.1462766D0
	777777;	644343;	720712	/ -0.3211819D0
	000000;	377631;	067426	/0.9992150D0
	.ENDC
/
/--------------------
/
/DCOS	/COMPUTE THE APPROX COSINE
/OF A D. P. FLT POINT VALUE IN
/THE FLTING ACCUMULATOR.
/
/METHOD: COS X=SIN (PI/2+X)
/
	.IFUND	%NDCOS
DCOS	0
	JMS	%FG.
	0
	JMS	DLAC
	XCT	.-2
	JMS	DFAD	/DOUBLE ADD PI/2
	.DSA	%DCOS1	/TO THE ARG AND
	JMS	%DSIN	/COMPUTE SIGN
	JMP*	DCOS	/ERROR EXIT
	ISZ	DCOS	/POINT TO
	JMP*	DCOS	/NORMAL EXIT
%DCOS1	1		/(PI/2) (1)
	311037		/       (2)
	552424		/       (3)
	.ENDC
/
/--------------------
/
/DSIN	/COMPUTE THE SINE
/
	.IFUND	%NDSIN
DSIN	0
	JMS	%FG.
	0
	JMS	DLAC
	XCT	.-2
	JMS	%DSIN
	JMP*	DSIN
	ISZ	DSIN	/POINT TO NORMAL
	JMP*	DSIN	/EXIT
	.ENDC
/
/--------------------
/
/COMPUTE SINE OF
/D P ARG IN THE FLOATING
/ACCUMULATOR.
/
	.IFUND	%NDCOS
	.IFUND	%NDSIN
	.IFUND	%NDTAN
%DSIN	0
	LAC	%FAC2	/GET SIGN WORD
	AND	(400000
	DAC	%DSIN1	/SAVE SIGN (ANS SIGN)
	LAC	%FAC2
	AND	(377777	/STRIP OFF SIGN BIT
	DAC	%FAC2	/PUT BACK ABSOLUTE VALUE
	JMS	DMPY	/MULTIPLY DOUBLE
	.DSA	%DSN20	/(2/PI * ARG)
	JMS	DDAC	/STORE DOUBLE
	.DSA	%DSN10	/(TEMP)
	JMS	%FUNF.	/UNFLOAT 
	JMP*	%DSIN	/ERROR EXIT
	DAC	%DSIN2	/SAVE FOR QUADRANT DETERMINATION
	JMS	%FLOT.	/FLOAT 
	JMS	DSBR	/REVERSE SUBTRACT DOUBLE
	.DSA	%DSN10	/(TEMP(ARG*2/PI)-ACC(INT ARG*2/PI) FOR FRACTN.
	LAC	%DSIN2	/GET INTEGER
	RCR		/BIT 17 TO LINK
	SNL		/IF SET COMPUTE 1-F (QUAD II OR IV)
	JMP	%DSINA	/OR TEST BIT 16
	JMS	DSBR	/REVERSE SUBTRACT DOUBLE
	.DSA	%DSIN3	/(1-F)
%DSINA	LAC	%DSIN2	/GET INTEGER AGAIN
	RTR		/BIT 16 TO LINK
	SNL		/IF SET-NEG ANS (QUAD III OR IV)
	JMP	%DSINB	/ COMPUTE THE POLY
	LAC	(400000	/GET ANS SIGN BIT
	XOR	%DSIN1	/CHANGE IT
	DAC	%DSIN1	/PUT IT BACK
%DSINB	JMS	DPOLY	/POLY EVALUATE
	.DSA	%DSINC	/NO OF COEFFS.
	LAC	%FAC2	/GET SIGN WORD OF ANS
	AND	(377777
	XOR	%DSIN1	/INSERT SIGN
	DAC	%FAC2
	ISZ	%DSIN	/NORMAL
	JMS	%FNOR.
	JMP*	%DSIN	/EXIT
/THE COEFFICIENTS ARE
%DSINC	777772
	777750; 351732; 333236 /0.54465285D-7
	777756; 761211; 460000 /-0.3595184353D-5
	777764; 250166; 553372 /0.16043839964D-3
	777771; 631322; 620016 /-0.4681752998D-2
	777775; 243153; 614634 /+0.796926260D-1
	000000; 645273; 634612 /-0.6459640975D0
	000001; 311037; 552422 /+1.57079632680D0
%DSIN1	0
%DSIN2	0
%DSN10	0
	0
	0
%DSN20	0		/(2/PI) (1)
	242763		/       (2)
	015566		/       (3)
%DSIN3	1		/FLOATING 1.0 (1)
	200000		/             (2)
	0		/             (3)
	.ENDC
	.ENDC
	.ENDC
/
/--------------------
/
/DEXP	/COMPUTE NATURAL EXPONEMTIAL
/OF A D. P. FLTING POINT NUMBER IN THE ACC.
/
	.IFUND	%NDEXP
DEXP	0
	JMS	%FG.
	0
	JMS	DLAC
	XCT	.-2
	LAC	%FAC2	/GET SIGN WORD
	DAC	%DEXP1	/SAVE SIGN OF ARG
	AND	(377777	/GET THE ABSOLUTE VALUE
	DAC	%FAC2
	JMS	DMPY
	.DSA	%DEX10	/ARG *LOG2 (E))
	JMS	DDAC	/STORE DOUBLE
	.DSA	%DEX20	/(TEMP)
	JMS	%FUNF.	/UNFLOAT PRODUCT
	JMP*	DEXP
	DAC	%DEX30	/STORE INT AS ANS EXPONENT
	ISZ	%DEX30	/BUMP EXPONENT FOR .5 SCALING
	JMS	%FLOT.	/FLOAT THE INTEGER
	JMS	DSBR	/REVERSE SUBTRACT DOUBLE
	.DSA	%DEX20	/(TEMP - FAC) FRACTION ONLY (F)
	JMS	DDAC	/STORE DOUBLE
	.DSA	%DEX20	/(F  INTO TEMP)
	LAW	-10	/SET ITERATION CNTR FOR 9 PASSES
	DAC	%DEXP2
	LAC	%DEXD	/GET ADDRESS OF  COEFFICIENTS
	DAC	%DEXP3	/SAVE IT AS A POINTER
	JMS	DLAC	/LOAD 
	.DSA	%DEXC	/C(9) AS TERM
%DEXPA	JMS	DMPY	/MULTIPLY DOUBLE
	.DSA	%DEX20	/(F*TERM)
%DEXPB	ISZ	%DEXP3	/POINT TO THE NEXT
	ISZ	%DEXP3
	ISZ	%DEXP3
	JMS	DFAD	/ADD DOUBLE
	XCT	%DEXP3	/C (PTR) AS NEW TERM
	ISZ	%DEXP2	/INCREMENT COUNTER
	JMP	%DEXPA	/NOT DONE CYCLE AGAIN
	JMS	DDAC	/DONE
	.DSA	%DEX20	/(RESULT TO TEMP)
	JMS	DMPY
	.DSA	%DEX20	/(RESULT*RESULT)
	JMS	DMPY
	.DSA	%DEX30	/(EXPONENT* FINAL RESULT)
	ISZ	DEXP	/NORMAL EXIT
	LAC	%DEXP1	/IF THE INPUT
	SMA		/ARG WAS .LT.0
	JMP*	DEXP	/GET ITS RECIPROCAL
	JMS	DDVR
	.DSA	%DEXP4	/BEFORE EXITING
	JMP*	DEXP
%DEXC	777745;	261767;	414110	/0.5180D-8
	777752;	200334;	225756	/0.119610D-6
	777756;	241410;	444530	/0.2406787D-5
	777762;	256607;	213216	/0.41667033D-4
	777766;	235452;	556272	/0.601133075D-3
	777771;	343260;	433536	/0.6938013677D-2
	777774;	365773;	677740	/0.6005662674D-1
	777777;	261344;	137700
%DEXP4     1;	200000;	000000	/1.0D0
%DEXD	.DSA	%DEXC
%DEXP1	0
%DEXP2	0			/COUNTER
%DEXP3	0			/POINTER
%DEX10	1			/LOG2 (E)  (1)
	270524			/	 (2)
	354514			/	 (3)
%DEX20	0
	0
	0
%DEX30	0
	200000
	0
	.ENDC
	.EOT
?*U*/