perm filename SAITRG[S,AIL] blob sn#044850 filedate 1973-07-14 generic text, type C, neo UTF8
COMMENT ⊗   VALID 00011 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002		TITLE	TRIG
C00004 00003	↑COSD:				ENTRY TO COSINE DEGREES ROUTINE
C00006 00004	S2:	SKIPGE	A		CHECK SIGN OF ORIGINAL ARG
C00008 00005	BEGIN	ASIN
C00010 00006	BEGIN	ACOS
C00011 00007	BEGIN	SQRT
C00014 00008	BEGIN	ATAN
C00015 00009	↑ATAN:	 			ENTRY TO ARCTANGENT ROUTINE
C00018 00010	BEGIN	ATAN2
C00020 00011	↑LOG:   ---------------------------------------------------------
C00021 ENDMK
C⊗;
	TITLE	TRIG
	INTERN	SIN,COS,SIND,COSD,ASIN,ACOS,SQRT,ATAN,ATAN2

BEGIN	SIN
;FLOATING POINT SINGLE PRECISION SINE AND COSINE FUNCTION
;IF THE ARGUMENT IS IN DEGREES, THE PROPER ENTRY POINTS ARE
;SIND AND COSD, WHILE IF THE ARGUMENT IS IN RADIANS, THE PROPER
;ENTRY POINTS ARE SIN AND COS.
;COSD CALLS SIND TO CALCULATE SIND(PI/2 + X)
;COS CALLS SIN TO CALCULATE SIN(PI/2 + X)
;SIND CALLS SIN AFTER A CONVERSION FROM DEGREES TO RADIANS

;THE ROUTINE CALCULATES SINES AFTER REDUCING THE ARGUMENT TO
;THE FIRST QUADRANT AND CHECKING THE OVERFLOW BITS TO DETERMINE
;THE QUADRANT OF THE ORIGINAL ARGUMENT
;000 - 1ST QUADRANT
;001 - 2ND QUADRANT, X=-(X-PI)
;010 - 3RD QUADRANT, X=-(X-PI)
;011 - 4TH QUADRANT, X=X-3*PI/2-PI/2
;THE ALGORITHM USES A MODIFIED TAYLOR SERIES TO CALCULATE
;THE SINE OF THE NORMALIZED ARGUMENT.
	A←	1
	B←	2
	C←	3
	P←	17
↑COSD:				;ENTRY TO COSINE DEGREES ROUTINE
	MOVE	A, -1(P)	;PICK UP THE ARGUMENT IN B
	FADR	A, [90.0]	;ADD 90 DEGREES
	FDVR	A, SCD1		;CONVERT TO RADIANS
	JRST	S1		;ENTER SINE ROUTINE
↑SIND:				;ENTRY TO SINE DEGREES ROUTINE
	MOVE	A, -1(P)	;PICK UP THE ARGUMENT IN B
	FDVR	A, SCD1		;CONVERT FROM DEGREES TO RADIANS
	JRST	S1		;ENTER SINE ROUTINE
↑COS:				;ENTRY TO COSINE RADIANS ROUTINE
	MOVE	A, -1(P)	;PICK UP THE ARGUMENT IN B
	FADR	A, PIOT		;ADD PI/2
	JRST	S1		;ENTER THE SINE ROUTINE
↑SIN:				;ENTRY TO SINE RADIANS ROUTINE
	MOVE	A, -1(P)	;PICK UP THE ARGUMENT IN A
S1:	MOVM	B,A		;GET ABSF OF ARGUMENT
	CAMG	B, SP2		;SINX = X IF X<2↑-10
	JRST	SUBEXIT
	FDV	B, PIOT		;DIVIDE X BY PI/2
	CAMG	B, [1.0]	;IS X/(PI/2) < 1.0?
	JRST	S2		;YES, ARG IN 1ST QUADRANT ALREADY
	MULI	B, 400		;NO, SEPARATE FRACTION AND EXP.
	ASH	C, -202(B)	;GET X MODULO 2PI
	MOVEI	B, 200		;PREPARE FLOATING FRACTION
	ROT	C, 3		;SAVE 3 BITS TO DETERMINE QUADRANT
	LSHC	B, 33		;ARGUMENT NOW IN RANGE (-1,1)
	FAD	B, [0]		;NORMALIZE THE ARGUMENT
	JUMPE	C, S2		;REDUCED TO FIRST QUAD IF BITS 00
	TLCE	C, 1000		;SUBTRACT 1.0  FROM ARG IF BITS ARE
	FSB	B, [1.0]	;01 OR 11
	TLCE	C, 3000		;CHECK FOR FIRST QUADRANT, 01
	TLNN	C, 3000		;CHECK FOR THIRD QUADRANT, 10
	MOVNS	B		;01,10
