TITLE	TAN	%EXTERNALREALFN TAN(%REAL X)

;THIS FUNCTION CALCULATES TAN(X) AS SIN(X)/COS(X)
;IF COS(X)=0 AN EVENT IS SIGNALLED

	SEARCH	IMPPRM
	ENTRY	TAN
	EXTERN	SIN
	EXTERN	COS
	TWOSEG	400000


TAN:	MOVEM	ARG1,ARGSAV		;SAVE THE ARGUMENT FOR USE LATER
	CALL	COS			;GET COS(X)
	JUMPE	AC,INFRES		;ERROR IF COS(ARG)=0
	MOVEM	AC,COSSAV		;SAVE THE COS FOR LATER
	MOVE	ARG1,ARGSAV		;GET THE ORIGINAL ARG BACK
	CALL	SIN			;GET SIN(X)
	FDVR	AC,COSSAV		;TAN=SIN/COS
	RETURN				;TAN RETURN

INFRES:	ERROR(1,1,AC,IMPSTR<INFINITE RESULT FROM TAN>)

	RELOC	0
ARGSAV:	0
COSSAV:	0

	PRGEND



TITLE	RAISE	%EXTERNALREALFN RAISE(%REAL NUMBER,TO POWER)

;THIS FUNCTION CALCULATES ARG1**ARG2 BY THE ALGORITHM:
;   A**B=EXP(B*LOG(A))
;THE CALLING SEQUENCE IS:
;	MOVE	ARG1,A
;	MOVE	ARG2,B
;	PUSHJ	P,RAISE
;THE RESULT IS RETURNED IN AC

	SEARCH	IMPPRM
	TWOSEG	400000
	ENTRY	RAISE
	EXTERN	LOG
	EXTERN	EXP

RAISE:	MOVEM	ARG2,EXPSAV		;SAVE THE POWER FOR LATER
	JUMPL	ARG1,ERR		;ERROR IF THE BASE IS<0
	JUMPE	ARG1,ZERARG		;ARG=0 IMPLIES RESULT=0
	CALL	LOG			;CALC LOG(A)
	FMPR	AC,EXPSAV		;GET B*LOG(A)
	MOVE	ARG1,AC
	CALL	EXP			;GET THE ANSWER IN AC
	RETURN				;RAISE RETURN

ERR:	ERROR(1,1,ARG1,IMPSTR<NEGATIVE BASE GIVEN TO RAISE>)

ZERARG:	MOVE	AC,ONE			;RESULT=1
	RETURN				;RAISE RETURN

	RELOC	0
EXPSAV:	0
ONE:	1.0			;CONSTANT 1

	PRGEND
TITLE	ACOS  %EXTERNALREALFN ACOS(%REAL X)
	SEARCH	IMPPRM
	ENTRY	ACOS
	EXTERN	SQR
	EXTERN	ATAN

;FLOATING POINT SINGLE PRECISION ARCCOSINE FUNCTION

;ACOS(X) IS CALCULATED IN THE FOLLOWING MANNER:
;	IF X > 0,	ACOS(X)=ATAN((SQRT(1-X^2))/X)
;	IF X < 0,	ACOS(X)=PI + ATAN((SQRT(1-X^2))/X)
;	IF X = 0,	ACOS(X)=PI/2

;THE RANGE OF DEFINITION FOR ACOS IS -1.0 TO +1.0.
;ARGUMENTS OUTSIDE OF THIS RANGE WILL CAUSE AN ERROR MESSAGE

;THE CALLING SEQUENCE FOR ACOS IS:
;	MOVE	ARG1,X
;	PUSHJ	P,ACOS


	TWOSEG	400000
ACOS:	MOVM	0,ARG1		;GET /ARG./ IN AC 0.
	CAMLE	0,ONE		;IS MAGNITUDE OF ARG. GT 1.0?
	JRST	TOOLRG		;YES, GO TO ERROR RETURN.
	JUMPE	0,ZERARG	;IF ARG=0, GO TO ZERARG.
	FMPR	0,0		;X^2 IN AC 0.
	JFCL			;SUPPRESS ERROR MESSAGE FROM OVTRAP, IF NECESSARY.
	MOVNS	0		;-X^2 IN AC 0.
	FAD	0,ONE		;1.0-X^2 IN AC 0.
	MOVEM	ARG1,ARGSAV	;SAVE THE ORIGINAL ARGUMENT FOR LATER
	MOVE	ARG1,0		;GET ARGUMENT FOR SQUARE ROOT
	CALL	SQR		;CALC	SQRT(1.0-X^2)
	IFN AC <MOVE 0,AC>	;GET RESULT IN 0 IF IT ISN'T THERE
	FDVR	0,ARGSAV	;(SQRT(1.0-X^2))/X IS IN AC 0
	JFCL			;SUPPRESS ERROR MESSAGE FROM OVTRAP, IF NECESSARY.
	MOVE	ARG1,0		;GET THE ARGUMENT FOR ATAN
	PUSHJ	P,ATAN		;FIND ATAN(SQRT(1.0-X^2)/X).
	SKIPGE	ARGSAV		;SKIP IF ORIGINAL ARG >= 0.
	FAD	0,PII		;ANSWER IS PI - ORIGINAL ANSWER.
	IFN AC<MOVE AC,0>	;IN CASE AC<>0
	RETURN			;AND RETURN FROM ACOS
