00100	TITLE	SAITRG
00200	
00250	;ALTERNATE PDP-10 MNEMONICS.
00300		OPDEF LAC[MOVE]↔OPDEF DAC[MOVEM]↔OPDEF GO[JRST]
00400		OPDEF LACM[MOVM]↔OPDEF LACN[MOVN]↔OPDEF DAP[HRRM]
00450	;SAIL CONVENTIONS.
00500		DEFINE SUBR(NAME){INTERN NAME↔NAME:}
00600		DEFINE ARG1<-1(17)>
00610		DEFINE ARG2<-2(17)>
00700		DEFINE POP1J<GO POP1J.>
00750		DEFINE POP2J<GO POP2J.>
00775	
00800		HALFPI:	201622077325 ;PI/2
00850		PI:	202622077325
00900	
01000	INTERN SIN,COS
01100	BEGIN SINCOS
01200		A←1 ↔ B←2 ↔ C←3
01300	↑COS:	SKIPA A,ARG1
01400	↑SIN:	SKIPA A,ARG1
01500		FADR  A,HALFPI			;COS(X) = SIN(X+π/2).
01600		MOVM B,A↔CAMG B,[17B5]↔POP1J	;FOR SMALL X, SIN(X)=X.
01700	
01800	;B ← (ABS(X)MODULO 2π)/HALFPI
01900	;C ← QUADRANT 0, 1, 2 OR 3.
02000		FDVR B,HALFPI
02100		LAC C,B↔FIX C,233000
02200		CAILE C,3↔GO[
02300		TRZ C,3↔FSC C,233
02400		FSBR B,C↔GO .-3]		;MODULO 2π.
02500		GO .+1(C)↔GO .+4↔JFCL↔GO[
02650		FSBRI B,(2.0)↔MOVNS B↔GO .+2]	;SIN(X+π)=SIN(-X)
02700		FSBRI B,(4.0)			;SIN(X+2π)=SIN(X)
02800		SKIPGE A↔MOVNS	B		;SIN(-X) = -SIN(X).
02900	
03000	;FOR -1 ≤ B ≤ +1 REPRESENTING -π/2 ≤ X ≤ +π/2,
03100	;COMPUTE SINE(X) APPROXIMATION BY TAYLOR SERIES.
03200		DAC B,C↔FMPR B,B	
03300		LAC A,[164475536722]↔FMP A,B
03400		FAD A,[606315546346]↔FMP A,B
03500		FAD A,[175506321276]↔FMP A,B
03600		FAD A,[577265210372]↔FMP A,B
03700		FAD A,HALFPI↔FMPR A,C↔POP1J
03800		LIT
03900	BEND
     

00100	;ACOS(X)= π/2 - ASIN(X).
00200	;GIVEN -1 ≤ X ≤ +1 RETURN 0 ≤ ACOS(X) ≤ +π.
00300	SUBR(ACOS)
00400	BEGIN ACOS
00500		PUSH 17,ARG1↔PUSHJ 17,ASIN
00600		MOVNS 1↔FADR 1,HALFPI↔POP1J
00700	BEND
00800	
00900	;ASIN(X)=ATAN(X/SQRT(1-X↑2)).
01000	;GIVEN -1 ≤ X ≤ +1 RETURN -π/2 ≤ ASIN(X) ≤ +π/2.
01100	SUBR(ASIN)
01200	BEGIN ASIN
01300		A←1 ↔ B←2
01400		LACN A,ARG1↔FMPR A,ARG1↔FADRI A,(1.0)
01500		JUMPE A,[		;WAS X EITHER -1.0 OR 1.0?
01600			LAC A,HALFPI
01700			SKIPGE ARG1
01800			MOVNS A↔POP1J]
01900		PUSH 17,A↔PUSHJ 17,SQRT
02000		LAC B,ARG1↔FDVR B,1↔DAC B,ARG1	;CALCULATE X/SQRT(1-X↑2)
02100		GO ATAN			;CALCULATE ATAN(SQRT(1-X↑2))
02200	BEND
02300	
02400	SUBR(LOG)
02500	BEGIN LOG
02600		MOVM ARG1↔SKIPE 1,0↔CAMN 0,[1.0]↔POP1J
02700		ASHC 0,-33↔ADDI 0,211000↔MOVSM 0,TMP1#
02800		MOVSI 0,(-128.5)↔FADM 0,TMP1
02900		ASH 1,-10↔TLC 1,200000↔FAD 1,[-0.70710678]
03000		LAC 0,1↔FAD 0,[1.4142135]↔FDV 1,0
03100		DAC 1,TMP2#↔FMP 1,1
03200		LAC 0,[0.59897864]↔FMP 0,1
03300		FAD 0,[0.96147063]↔FMP 0,1
03400		FAD 0,[2.88539120]↔FMP 0,TMP2↔FAD 0,TMP1
03500		FMP 0,[0.69314718]↔LAC 1,0↔POP1J
03550		LIT↔VAR
03600	BEND
     