S2:	SKIPGE	A		;CHECK SIGN OF ORIGINAL ARG
	MOVNS	B		;SIN(-X) = -SIN(X)
	MOVEM	B, C 		;STORE REDUCED ARGUMENT
	FMPR	B, B		;CALCULATE X↑2
	MOVE	A, SC9		;GET FIRST CONSTANT
	FMP	A, B		;MULTIPLY BY X↑2
	FAD	A, SC7		;ADD IN NEXT CONSTANT
	FMP	A, B		;MULTIPLY BY X↑2
	FAD	A, SC5		;ADD IN NEXT CONSTANT
	FMP	A, B		;MULTIPLY BY X↑2
	FAD	A, SC3		;ADD IN NEXT CONSTANT
	FMP	A, B		;MULTIPLY BY X↑2
	FAD	A, PIOT		;ADD IN LAST CONSTANT
S2B:	FMPR	A, C 		;MULTIPLY BY X
	JRST	SUBEXIT
SC3:	577265210372		;-0.64596371106
SC5:	175506321276		;0.07968967928
SC7:	606315546346		;0.00467376557
SC9:	164475536722		;0.00015148419
SP2:	170000000000		;2**-10
SCD1:	206712273406
PIOT:	201622077325		;PI/2
	BEND			;SINE
BEGIN	ASIN
;FLOATING POINT SINGLE PRECISION ARCSINE FUNCTION
;THE ARCSINE IS CALCULATED WITH THE FOLLOWING ALGORITHM:
;	ASIN(X) = ATAN(X/SQRT(1-X↑2))
;THE RANGE OF DEFINITION FOR ASIN IS (-1.0,1.0)
;OTHER ARGUMENTS WILL GENERATE MEANINGLESS RESULTS
	A←	1
	B←	2
	P←	17
↑ASIN:				;ENTRY TO ASIN ROUTINE
	MOVN	A, -1(P)	;SAVE THE NEGATIVE OF ARG
	FMP	A, -1(P)	;CALCULATE -(X↑2)
	FAD	A, [1.0]	;CALCULATE 1-(X↑2)
	JUMPE	A, ASIN1	;WAS X EITHER -1.0 OR 1.0?
	PUSH	P, A		;NO, CALCULATE SQRT(1-X↑2)
	PUSHJ   P, SQRT		;ADDRESS OF ARGUMENT OF SQRT
	MOVE	B, -1(P)	;GET THE ARGUMENT BACK AGAIN
	FDVM	B, A		;CALCULATE X/SQRT(1-X↑2)
	PUSH	P, A		;CALCULATE ATAN(SQRT(1-X↑2))
	PUSHJ   P, ATAN		;ADDRESS FOR ATAN ROUTINE
	JRST	SUBEXIT
ASIN1:	MOVE	A, PIOT		;ANSWER IS EITHER PI/2 OR-PI/2
	SKIPGE	-1(P)		;WAS ORIGINAL ARGUMENT POSITIVE?
	MOVNS	A		;NO, GET -PI/2
	JRST	SUBEXIT
PIOT:	201622077325		;PI/2
	BEND			;ASIN
BEGIN	ACOS
;FLOATING POINT SINGLE PRECISION ARC-COSINE FUNCTION
;ACOS(X) IS CALCULATED AS
;	ACOS(X)=(PI)/2 - ASIN(X)
;THE RANGE OF DEFINITION FOR ACOS IS (-1.0,1.0)
;ANY OTHER ARGUMENT WILL PRODUCE AN UNDEFINED ANSWER
	A←	1
	P←	17
↑ACOS:				;ENTRY TO SUBROUTINE
	MOVE	A, -1(P)	;GET ARGUMENT IN A
	PUSH    P, A    	;CALL ASIN ROUTINE
	PUSHJ   P, ASIN		;ADDRESS OF ARGUMENT
	MOVNS	A		;GET -(ASIN(X))
	FAD	A, PIOT		;CALCULATE (PI/2)-ASIN(X)
	JRST	SUBEXIT
PIOT:	201622077325
	BEND			;ACOS