ZERARG:	MOVE	AC,PI2		;ANSWER IS PI/2.
	RETURN			;ACOS RETURN
TOOLRG:	ERROR(6,2,0,IMPSTR<ILLEGAL ARG TO ACOS>)
				;CAUSE A LIBRARY FUNCTION EVENT

	RELOC	0
ONE:	201400000000		;1.0
PI2:	201622077325
PII:	202622077325
ARGSAV:	0
	PRGEND



TITLE	ASIN	%EXTERNALREALFN ASIN(%REAL X)
;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 CAUSE AN EVENT

	A=	0
	B=	1

	SEARCH	IMPPRM
	TWOSEG	400000
	ENTRY	ASIN
	EXTERN	ATAN
	EXTERN	SQR

ASIN:	MOVM	B,ARG1		;GET MAGNITUDE OF ARG. IN B
	MOVEM	ARG1,ASIN2	;SAVE THE ARGUMENT FOR FURTHER USE
	CAMLE	B,ONE		;IS THE MAGNITUDE OF THE ARG. LE 1.0?
	JRST	TOOLRG		;NO, GO TO ERROR RETURN.
	MOVN	A,ARG1		;GET THE NEGATIVE OF ARG
	FMP	A,ARG1		;CALCULATE -(X^2)
	JFCL			;SUPPRESS ERROR MESSAGE FROM OVTRAP, IF NECESSARY.
	FAD	A, ONE		;CALCULATE 1-(X^2)
	JUMPE	A, ASIN1	;WAS X EITHER -1.0 OR 1.0?
	MOVE	ARG1,A		;SET UP FOR CALL TO SQR
	CALL	SQR		;CALCULATE SQRT(1-X^2)
	MOVE	B,ASIN2		;GET THE ARGUMENT BACK AGAIN
	FDV	B, A		;CALCULATE X/SQRT(1-X^2)
	MOVEM	B,ARG1		;THEN
	CALL	ATAN		;CALCULATE ATAN(X/SQR(1-X^2))
	RETURN			;ASIN RETURN

ASIN1:	MOVE	AC,PIOT	;ANSWER IS EITHER PI/2 OR-PI/2
	SKIPG	ASIN2		;WAS ORIGINAL ARGUMENT POSITIVE?
	MOVNS	AC		;NO, GET -PI/2
	RETURN			;ASIN RETURN

TOOLRG:	ERROR(6,2,0,IMPSTR<ILLEGAL ARG TO ASIN>)

	RELOC 0			;NOW THE DATA SEGMENT

ASIN2:	0			;STORAGE FOR ARGUMENT
PIOT:	201622077325		;PI/2
ONE:	1.0
	PRGEND
TITLE	ATAN	%EXTERNALREALFN ATAN(%REAL X)

;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)

;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
;	MOVE	ARG1,X
;	PUSHJ	P,ATAN
;THE ANSWER IS RETURNED IN ACCUMULATOR AC

	SEARCH	IMPPRM
	TWOSEG	400000
	ENTRY	ATAN

	A=	0
	B=	1
	C=	2
	D=	3

ATAN:	MOVE	A, ARG1	;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	AT5		;ATAN(X) = X
	MOVEM	D, D1		;SAVE ACCUMULATOR D
	HLLO	D, A		;SAVE SIGN, SET RH(D) = -1
	CAML	B, A2		;IF A>2^33, THEN RETURN WITH
	JRST	AT4		;ATAN(X) = PI/2
	MOVEM	C, C1		;SAVE ACCUMULATOR C
	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, C3		;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, C3		;MULTIPLY BY ORIGINAL ARGUMENT
	TRNE	D, -1		;CHECK .G. 1.0 INDICATOR
	FSB	A, PIOT		;ATAN(A) = -(ATAN(1/A)-PI/2)
	SKIPA	C, C1		;RESTORE ACCUMULATOR C AND SKIP
