perm filename TRIGS[X,AIL] blob
sn#123322 filedate 1974-11-26 generic text, type T, neo UTF8
COMMENT ⊗ VALID 00014 PAGES VERSION 17-1(3)
RECORD PAGE DESCRIPTION
00001 00001
00003 00002 HISTORY
00004 00003 NUMERICAL ROUTINES FOR SAIL
00005 00004 ROUTINES FOR HANDLING UNDER/OVER FLOW.
00011 00005 FLOATING POINT SINGLE PRECISION SINE AND COSINE FUNCTION
00017 00006 FLOATING POINT SINGLE PRECISION SQUARE ROOT FUNCTION
00021 00007 PSEUDO RANDOM NUMBER GENERATOR AND INITIALIZING ROUTINE
00024 00008 FLOATING POINT SINGLE PRECISION ARCTANGENT FUNCTION
00029 00009 FLOATING POINT SINGLE PRECISION ARCSINE FUNCTION
00032 00010 FLOATING POINT SINGLE PRECISION ARCCOSINE FUNCTION
00035 00011 FLOATING POINT SINGLE PRECISION ARCTANGENT OF TWO ARGUMENTS
00039 00012 FLOATING POINT SINGLE PRECISION HYPERBOLIC TANGENT ROUTINE
00043 00013 FLOATING POINT SINGLE PRECISION HYPERBOLIC COSINE FUNCTION.
00046 00014 FLOATING POINT SINGLE PRECISION HYPERBOLIC SINE FUNCTION.
00050 ENDMK
⊗;
COMMENT ⊗HISTORY
AUTHOR,REASON
021 102100000003 ⊗;
COMMENT ⊗
VERSION 17-1(3) 12-12-73 BY RHT & RFS BUG #PW# NEED SEVERAL APR BIT ENABLINGS
VERSION 17-1(2) 12-12-73
VERSION 17-1(1) 12-11-73 BY JRL BUG #PV# KEEP STACK HEIGHT CONSISTENT WITHIN ATAN
⊗;
SUBTTL NUMERICAL ROUTINES FOR SAIL
COMMENT ⊗
This set of routines requires that some pieces be in the low
segment. We have decided (for this and other reasons) that
it shall not be part of the standard SAIL upper segment.
As a result, the switches UP, LOW, NOUP and NOLOW do not appear
in here; this file is just not included in the second segment
assembly.
When the file is used to make a library, the parameters to the
COMPIL macro (and RENSW, set before all) determine where the code
will lie.
⊗
;ROUTINES FOR HANDLING UNDER/OVER FLOW.
NOUP <
COMPIL(TGI,<TRIGINI>,<JOBTPC,JOBAPR,OVPCWD>
,<TRIG ROUTINE INTERRUPT HANDLER>,<.RSEED>,INHIBIT)
BEGIN UNDER
OV←←400000
FOV←←40000
ZDV←←40
FXU←←100
IFE ALWAYS,< ;MORE EXTERNALS, CONDITIONALLY ASSEMBLED.
EXPO <
EXTERNAL INTMAP,ENABLE
>;EXPO
EXTERNAL APRACS,GOGTAB
NOEXPO <
EXTERNAL INTTBL
>;NOEXPO
>
;The underflow processing is turned
;on with the call TRIGINIT( <rout-name> ). This uses the
;INTMAP and ENABLE SAIL functions to set up interrupt
;traps (export); at Stanford, it turns on the APRENB system
;itself.
;If an interrupts happens, and it is
;not followed by a JFCL (indicating that overflow from this
;spot is ok), then the <rout-name> is called.
;It may peer around at things, change JOBTPC, etc., just
;as documented for any interrupt routine. When it returns,
;the interrupt is dismissed.
;IF THIS CODE IS LOADED, INITIALIZE IT!
TRGIN: 0
XWD 000000,TRGIQ ;SPECIAL INIT IN FIRST SYS PHASE
0 ;ONLY 1 ROUTINE TO CALL.
LINK %INLNK,TRGIN
TRGIQ: PUSH P,[0] ;NO ROUTINE.
PUSHJ P,TRIGINI ;INITIALIZATION
POPJ P, ;RETURN.
HERE(TRIGINIT)
JFCL 17,.+1 ;CLEAR NUMERIC FLAGS.
EXPO <
PUSH P,[=32] ;ARM FOR FLOATING OVERFLOWS
PUSH P,[FLTOV] ;ROUTINE.
PUSH P,[0] ;NOT DEFERRED
PUSHJ P,INTMAP ;SET IT ALL UP.
PUSH P,[=32]
PUSHJ P,ENABLE ;AND ENABLE IT.
;;#PW# RHT SINCE APR DISPATCHER ASSUMES ONLY 1 BIT & DEC MAY SET SEVERAL
PUSH P,[=29] ;ARM FOR regular, too.
PUSH P,[FLTOV] ;ROUTINE.
PUSH P,[0] ;NOT DEFERRED
PUSHJ P,INTMAP ;SET IT ALL UP.
PUSH P,[=29]
PUSHJ P,ENABLE ;AND ENABLE IT.
;;#PW#
>;EXPO
NOEXPO <
MOVE USER,GOGTAB ;LOOK TO SEE IF SET UP FOR
SKIPE DISPAT(USER) ;INTERRUPTS YET; IF NOT,
JRST .+3
PUSH P,[=128] ;ENABLE THEM NOW.
PUSHJ P,INTTBL ;...
MOVEI A,10 ;TURN ON ALL OVERFLOWS.
CALL6 (A,APRENB)
MOVEI A,FLTOV
MOVEM A,JOBAPR ;...
>;NOEXPO
POP P,1 ;RETURN ADDRESS
POP P,OVROUT ;USER'S ROUTINE.
JRST (1) ;RETURN.
FLTOV: ;HERE WHEN AN INTERRUPT HAPPENS
NOEXPO <
MOVEM 1,SAVP ;GET AN ACCUMULATOR.
>;NOEXPO
MOVE 1,JOBTPC ;COME HERE WITH AC'S SET UP.
MOVEM 1,OVPCWD ;SAVE FOR LOOKING.
TLNN 1,FXU ;WHAT KIND OF INTERRUPT
JRST NOFX ;NOT FLOATING UNDERFLOW.
MOVE 1,-1(1) ;GET OPCODE.
TLNN 1,40000 ;CHECK FOR FSC
TLZ 1,2000
DPB 1,[POINT 29,SEW,35] ;CHANGE INSTRUCTION
EXPO < ;GET BACK INTRRUPT AC'S
PUSH P,16
MOVEM P,SAVP
MOVE 17,[XWD APRACS,0]
BLT 17,17 ;RESTORE AC'S
>;EXPO
NOEXPO <
MOVE 1,SAVP ;GET BACK SAVED AC.
>;NOEXPO
SEW: SETZ 0,0 ;MODIFIED!!!!
NOEXPO <
MOVEM 1,SAVP ;SAVE AGAIN.
>;NOEXPO
EXPO <
MOVEM 17,APRACS+17
MOVEI 17,APRACS
BLT 17,APRACS+16
MOVE P,SAVP
POP P,16
>;EXPO
NOFX: HLRZ 1,@JOBTPC ;CHECK IF NEXT INSTR IS JFCL
ANDCMI 1,777
CAIE 1,(<JFCL>) ;CHECK.
JRST USRRT ;GO TO USER ROUTINE.
HRRZ 1,@JOBTPC ;GET EFFECTIVE ADDRESS
HRRM 1,JOBTPC
RET: MOVSI 1,OV+FOV+FXU+ZDV
ANDCAM 1,JOBTPC
NOEXPO <
MOVE 1,SAVP ;RESTORE THE ACCUMULATOR
JRST 2,@JOBTPC
>;NOEXPO
EXPO <
POPJ P, ;RETURN TO INTERRUPT HANDLER.
>;EXPO
USRRT: SKIPN OVROUT
JRST RET ;NO USER ROUTINE.
NOEXPO <
MOVEI 1,APRACS ;SAVE AC'S HERE FOR USER TO SEE
BLT 1,APRACS+17
MOVE 1,SAVP
MOVEM 1,APRACS+1
MOVE USER,GOGTAB
MOVE P,IPDP(USER)
MOVE SP,ISPDP(USER)
>;NOEXPO
PUSHJ P,@OVROUT ;CALL USER'S ROUTINE.
NOEXPO <
MOVE 1,APRACS+1 ;RESTORE AC'S
MOVEM 1,SAVP
MOVSI 1,APRACS
BLT 1,17 ;AND RESTORE ALL OTHERS
>;NOEXPO
JRST RET
SAVP: 0
OVROUT: 0
BEND UNDER
.RSEED: =524287
ENDCOM (TGI)
>;NOUP
NOLOW < ;REST IS REENTRANT
COMPIL(TG1,<SIN$,COS$,SIND$,COSD$,SQRT$,RAN$>,<TRIGINI,X22,.RSEED>,<SIN, SQRT ROUTINES>)
;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.
;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:
; PUSH P,ARG
; PUSHJ P,SIN$ (OR COS$,SIND$, OR COSD$)
;THE ANSWER IS RETURNED IN ACCUMULATOR 1
BEGIN SIN$
INTEGER SAVEAA,SAVEBB,SAVECC
AA←13
BB←14
CC←15
HERE(COSD$) ;ENTRY TO COSINE DEGREES ROUTINE.
MOVEM BB,SAVEBB ;SAVE ACCUMULATORS
MOVE BB,-1(P) ;PICK UP THE ARG.
FADR BB,CD1 ;ADD 90 DEGREES.
JRST CONVER ;CONVERT TO RADIANS
;THEN GO TO SIN ROUTINE
HERE(SIND$) ;ENTRY TO SINE DEGREES ROUTINE.
MOVEM BB,SAVEBB ;SAVE ACCUMULATORS
MOVE BB,-1(P) ;PICK UP THE ARG.
CONVER: FDVR BB,SCD1 ;CONVERT TO RADIANS
JFOV .+1 ;SUPPRESS ERROR MESSAGE ON UNDERFLOW.
;SPECIAL INTERRUPT CODE WILL SET
; BB TO 0 ON UNDERFLOW
JRST S1 ;ENTER SINE ROUTINE.
HERE(COS$) ;ENTRY TO COSINE RADIANS ROUTINE.
MOVEM BB,SAVEBB ;SAVE ACCUMULATORS
MOVE BB,-1(P) ;PICK UP THE ARG.
FADR BB,PIOT ;ADD PI/2.
JRST S1 ;ENTER SINE ROUTINE.
HERE(SIN$) ;ENTRY TO SINE RADIANS ROUTINE.
MOVEM BB,SAVEBB ;SAVE ACCUMULATORS
MOVE BB,-1(P) ;PICK UP THE ARG.
S1: MOVEM BB,-1(P) ;SAVE THE ARG.
MOVEM AA,SAVEAA ;SAVE ACCUMULATORS
MOVEM CC,SAVECC
MOVMS BB ;GET ABS OF ARG.
CAMG BB,SP2 ;SIN(X)=X IF X LESS THAN 2↑-9.
JRST S3A ;EXIT WITH ARG. IN B.
FDV BB,PIOT ;DIVIDE X BY PI/2.
CAMG BB,ONE ;IS X/(PI/2) LESS THAN 1.0 ?
JRST S2 ;YES,ARG IN 1ST QUADRANT ALREADY.
MULI BB,400 ;NO,SEPARATE FRACTION AND EXP.
ASH CC,-202(BB) ;GET X MODULO 2PI.
JFOV .+1 ;SUPRESS ERROR MESSAGE FROM OVTRAP.
;SPECIAL INTERRUPT CODE WILL
;RETURN WITHOUT ATTEMPTING A
;FIXUP
MOVEI BB,200 ;PREPARE FLOATING FRACTION.
ROT CC,3 ;SAVE THREE BITS TO DETERMINE QUADRANT.
LSHC BB,33 ;ARGUMENT NOW IN THE RANGE (-1,1).
FAD BB,SP3 ;NORMALIZE THE ARGUMENT.
JUMPE CC,S2 ;REDUCED TO 1ST QUAD IF BITS 000.
TLCE CC,1000 ;SUBTRACT 1.0 FROM ARG IF BITS ARE
FSB BB,ONE ;001 OR 011.
TLCE CC,3000 ;CHECK FOR FIRST QUADRANT, 001.
TLNN CC,3000 ;CHECK FOR THIRD QUADRANT, 010.
MOVNS BB ;001,010.
S2: SKIPGE -1(P) ;CHECK SIGN OF ORIGINAL ARG.
MOVNS BB ;SIN(-X)=-SIN(X).
MOVEM BB,-1(P) ;STORE REDUCED ARG.
FMPR BB,BB ;CALCULATE X↑X
MOVE AA,SC9 ;GET 1ST CONSTANT.
FMP AA,BB ;MULTIPLY BY X↑2
FAD AA,SC7 ;ADD IN NEXT CONSTANT.
FMP AA,BB ;MULTIPLY BY X↑2.
FAD AA,SC5 ;ADD IN NEXT CONSTANT.
FMP AA,BB ;MULTIPLY BY X↑2.
FAD AA,SC3 ;ADD IN NEXT CONSTANT.
FMP AA,BB ;MULTIPLY BY X↑2.
FAD AA,PIOT ;ADD IN LAST CONSTANT.
S2B: FMPR AA,-1(P) ;MULTIPLY BY X.
SKIPA 1,AA ;ANSWER IN 1
S3A: MOVE 1,-1(P) ;ANSWER IN 1.
MOVE AA,SAVEAA
MOVE BB,SAVEBB
MOVE CC,SAVECC
SUB P,X22
JRST @2(P)
SC3: 577265210372
SC5: 175506321276
SC7: 606315546346
SC9: 164475536722
SP2: 170000000000
SP3: 0
CD1: 90.0
SCD1: 206712273406
PIOT: 201622077325
ONE: 1.0
BEND SIN$
;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 LESS THAN F LESS THAN 1
;SQRT(X) IS THEN CALCULATED AS (SQRT(F))*(2**B)
;SQRT(F) IS CALCULATED BY A LINEAR APPROXIMATION, THE NATURE
;OF WHICH DEPENDS ON WHETHER 1/4 LESS THAN F LESS THAN 1/2 OR 1/2 LESS THAN F LESS THAN 1,
;FOLLOWED BY TWO ITERATIONS OF NEWTON'S METHOD.
;THE CALLING SEQUENCE FOR THE SQUARE ROOT IS AS FOLLOWS:
; PUSH 17,ARG
; PUSHJ 17,$SQRT
;THE ANSWER IS RETURNED IN ACCUMULATOR 1.
BEGIN SQRT$
INTEGER SAVEAA,SAVEBB,SAVECC
AA←13
BB←14
F←15
HERE(SQRT$) ;ENTRY TO SQUARE ROOT ROUTINE
MOVEM AA,SAVEAA
MOVEM BB,SAVEBB
MOVEM F,SAVECC
SKIPG BB,-1(P) ;PICK UP ARG. CHECK IF GREATER THAN 0
JRST SQRT4 ;NO, HANDLE NON-POSITIVE ARGUMENT
MOVEI AA,0 ;GET EXPONENT TO AA
LSHC AA,=9
SUBI AA,201 ;GET TRUE EXPONENT + 1
ROT AA,-1 ;DIVIDE BY 2
;AA HAS SCAL FACTOR.
JUMPL AA,SQRT3 ;JUMP IF FRACTION GREATER THAN .5
;FRACTION LESS THAN .5
LSH BB,=-9 ;RESTORE POSITION OF FRACTION IN BB
FSC BB,177 ;AND FIX UP EXPONENT .25 LESS THAN F LESS THAN .5
MOVEM BB,F ;SAVE FRACTION
;COMPUTE LINEAR APPROX #1
FMPRI BB,200640
FADRI BB,177465
SQRT1: MOVE 1,F ;1ST ITERATION OF NEWTON
FDV 1,BB ; F/APPROX
FAD BB,1 ; APPROX + F/APPROX
FSC BB,-1 ; .5*( APPROX + F/APPROX)
MOVE 1,F ;2ND ITERATION OF NEWTON
FDV 1,BB ; F/APPROX
FADR 1,BB ; APPROX + F/APPROX
FSC 1,(AA) ;HALVE AND SCALE EXPONENT
EXIT: SUB P,X22
MOVE AA,SAVEAA
MOVE BB,SAVEBB
MOVE F,SAVECC
JRST @2(P)
;HERE ON F GREATER THAN= .5
SQRT3: LSH BB,=-9 ;RESTORE POSITION OF FRACTION IN BB
FSC BB,200 ;AND FIX UP EXPONENT .5 LESS THAN= F LESS THAN 1
MOVEM BB,F ;SAVE FRACTION
;COMPUTE LINEAR APPROX #2
FMPRI BB,200450
FADRI BB,177660
JRST SQRT1 ;NOW GO ITERATE
SQRT4: JUMPE BB,ZERO
ERR <SQRT: Negative argument - 0 returned>,1
ZERO: MOVEI 1,0 ;HERE ON NON-POSITIVE ARG. RETURN ZERO
JRST EXIT
BEND SQRT$
;PSEUDO RANDOM NUMBER GENERATOR AND INITIALIZING ROUTINE
;METHOD SUGGESTED BY D. H. LEHMER
;CALLING SEQUENCE FOR FUNCTION RAN:
;PUSH 17,ARG
;PUSHJ 17,$RAN
;IF ARG NEQ 0, ARG IS USED AS PREVIOUS RANDOM NO.
;(I.E. AS THE STARTING VALUE), OTHERWISE THE PREVIOUS
;VALUE (WHICH WAS STORED LOCALLY) IS USED.
;ANSWER IS RETURNED IN ACCUMULATOR 1 AS A SINGLE
;PRECISION FLOATING POINT NUMBER IN THE RANGE
;0 LESS THAN X LESS THAN 1
BEGIN RAN$
INTEGER SAVEAA,SAVEBB
AA←13
BB←14
INTERNAL RAN$
RAN$:
MOVEM AA,SAVEAA
MOVEM BB,SAVEBB
MOVE AA,-1(P) ;IF ARG = 0 THEN
JUMPE AA,R1 ;USE PREVIOUS RANDOM NO.
TLZ AA,760000 ;OTHERWISE MASK 5 BITS
MOVEM AA,.RSEED ;AND STORE NEW NO.
R1: MOVE AA,K ;GET K [14**29(MOD2**31 -1)]
MUL AA,.RSEED ;MULTIPLY WITH LAST RANDOM NUMBER
ASHC AA,4 ;SEPARATE RESULT IN TWO 31 BIT WORDS
LSH BB,-4
ADD AA,BB ;ADD THEM TOGETHER
TLZE AA,760000 ;SKIP IF RESULT LESS THAN 31 BITS
ADDI AA,1
MOVEM AA,.RSEED ;STORE NEW RN IN INTEGER MODE
HLRZ 1,AA ;CONVERT TO FP IN TWO STEPS IN
FSC 1,216 ;ORDER TO LOOSE NO LOW ORDER
HRLI AA,0 ;BITS
FSC AA,174
FAD 1,AA
MOVE AA,SAVEAA
MOVE BB,SAVEBB
SUB P,X22
JRST @2(P) ;EXIT
K: =630360016 ;14**29(MOD 2**31 -1)
;.RSEED:=524287 ;STARTING VALUE
BEND RAN$
ENDCOM(TG1)
;; #PV# ! (1 OF 3) INCLUDE X11 IN COMPIL
COMPIL(TG2,<ATAN$,ASIN$,ACOS$,ATAN2$>,<OVPCWD,TRIGINI,X33,SQRT$,X22,X11>
,<ARC-TRIG ROUTINES>)
;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 LESS THAN X LESS THAN= 1
;IF X GREATER THAN 1, THEN ATAN(X) = PI/2 - ATAN(1/X)
;AC DD IS USED INTERNALLY TO KEEP TRACK OF CASES
;IF X GREATER THAN 1, THEN RH(DD) =-1, AND LH(DD) = -SGN(X)
;IF X LESS THAN 1, THEN RH(DD) = 0, AND LH(DD) = SGN(X)
;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
; PUSH 17,ARG
; PUSHJ 17,$ATAN
;THE ANSWER IS RETURNED IN ACCUMULATOR 1
BEGIN ATAN$
INTEGER SAVEBB,SAVECC,SAVEDD
AA←1
BB←13
CC←14
DD←15
HERE(ATAN$) ;ENTRY TO ARCTANGENT ROUTINE
MOVEM BB,SAVEBB
MOVEM CC,SAVECC
MOVEM DD,SAVEDD
MOVE AA,-1(P) ;PICK UP THE ARGUMENT IN AA
ATAN1: MOVM BB, AA ;GET ABSF OF ARGUMENT
CAMG BB, A1 ;IF X LESS THAN 2↑-33, THEN RETURN WITH...
JRST EXIT ;ATAN(X) = X
HLLO DD, AA ;SAVE SIGN, SET RH(DD) = -1
CAML BB, A2 ;IF AA GREATER THAN 2↑33, THEN RETURN WITH
JRST AT4 ;ATAN(X) = PI/2
MOVSI CC, (1.0) ;FORM 1.0 IN CC
CAMG BB, CC ;IS ABSF(X) GREATER THAN 1.0?
TRZA DD, -1 ;IF BB .LE. 1.0, THEN RH(DD) = 0
FDVM CC, BB ;BB IS REPLACED BY 1.0/BB
TLC DD, (DD) ;XOR SIGN WITH .G. 1.0 INDICATOR
PUSH P,BB ;SAVE THE ARGUMENT
FMP BB, BB ;GET BB↑2
MOVE CC, KB3 ;PICK UP A CONSTANT
FAD CC, BB ;ADD BB↑2
MOVE AA, KA3 ;ADD IN NEXT CONSTANT
FDVM AA, CC ;FORM -A3/(B↑2 + B3)
FAD CC, BB ;ADD BB↑2 TO PARTIAL SUM
FAD CC, KB2 ;ADD B2 TO PARTIAL SUM
MOVE AA, KA2 ;PICK UP -A2
FDVM AA, CC ;DIVIDE PARTIAL SUM BY -A2
FAD CC, BB ;ADD BB↑2 TO PARTIAL SUM
FAD CC, KB1 ;ADD B1 TO PARTIAL SUM
MOVE AA, KA1 ;PICK UP A1
FDV AA, CC ;DIVIDE PARTIAL SUM BY A1
FAD AA, KB0 ;ADD B0
FMP AA,(P) ;MULTIPLY BY ORIGINAL ARGUMENT
TRNE DD, -1 ;CHECK .G. 1.0 INDICATOR
FSB AA, PIOT ;ATAN(AA) = -(ATAN(1/AA)-PI/2)
;; #PV# ! (2 OF 3) JRL BOTH PATHS TO AT4 SHOULD LEAVE STACK SAME HEIGHT
SUB P,X11
SKIPA 0,0 ;SKIP
AT4: MOVE AA, PIOT ;GET PI/2 AS ANSWER
SKIPGE DD ;LH(DD) = -SGN(BB) IF BB GREATER THAN 1.0
MOVNS AA ;NEGATE ANSWER
;; #PV# ! (3 OF 3) USED TO BE SUB P,X33 BUT NOW STACK DECREMENTED EARLIER
EXIT: MOVE BB,SAVEBB
MOVE CC,SAVECC
MOVE DD,SAVEDD
SUB P,X22
JRST @2(P)
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$
;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 ERROR MESSAGE TO BE
;TYPED AND AN ANSWER OF ZERO TO BE RETURNED.
;CALLING SEQUENCE:
; PUSH P,ARG
; PUSHJ P,$ASIN
;
;RESULT RETURNED IN AC 1.
;
;
BEGIN ASIN$
INTEGER SAVEAA
AA←13
BB←1
HERE(ASIN$) ;ENTRY TO ASIN ROUTINE
MOVE AA,SAVEAA
MOVM BB,-1(P) ;GET MAGNITUDE OF ARG. IN BB
CAMLE BB,ONE ;IS THE MAGNITUDE OF THE ARG. LE 1.0?
JRST TOOLRG ;NO, GO TO ERROR RETURN.
MOVN AA,-1(P) ;GET THE NEGATIVE OF ARG
FMP AA,-1(P) ;CALCULATE -(X↑2)
JFOV .+1 ;SUPPRESS ERROR MESSAGE FROM OVTRAP
;ON UNDERFLOW, THE SPECIAL INTERRUPT
;CODE SETS AA TO 0
FAD AA, ONE ;CALCULATE 1-(X↑2)
JUMPE AA, ASIN1 ;WAS X EITHER -1.0 OR 1.0?
PUSH P,AA ;NO,
PUSHJ P,SQRT$ ;CALCULATE SQRT(1-X↑2)
MOVE AA,-1(P) ;GET THE ARGUMENT BACK AGAIN
FDV AA,BB ;CALCULATE X/SQRT(1-X↑2)
PUSH P,AA ;THEN
PUSHJ P,ATAN$ ;CALCULATE ATAN(X/SQRT(1-X↑2)),.
EXIT: MOVE AA,SAVEAA
SUB P,X22
JRST @2(P)
TOOLRG: ERR <ASIN: Argument mangitude greater than 1.0; 0 returned>,1
TDZA BB,BB
ASIN1: MOVE BB, PIOT ;ANSWER IS EITHER PI/2 OR-PI/2
SKIPG -1(P) ;WAS ORIGINAL ARGUMENT POSITIVE?
MOVNS BB ;NO, GET -PI/2
JRST EXIT
PIOT: 201622077325 ;PI/2
ONE: 1.0
BEND ASIN$
;FLOATING POINT SINGLE PRECISION ARCCOSINE FUNCTION
;ACOS(X) IS CALCULATED IN THE FOLLOWING MANNER:
; IF X GREATER THAN 0, ACOS(X)=ATAN((SQRT(1-X↑2))/X)
; IF X LESS THAN 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
;TO BE TYPED AND WILL RETURN AN ANSWER OF ZERO.
;THE CALLING SEQUENCE FOR ACOS IS:
; PUSH 17,ARG
; PUSHJ 17,$ACOS
;THE RESULT IS RETURNED IN AC 1
BEGIN ACOS$
HERE(ACOS$) ;ENTRY TO ACOS ROUTINE.
MOVM 1,-1(P) ;GET /ARG./ IN AC 1.
CAMLE 1,ONE ;IS MAGNITUDE OF ARG. GT 1.0?
JRST TOOLRG ;YES, GO TO ERROR RETURN.
JUMPE 1,ZERARG ;IF ARG=0, GO TO ZERARG.
FMPR 1,1 ;X↑2 IN AC 1.
JFOV .+1 ;SUPPRESS ERROR MESSAGE FROM OVTRAP
;ON UNDERFLOW THE SPECIAL
;INTERRUPT CODE WILL SET
;AC 1 TO 0
MOVNS 1 ;-X↑2 IN AC 1.
FAD 1,ONE ;1.0-X↑2 IN AC 1.
PUSH P,1
PUSHJ P,SQRT$ ;CALC. $SQRT(1.0-X↑2).
FDVR 1,-1(P) ;($SQRT(1.0-X↑2))/X IS IN AC 1.
JFOV .+1 ;SUPPRESS ERROR MESSAGE FROM OVTRAP
;ON UNDERFLOW THE SPECIAL INTERRUPT
;CODE WILL SET AC 1 TO 0
;ON OVERFLOW AC 1 WILL BE SET
;TO LARGEST MAGNITUDE WITH
;CORRECT SIGN
PUSH P,1
PUSHJ P,ATAN$ ;FIND $ATAN($SQRT(1.0-X↑2)/X).
SKIPL -1(P) ;SKIP IF ORIGINAL ARG LESS THAN 0.
JRST EXIT ;RETURN.
FAD 1,PII ;ANSWER IS PI + ANSWER IN AC 1.
EXIT: SUB P,X22
JRST @2(P)
TOOLRG: ERR <ACOS: Argument magnitude greater than 1.0; 0 returned>,1
TDZA 1,1 ;RETURN ZERO.
ZERARG: MOVE 1,PI2 ;ANSWER IS PI/2
JRST EXIT
ONE: 1.0
PI2: 201622077325
PII: 202622077325
BEND ACOS$
;FLOATING POINT SINGLE PRECISION ARCTANGENT OF TWO ARGUMENTS
;---------------------------------------------------------
;RETURNS ARCTANGENT OF A/B
;IF ARGUMENT IS IN 2ND QUADRANT, ATAN2(A/B) = PI + ATAN(A/B)
;IF ARGUMENT IS IN 3RD QUADRANT, ATAN2(A/B) = ATAN(A/B) - PI
;IF A/B OVERFLOWS (OR DIVIDE CHECK), THEN RETURN
; +PI/2 IF A GREATER THAN= 0, AND
; -PI/2 IF A LESS THAN 0.
;IF A/B UNDERFLOWS, THEN RETURN
; 0 IF B GREATER THAN= 0, AND
; +PI IF B LESS THAN 0 AND A GREATER THAN= 0,
; -PI IF B LESS THAN 0 AND A LESS THAN 0.
;THERE IS NO RESTRICTION ON THE ARGUMENTS
;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
; PUSH 17,ARG1
; PUSH 17,ARG2
; PUSHJ 17,$ATAN2
;THE ANSWER IS RETURNED IN ACCUMULATOR 1.
BEGIN ATAN2$
INTEGER SAVEAA
AA←13
BB←1
HERE(ATAN2$) ;ENTRY POINT TO ATAN2 ROUTINE
MOVEM AA,SAVEAA
JFOV .+1 ;CLEAR FLAGS FOR GOOD MEASURE (JAM)
MOVE AA,-2(P) ;PICK UP FIRST ARGUMENT
MOVE BB,-1(P) ;PICK UP SECOND ARGUMENT
FDVR AA, BB ;FORM AA/BB
JFOV OVUNFO ;EXTRA JFCL BECAUSE OF FDV HARDWARE BUG
JFOV OVUNFO ;SUPPRESS ERROR MESSAGE FROM
;OVTRAP IF NECESSARY AND GO TO
;OVUNFO IF AN EXCEPTION OCCURRED
;SPECIAL INTERRUPT CODE SETS
;OVPCWD AND RETURNS BEST VALUE
;IT CAN IN A.
PUSH P,AA ;CALCULATE ATAN(AA/BB)
PUSHJ P,ATAN$
SKIPL -1(P) ;IF BB GREATER THAN 0, SGN(ATAN2)=SGN(A)
JRST EXIT ;EXIT
JUMPGE BB, ATAN2A ;IS BB POSITIVE?
FADR BB, PII ;NO, SECOND QUADRANT, ADD PI
EXIT: MOVE AA,SAVEAA
SUB P,X33
JRST @3(P)
ATAN2A: FSBR BB, PII ;YES,3RD QUADRANT, SUBTRACT PI
JRST EXIT ;EXIT
OVUNFO: SKIPN AA,OVPCWD ;PICK UP FLAGS.
JSP AA,.+1
TLNE AA,100 ;SKIP IF OVERFLOW
JRST UNDER
MOVE BB,HALFPI ;ANSWER TO PI OVER 2
SKIPGE -2(P) ;SKIP IF ANS IS TO BE +
MOVNS BB
JRST EXIT
UNDER: JUMPL BB,BNEG
MOVEI BB,0
JRST EXIT ;RETURN 0
BNEG: MOVE BB,PII
SKIPGE -2(P)
MOVNS BB
JRST EXIT
PII: 202622077325 ;PII
HALFPI: 201622077325 ;PII/2
BEND ATAN2$
ENDCOM(TG2)
COMPIL(TG3,<TANH$,COSH$,SINH$>,<EXP$,X22,TRIGINI>,<HYPERBOLIC FUNCTIONS>)
;FLOATING POINT SINGLE PRECISION HYPERBOLIC TANGENT ROUTINE
;---------------------------------------------------------
;THIS ROUTINE CALCULATES THE TANH BY THE FOLLOWING ALGORITHM:
;IF ABSF(X) LESS THAN .00034, THEN TANH(X) = X
;IF ABSF(X) GREATER THAN 12.0, THEN TANH(X) = 1.0*SIGN(X)
;IF 0.17 LESS THAN= X LESS THAN 12.0, THEN TANH IS CALCULATED AS
; TANH(X) = 1.0 - 2(1.0 + EXP(2*X))**-1
;IF .00034 LESS THAN= X LESS THAN 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:
; PUSH 17,ARG
; PUSHJ 17,$TANH
;THE ANSWER IS RETURNED IN ACCUMULATOR 1
BEGIN TANH$
INTEGER SAVEBB,SAVTM1,SAVTM2
AA← 1
BB←13
TM1←14
TM2←15
HERE(TANH$) ;ENTRY TO TANH ROUTINE
MOVEM BB,SAVEBB
MOVEM TM1,SAVTM1
MOVEM TM2,SAVTM2
MOVE AA,-1(P) ;PICK UP THE ARGUMENT
MOVM BB, AA ;GET ABSF(ARGUMENT)
CAMGE BB, TA ;RETURN TANH(X)=X IF
JRST EXIT ;ABSF(X) .LE. .00034
CAMLE BB, T2 ;RETURN TANH(X) = 1.0*SIGN(X) IF
JRST TH5 ;ARGUMENT GREATER THAN 12.0
CAMGE BB, T3 ;USE RATIONAL APPROXIMATION IF
JRST TH3 ;ARGUMENT IS LESS THAN 0.17
FMPRI BB,202400 ;GET 2*ARG.
PUSH P,BB ;CALCULATE EXP(2X)
PUSHJ P,EXP$
MOVSI BB, (1.0) ;FORM 1.0
FAD AA, BB ;1 + EXP(2X)
FDVM BB, AA ;(1 + EXP(2X))**-1
FMPRI AA,202400 ;2*(1 + EXP(2X))**-1
FSBRM BB, AA ;1 - 2*(1 + EXP(2X))**-1
SKIPGE -1(P) ;SKIP AHEAD IF ARG WAS GREATER THAN= 0.
MOVNS AA ;OTHERWISE,NEGATE THE ANSWER.
EXIT: MOVE BB,SAVEBB
MOVE TM1,SAVTM1
MOVE TM2,SAVTM2
SUB P,X22
JRST @2(P)
TH3: FMP AA, T7 ;FORM 4*X*LOG(E) BASE 2
MOVEM AA, TM1 ;SAVE IT IN TM1
FMP AA, AA ;SQUARE IT
MOVEM AA, TM2 ;SAVE IT
FAD AA, T4 ;FORM F↑2 + T4
MOVE BB, T5 ;GET T5 IN ACCUMULATOR BB
FDV BB, AA ;T5/(F↑2 + T4)
FAD BB, T6 ;T6 + T5/(F↑2 + T4)
FMP BB, TM2 ;MULTIPLY BY F↑2
FAD BB, T7 ;ADD T7 (4*LOG(E) BASE 2)
MOVE AA, TM1 ;GET F IN ACCUMULATOR AA
TH5: FDV AA, BB ;DIVIDE F BY PARTIAL SUM
JRST EXIT ;EXIT
TA: 165544410070 ;0.00034
T2: 204600000000 ;12.0
T3: 176534121727 ;0.17
T4: 211535527022 ;349.6699888
T5: 204704333567 ;14.1384514018
T6: 173433723376 ;0.01732867951
T7: 203561250731 ;5.7707801636
BEND TANH$
;FLOATING POINT SINGLE PRECISION HYPERBOLIC COSINE FUNCTION.
;---------------------------------------------------------
;COSH(X) IS CALCULATED AS FOLLOWS:
; IF ABS(X) LESS THAN= 88.029,
; COSH(X) = 1/2(EXP(X) + 1.0/EXP(X))
; IF ABS(X) GREATER THAN 88.029 AND (ABS(X)-LN(2)) LESS THAN= 88.029,
; COSH(X) = EXP(ABS(X)-LN(2))
; IF (ABS(X)-LN(2)) GREATER THAN 88.029,
; COSH(X)=377777777777
; AND AN ERROR MESSAGE IS RETURNED.
;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
; PUSH 17,ARG
; PUSHJ 17,$COSH
;THE ANSWER IS RETURNED IN AC 1.
BEGIN COSH$
INTEGER SAVEAA
AA←13
HERE(COSH$) ;ENTRY TO HYPERBOLIC COSINE ROUTINE.
MOVEM AA,SAVEAA
MOVE 1,-1(P) ;PICK UP THE ARGUMENT.
MOVM AA,1 ;PUT ABS(X) IN AA
CAMLE 2,EIGHT8 ;IF ABS(X) GREATER THAN 88.029,
JRST OV88 ;GO TO OV88.
PUSH P,AA ;O'E, CALC. EXP(ABS(X))
PUSHJ P,EXP$
MOVSI AA,(1.0) ;PUT 1.0 IN AA
FDVR AA,1 ;CALC. 1.0/EXP(ABS(X)).
FADR 1,AA ;CALC. EXP(ABS(X)) + EXP(-ABS(X)).
FDVRI 1,202400 ;DIVIDE THIS BY 2.0.
EXIT: MOVE AA,SAVEAA
SUB P,X22 ;RETURN.
JRST @2(P)
OV88: FSBR AA,LN2BE ;FORM ABS(X)-LN(2).
CAMG AA,EIGHT8 ;OVERFLOW?
JRST EXPP ;NO,GO AHEAD.
ERR <COSH: Result too large - largest positive number returned>,1
HRLOI 1,377777 ;ANSWER = +INFINITY.
JRST EXIT ;RETURN
EXPP: PUSH P,AA ;CALC. EXP(ABS(X)-LN(2)).
PUSHJ P,EXP$
JRST EXIT ;RETURN.
EIGHT8: 207540074636 ;88.029
LN2BE: 200542710300 ;LOG(2) BASE E.
BEND COSH$
;FLOATING POINT SINGLE PRECISION HYPERBOLIC SINE FUNCTION.
;---------------------------------------------------------
;SINH IS CALCULATED AS FOLLOWS:
; IF ABS(X) GREATER THAN 88.029,
; SINH(X)=(EXP[ABS(X)-LN(2)])*SIGN(X)
; IF ABS(X) LESS THAN= 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:
; PUSH 17,ARG
; PUSHJ 17,$SINH
;THE ANSWER IS RETURNED IN AC 1.
BEGIN SINH$
INTEGER SAVEAA,SAVEBB,SAVSX2
AA←13
BB←14
SX2←15
HERE(SINH$) ;ENTRY TO HYPERBOLIC SINE ROUTINE.
MOVEM AA,SAVEAA
MOVEM BB,SAVEBB
MOVEM SX2,SAVSX2
MOVE 1,-1(P) ;PICK UP THE ARG.
MOVM AA,1 ;GET MAGNITUDE OF ARG IN AA
CAMLE AA,EIGHT8 ;IF ABS(X) GREATER THAN 88.029,
JRST OV88 ;THEN GO TO OV88.
CAMG AA,ONE10T ;IF ABS(X) LESS THAN= 0.10,
JRST SERIES ;THEN GO TO SERIES.
PUSH P,AA ;CALCULATE EXP(ABS(X)).
PUSHJ P,EXP$ ;ABS(X) IS IN AA
HRLZI BB,576400 ;PUT -1.0 IN BB
FDVR BB,1 ;CALC. -EXP(-ABS(X)).
FADR 1,BB ;CALC. EXP(ABS(X))-EXP(-ABS(X)).
FDVRI 1,202400 ;CALC. THIS/2.0
SKIPGE -1(P) ;ANSWER IS POSITIVE.
MOVNS 1,1 ;ANSWER IS NEGATIVE.
EXIT: MOVE AA,SAVEAA
MOVE BB,SAVEBB
MOVE SX2,SAVSX2
SUB P,X22
JRST @2(P)
SERIES: FMPR AA,AA ;CALC. X↑2.
JFOV .+1 ;SUPPRESS ERROR MESSAGE FROM OVTRAP.
;ON UNDERFLOW, SPECIAL INTERRUPT
;CODE RETURNS 0.
MOVEM AA,SX2 ;SAVE X↑2 IN SX2.
FDVR 2,ONE120 ;CALC.X↑2/120
JFOV .+1 ;SUPPRESS ERROR MESSAGE FROM OVTRAP.
;ON UNDERFLOW, SPECIAL INTERRUPT
;CODE RETURNS 0.
FADR AA,ONESIX ;CALC. (X↑2/120)+1/6
FMPR AA,SX2 ;MULTIPLY IT BY X↑2.
JFOV .+1 ;SUPPRESS ERROR MESSAGE FROM OVTRAP.
;ON UNDERFLOW SPECIAL INTERRUPT
;CODE RETURNS 0.
FADRI AA,(1.0) ;ADD 1.0.
FMPR 1,AA ;MULTIPLY BY X.
JRST EXIT ;RETURN.
OV88: FSBR AA,LN2BE ;CALC.ABS(X)-LN(2)
CAMG AA,EIGHT8 ;OVERFLOW?
JRST EXPP ;NO,GO TO CALC.
ERR <SINH: Result too large - largest positive number returned>,1
HRLOI 1,377777 ;SET ANS.=INFINITY.
JRST EXPP+2 ;GO TO SET SIGN OF ANS.
EXPP: PUSH P,AA ;CALC. EXP
PUSHJ P,EXP$
SKIPGE -1(P) ;RETURN ANS. GREATER THAN 0 IF X GREATER THAN 0.
MOVNS 1,1 ;O'E, ANS. LESS THAN 0.
JRST EXIT ;RETURN.
LN2BE: 200542710300 ;LN(2)
EIGHT8: 207540074636 ;88.029
ONE10T: 0.10
ONE120: 207740000000 ;120.0
ONESIX: 0.16666667
BEND SINH$
ENDCOM(TG3)
>;NOLOW