BEGIN	SQRT
;FLOATING POINT SINGLE PRECISION SQUARE ROOT FUNCTION
;THE SQUARE ROOT OF THE ABSOLUTE VALUE OF THE ARGUMENT IS
;CALCULATED. THE ARGUMENT IS WRITTEN IN THE FORM
;	X=	F*(2**2B)	WHERE 0<F<1
;SQRT(X) IS THEN CALCULATED AS (SQRT(X))*(2**B)
;SQRT(F) IS CALCULATED BY A LINEAR APPROXIMATION, THE NATURE
;OF WHICH DEPENDS ON WHETHER 1/4 < F < 1/2 OR 1/2 < F < 1,
;FOLLOWED BY TWO ITERATIONS OF NEWTON'S METHOD.
	A←	1
	B←	2
	C←	3
	D←	4
	P←	17
↑SQRT:	 			;ENTRY TO SQUARE ROOT ROUTINE
	MOVE	B,-1(P)		;PICK UP THE ARGUMENT IN B
SQRT1:	JUMPLE	B,[SETZ A, ↔ JRST SUBEXIT]
	ASHC	B, -33		;PUT EXPONENT IN B, FRACTION IN C
	SUBI	B, 201		;SUBTRACT 201 FROM EXPONENT
	ROT	B, -1		;CUT EXP IN HALF, SAVE ODD BIT
	HRRZM	B, D		;SAVE FOR FUTURE SCALING OF ANS
	LSH	B, -43		;GET BIT SAVED BY PREVIOUS INST.
	ASH	C, -10		;PUT FRACTION IN PROPER POSITION
	FSC	C, 177(B)	;PUT EXPONENT OF FRACT TO -1 OR 0
	MOVEM	C, A		;SAVE IT. 1/4 < F < 1
	FMP	C, S1(B)	;LINEAR FIRST APPROX,DEPENDS ON
	FAD	C, S2(B)	;WHETHER 1/4<F<1/2 OR 1/2<F<1.
	MOVE	B, A		;START NEWTONS METHOD WITH FRAC
	FDV	B, C		;CALCULATE X(0)/X(1)
	FAD	C, B		;X(1) + X(0)/X(1)
	FSC	C, -1		;1/2(X(1) + X(0)/X(1))
	FDV	A, C		;X(0)/X(2)
	FADR	A, C		;X(2) + X(0)/X(2)
	FSC	A,@D		;SCALE ANSWER FOR NEWTON AND EXP
↑SUBEXIT:	SUB	P,[XWD 2,2]
		JRST	@2(P)		;EXIT
S1:	0.8125			;CONSTANT, USED IF 1/4<FRAC<1/2
	0.578125		;CONSTANT, USED IF 1/2<FRAC<1
S2:	0.302734		;CONSTANT, USED IF 1/4<FRAC<1/2
	0.421875		;CONSTANT, USED IF 1/2<FRAC<1
	BEND			; SQRT
BEGIN	ATAN
;FLOATING POINT SINGLE PRECISION ARCTANGENT FUNCTION
;ATAN(X) = X(B0+A1(Z+B1-A2(Z+B2-A3(Z+B3)**-1)**-1)**-1)
;WHERE Z=X↑2, IF 0<X<=1
;IF X>1, THEN ATAN(X) = PI/2 - ATAN(1/X)
;IF X>1, THEN RH(D) =-1, AND LH(D) = -SGN(X)
;IF X<1, THEN RH(D) = 0, AND LH(D) =  SGN(X)
	A←	1
	B←	2
	C←	3
	D←	4
	E←	5
	P←	17
↑ATAN:	 			;ENTRY TO ARCTANGENT ROUTINE
	MOVE	A,-1(P)		;PICK UP THE ARGUMENT IN A