AT4:	MOVE	A, PIOT		;GET PI/2 AS ANSWER
	SKIPGE	D		;LH(D) = -SGN(B) IF B>1.0
	MOVNS	A		;NEGATE ANSWER
	IFN AC<MOVE AC,0>	;GET ANSWER IN AC IF IT ISN'T
	MOVE	D, D1		;RESTORE ACCUMULATOR
AT5:	RETURN			;ATAN TERURN

	RELOC	0		;DATA SEGMENT

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
C1:	0
C3:	0
D1:	0
PIOT:	201622077325		;PI/2
	PRGEND
TITLE	SINH	%EXTERNALREALFN SINH(%REAL X)

;FLOATING POINT SINGLE PRECISION HYPERBOLIC SINE FUNCTION.

;SINH IS CALCULATED AS FOLLOWS:
;	IF ABS(X)>88.029,
;		SINH(X)=(EXP[ABS(X)-LN(2)])*SIGN(X)
;	IF ABS(X)<=0.10,
;		SINH(X)=X+(X**3)/6+(X**5)/120
;	FOR ALL OTHER VALUES OF X,
;		SINH(X)=1/2[EXP(X)-1/EXP(X)]

;THE CALLING SEQUENCE IS:
;	MOVE	ARG1,X
;	PUSHJ	P,SINH

;THE ANSWER IS RETURNED IN AC.

	SEARCH	IMPPRM
	ENTRY	SINH
	EXTERN	EXP
	TWOSEG	400000

SINH:	MOVE	0,ARG1		;PICK UP THE ARG.
	MOVEM	2,SAVE2		;SAVE AC 2.
	MOVEM	0,ARGTMP	;SAVE ARG.
	MOVM	2,0		;GET MAGNITUDE OF ARG IN AC 2.
	CAMLE	2,EIGHT8	;IF ABS(X)>88.029,
	JRST	OV88		;THEN GO TO OV88.
	CAMG	2,ONE10T	;IF ABS(X)<=0.10,
	JRST	SERIES		;THEN GO TO SERIES.
	MOVE	ARG1,2		;PREPARE FOR EXP
	CALL	EXP		;CALCULATE EXP(ABS(X)).
	IFN AC<MOVE 0,AC>
	HRLZI	1,576400	;PUT -1.0 IN AC 1.
	FDVR	1,0		;CALC. -EXP(-ABS(X)).
	FADR	0,1		;CALC. EXP(ABS(X))-EXP(-ABS(X)).
	FDVRI	0,202400	;CALC. THIS/2.0
	SKIPGE	ARGTMP		;ANSWER IS POSITIVE.
	MOVNS	0,0		;ANSWER IS NEGATIVE.
	MOVE	2,SAVE2		;RESTORE AC 2.
	IFN AC<MOVE AC,0>	;COND ARG SELECT
	RETURN			;SINH RETURN

SERIES:	FMPR	2,2		;CALC. X^2.
	JFCL			;SUPPRESS ERROR MESSAGE FROM OVTRAP.
	MOVEM	2,SX2		;SAVE X^2 IN SX2.
	FDVR	2,ONE120	;CALC.X^2/120
	JFCL			;SUPPRESS ERROR MESSAGE FROM OVTRAP.
	FADR	2,ONESIX	;CALC. (X^2/120)+1/6
	FMPR	2,SX2		;MULTIPLY IT BY X^2.
	JFCL			;SUPPRESS ERROR MESSAGE FROM OVTRAP.
	FADRI	2,201400	;ADD 1.0.
	FMPR	0,2		;MULTIPLY BY X.
	MOVE	2,SAVE2		;RESTORE AC 2.
	IFN AC<MOVE AC,0>	;COND RESULT SELECT
	RETURN			;SINH RETURN

OV88:	FSBR	2,LN2BE		;CALC.ABS(X)-LN(2)
	CAMG	2,EIGHT8	;OVERFLOW?
	JRST	EXPP		;NO,GO TO CALC.
	ERROR(1,1,0,IMPSTR<ILLEGAL SINH ARG>)

EXPP:	MOVE	ARG1,2			;GET READY FOR EXP
	CALL	EXP			;CALCULATE EXP
	IFN AC<MOVE 0,AC>
EXPP1:	SKIPGE	ARGTMP		;RETURN ANS. >0 IF X>0.
	MOVNS	0		;O'E, ANS. <0.
	MOVE	2,SAVE2		;RESTORE AC 2.
	IFN AC <MOVE AC,0>	;RESULT IN AC IF NOT THERE
	RETURN			;AND RETURN

	RELOC	0
