00100	TITLE LSPTRG - LISP TRIG FUNCTIONS - BGB - 7 APRIL 1973.
00200		EXTERN NUMVAL,MAKNUM
00300	
00400	;ALTERNATE PDP-10 MNEMONICS.
00500		OPDEF LAC[MOVE]↔OPDEF DAC[MOVEM]↔OPDEF GO[JRST]
00600		OPDEF LACM[MOVM]↔OPDEF LACN[MOVN]↔OPDEF DAP[HRRM]
00700	;LISP CONVENTIONS.
00800		DEFINE SUBR(NAME){INTERN NAME↔NAME:}
00900		DEFINE ARG1<1>
01000		DEFINE ARG2<2>
01100		DEFINE POP1J<GO POP1J.>
01200		DEFINE POP2J<GO POP2J.>
01300	
01400		HALFPI:	201622077325 ;PI/2
01500		PI:	202622077325
01600	
01700	INTERN SIN,COS
01800	BEGIN SINCOS
01900		A←1 ↔ B←2 ↔ C←3
02000	↑COS:	PUSHJ 14,NUMVAL
02100		HRLI 2,(<MOVEI 2,>)↔CAME 2,MAKNUM+2↔GO[FSC 1,233
02200		FMPR 1,[0.017453292]↔GO .+1]
02300		FADR A,HALFPI↔GO SIN+4	;COS(X) = SIN(X+π/2).
02400	↑SIN:	PUSHJ 14,NUMVAL
02500		HRLI 2,(<MOVEI 2,>)↔CAME 2,MAKNUM+2↔GO[FSC 1,233
02600		FMPR 1,[0.017453292]↔GO .+1]
02700		MOVM B,A↔CAMG B,[17B5]↔POP1J	;FOR SMALL X, SIN(X)=X.
02800	
02900	;B ← (ABS(X)MODULO 2π)/HALFPI
03000	;C ← QUADRANT 0, 1, 2 OR 3.
03100		FDVR B,HALFPI
03200		LAC C,B↔FIX C,233000
03300		CAILE C,3↔GO[
03400		TRZ C,3↔FSC C,233
03500		FSBR B,C↔GO .-3]		;MODULO 2π.
03600		GO .+1(C)↔GO .+4↔JFCL↔GO[
03700		FSBRI B,(2.0)↔MOVNS B↔GO .+2]	;SIN(X+π)=SIN(-X)
03800		FSBRI B,(4.0)			;SIN(X+2π)=SIN(X)
03900		SKIPGE A↔MOVNS	B		;SIN(-X) = -SIN(X).
04000	
04100	;FOR -1 ≤ B ≤ +1 REPRESENTING -π/2 ≤ X ≤ +π/2,
04200	;COMPUTE SINE(X) APPROXIMATION BY TAYLOR SERIES.
04300		DAC B,C↔FMPR B,B	
04400		LAC A,[164475536722]↔FMP A,B
04500		FAD A,[606315546346]↔FMP A,B
04600		FAD A,[175506321276]↔FMP A,B
04700		FAD A,[577265210372]↔FMP A,B
04800		FAD A,HALFPI↔FMPR A,C↔POP1J
04900		LIT
05000	BEND
     

00100	;ACOS(X)= π/2 - ASIN(X).
00200	;GIVEN -1 ≤ X ≤ +1 RETURN 0 ≤ ACOS(X) ≤ +π.
00300	SUBR(ACOS)
00400	BEGIN ACOS
00500		PUSHJ 14,ASIN
00600		PUSHJ 14,NUMVAL
00700		MOVNS 1↔FADR 1,HALFPI
00800		POP1J
00900	BEND
01000	
01100	;ASIN(X)=ATAN(X/SQRT(1-X↑2)).
01200	;GIVEN -1 ≤ X ≤ +1 RETURN -π/2 ≤ ASIN(X) ≤ +π/2.
01300	SUBR(ASIN)
01400	BEGIN ASIN
01500		A←1 ↔ B←2
01600		PUSH 14,NUMVAL
01700		LACN A,ARG1↔FMPR A,ARG1↔FADRI A,(1.0)
01800		JUMPE A,[		;WAS X EITHER -1.0 OR 1.0?
01900			LAC A,HALFPI
02000			SKIPGE ARG1
02100			MOVNS A↔POP1J]
02200		PUSHJ 14,SQRT↔PUSHJ 14,NUMVAL
02300		LAC B,ARG1↔FDVR B,1↔DAC B,ARG1	;CALCULATE X/SQRT(1-X↑2)
02400		GO ATAN+1		;CALCULATE ATAN(SQRT(1-X↑2))
02500	BEND
02600	
     

00100	SUBR(LOG)
00200	BEGIN LOG
00300		PUSHJ 14,NUMVAL
00400		HRLI 2,(<MOVEI 2,>)↔CAME 2,MAKNUM+2↔FSC 1,233
00500		MOVM 2,ARG1↔SKIPE 1,2↔CAMN 2,[1.0]↔POP1J
00600		ASHC 2,-33↔DAC 3,1
00700		ADDI 2,211000↔MOVSM 2,TMP1#
00800		MOVSI 2,(-128.5)↔FADM 2,TMP1
00900		ASH 1,-10↔TLC 1,200000↔FAD 1,[-0.70710678]
01000		LAC 2,1↔FAD 2,[1.4142135]↔FDV 1,2
01100		DAC 1,TMP2#↔FMP 1,1
01200		LAC 2,[0.59897864]↔FMP 2,1
01300		FAD 2,[0.96147063]↔FMP 2,1
01400		FAD 2,[2.88539120]↔FMP 2,TMP2↔FAD 2,TMP1
01500		FMP 2,[0.69314718]↔LAC 1,2↔POP1J
01600		LIT↔VAR
01700	BEND
     

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

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