00100	SUBR(SQRT)
00200	BEGIN	SQRT
00300		A←1 ↔ B←2 ↔ C←3
00400		LACM B,ARG1↔JUMPE B,POP1J.
00500	
00600	;LET X=F*(2↑2B) WHERE 0.25<F<1.00 THEN SQRT(X)=SQRT(F)*(2↑B).
00700		ASHC B,-=27↔SUBI B,201	;GET EXPONENT IN B, FRACTION IN C.
00800		ROT B,-1		;CUT EXP IN HALF, SAVE ODD BIT
00900		DAP B,L↔LSH B,-=35	;USE THAT ODD BIT.
01000		ASH C,-10↔FSC C,177(B)	;0.25 < FRACTION < 1.00
01100	
01200	;LINEAR APPROXIMATION TO SQRT(F).
01300		DAC C,A
01400		FMP C,[0.8125↔0.578125](B)
01500		FAD C,[0.302734↔0.421875](B)
01600	
01700	;TWO ITERATIONS OF NEWTON'S METHOD.
01800		LAC B,A
01900		FDV B,C↔FAD C,B↔FSC C,-1
02000		FDV A,C↔FADR A,C↔L: FSC A,0
02100	↑POP1J.: SUB 17,[XWD 2,2]
02200		GO @2(17)
02300		LIT
02400	BEND
     

00001	;ATAN(X) = X*(B0+A1 / (Z+B1-A2 / (Z+B2-A3 / (Z+B3))) )
00020	;WHERE Z=X↑2, IF 0<X<=1
00039	;IF X>1, THEN ATAN(X) = PI/2 - ATAN(1/X)
00058	;IF X>1, THEN RH(D) =-1, AND LH(D) = -SGN(X)
00077	;IF X<1, THEN RH(D) = 0, AND LH(D) =  SGN(X)
00100	SUBR(ATAN)
00200	BEGIN ATAN
00300		A←1 ↔ B←2 ↔ C←3 ↔ D←4 ↔ E←5 ↔ P←17
01000		LAC	A,ARG1		;PICK UP THE ARGUMENT IN A
01100	ATAN1:	LACM	B, A		;GET ABSF OF ARGUMENT
01200		CAMG	B, A1		;IF X<2↑-33, THEN RETURN WITH...
01300		POP1J		;ATAN(X) = X
01400		HLLO	D, A		;SAVE SIGN, SET RH(D) = -1
01500		CAML	B, A2		;IF A>2↑33, THEN RETURN WITH
01600		GO[LAC A,HALFPI ↔POP1J];	ATAN(X) = PI/2
01700		MOVSI	C, 201400	;FORM 1.0 IN C
01800		CAMG	B, C		;IS ABSF(X)>1.0?
01900		TRZA	D, -1		;IF B ≤ 1.0, THEN RH(D) = 0
02000		FDVM	C, B		;B IS REPLACED BY 1.0/B
02100		TLC	D, (D)		;XOR SIGN WITH > 1.0 INDICATOR
02200	
02300		DAC B,E↔FMP B,B
02400		LAC C,B↔FAD C,KB3↔LAC A,KA3↔FDVM A,C
02500		FAD C,B↔FAD C,KB2↔LAC A,KA2↔FDVM A,C
02600		FAD C,B↔FAD C,KB1↔LAC A,KA1↔FDV  A,C
02700		FAD A,KB0↔FMP A,E
02800	
02900		TRNE	D, -1		;CHECK > 1.0 INDICATOR
03000		FSB	A, HALFPI		;ATAN(A) = -(ATAN(1/A)-PI/2)
03100		SKIPGE	D		;LH(D) = -SGN(B) IF B>1.0
03200		MOVNS A		;NEGATE ANSWER
03300		POP1J		;EXIT
03400	A1:	145000000000		;2↑-33
03500	A2:	233000000000		;2↑33
03600	
03700	KB0:	176545543401		;0.1746554388
03800	KB1:	203660615617		;6.762139240
03900	KB2:	202650373270		;3.316335425
04000	KB3:	201562663021		;1.448631538
04100	
04200	KA1:	202732621643		;3.709256262
04300	KA2:	574071125540		;-7.106760045
04400	KA3:	600360700773		;-0.2647686202
04500	BEND
     

00100	; OMEGA ← ATAN2(Y,X).
00200	SUBR(ATAN2)
00300	BEGIN	ATAN2
00400		Y←1 ↔ X←2
00500		LACM Y,ARG2↔LACM X,ARG1
00600		CAML Y,X↔GO L1
00700	
00800	;HORIZONTAL TO π/2; ABS(Y) < ABS(X).
00900		LAC  Y,ARG2↔FDVR Y,ARG1
01000		PUSH 17,Y↔PUSHJ 17,ATAN	; ARCTAN(Y/X)
01100		SKIPL ARG1↔POP2J	;1ST & 2ND QUADRANTS.
01200		JUMPGE Y,[
01300		FSBR Y,PI.↔POP2J]	;3RD QUADRANT.
01400		FADR Y,PI.↔POP2J	;2ND QUADRANT.
01500	
01600	;VERTICAL TO π/2; ABS(X) < ABS(Y).
01700	L1:	LACN X,ARG1↔FDVR X,ARG2
01800		PUSH 17,X↔PUSHJ 17,ATAN	;ARCTAN(X/Y)
01900		SKIPG ARG2↔GO[
02000		FSB Y,PI2↔POP2J]
02100		FADR	Y,PI2
02200	POP2J.:	SUB 17,[XWD 3,3]↔GO @3(17)
02300	
02400	PI.:	202622077325
02500	PI2:	201622077325
02600	BEND
02700	
02800	END
02900	SAITRG.FAI - EOF.