SAVE2:	0
LN2BE:	200542710300		;LN(2)
EIGHT8:	207540074636		;88.029
ARGTMP:	0
ONE10T:	0.10
SX2:	0
ONE120:	207740000000		;120.0
ONESIX:	0.16666667
	PRGEND
TITLE	COSH	%EXTERNALREALFN COSH(%REAL X)

;FLOATING POINT SINGLE PRECISION HYPERBOLIC COSINE FUNCTION.

;COSH(X) IS CALCULATED AS FOLLOWS:
;	IF ABS(X) <= 88.029,
;		COSH(X) = 1/2(EXP(X) + 1.0/EXP(X))
;	IF ABS(X) > 88.029 AND (ABS(X)-LN(2)) <= 88.029,
;		COSH(X) = EXP(ABS(X)-LN(2))

;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
;	MOVE	ARG1,X
;	PUSHJ	P,COSH
;THE ANSWER IS RETURNED IN AC.

	SEARCH	IMPPRM
	ENTRY	COSH
	EXTERN	EXP
	TWOSEG	400000

COSH:	MOVE	0,ARG1		;PICK UP THE ARGUMENT.
	MOVEM	2,SAVE2		;SAVE AC 2.
	MOVM	2,0		;PUT ABS(X) IN AC 2.
	CAMLE	2,EIGHT8	;IF ABS(X) > 88.029,
	JRST	OV88		;GO TO OV88.
	MOVE	ARG1,2
	CALL	EXP
	IFN AC<MOVE 0,AC>	;CALCULATE EXP(ABS(X))
	MOVSI	2,201400	;PUT 1.0 IN AC 2.
	FDVR	2,0		;CALC. 1.0/EXP(ABS(X)).
	FADR	0,2		;CALC. EXP(ABS(X)) + EXP(-ABS(X)).
	FDVRI	0,202400	;DIVIDE THIS BY 2.0.
	MOVE	2,SAVE2		;RESTORE AC 2.
	IFN AC<MOVE AC,0>
	RETURN			;COSH RETURN

OV88:	FSBR	2,LN2BE		;FORM ABS(X)-LN(2).
	CAMG	2,EIGHT8	;OVERFLOW?
	JRST	EXPP		;NO,GO AHEAD.
	ERROR(1,1,0,IMPSTR<ARG TOO LARGE FOR COSH>)

EXPP:	MOVE	ARG1,2
	CALL	EXP
	MOVE	2,SAVE2		;RESTORE 
AC 2.
	RETURN			;COSH RETURN

	RELOC	0
SAVE2:	0
EIGHT8:	207540074636	;88.029
LN2BE:	200542710300	;LOG(2) BASE E.

	PRGEND
TITLE	TANH	%EXTERNALREALFN TANH(%REAL X)

;FLOATING POINT SINGLE PRECISION HYPERBOLIC TANGENT ROUTINE

;THIS ROUTINE CALCULATES THE TANH BY THE FOLLOWING ALGORITHM:
;IF ABSF(X) <.00034, THEN TANH(X) = X
;IF ABSF(X) >12.0, THEN TANH(X) = 1.0*SIGN(X)
;IF 0.17 <= X < 12.0, THEN TANH IS C
ALCULATED AS
;	TANH(X) = 1.0 - 2(1.0 + EXP(2*X))**-1
;IF .00034 <= X < 0.17, THEN TANH IS CALCULATED AS
;TANH(X) = F(A+F^2(B+C(D+F^2)**-1))**-1
;WHERE X = 4*LOG(E)  (BASE 2)

;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
;	MOVE	ARG1,X
;	PUSHJ	P,TANH
;THE ANSWER IS RETURNED IN ACCUMULATOR AC

	SEARCH	IMPPRM
	ENTRY	TANH
	EXTERN	EXP
	TWOSEG	400000

	A=	0
	B=	1

TANH:	MOVE	A, ARG1	;PICK UP THE ARGUMENT
	MOVM	B, A		;GET ABSF(ARGUMENT)
	CAMGE	B, KT1		;RETURN TANH(X)=X IF 
	JRST	TH2		;ABSF(X) .LE. .00034
	CAMLE	B, KT2		;RETURN TANH(X) = 1.0 IF
	JRST	TH5		;ARGUMENT GREATER THAN 12.0
	CAMGE	B, KT3		;USE RATIONAL APPROXIMATION IF
	JRST	TH3		;ARGUMENT IS LESS THAN 0.17
	MOVEM	A,SAVEA		;SAVE ARG. 
	FMPRI	B,202400	;GET 2*ARG.
	MOVE	ARG1,B		;SET UP FOR EXP
	CALL	EXP
	IFN AC<MOVE 0,AC>	;COND RESULT SELECT
	MOVSI	B, 201400	;FORM 1.0
	FAD	A, B		;1 + EXP(2X)
	FDVM	B, A		;(1 + EXP(2X))**-1
	FMPRI	A,202400	;2*(1 + EXP(2X))**-1
	FSBRM	B, A		;1 - 2*(1 + EXP(2X))**-1
	SKIPGE	SAVEA		;SKIP AHEAD IF ARG WAS >=0.
	MOVNS	A		;OTHERWISE,NEGATE THE ANSWER.