ATAN1:	MOVM	B, A		;GET ABSF OF ARGUMENT
	CAMG	B, A1		;IF X<2↑-33, THEN RETURN WITH...
	JRST	SUBEXIT		;ATAN(X) = X
	HLLO	D, A		;SAVE SIGN, SET RH(D) = -1
	CAML	B, A2		;IF A>2↑33, THEN RETURN WITH
	JRST	[MOVE A,PIOT ↔ JRST SUBEXIT];	ATAN(X) = PI/2
	MOVSI	C, 201400	;FORM 1.0 IN C
	CAMG	B, C		;IS ABSF(X)>1.0?
	TRZA	D, -1		;IF B .LE. 1.0, THEN RH(D) = 0
	FDVM	C, B		;B IS REPLACED BY 1.0/B
	TLC	D, (D)		;XOR SIGN WITH .G. 1.0 INDICATOR
	MOVEM	B, E 		;SAVE THE ARGUMENT
	FMP	B, B		;GET B↑2
	MOVE	C, KB3		;PICK UP A CONSTANT
	FAD	C, B		;ADD B↑2
	MOVE	A, KA3		;ADD IN NEXT CONSTANT
	FDVM	A, C		;FORM -A3/(B↑2 + B3)
	FAD	C, B		;ADD B↑2 TO PARTIAL SUM
	FAD	C, KB2		;ADD B2 TO PARTIAL SUM
	MOVE	A, KA2		;PICK UP -A2
	FDVM	A, C		;DIVIDE PARTIAL SUM BY -A2
	FAD	C, B		;ADD B↑2 TO PARTIAL SUM
	FAD	C, KB1		;ADD  B1 TO PARTIAL SUM
	MOVE	A, KA1		;PICK UP A1
	FDV	A, C		;DIVIDE PARTIAL SUM BY A1
	FAD	A, KB0		;ADD B0
	FMP	A, E 		;MULTIPLY BY ORIGINAL ARGUMENT
	TRNE	D, -1		;CHECK .G. 1.0 INDICATOR
	FSB	A, PIOT		;ATAN(A) = -(ATAN(1/A)-PI/2)
	SKIPGE	D		;LH(D) = -SGN(B) IF B>1.0
	MOVNS	A		;NEGATE ANSWER
	JRST	SUBEXIT		;EXIT
A1:	145000000000		;2**-33
A2:	233000000000		;2**33
KB0:	176545543401		;0.1746554388
KB1:	203660615617		;6.762139240
KB2:	202650373270		;3.316335425
KB3:	201562663021		;1.448631538
KA1:	202732621643		;3.709256262
KA2:	574071125540		;-7.106760045
KA3:	600360700773		;-0.2647686202
PIOT:	201622077325		;PI/2
	BEND			;ATAN
BEGIN	ATAN2
;FLOATING POINT ARCTANGENT OF TWO ARGUMENTS
;RETURNS ARCTANGENT OF A/B
;THERE IS NO RESTRICTION ON THE ARGUMENTS
;THE ANSWER IS RETURNED IN ACCUMULATOR A.
	A←	1
	B←	2
	P←	17
↑ATAN2:	 			;ENTRY POINT TO ATAN2 ROUTINE
	MOVM    A,-2(P)
	MOVM    B,-1(P)
	CAML    A,B
	JRST    FOO
	MOVE	A, -2(P)		;PICK UP FIRST ARGUMENT
	FDVR	A, -1(P)	;FORM A/B
	PUSH	P, A   		;CALCULATE ATAN (A/B)
	PUSHJ	P, ATAN
	SKIPL	-1(P)		;IF B>0, SGN(ATAN2)=SGN(A)
	JRST	EXIT2		;EXIT
	JUMPGE	A, ATAN2A	;IS A POSITIVE?
	FADR	A, PI.		;NO, SECOND QUADRANT, ADD PI
	JRST	EXIT2		;EXIT
ATAN2A:	FSB	A, PI.		;YES,3RD QUADRANT, SUBTRACT PI
	JRST	EXIT2		;EXIT
FOO:	MOVN    A,-1(P)
	FDVR    A,-2(P)
	PUSH	P,A
	PUSHJ	P,ATAN
	SKIPG   -2(P)
	JRST	XYZ
	FADR	A,PI2
	JRST	EXIT2
XYZ:	FSB	A,PI2
EXIT2:	SUB	P,[XWD 3,3]
	JRST	@3(P)
PI.:	202622077325
PI2:	201622077325
	BEND			;ATAN2
↑LOG: ;  ---------------------------------------------------------
INTERN LOG
BEGIN LOG
	OPDEF LAC[MOVE]↔OPDEF DAC[MOVEM]
	MOVM -1(17)↔SKIPE 1,0↔CAMN 0,[1.0]↔JRST SUBEXIT
	ASHC 0,-33↔ADDI 0,211000↔MOVSM 0,TMP1#
	MOVSI 0,(-128.5)↔FADM 0,TMP1
	ASH 1,-10↔TLC 1,200000↔FAD 1,[-0.70710678]
	LAC 0,1↔FAD 0,[1.4142135]↔FDV 1,0
	DAC 1,TMP2#↔FMP 1,1
	LAC 0,[0.59897864]↔FMP 0,1
	FAD 0,[0.96147063]↔FMP 0,1
	FAD 0,[2.88539120]↔FMP 0,TMP2↔FAD 0,TMP1
	FMP 0,[0.69314718]↔LAC 1,0↔JRST SUBEXIT
	LIT↔VAR
BEND;-------------------------------------------------------------
END