TH2:	IFN AC<MOVE AC,0>
	RETURN			;TANH RETURN

TH3:	FMP	A, KT7		;FORM 4*X*LOG(E) BASE 2
	MOVEM	A, TM1		;SAVE IT IN TM1
	FMP	A, A		;SQUARE IT
	MOVEM	A, TM2		;SAVE IT
	FAD	A, KT4		;FORM F^2 + T4
	MOVE	B, KT5		;GET T5 IN ACCUMULATOR B
	FDV	B, A		;KT5/(F^2 + KT4)
	FAD	B, KT6		;KT6 + KT5/(F^2 + KT4)
	FMP	B, TM2		;MULTIPLY BY F^2
	FAD	B, KT7		;ADD T7 (4*LOG(E) BASE 2)
	MOVE	A, TM1		;GET F IN ACCUMULATOR A
TH5:	FDV	A, B		;DIVIDE F BY PARTIAL SUM
	IFN AC<MOVE AC,0>
	RETURN			;TANH RETURN

	RELOC	0
KT1:	165544410070		;0.00034
KT2:	204600000000		;12.0
KT3:	176534121727		;0.17
KT4:	211535527022		;349.6699888
KT5:	204704333567		;14.1384514018
KT6:	173433723376		;0.01732867951
KT7:	203561250731		;5.7707801636

TM1:	0
TM2:	0
SAVEA:	0
	PRGEND
TITLE	LGT	%EXTERNALREALFN LGT(%REAL X)

;SINGLE PRECISION LOG TO THE BASE 10 FUNCTION
;LOG10(X) IS CALCULATED AS LOG10(E)*LOGE(X)

;THE CALLING SEQUENCE IS:
;	MOVE	ARG1,X
;	PUSHJ	P,LGT
;THE ANSWER IS RETURNED IN ACCUMULATOR AC

	A=	0
	B=	1

	SEARCH	IMPPRM
	TWOSEG	400000
	ENTRY	LGT
	EXTERN	LOG

LGT:	JUMPE	ARG1,LZERO	;CHECK FOR ZERO ARG.
	CALL	LOG		;CALC LOG(X)
	FMPR	AC,LOG10A	;MULTIPLY IT BY LOG10(E).
	RETURN			;LGT RETURN

LZERO:	ERROR(1,1,0,IMPSTR<ZERO ARG FOR LGT>)

	RELOC	0
LOG10A:	177674557305

	PRGEND
TITLE	LOG	%EXTERNALREALFN LOG(%REAL X)

;FLOATING POINT SINGLE PRECISION LOGARITHM FUNCTION
;LOG(ABSF(X)) IS CALCULATED BY THE SUBROUTINE, AND AN

;FOR LOGE(X), THE ALGORITHM IS:
;	LOGE(X) = (I + LOG2(F))*LOGE(2)
;	WHERE X = (F/2)*2^(I+1), AND LOG2(F) IS GIVEN BY
;	LOG2(F) = C1*Z + C3*Z^3 + C5*Z^5 - 1/2
;	AND Z = (F-SQRT(2))/(F+SQRT(2))

;THE CALLING SEQUENCE FOR THE ROUTINE IS AS FOLLOWS:
;	MOVE	ARG1,X
;	PUSHJ	P,LOG
;THE ANSWER IS RETURNED IN ACCUMULATOR AC

	A=0
	B=1


	ENTRY	LOG
	TWOSEG	400000
	SEARCH	IMPPRM

LOG:	MOVE	A, ARG1		;GET ABSF(X)
	JUMPG	A,ALOGOK	;ARG IS GREATER THAN 0
	JUMPE	A, LZERO	;CHECK FOR ZERO ARGUMENT
	ERROR(1,1,0,IMPSTR<NEGATIVE ARGUMENT TO LOG>)

ALOGOK:	CAMN	A, ONE		;CHECK FOR 1.0 ARGUMENT
	JRST	ZERANS		;IT IS 1.0 RETURN ZERO ANS.
	ASHC	A, -33		;SEPARATE FRACTION FROM EXPONENT
	ADDI	A, 211000	;FLOAT THE EXPONENT AND MULT. BY 2
	MOVSM	A, LS		;NUMBER NOW IN CORRECT FL. FORMAT
	MOVSI	A, 567377	;SET UP -401.0 IN A
	FADM	A, LS		;SUBTRACT 401 FROM EXP.*2
	ASH	B, -10		;SHIFT FRACTION FOR FLOATING
	TLC	B, 200000	;FLOAT THE FRACTION PART
	FAD	B, L1		;B = B-SQRT(2.0)/2.0
	MOVE	A, B		;PUT RESULTS IN A
	FAD	A, L2		;A = A+SQRT(2.0)
	FDV	B, A		;B = B/A
	MOVEM	B, LZ		;STORE NEW VARIABLE IN LZ
	FMP	B, B		;CALCULATE Z^2
	MOVE	A, L3		;PICK UP FIRST CONSTANT
	FMP	A, B		;MULTIPLY BY Z^2
	FAD	A, L4		;ADD IN NEXT CONSTANT
	FMP	A, B		;MULTIPLY BY Z^2
	FAD	A, L5		;ADD IN NEXT CONSTANT
	FMP	A, LZ		;MULTIPLY BY Z
	FAD	A, LS		;ADD IN EXPONENT TO FORM LOG2(X)
	FMP	A, L7		;MULTIPLY TO FORM LOGE(X)
	IFN AC<MOVE AC,0>	;GET RESULT
	RETURN			;LOG RETURN

LZERO:	ERROR(1,1,0,IMPSTR<ZERO ARG TO LOG>)

ZERANS:	MOVEI	AC, 0		;MAKE ANSWER ZERO
	RETURN			;LOG RETURN


	RELOC	0
ONE:	201400000000
L1:	577225754146		;-0.707106781187
L2:	201552023632		;1.414213562374
L3:	200462532521		;0.5989786496
L4:	200754213604		;0.9614706323
L5:	202561251002		;2.8853912903
L7:	200542710300		;0.69314718056

LS:	0
LZ:	0

	PRGEND
TITLE	COS	%EXTERNALREALFN COS(%REAL X)

;FLOATING POINT SINGLE PRECISION COSINE FUNCTION
;COS CALLS SIN TO CALCULATE SIN(PI/2+X)

;THE ROUTINE IS CALLED AS FOLLOWS
;	MOVE	ARG1,X
;	PUSHJ	P,COS
;THE ANSWER IS RETURNED IN ACCUMULATOR AC

	SEARCH	IMPPRM
	ENTRY	COS
	EXTERN	SIN
	TWOSEG	400000

COS:	FADR	ARG1,PIOT	;ADD PI/2.
	CALL	SIN		;CALCULATE SIN(PI/2+X)
	RETURN			;COS RETURN
	RELOC 0
PIOT:	201622077325		;PI/2
	PRGEND

TITLE	SIN	%EXTERNALREALFN SIN(%REAL X)

;FLOATING POINT SINGLE PRECISION SINE FUNCTION


;THIS 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
;010 - 3RD QUADRANT
;011 - 4TH QUADRANT
;THE ALGORITHM USES A MODIFIED TAYLOR SERIES TO CALCULATE 
;THE SINE OF THE NORMALIZED ARGUMENT.

;THE ROUTINES ARE CALLED IN THE FOLLOWING MANNER:
;	MOVE	ARG1,X
;	PUSHJ	P,SIN
;THE RESULT IS RETURNED IN ACCUMULATOR AC

	SEARCH	IMPPRM
	TWOSEG	400000
	ENTRY	SIN

	A=0
	B=1
	C=2
	D=3

SIN:	MOVE	B,ARG1		;PICK UP THE ARG.
S1:	MOVEM	B,SX		;SAVE THE ARG.
	MOVMS	B		;GET ABS OF ARG.
	CAMG	B,SP2		;SIN(X)=X IF X<2^-9.
	JRST	S3A		;EXIT WITH ARG. IN A.
	MOVEM	C,SC		;SAVE AC C.
	FDV	B,PIOT		;DIVIDE X BY PI/2.
	CAMG	B,ONE		;IS X/(PI/2) < 1.0 ?
	JRST	S2		;YES,ARG IN 1ST QUADRANT ALREADY.
	MULI	B,400		;NO,SEPARATE FRACTION AND EXP.
	LSH	C,-202(B)	;GET X MODULO 2PI.
	TLZ	C,(1B0)		;SUPRESS ERROR MESSAGE FROM OVTRAP.
	MOVEI	B,200		;PREPARE FLOATING FRACTION.
	ROT	C,3		;SAVE THREE BITS TO DETERMINE QUADRANT.
	LSHC	B,33		;ARGUMENT NOW IN THE RANGE (-1,1).
	FAD	B,SP3		;NORMALIZE THE ARGUMENT.
	JUMPE	C,S2		;REDUCED TO 1ST QUAD IF BITS 000.
	TLCE	C,1000		;SUBTRACT 1.0 FROM ARG IF BITS ARE
	FSB	B,ONE		;001 OR 011.
	TLCE	C,3000		;CHECK FOR FIRST QUADRANT, 001.
	TLNN	C,3000		;CHECK FOR THIRD QUADRANT, 010.
	MOVNS	B		;001,010.
S2:	SKIPGE	SX		;CHECK SIGN OF ORIGINAL ARG.
	MOVNS	B		;SIN(-X)=-SIN(X).
	MOVEM	B,SX		;STORE REDUCED ARG.
	FMPR	B,B		;CALCULATE X^X
	MOVE	A,SC9		;GET 1ST 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,SX		;MULTIPLY BY X.
	SKIPA	C,SC		;RESTORE AC C.
S3A:	MOVE	A,SX		;ANSWER IN X.
	IFN AC<MOVE AC,0>	;RESULT IN AC
	RETURN			;SIN RETURN

	RELOC	0
SC3:	577265210372
SC5:	175506321276
SC7:	606315546346
SC9:	164475536722

SP2:	170000000000
SP3:	0
SX:	0
CD1:	90.0
SCD1:	206712273406
PIOT:	201622077325
SC:	0
ONE:	1.0
	PRGEND
TITLE	SQR	%EXTERNALREALFN SQR(%REAL X)

;FLOATING POINT SINGLE PRECISION SQUARE ROOT FUNCTION
;THE SQUARE ROOT OF THE ABSOLUTE VALUE OF THE ARGUMENT IS
;CALCULATED.  THE FIRST GUESS IS CALCULATED TO BE OPTIMUM
;FOR NUMBERS BETWEEN 1/2 AND 2
;FOLLOWED BY TWO ITERATIONS OF NEWTON'S METHOD.

;THE CALLING SEQUENCE FOR THE SQUARE ROOT IS AS FOLLOWS:
;	MOVE	ARG1,X
;	PUSHJ	P,SQR
;THE ANSWER IS RETURNED IN ACCUMULATOR AC.

	SEARCH	IMPPRM
	TWOSEG	400000
	ENTRY	SQR

	A=0
	B=1

SQR:	SKIPG	B,ARG1		;PICK UP ARG. CHECK IF > 0
	JRST	SQRTLE		;NO, HANDLE NON-POSITIVE ARGUMENT

SQRTP:	MOVE	A,B		;SAVE NUMBER
	LSH	B,-1		;DIVIDE EXP BY 2
	TLZE	B,400		;WAS EXPONENT ODD?
	JRST	SQRT2		;YES

;HERE WHEN EXPONENT WAS EVEN.  B CONTAINS AN UNNORMALIZED FLOATING
; POINT NUMBER, THE FRACTION PART OF WHICH IS 1/2 THE FRACTION OF
; THE ARGUMENT.  OUR INITIAL GUESS IS MADE BY A LINEAR APPROXIMATION
; USING Y0 = SE (X + C), WHERE SE AND C ARE CONSTANTS USED FOR
; EVEN EXPONENTS IN X.

	ADD	B,CONST1		;COMPUTE LINEAR APPROXIMATION
	FMPRI	B,301454	;RESCALE EXPONENT
	JRST	SQRT3

;HERE WITH ODD EXPONENT, USE Y0 = SO * (X+C).

SQRT2:	ADD	B,CONST1		;LINEAR APPROXIMATION
	FMPRI	B,301650		;RESCALE EXPONENT

SQRT3:	FDV	A,B		;ORIGINAL / INITIAL GUESS
	FAD	B,A		;AVERAGE THEM
	FSC	B,-1

	MOVM	A,ARG1		;GET ORIGINAL NUMBER
	FDV	A,B		;SECOND ITERATION
	FADR	A,B
	FSC	A,-1		;AVERAGE THIRD GUESS WITH SECOND

	IFN AC<MOVE AC,0>	;GET RESULT IN AC IF IT ISN'T
	RETURN			;SQR RETURN

SQRTLE:	JUMPE	B,ZERO
	ERROR(6,2,0,IMPSTR<ILLEGAL ARG TO SQR>)

ZERO:	MOVEI	AC,0		;HERE ON NON-POSITIVE ARG. RETURN ZERO
	RETURN			;SQRT RETURN


	RELOC	0		;DATA SEGMENT

CONST1:	XWD	267,607000		;CONSTANT FOR 1ST APPROX
	PRGEND
TITLE	EXP	%EXTERNALREALFN EXP(%REAL X)

;FLOATING POINT SINGLE PRECISION EXPONENTIAL FUNCTION
;IF X<=-89.415..., THE PROGRAM RETURNS ZERO AS THE ANSWER
;IF X>= 88.029   AN OVERFLOW EVENT IS SIGNALLED
;THE RANGE OF THE ARGUMENT IS REDUCED AS FOLLOWS:
;EXP(X) = 2**(X*LOG(E)BASE2) = 2**(M+F)
;WHERE M IS AN INTEGER AND F IS A FRACTION
;2**M IS CALCULATED BY ALGEBRAICALLY ADDING M TO THE EXPONENT
;OF THE RESULT OF 2**F. 2**F IS CALCULATED AS

;2**F = 2(0.5+F(A+B*F^2 - F-C(F^2 + D)**-1)**-1

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
;	MOVE	ARG1,X
;	PUSHJ	P,EXP
;THE ANSWER IS RETURNED IN ACCUMULATOR AC

	A=	0
	B=	1
	C=	2
	D=	3




	SEARCH	IMPPRM
	ENTRY	EXP
	TWOSEG	400000


EXP:	MOVE	B,ARG1		;PICK UP THE ARGUMENT IN B
	CAMGE	B,E77		;IS EXP. < -89.41...?
	JRST	OUT2		;YES, GO TO EXIT.
	CAMG	B,E7		;IS EXP. > +88.029...?
	JRST	EXP1		;GO TO STANDARD ALGORITHM.
	HALT			;!*** TEMP ***!
	HRLOI	A, 377777	;GET LARGEST FLOATING NUMBER
	RETURN		;RETURN

OUT2:	MOVEI	AC,0		;EXPONENT <-88 IMPLIES RESULT=0
	RETURN			;EXP RETURN

EXP1:	MOVEM	C, ES1		;SAVE ACCUMULATOR C
	MOVEM	D, ES3		;SAVE ACCUMULATOR D
	SETZM	ES2		;INITIALIZE ES2
	MULI	B, 400		;SEPARATE FRACTION AND EXPONENT
	TSC	B, B		;GET A POSITIVE EXPONENT
	MUL	C, E5		;FIXED POINT MULTIPLY BY LOG2(E)
	ASHC	C, -242(B)	;SEPARATE FRACTION AND INTEGER
	AOSG	C		;ALGORITHM CALLS FOR MULT. BY 2
	AOS	C		;ADJUST IF FRACTION WAS NEGATIVE
	HRRM	C, TEMP		;SAVE FOR FUTURE SCALING
	JUMPG	D,ASHH		;GO AHEAD IF ARG > 0.
	TRNN	D,377		;ARE ALL THESE BITS 0?
	JRST	ASHH		;YES, GO AHEAD.
	ADDI	D,200		;NO, FIX UP.
ASHH:	ASH	D, -10		;MAKE ROOM FOR EXPONENT
	TLC	D, 200000	;PUT 200 IN EXPONENT BITS
	FADB	D, ES2		;NORMALIZE, RESULTS TO D AND ES2
	FMP	D, D		;FORM X^2
	MOVE	A, E2		;GET FIRST CONSTANT
	FMP	A, D		;E2*X^2 IN A
	FAD	D, E4		;ADD E4 TO RESULTS IN D
	MOVE	B, E3		;PICK UP E3
	FDV	B, D		;CALCULATE E3/(F^2 + E4)
	FSB	A, B		;E2*F^2-E3(F^2 + E4)**-1
	MOVE	C, ES2		;GET F AGAIN
	FSB	A, C		;SUBTRACT FROM PARTIAL SUM
	FAD	A, E1		;ADD IN E1
	FDVM	C, A		;DIVIDE BY F
	FAD	A, E6		;ADD 0.5
EX1:	MOVE	B,TEMP		;GET THE SCALE FACTOR
	FSC	A,(B)		;SCALE THE RESULTS
	MOVE	C, ES1		;RESTORE ACCUMULATOR C
	MOVE	D, ES3		;RESTORE ACCUMULATOR D
	IFN AC<MOVE AC,A>	;GET THE ANSWER IF ITS NOT THERE ALREADY
	RETURN			;AND RETURN

	RELOC	0
E1:	204476430062		;9.95459578
E2:	174433723400		;0.03465735903
E3:	212464770715		;617.97226953
E4:	207535527022		;87.417497202
E5:	270524354513		;LOG(E), BASE 2
E6:	0.5
E7:	207540074636		;88.029...
E77:	570232254037		;-89.415986
ES1:	0
ES2:	0
ES3:	0
TEMP:	0

	PRGEND
	END