Trailing-Edge
-
PDP-10 Archives
-
BB-4157F-BM_1983
-
fortran/ots-debugger/mthgdb.mac
There are 9 other files named mthgdb.mac in the archive. Click here to see a list.
SEARCH MTHPRM
TV MTHGDB GFLOAT DOUBLE PRECISION ROUTINES ,7(3215)
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SUBTTL REVISION HISTORY
COMMENT \
***** Begin Revision History *****
1100 CKS --------
Creation.
1464 DAW -------
Put LERR's in new format.
1673 CDM 9-Sept-81
Added functions (GDIM, GINT, GNINT, GPROD IGNINT).
1721 CDM 15-Sept-81
Added fix to GDIM to prevent overflow for GDIM(-INF,+INF).
1724 BL 17-Sep-81
Clean up error messages.
***** Begin Mathlib *****
3055 CKS/JLC 17-Mar-82
Fix underflow/overflow error message in GEXP2.
3200 JLC 10-May-82
Mathlib integration. Change all LERRs to $LCALLs. Change
TWOSEG/RELOC to SEGMENT macros.
3201 RVM 20-May-82
Fix problem that made an argument of plus or minus one
to GASIN and GACOS illegal.
3206 AHM 8-Jun-82
Remove CKS's older version of GINT. in favor of CDM's newer
version.
3215 RVM 20-Aug-82
Change a CAIGE to a CAIG in GINT. Numbers less than 1.00D0
were not being set to 0.00D0.
***** End Revision History *****
\
PRGEND
TITLE GACOS ARC SINE AND ARC COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 19, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GACOS
EXTERN GACOS.
GACOS=GACOS.
PRGEND
TITLE GASIN ARC SINE AND ARC COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 19, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GASIN
EXTERN GASIN.
GASIN=GASIN.
PRGEND
TITLE GASIN. ARC SINE AND ARC COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 19, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DASIN(X) AND DACOS(X) ARE CALCULATED AS FOLLOWS
; LET R(G) = (G*(RP1+G*(RP2+G*(RP3+G*(RP4+G*RP5)))))/(Q0+G*(Q1
; +G*(Q2+G*(Q3+G*(Q4+G)))))
; (RP1,RP2,RP3,RP4,RP5,Q0,Q1,Q2,Q3,Q4, AND Q5 ARE GIVEN BELOW)
; LET S = Y + Y*R(G)
; FOR SITUATIONS 1. AND 3. BELOW,
; G = Y**2, WHERE Y = ABS(X)
; FOR SITUATIONS 2. AND 4., G = (1.0 - ABS(X))/2.0
; AND Y = -2.0*SQRT(G)
; LET W = ABS(X) AND TERM THE RESULT T.
; THEN, CONSIDER THE FOUR SITUATIONS:
; 1. DASIN, FOR W <= 0.5. T = S, AND T IS NEGATED
; FOR NEGATIVE X
; 2. DASIN, FOR W > 0.5. T = PI/2 + S, AND T IS NEGATED
; FOR NEGATIVE X
; 3. DACOS, FOR W <= 0.5. T = PI/2 - S IF X > 0
; = PI/2 + S IF X < 0
; 4. DACOS FOR W > 0.5. T = -S IF X > 0
; = PI + S IF X < 0.
;THE RANGE OF DEFINITION FOR DASIN/DACOS IS ABS(X) <= 1.0. AND ERROR MESSAGES
; WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE. DASIN/DACOS WILL BE SET
; TO + MACHINE INFINITY.
;REQUIRED (CALLED) ROUTINES: DSQRT
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; MOVEI L,ARG
; PUSHJ P,GASIN
; OR
; PUSHJ P,GACOS
;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (GACOS,.) ;ENTRY TO GACOS ROUTINE
HRRZI T0,2 ;SET FLAG TO TWO
MOVEM T0,FLAG ;MOVE FLAG TO MEMORY
JRST GETX ;GO TO GETX
HELLO (GASIN,.) ;ENTRY TO GASIN ROUTINE
SETZM FLAG ;SET FLAG TO ZEROES
GETX: DMOVE T0,@(L) ;OBTAIN X
JUMPGE T0,XPOS ;IF X IS NEGATIVE
DMOVN T0,T0 ;Y = -X
XPOS: CAMGE T0,HALF ;IF Y IS .LT. 1/2
JRST LE ;GO TO LE
CAME T0,HALF ;IF Y IS .GT. 1/2
JRST GT ;GO TO GT
JUMPE T1,LE ;
GT: CAMLE T0,ONE ;IF Y IS .GT. ONE
JRST MSTK ;GO TO MSTK
CAME T0,ONE ;IF Y IS .LT. ONE
JRST ALG ;GO TO ALG
JUMPE T1,ALG ;[3201] CHECK SECOND WORD
MSTK: $LCALL AOI
;LERR (LIB,%,<DASIN or DACOS: ABS(arg) .GT. 1.0; result = +infinity>)
HRLOI T0,377777 ;RESULT =
SETO T1, ;+MACHINE INFINITY
GOODBY (1) ;RETURN
ALG: PUSH P,T2 ;SAVE ACCUMULATORS
PUSH P,T3
PUSH P,T4
PUSH P,T5
MOVN T5,FLAG ;GET A COPY OF - FLAG
HRRZI T4,2 ;SET T4 TO 2
ADD T5,T4 ;I = 2-FLAG
HRLZI T2,200140 ;SET T2,T3 TO
SETZ T3, ;ONE
GFSB T2,T0 ;G = 1-Y
EXTEND T2,[GFSC -1] ;* 1/2
DMOVEM T2,TEMP ;SAVE A COPY OF G
FUNCT GSQRT.,<TEMP> ;Y = GSQRT(G)
EXTEND T0,[GFSC 1] ;* 2
DMOVN T0,T0 ;Y = -Y
MOVEM T5,I ;MOVE I TO MEMORY
JRST GOTY ;GO TO GOTY
LE: PUSH P,T2 ;SAVE ACCUMULATORS
PUSH P,T3
PUSH P,T4
PUSH P,T5
MOVE T5,FLAG ;I = FLAG
MOVEM T5,I ;MOVE I TO MEMORY
CAMGE T0,TWOM30 ;IF Y < TWOM30
JRST GOTRES ;GO TO GOTRES
DMOVE T2,T0 ;SAVE A COPY OF Y
GFMP T2,T2 ;G = Y*Y
GOTY: DMOVE T4,T2 ;SAVE A COPY OF G
GFAD T4,Q4 ;XDEN = G + Q4
GFMP T4,T2 ;*G
GFAD T4,Q3 ;+Q3
GFMP T4,T2 ;*G
GFAD T4,Q2 ;+Q2
GFMP T4,T2 ;*G
GFAD T4,Q1 ;+Q1
GFMP T4,T2 ;*G
GFAD T4,Q0 ;+Q0
DMOVEM T2,TEMP ;SAVE A COPY OF G
GFMP T2,RP5 ;XNUM = G*RP5
GFAD T2,RP4 ;+RP4
GFMP T2,TEMP ;*G
GFAD T2,RP3 ;+RP3
GFMP T2,TEMP ;*G
GFAD T2,RP2 ;+RP2
GFMP T2,TEMP ;*G
GFAD T2,RP1 ;+RP1
GFMP T2,TEMP ;*G
GFDV T2,T4 ;RESULT = XNUM/XDEN
GFMP T2,T0 ;*Y
JRST GETI ;GO TO GETI
GOTRES: DMOVE T2,AHI ;ZERO T2
GETI: MOVE T5,I ;GET A COPY OF I
SKIPE FLAG ;IF FLAG IS .NE. 0
JRST NE ;GO TO NE
GFAD T0,ALO(T5)
GFAD T0,T2
GFAD T0,AHI(T5)
SKIPGE @(L) ;IF X IS NEGATIVE
DMOVN T0,T0 ;RESULT = -RESULT
POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4
POP P,T3
POP P,T2
GOODBY (1) ;RETURN
NE: SKIPGE @(L) ;IF X IS NEGATIVE
JRST NEGX ;GO TO NEGX
DMOVN T0,T0 ;RESULT = -RESULT
GFAD T0,ALO(T5)
GFSB T0,T2
GFAD T0,AHI(T5)
POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4
POP P,T3
POP P,T2
GOODBY (1) ;RETURN
NEGX: GFAD T0,BLO(T5)
GFAD T0,T2
GFAD T0,BHI(T5)
POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4
POP P,T3
POP P,T2
GOODBY (1) ;RETURN
RP1: DOUBLE 577211206522,276137263562 ;-.27368494524164255994D+2
RP2: DOUBLE 200671152471,260255221215 ;.57208227877891731407D+2
RP3: DOUBLE 577130237232,262622254077 ;-.39688862997504877339D+2
RP4: DOUBLE 200450470273,047303444713 ;.10152522233806463645D+2
RP5: DOUBLE 577723321022,120636063337 ;-.69674573447350646411D+0
Q0: DOUBLE 576726744776,016507406363 ;-.16421096714498560795D+3
Q1: DOUBLE 201164111170,200753410270 ;.41714430248260412556D+3
Q2: DOUBLE 576620210610,035225367030 ;-.38186303361750149284D+3
Q3: DOUBLE 201045571744,262621615746 ;.15095270841030604719D+3
Q4: DOUBLE 577220264274,210147474061 ;-.23823859153670238830D+2
AHI: DOUBLE 0,0
A1HI: DOUBLE 200162207732,242102643021 ;PI/2
BHI: DOUBLE 200262207732,242102643021 ;PI
B1HI: DOUBLE 200162207732,242102643021 ;PI/2
ALO: DOUBLE 0,0 ;
A1LO: DOUBLE 170551423063,024270033407 ;NEXT 59 BITS OF PI/2
BLO: DOUBLE 170651423063,024270033407 ;NEXT 59 BITS OF PI
B1LO: DOUBLE 170551423063,024270033407 ;NEXT 59 BITS OF PI/2
TWOM30: 174340000000 ;HIGH ORDER PART OF 2**(-30)
ONE: 200140000000 ;1.0
HALF: 200040000000 ;1/2
SEGMENT DATA
I: 0
FLAG: 0
TEMP: DOUBLE 0,0
PRGEND
TITLE GATAN2 ARC TAN FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MAY 9, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GATAN2
EXTERN GATN2.
GATAN2=GATN2.
PRGEND
TITLE GATAN ARC TAN FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL Chris Smith/CKS
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GATAN
EXTERN GATAN.
GATAN=GATAN.
PRGEND
TITLE GATAN. ARC TAN FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DATAN(X) is computed as follows:
;
;If X < 0, compute DATAN(|X|) below, then DATAN(X) = -DATAN(-X).
;
;If X > 0, use the identity
;
; DATAN(X) = DATAN(XHI) + DATAN(Z)
; Z = (X - XHI) / (1 + X*XHI)
;
;where XHI is chosen so that |Z| <= tan(pi/32).
;
;XHI is chosen to be exactly representable as a double precision number,
;and so Z can be calculated without loss of significance.
;
;DATAN(XHI) is found by table lookup. It is stored as ATANHI + ATANLO to
;provide guard bits for the final addition to DATAN(Z).
;
;DATAN(Z) is evaluated by means of a polynomial approximation from Hart et al.
;(formula 4904).
;
;If X < tan(pi/32), DATAN(X) = DATAN(Z).
;If X > 1/tan(pi/32), DATAN(X) = pi/2 - DATAN(1/X).
;
;If tan(pi/32) < X < 1/tan(pi/32), an appropriate XHI is obtained by indexing
;into a table. The table tells which XHI to use for various ranges of X.
;The index into the table is formed from the low 3 exponent bits and the high
;3 fraction bits of X.
SEARCH MTHPRM
SEGMENT CODE
SALL
T6=6 ;MORE AC NAMES
T7=7
T8=10
HELLO (GATAN,.) ;DATAN ENTRY
PUSH P,T2 ;SAVE REGISTERS
PUSH P,T3
PUSH P,T4
PUSH P,T5
PUSH P,T6
DMOVE T0,@(L) ;GET ARGUMENT X
MOVEM T0,SGNFLG ;SAVE ARGUMENT SIGN FOR RESULT
CAIGE T0,0 ;GET |X|
DMOVN T0,T0
CAML T0,MAXX ;IS X LARGE?
JRST LARGEX ;YES, GO USE ATAN(X) = PI/2 - ATAN(1/X)
SETZ T2, ;T2 WILL GET OFFSET INTO XHI TABLES
CAMG T0,MINX ;IS X SMALL ENOUGH THAT NO ARG REDUCTION IS REQUIRED?
JRST CALC ;YES, JOIN CALCULATION BELOW
MOVE T3,T0 ;GET HIGH WORD OF X
LSHC T2,^D12 ;GET EXPONENT, SHIFT HIGH FRACTION BIT
;(ALWAYS 1) INTO SIGN BIT OF T3
ASHC T2,3 ;GET THREE FRACTION BITS, LEAVING
;THE 1 BEHIND
HLRZ T2,OFFSET-20000+24(T2) ;GET OFFSET INTO XHI TABLES
DMOVE T3,T0 ;GET A COPY OF X
GFSB T0,XHI(T2) ;GET X-XHI
GFMP T3,XHI(T2) ; X*XHI
GFAD T3,ONE ; 1 + X*XHI
GFDV T0,T3 ; (X-XHI) / (1 + X*XHI)
;Here T0 has the reduced argument Z with |Z| <= tan(pi/32). T2 has the
;index into the ATAN(XHI) tables ATANHI and ATANLO. SGNFLG is set negative
;if the result should be negated.
CALC: DMOVE T3,T0 ;GET |Z|
CAIGE T3,0
DMOVN T3,T3
CAMG T3,EPS ;IS IT SMALL ENOUGH THAT ATAN(Z) = Z?
JRST SMALLX ;YES, BE FAST, AVOID UNDERFLOW
GFMP T3,T3 ;GET Z**2
DMOVE T5,P06 ;GET P(Z**2)
GFMP T5,T3
GFAD T5,P05
GFMP T5,T3
GFAD T5,P04
GFMP T5,T3
GFAD T5,P03
GFMP T5,T3
GFAD T5,P02
GFMP T5,T3
GFAD T5,P01
GFMP T3,T5
GFMP T3,T0 ; * Z
GFAD T0,T3 ; + Z = ATAN(Z)
SMALLX: GFAD T0,ATANLO(T2) ; + ATAN(XHI) LOW
GFAD T0,ATANHI(T2) ; + ATAN(XHI) HI = DATAN(X)
SKIPG SGNFLG ;ATTACH SIGN TO RESULT
DMOVN T0,T0
RET: POP P,T6 ;RETURN
POP P,T5
POP P,T4
POP P,T3
POP P,T2
POPJ P,
LARGEX: DMOVE T3,T0 ;GET -1/X
DMOVN T0,ONE
GFDV T0,T3
MOVEI T2,PI2OFFS ;GET OFFSET OF PI/2
JRST CALC ;GO COMPUTE PI/2 + ATAN(-1/X)
SUBTTL DATAN2 (Y,X)
;To compute DATAN2(Y,X), let U = |Y| and V = |X|, and compute DATAN(U/V).
;Then find DATAN2(Y,X) based on the signs of Y and X as follows:
;
; X Y DATAN2(Y/X)
;
; pos pos DATAN(U/V)
; pos neg -DATAN(U/V)
; neg pos -(DATAN(U/V) - pi)
; neg neg DATAN(U/V) - pi
;
;The add of -pi is combined with the add of DATAN(XHI) which is the last step
;of the DATAN algorithm.
;
;The reduced argument for DATAN is Z = (U/V - XHI) / (1 + U/V * XHI).
;This is rewritten as (U - V*XHI) / (V + U*XHI). To accurately calculate the
;numerator, find VHI and VLO with
;
; V = VHI + VLO
; VHI has at most 27 significant bits
; VLO has at most 35 significant bits
;
;and choose XHI with at most 13 significant bits. Then VHI*XHI and VLO*XHI can
;be exactly represented as double precision numbers, and the numerator is
;
; U - V*XHI = (U - VHI*XHI) - VLO*XHI
HELLO (GATN2.,GATAN2) ;DATAN2 ENTRY
PUSH P,T2 ;SAVE REGISTERS
PUSH P,T3
PUSH P,T4
PUSH P,T5
PUSH P,T6
DMOVE T0,@0(L) ;GET Y
GFDV T0,@1(L) ;GET Y/X
JFCL EXCEP ;OVERFLOW AND UNDERFLOW CAN OCCUR
DMOVEM T0,SGNFLG ;RESULT SHOULD BE MULTIPLIED BY SGN(Y/X)
DMOVE T3,T0 ;GET |Y/X|, ATAN ARG
CAIGE T3,0
DMOVN T3,T3
CAML T3,MAXX ;IS |Y/X| LARGE?
JRST LARGE2 ;YES, GO USE ATAN(Y/X) = PI/2 - ATAN(X/Y)
SETZ T2, ;GET OFFSET INTO ATAN TABLES
CAMG T3,MINX ;IS |Y/X| SMALL ENOUGH TO USE POLYNOMIAL DIRECTLY?
JRST [DMOVE T0,T3 ;YES, DO SO
JRST SMALL2]
LSHC T2,^D12 ;GET INDEX INTO OFFSET TABLES
ASHC T2,3
HLRZ T2,OFFSET-20000+24(T2) ;GET INDEX INTO XHI TABLES
PUSH P,T7 ;SAVE MORE REGISTERS
PUSH P,T8
DMOVE T0,@0(L) ;GET |Y| = U
CAIGE T0,0
DMOVN T0,0
DMOVEM T0,USAVE ;SAVE FOR LATER
DMOVE T3,@1(L) ;GET |X| = V
CAIGE T3,0
DMOVN T3,T3
MOVE T5,T3 ;GET A COPY OF V
DMOVE T7,T3 ;GET ANOTHER
SETZ T6, ;GET HIGH 27 BITS OF V = VHI
GFSB T7,T5 ;GET LOW 35 BITS OF V = VLO
GFMP T5,XHI(T2) ;GET V*XHI = VHI * XHI
GFMP T7,XHI(T2) ; + VLO * XHI
GFSB T0,T5 ;GET (U - VHI*XHI)
GFSB T0,T7 ; - VLO*XHI
DMOVE T5,USAVE ;GET U BACK
GFMP T5,XHI(T2) ;GET U * XHI
GFAD T3,T5 ;GET V + U*XHI
GFDV T0,T3 ;GET (U - V*XHI) / (V + U*XHI)
POP P,T8 ;RESTORE REGISTERS
POP P,T7
SMALL2: SKIPGE @1(L) ;IF SECOND ARG (X) IS NEGATIVE
ADDI T2,MPIOFFS ; ADD -PI TO RESULT
JRST CALC ;GO GET ATAN AND RETURN
LARGE2: DMOVN T0,@1(L) ;GET -X/Y
GFDV T0,@0(L)
MOVMM T0,SGNFLG ;SET SGNFLG POSITIVE
MOVEI T2,PI2OFFS ;ADD PI/2 TO RESULT IF FIRST ARG (Y) IS POSITIVE
SKIPGE @0(L)
MOVEI T2,MPI2OFFS ;ADD -PI/2 TO RESULT IF FIRST ARG (Y) IS NEGATIVE
JRST CALC ;GO COMPUTE +/- PI/2 + ATAN(-X/Y)
EXCEP: SKIPN @1(L) ;CHECK FOR DIVIDE BY 0
JRST DIVCHK ;IF DIVIDE CHECK, GO CHECK FOR ATAN2(0,0)
JUMPN T0,OVER ;IF OVERFLOW, RESULT IS +/- PI/2
SKIPL @1(L) ;IF UNDERFLOW, CHECK SECOND ARGUMENT
$LCALL RUN,RET
; LERR (LIB,%,<DATAN2: result underflow>,,RET)
;IF SECOND ARG (X) POSITIVE, RESULT UNDERFLOWS
DMOVN T0,MPI ;ELSE RESULT IS PI WITH SIGN OF FIRST ARG
JRST YSIGN ;GO ATTACH SIGN AND RETURN
DIVCHK: SKIPE @0(L) ;CHECK FOR BOTH ARGS 0
JRST OVER ;ATAN2(NONZERO,0) IS SAME AS OVERFLOW
$LCALL BAZ
;LERR (LIB,%,<DATAN2: both arguments are zero, result=0.0>)
SETZB T0,T1 ;RETURN 0
JRST RET
OVER: DMOVE T0,PI2 ;OVERFLOW, RESULT IS PI/2 WITH SIGN OF FIRST ARG
YSIGN: SKIPGE @0(L) ;ATTACH SIGN OF FIRST ARGUMENT (Y)
DMOVN T0,T0
JRST RET ;RETURN
SUBTTL TABLES
;This table is indexed by the low 3 exponent bits and the high 3 fraction bits
;of X, where MINX < X < MAXX. It gives the offset into XHI, ATANHI, and ATANLO
;of a suitable XHI.
DEFINE OFFS (X) <
IRP X,<XWD 2*X,X>
>
RADIX 10
OFFSET: OFFS <1,1,1,2,2,3,3,4,4,4,5,5,5,6,6,7,7,7,8,8,8,9,9,9>
OFFS <10,10,10,10,11,11,11,12,12,12,12,12,13,13,13,13>
OFFS <13,13,14,14,14,14,14,14,14,14,14,14,14,14,14>
RADIX 8
PI2OFFS=2*^D15 ;OFFSET INTO ATANHI AND ATANLO OF PI/2
MPIOFFS=2*^D16 ;OFFSET OF -PI
MPI2OFFS=2*^D31 ;OFFSET OF -PI/2
XHI: EXP 000000000000,0 ; .0000000000000000000 not used
EXP 177565774000,0 ; .1054534912109375000 X1
EXP 177640774000,0 ; .1288757324218750000 X2
EXP 177647774000,0 ; .1562194824218750000
EXP 177661764000,0 ; .1952209472656250000
EXP 177677740000,0 ; .2497558593750000000
EXP 177747754000,0 ; .3121948242187500000
EXP 177761720000,0 ; .3898925781250000000
EXP 177777630000,0 ; .4984130859375000000
EXP 200051574000,0 ; .6522216796875000000
EXP 200067404000,0 ; .8673095703125000000
EXP 200145344000,0 ; 1.170166015625000000
EXP 200164510000,0 ; 1.645019531250000000
EXP 200251104000,0 ; 2.570800781250000000
EXP 200352700000,0 ; 5.359375000000000000 X14
ATANHI: EXP 000000000000,000000000000 ; .0000000000000000000 ATAN(0)
EXP 177565626152,025171771712 ; .1050651824695016828 ATAN(X1)
EXP 177640637315,230510354743 ; .1281692623534198424
EXP 177647527650,022607523014 ; .1549669515088296697
EXP 177661266130,275334150302 ; .1927961221331547523
EXP 177676517562,252343141216 ; .2447488705187413410
EXP 177746567507,360775545336 ; .3026068193134274480
EXP 177757453662,234657557072 ; .3717628303965550499
EXP 177773136266,273706162457 ; .4623772720674932798
EXP 200044771623,334166732152 ; .5779354489017924541
EXP 200055563263,207236241201 ; .7144577222199199901
EXP 200067214043,155315244501 ; .8636495725719767220
EXP 200140622720,225243500416 ; 1.024591515455200379
EXP 200146311711,011245646124 ; 1.199822549393342833
EXP 200154271467,254276104335 ; 1.386328658261967262 ATAN(X14)
PI2: EXP 200162207732,242102643022 ; 1.570796326794896620 PI/2
MPI: EXP 577515570045,135675134756 ;-3.141592653589793241 -PI
EXP 577517324610,256420774615 ;-3.036527471120291553 -PI + ATAN(X1)
EXP 577517622022,067321553615 ;-3.013423391236373393 -PI + ATAN(X2)
EXP 577520155437,337025522117 ;-2.986625702080963569
EXP 577520643352,351552743372 ;-2.948796531456638489
EXP 577521515034,210413303027 ;-2.896843783071051899
EXP 577522447016,133774711512 ;-2.838985834276365791
EXP 577523535433,261363112666 ;-2.769829823193238186
EXP 577525103674,065265753224 ;-2.679215381522299960
EXP 577526766412,124732723411 ;-2.563657204688000783
EXP 577531124722,077544605217 ;-2.427134931369873246
EXP 577533433056,071160406077 ;-2.277943081017816514
EXP 577536101415,250416775165 ;-2.117001138134592862
EXP 577601672023,305040140060 ;-1.941770104196450408
EXP 577607651602,150070376271 ;-1.755263995327825979 -PI + ATAN(X14)
EXP 577615570045,135675134756 ;-1.570796326794896620 -PI/2
ATANLO: EXP 000000000000,000000000000 ; .0000000000000000000
EXP 607611362655,250215702700 ;-.9237013676958846587E-19
EXP 170143152717,266776243666 ; .5964602757438210869E-19
EXP 170154077372,225515747423 ; .7474896823762959724E-19
EXP 607735160270,370510376242 ;-.2946026702232522009E-19
EXP 610004420274,014224300071 ;-.2518569148307390217E-19
EXP 170347024456,077470314101 ; .2645467902793737610E-18
EXP 607510312505,015573123116 ;-.1883944550948088893E-18
EXP 170254520527,135456552170 ; .1513056980995441703E-18
EXP 170452555761,155354735174 ; .5788933265131145596E-18
EXP 170453236423,212017737506 ; .5869551379299690115E-18
EXP 607321026727,320125051737 ;-.6363620490158901181E-18
EXP 170443623615,021735205127 ; .4850262996822095826E-18
EXP 607222115425,120272454763 ;-.1242727478712024599E-17
EXP 607226174644,043107635453 ;-.1131804334593366021E-17
EXP 607223046146,050560067016 ;-.1217705177797396539E-17
EXP 170654731631,327217710762 ; .2435410355594793078E-17
EXP 607123161307,104424247040 ;-.2427449340111014898E-17
EXP 607106015110,124777656057 ;-.3142794913755447875E-17
EXP 170471157106,257651240651 ; .7754358478556155783E-18
EXP 170674303434,273125014744 ; .3273311826560871401E-17
EXP 170570727666,236602064744 ; .1542862926123315625E-17
EXP 170543470577,076355704763 ; .9652336698973597403E-18
EXP 607111346356,050107456126 ;-.2957154527430437100E-17
EXP 170577335536,232205477162 ; .1719354315705933698E-17
EXP 607336325130,312454401102 ;-.4551432698457065547E-18
EXP 607127601336,271623700703 ;-.2181804934405659197E-17
EXP 607101137417,313245123351 ;-.3405122121351518328E-17
EXP 170665676575,033607152207 ; .2920436655277002656E-17
EXP 170554001110,376732276727 ; .1192682876882768478E-17
EXP 170560060327,321547457416 ; .1303606021001427053E-17
EXP 170554731631,327217710762 ; .1217705177797396539E-17
;COEFFICIENTS OF APPROXIMATION POLYNOMIAL ATAN(X) = X*P(X**2)
P06: EXP 177546176241,173026455171 ; .074700604980000000
P05: EXP 600221360346,262546426550 ;-.090879628821850000
P04: EXP 177570707036,320713526601 ; .111110916853003200
P03: EXP 600133333333,171014105013 ;-.142857142198848259
P02: EXP 177663146314,314620200167 ; .199999999998937080
P01: EXP 600025252525,125252526621 ;-.333333333333332690
ONE: EXP 200140000000,000000000000 ; 1.00000000000000000
EPS: EXP 174372113547 ;LARGEST X WITH DATAN(X) = X (HIGH WORD)
;(IE, ALL DOUBLE PRECISION X WITH HIGH
; WORD OF X <= EPS HAVE DATAN(X)=X).
MINX: EXP 177562332734 ;TAN(PI/32) (HIGH WORD, ROUNDED DOWN)
MAXX: EXP 200450471543 ;1/TAN(PI/32) (HIGH WORD, ROUNDED UP)
SEGMENT DATA
SGNFLG: BLOCK 1 ;SIGN TO BE ATTACHED TO RESULT
USAVE: BLOCK 2 ;TEMP STORAGE FOR U
PRGEND
TITLE GCOSH HYPERBOLIC COSINE FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 7, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GCOSH
EXTERN GCOSH.
GCOSH=GCOSH.
PRGEND
TITLE GCOSH. HYPERBOLIC COSINE FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 7, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DCOSH(X) IS CALCULATED AS FOLLOWS
; LET V BE APPROXIMATELY 2 SO THAT LNV AND ABS(X)+LN V
; CAN BE EXACTLY REPRESENTED WHEN X IS EXACTLY REPRESENTABLE.
; THEN, LETTING Y = ABS(X), FOR
; Y <= 709.089565, DCOSH = (EXP(Y)+EXP(-Y))/2,
; 709.089565 <= Y < 1024 * LN(2)
; DCOSH = (V/2)*EXP(Y - LN V),
; Y >= 1024 * LN(2), DCOSH = +MACHINE INFINITY
;THE RANGE OF DEFINITION FOR DCOSH IS ABS(X) < 1024 * LN(2) AND ERROR MESSAGES
; WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE. DCOSH WILL BE SET
; TO + MACHINE INFINITY.
;REQUIRED (CALLED) ROUTINES: GEXP
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; MOVEI L,ARG
; PUSHJ P,GCOSH
;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (GCOSH,.) ;ENTRY FOR GCOSH ROUTINE
DMOVE T0,@(L) ;OBTAIN X
JUMPGE T0,DCSH ;IF X IS NEGATIVE
DMOVN T0,T0 ;Y = -X
DCSH: CAMGE T0,TWENT2 ;IF ABS(X) < 22, GO TO
JRST ALG2 ; FULL CALCULATION.
CAMG T0,BIGX ;IF Y IS .LE. BIGX
JRST EASY ;THEN GO TO EASY
CAMGE T0,YMAXHI ;IF HI OF Y < YMAXHI
JRST GETW ; GO TO GETW
CAME T0,YMAXHI ;IF HI OF Y > YMAXHI
JRST MSTK ;GO TO MSTK
CAMGE T1,YMAXLO ;HI OF Y = YMAXHI
JRST GETW ; IF LO OF Y .LE. YMAXLO, GO TO GETW
MSTK: HRLOI T0,377777 ;SET RESULT TO
SETO T1, ;MACHINE INFINITY
$LCALL ROV
;LERR (LIB,%,<DCOSH: result overflow>)
GOODBY (1) ;RETURN
GETW: GFAD T0,LNV ;W = Y-LNV
DMOVEM T0,TEMP ;MOVE W TO A TEMP REGISTER
FUNCT GEXP.,<TEMP> ;Z = EXP(W)
DMOVEM T0,TEMP ;SAVE A COPY OF Z
GFMP T0,CON1 ;RESULT = Z*CON1
GFAD T0,TEMP ;+Z
JFCL MSTK
GOODBY (1) ;RETURN
EASY: DMOVEM T0,TEMP ;MOVE Y TO TEMP
FUNCT GEXP.,<TEMP> ;GET ITS EXP
EXTEND T0,[GFSC -1] ;DIV BY 2
GOODBY(1) ;RETURN
ALG2: CAMG T0,TM30 ;IF ABS(X) .LT. 2**(-30) THE
JRST TINY ; RESULT IS 1.
PUSH P,T2 ;SAVE MORE ACCUMULATORS
PUSH P,T3
PUSH P,T4
PUSH P,T5
DMOVEM T0,TEMP ;MOVE Y TO A TEMPORARY REGISTER
FUNCT GEXP.,<TEMP> ;Z = EXP(Y)
DMOVE T2,T0 ;SAVE A COPY OF Z
HRLZI T4,200140 ;SET T4 TO 1.0
SETZ T5, ;CLEAR SECOND WORD
GFDV T4,T2 ;1/Z
GFAD T0,T4 ;+Z
EXTEND T0,[GFSC -1] ;*1/2
POP P,T5
POP P,T4 ;RESTORE ACCUMULATORS
POP P,T3
POP P,T2
GOODBY (1) ;RETURN
TINY: HRLZI T0,200140 ;RESULT IS 1
SETZ T1,
GOODBY (1)
ONE: 200140000000 ;1.0
BIGX: 201254242673 ;709.089565
YMAXHI: 201254271027 ;709.782712893383998706
YMAXLO: 367643475715 ;
LNV: DOUBLE 577723506750,010134300000 ;-.693147180559947174D0
CON1: DOUBLE 172041452000,000000000000 ;.186417721E-14
TWENT2: 200554000000 ;22
TM30: 174340000000 ;2**(-30)
SEGMENT DATA
TEMP: DOUBLE 0,0
PRGEND
TITLE GEXP EXPONENTIAL FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL.INC APRIL 3, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GEXP
EXTERN GEXP.
GEXP=GEXP.
PRGEND
TITLE GEXP. EXPONENTIAL FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC.APRIL 3, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DEXP(X) IS CALCULATED AS FOLLOWS
; IF X <= -710.475860073943942, EXP = 0
; IF X > 709.089565712824051, EXP = +MACHINE INFINITY
; IF X = 0.0, EXP = 1
; OTHERWISE,
; THE ARGUMENT REDUCTION IS:
; LET X1 = [X], THE GREATEST INTEGER IN X
; X2 = X - X1
; N = THE NEAREST INTEGER TO X/LN(2)
; THE REDUCED ARGUMENT IS G = ((X1 - N*C1)+X2)+N*C2
; WHERE C1 = .543 (OCTAL),
; AND C2 IS GIVEN BELOW
; THE CALCULATION IS:
; EXP = R(G)*2**(N+1)
; WHERE R(G) = 0.5 + G*P/(Q - G*P)
; P = ((((P2*G**2)+P1)*G**2)+P0)*G**2
; Q = (((((Q3*G**2)+Q2)*G**2)+Q1)*G**2)+Q0
; P0, P1, P2, Q0, Q1, Q2, AND Q3 ARE GIVEN BELOW AS
; XP0, XP1, XP2, XQ0, XQ1, XQ2, AND XQ3 .
;THE RANGE OF DEFINITION FOR DEXP IS GIVEN ABOVE, AND ERROR MESSAGES
; WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE
;REQUIRED (CALLED) ROUTINES: NONE
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; MOVEI L,ARG
; PUSHJ P,GEXP
;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0,
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1.
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (GEXP,.) ;ENTRY TO GEXP ROUTINE
DMOVE T0,@(L) ;GET DP ARGUMENT
CAMLE T0,SMLXHI ;IF HI OF X > SMLXHI
JRST NXTCHK ; GO TO NXTCHK
CAME T0,SMLXHI ;IF HI OF X < SMLXHI
JRST OUT2 ; GO TO OUT2
CAMGE T1,SMLXLO ;HI PART = SMLXHI,
JRST OUT2 ; IF LO OF X < SMLXLO, OUT2
NXTCHK: CAMGE T0,BIGXHI ;IF HI OF X < BIGXHI
JRST EXP1 ; GO TO EXP1
CAME T0,BIGXHI ;IF HI OF X > BIGXHI
JRST MSTK ; GO TO MSTK
CAMG T1,BIGXLO ;IF LO OF X .LE. BIGXLO,
JRST EXP1 ; GO TO EXP1
MSTK: $LCALL ROV
;LERR (LIB,%,<DEXP: result overflow>)
HRLOI T0,377777 ;DEXP = +MACHINE INFINITY
HRLOI T1,377777
GOODBY (1) ;RETURN
OUT2: $LCALL RUN
;LERR (LIB,%,<DEXP: result underflow>)
MOVEI T0,0 ;EXP = 0
MOVEI T1,0
GOODBY (1) ;RETURN
EXP1: JUMPE T0,[MOVSI T0,200140 ;RETURN 1.0 FOR ARG OF ZERO
GOODBY (1)] ;EXIT
PUSH P,T2 ;SAVE ACCUMULATORS
PUSH P,T3
PUSH P,T4
PUSH P,T5
MOVM T2,T0 ;GET ABS(ARGHI)
CAML T2,LN2OV2 ;IF .GE. (LN(2))/2
JRST REDUCE ; MUST REDUCE ARGUMENT.
SETZ T4,
MOVEM T4,SAVEN ;SET N TO ZERO.
DMOVE T4,T0 ;GET COPY OF ARG.
JRST MERGE ;MERGE WITH MAIN FLOW.
REDUCE: DMOVE T2,T0 ;GET COPY OF ARG
GFMP T2,RNDLN2 ;ARG/LN2
EXTEND T2,[GFIXR T2] ;NEAREST INTEGER = N
EXTEND T2,[GFLTR T2] ;FLOAT N
MOVEM T2,SAVEN ;SAVEN
GFMP T2,C1 ;N*C1
GFAD T0,T2 ;X + N*C1
MOVE T2,SAVEN ;RETRIEVE N
MOVEI T3,0 ;ZERO SECOND WORD
GFMP T2,C2 ;FORM N*C2
GFAD T0,T2 ;N*C2 + (N*C1 + X)
DMOVE T4,T0 ;SAVE G
MOVM T2,T4 ;GET ABS(G)
MERGE: CAML T2,TWOM30 ;IF REDUCED ARG IS>= 2**-30
JRST APPRX ; GO TO APPRX
GFAD T0,ONE ;R(G) = 1. + G
EXTEND T0,[GFSC -1] ;*.5
JRST BRNCH ;GO TO BRNCH
APPRX: GFMP T4,T4 ;Z = G*G
DMOVE T2,T4 ;SAVE Z
GFMP T4,XP2 ;Z*XP2
GFAD T4,XP1 ;+XP1
GFMP T4,T2 ;* Z
GFAD T4,XP0 ;+ XP0
GFMP T0,T4 ;* G
DMOVE T4,T2 ;SAVE Z
GFMP T4,XQ3 ;XQ3*Z
GFAD T4,XQ2 ;+XQ2
GFMP T4,T2 ;*Z
GFAD T4,XQ1 ;+ XQ1
GFMP T4,T2 ;* Z
GFAD T4,XQ0 ; + XQ0
GFSB T4,T0 ; XQ - G*XP
GFDV T0,T4 ;(G*XP)/(XQ-G*XP)
GFAD T0,XQ0 ; + .5
BRNCH: MOVE T4,SAVEN ;RETRIEVE N
SETZ T5, ;T5=0
EXTEND T4,[GFIX T4]
ADDI T4,1 ;N = N+1
EXTEND T0,[GFSC 0(T4)] ;ADD N TO THE EXPONENT
POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4
POP P,T3
POP P,T2
RET: GOODBY (1) ; RETURN
SMLXHI: 576523460613 ;-710.475860073943942
SMLXLO: 202140360224 ;
BIGXHI: 201254242673 ;709.089565712824051
BIGXLO: 161647554056 ;
RNDLN2: DOUBLE 200156125073,051270137606 ;1.44269504088896341 = 1/LN2
ONE: DOUBLE 200140000000,000000000000 ;1.0D0
TWOM30: DOUBLE 174340000000,000000000000 ;2**-30
C1: DOUBLE 577723500000,000000000000 ;-0.693359375D0
C2: DOUBLE 176467500202,343020625033 ;2.12194440054690583D-4
XP0: DOUBLE 177740000000,000000000000 ;0.250D0
XP1: DOUBLE 177176035137,221241124545 ;0.757531801594227767D-2
XP2: DOUBLE 176241055041,127151103610 ;0.315551927656846464D-4
XQ0: DOUBLE 200040000000,000000000000 ;0.5D0
XQ1: DOUBLE 177472134502,216775477655 ;0.568173026985512218D-1
XQ2: DOUBLE 176651274142,341215233732 ;0.631218943743985036D-3
XQ3: DOUBLE 175462315430,147120212512 ;0.751040283998700461D-6
LN2OV2: DOUBLE 177754271027,367643475715
SEGMENT DATA
SAVEN: 0
PRGEND
TITLE GEXP2. DOUBLE ** INTEGER EXPONENTIATION
SUBTTL MARY PAYNE/MHP/CKS
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
SEGMENT CODE
SALL
; Mary Payne, July 30, 1982
HELLO (GEXP2,.)
DMOVEM T2,SAVE2 ;Save T2-T4
MOVEM T4,SAVE4
DMOVE T0,ONE ;Floating 1 to T0-T1
MOVM T2,@1(L) ;|exponent| to T2
JUMPE T2,EXP0 ;Exponent = 0 is special
DMOVE T3,@0(L) ;Base to T3-T4.
JUMPN T3,STEP1 ;If base not 0 go to main flow
JRST BASE0 ;Else to special code
LOOP: GFMP T3,T3 ;Square current result
JOV OVER2 ;Over/underflow possible
STEP1: TRNE T2,1 ;If exponent is odd
GFMP T0,T3 ; update current result
JOV OVER ; Branch on over/underflow
LSH T2,-1 ;Discard low bit of exponent
JUMPN T2,LOOP ;Iterate if not 0
SKIPL @1(L) ;If exponent > 0
JRST RET ; return
DMOVE T2,T0 ;Copy result
DMOVE T0,ONE ;Get reciprocal of
GFDV T0,T2 ;result; Underflow impossible
JOV OVMSG ; On overflow get message
JRST RET ;Else return
EXP0: SKIPE @0(L) ;Exponent 0. If base not
JRST RET ; 0, result is 1. Return
$LCALL ZZZ ;zero**zero is indeterminate, result=zero
SETZB T0,T1 ;Store 0
JRST RET
BASE0: SKIPL @1(L) ;If exponent > 0
JRST ZERO ; result is 0
$LCALL ROV ;Else overflow
HRLOI T0,377777 ;Store +biggest
HRLOI T1,377777
JRST RET
ZERO: SETZB T0,T1 ;Result is 0
JRST RET ;Return
;
;The following block of code deals with over/underflow in the
;square operation at LOOP:. Note that the "exponent" cannot be
;0 -- LOOP: is entered only if T2 is not 0. Moreover, if T2 is
;not 1 subsequent operations will aggravate the over/underflow
;condition in such a way that both the result of the iteration
;and its reciprocal will have the same exception as currently
;indicated. If, however, T2 = 1, and the square overflowed, it
;is possible that its reciprocal will be in range. We therefore
;complete the current pass through the loop, and if the LSH of T2
;makes it zero, we join the handling at OVER: for overflow/underflow
;on the MUL of T0 by T3. Note that no exception can occur on the
;MUL of T0 by a wrapped over/underflow of the square, so that the
;exception flags will still be valid after this step.
;
OVER2: GFMP T0,T3 ;No over/underflow. Hence flags
; from square of T3 still valid
LSH T2,-1 ;Discard low bit of exponent
JUMPE T2,OVER ;If T2 = 0, T0 has wrapped final
; result or its reciprocal
; which may be in range
;Final product surely
JSP T2,.+1 ;over/underflows. Get exception flags
TLNE T2,(1B11) ;If underflow flag set, reciprocal
JRST UNDER ; overflows. Go test sign of exponent
SKIPL @1(L) ;For overflow, if exponent > 0
JRST UNDMSG ; final result underflows.
JRST OVMSG ;Else reciprocal gives overflow
;
;The rest of the code handles over/underflow on the product of
;T0 by T3 and calculation of the reciprocal, if this is done.
;
OVER: JSP T2,.+1 ;Get exception flags
TLNE T2,(1B11) ;If underflow flag set
JRST UNDER ; underflow on product
SKIPL @1(L) ;Else, overflow on result if
JRST OVMSG ; exponent > 0. Get message
DMOVE T2,T0 ;Copy result
DMOVE T0,ONE ;For exponent < 0, get
GFDV T0,T2 ;reciprocal of wrapped overflow
JOV RET ;Underflow impossible; overflow
; compensates previous overflow
JRST UNDMSG ;Else, get underflow message
UNDER: SKIPL @1(L) ;Product underflowed. If exponent
JRST UNDMSG ; >/= 0, result underflows
;Else reciprocal overflows
OVMSG: $LCALL ROV ;Result overflow
JUMPL T0,NEGOV ;If result > 0
HRLOI T0,377777 ;Store +BIGGEST
HRLOI T1,377777
JRST RET ; and return
NEGOV: MOVSI T0,400000 ;If result < 0, store -BIGGEST
MOVEI T1,1
JRST RET ; and return
UNDMSG: $LCALL RUN ;Result underflow
SETZ T0,
RET: DMOVE T2,SAVE2 ;Restore T2-T4
MOVE T4,SAVE4
POPJ P,
ONE: EXP 200140000000,0 ;G-floating 1.0
SEGMENT DATA
SAVE2: BLOCK 2 ;Temp for T2-T3
SAVE4: BLOCK 1 ;Temp for T4
PRGEND
TITLE GEXP3. POWER FUNCTION
; (DOUBLE PRECISION)
SUBTTL IMSL, INC. JUNE 4, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;GEXP3 CALCULATES X**Y WHERE X AND Y ARE FLOATING POINT VALUES IN THE
;FOLLOWING RANGES
; 0.0 < X < + MACHINE INFINITY (X MAY EQUAL 0.0 IF Y > 0.0 AND
; X MAY BE LESS THAN 0.0 IF Y IS AN INTEGER. IF X IS < 0.0
; AND Y IS NOT AN INTEGER, A WARNING ERROR IS ISSUED AND
; DABS(X)**Y IS CALCULATED.)
; -1024.9375 <= FLOAT(INT((Y*LOG2(X))*16))/16 < 1023.0
; X**Y IS CALCULATED AS 2**W WHERE W = Y*LOG2(X). LOG2(X) IS
; CALCULATED AS FOLLOWS;
; X = F*(2**M), 1/2 <= F < 1. PICK P SUCH THAT P IS AN ODD
; INTEGER < 16 AND LET A = 2**(-P/16). NOW X = ((2**M)*A)*(F/A)
; LOG2(X) = M+LOG2(A) + LOG2(F/A) OR
; LOG2(X) = M-(P/16) + LOG2(F/A) .
; LET U1 = M-(P/16) AND
; U2 = LOG2(F/A) = LOG2((1+S)/(1-S)).
; THEN LOG2(X) = U1 + U2.
; AND S = (F-A)/(F+A). A RATIONAL
; APPROXIMATION IS USED TO EVALUATE U2. U1 AND U2 ARE THEN
; USED TO DETERMINE W1 AND W2 WHERE W=W1+W2
; AND W1 = FLOAT(INT(W*16.0))/16.0. FINALLY Z=X**Y=2**W
; IS RECONSTRUCTED AS Z = (2**W1) * (2**W2) WHERE
; W1 = M1-P1/16 AND 2**W2 IS EVALUATED FROM ANOTHER RATIONAL
; FUNCTION.
;
;THE RANGE OF DEFINITION FOR GEXP3 IS GIVEN ABOVE, AND ERROR MESSAGES
; WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE
;REQUIRED (CALLED) ROUTINES: NONE
;REGISTERS T2, T3, T4, T5, P1, P2, P3, AND P4 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; MOVEI L,ARG
; PUSHJ P,EXP3
;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN T1
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (GEXP3,.) ;ENTRY TO GEXP3. ROUTINE
PUSH P,T2 ;SAVE ACCUMULATORS
PUSH P,T3
DMOVE T0,@(L) ;GET THE BASE
DMOVE T2,@1(L) ;GET THE EXPONENT
SETZM IFLAG ;SET INTEGER EXP FLAG TO 0
JUMPG T0,XOK ;IF X IS NOT GT 0
JUMPE T0,X0 ;AND IF X IS NOT = 0
PUSH P,T4 ;SAVE ACCUMULATORS
PUSH P,T5
PUSH P,P1
PUSH P,P2
DMOVE T4,T2 ;GET A COPY OF Y
JUMPGE T4,YP ;IF Y IS NEGATIVE
DMOVN T4,T4 ;NEGATE IT
YP: MOVE P2,T4 ;GET HIGH ORDER OF Y
SETZ P1, ;SET P1 TO ZERO
LSHC P1,14 ;GET EXPONENT OF Y
SUBI P1,1765 ;GET SHIFTING FACTOR
LSHC T4,(P1) ;SHIFT OFF EXP AND PART OF INTEGER
TLNE T4,400000 ;IF Y IS ODD
SETOM IFLAG ;SET IFLAG TO ONES
LSHC T4,1 ;SHIFT OFF LAST OF INTEGER PART
DMOVN T0,T0 ;NEGATE X
JUMPN T5,MSTK ;IF LO OF Y .NE. 0, GO TO MSTK
JUMPE T4,YINT ;IF HI OF Y = 0, GO TO YINT
MSTK: $LCALL NNA
;LERR (LIB,%,<GEXP3: negative ** non-integer; ABS(base) used instead of base>)
JRST CONT ;GO TO CONT
X0: JUMPG T2,RET1 ;IF Y IS NOT GT 0
JUMPE T2,ZERZER ;SEPARATE OUT 0**0 CASE
$LCALL ZNI
;LERR (LIB,%,<GEXP3: base is 0 and exponent is le 0; result = infinity>)
HRLOI T0,377777 ;RESULT = INFINITY
HRLOI T1,377777
POPRET: POP P,T3 ;RESTORE ACCUMULATORS
POP P,T2
GOODBY (1) ;RETURN
ZERZER: $LCALL ZZZ ;"0**0 undefined, result = 0"
JRST POPRET ;GO RETURN
XOK: PUSH P,T4 ;SAVE ACCUMULATORS
PUSH P,T5
PUSH P,P1
PUSH P,P2
YINT: JUMPN T3,CONT ;IF LO OF Y IS NOT 0, GO TO CONT
JUMPN T2,YONE ;IF Y .NE. 0 GO TO YONE, OTHERWISE
DMOVE T0,A1 ;SET RESULT TO 1.0
JRST RET2 ;GO TO RET2
YONE: CAMN T2,A1 ;IF Y = 1.0
JRST RET2 ;GO TO RET2
CONT: PUSH P,P3
PUSH P,P4
MOVE T4,T0 ;OBTAIN THE EXPONENT
ASH T4,-30 ;SHIFT MANTISSA OFF
SUBI T4,2000 ;SUBTRACT 1024 FROM EXPONENT
AND T0,MASK1 ;EXTRACT MANTISSA
IOR T0,MASK2 ;SET EXPONENT TO 0
;THE FOLLOWING TESTS DETERMINE P
MOVEI T5,2 ;NP = 1
CAMLE T0,A1+20 ;IF F GT A1(9)
JRST NXT2 ;THEN GO TO NXT2
CAME T0,A1+20 ;
JRST OK1 ;
CAMLE T1,A1+21 ;
JRST NXT2
OK1: ADDI T5,20 ;
NXT2: CAMLE T0,A1+6(T5)
JRST NXT4 ;
CAME T0,A1+6(T5) ;
JRST OK2 ;
CAMLE T1,A1+7(T5) ;
JRST NXT4 ;
OK2: ADDI T5,10 ;
NXT4: CAMLE T0,A1+2(T5) ;
JRST NXT3 ;
CAME T0,A1+2(T5)
JRST OK3 ;
CAMLE T1,A1+3(T5) ;
JRST NXT3 ;
OK3: ADDI T5,4 ;
NXT3: MOVEM T5,PTEMP ;SAVE A COPY OF P
DMOVE P1,T0 ;SAVE A COPY OF F
GFAD P1,A1(T5) ;F+A1(2*P+2)
GFSB T0,A1(T5) ;F-A1(2*P+2)
SUBI T5,2 ;FORM [(2*P+2)-
ASH T5,-1 ;2/2]
GFSB T0,A2(T5) ;Z=(F-A1(P+1))-A2((P+1)/2)
GFDV T0,P1 ;/(F+A1(P+1))
EXTEND T0,[GFSC 1] ;
DMOVE P1,T0 ;SAVE A COPY OF Z
GFMP P1,P1 ;FORM Z**2
DMOVE P3,P1 ;SAVE A COPY OF Z**2
GFMP P3,RP4 ;R(Z)=RP4*Z**2
GFAD P3,RP3 ;+RP3
GFMP P3,P1 ;*Z**2
GFAD P3,RP2 ;+RP2
GFMP P3,P1 ;*Z**2
GFAD P3,RP1 ;+RP1
GFMP P3,P1 ;*Z**2
GFMP P3,T0 ;*Z
GFAD P3,T0 ;+Z
GFMP P3,C ;U2 = R(Z) * C
ASH T4,4 ;U1 = M*16
MOVE T5,PTEMP ;GET A COPY OF P
ASH T5,-1
SUB T4,T5 ;-P
FLTR T4,T4 ;FLOAT IT
SETZ T5, ;
EXTEND T4,[GDBLE T4] ;
EXTEND T4,[GFSC -4] ;
DMOVE T0,T2 ;SAVE A COPY OF Y
JUMPGE T2,YPOS ;IF Y IS NEGATIVE
DMOVN T0,T0 ;NEGATE IT
YPOS: MOVE P1,T0 ;GET COPY OF HIGH ORDER OF Y
ASH P1,-30 ;SHIFT OFF MANTISSA
CAIGE P1,2024 ;IF THE EXPONENT IS LESS THAN 2024
JRST SMALLY ;GO TO SMALLY
CAIE P1,2024 ;IF THE EXPONENT IS > 2024
JRST LARGEY ;GO TO LARGEY
SETZ T1, ;ZERO SECOND WORD
JRST SGNCHK ;GO TO SGNCHK
SMALLY: SUBI P1,1775 ;GET SHIFTING FACTOR
SETZ T1, ;SET Y1 TO 0
JUMPGE P1,GETY1 ;IF P1 IS .GE. 0
SETZ T0, ;THEN GO TO GETY1
JRST GETY2 ;
GETY1: AND T0,MSK3(P1) ;GET Y1
JRST SGNCHK ;GO TO CHECK SIGN
LARGEY: CAIL P1,2067 ;IF EXPONENT IS .GE. 2067
JRST SGNCHK ;GO TO SGNCHK
SUBI P1,2025 ;GET SHIFTING FACTOR
AND T1,REDMSK(P1) ;GET LOW ORDER PART OF Y1
SGNCHK: JUMPGE T2,GETY2 ;IF Y IS NEGATIVE
DMOVN T0,T0 ;NEGATE Y1
GETY2: DMOVE P1,T2 ;SAVE A COPY OF Y
GFSB T2,T0 ;Y2 = Y-Y1
GFMP T2,T4 ;FORM Y2*U1
JFCL
GFMP P1,P3 ;FORM Y*U2
JFCL
GFAD T2,P1 ;W = Y*U2+Y2*U1
GFMP T4,T0 ;FORM U1*Y1
JFCL
DMOVE T0,T2 ;RECONSTRUCT W
GFAD T0,T4
JFCL
CAMG T0,BIGW ;IF W IS NOT TOO BIG
JRST WOK ;GO TO WOK
OVFL: $LCALL ROV
;LERR (LIB,%,<GEXP3: result overflow>)
HRLOI T0,377777 ; RESULT = INFINITY
HRLOI T1,377777
JRST RET3 ;GO TO RET3
WOK: CAML T0,SMALLW ;IF W IS NOT TOO SMALL
JRST WOK2 ;THEN PROCEED
UNFL: $LCALL RUN
;LERR (LIB,%,<GEXP3: result underflow>)
SETZ T0, ; RESULT = 0
SETZ T1,
POP P,P4 ;RESTORE ACCUMULATORS
POP P,P3
POP P,P2
POP P,P1
POP P,T5
POP P,T4
POP P,T3
POP P,T2
GOODBY (1) ;RETURN
WOK2: SETZ P2, ;ZERO P2
MOVM P1,T2 ;GET HIGH ORDER OF W
MOVE P3,P1 ;
ASH P3,-30 ;SHIFT OFF MANTISSA
SUBI P3,1775 ;GET SHIFTING FACTOR
JUMPGE P3,GETW1 ;IF P3 .GE. 0
;THEN GO TO GETW1
SETZ P1, ;OTHERWISE, SET W1 = 0
JRST GETW2 ;GO TO GETW2
GETW1: AND P1,MSK3(P3) ;GETW1
JUMPG T2,GETW2 ;IF W IS NEGATIVE
DMOVN P1,P1 ;NEGATE W1
GETW2: GFSB T2,P1 ;W2 = W-W1
GFAD T4,P1 ;W = W1+U1*Y1
MOVM T0,T4 ;SAVE A COPY OF ABS(W)
MOVE P1,T0 ;
SETZ T1, ;ZERO T1
ASH P1,-30 ;SHIFT OFF MANTISSA
SUBI P1,1775 ;GET SHIFTING FACTOR
JUMPGE P1,GTW1 ;IF P1 IS .GE. 0
;THEN GO TO GTW1
SETZ T0, ;OTHERWISE, SET W1 TO 0
JRST GTW2 ;GO TO GET W2
GTW1: AND T0,MSK3(P1) ;GET W1
JUMPGE T4,GTW2 ;IF W IS NEGATIVE
DMOVN T0,T0 ;NEGATE W1
GTW2: GFSB T4,T0 ;FORM W-W1
GFAD T2,T4 ;W2 = W2+(W-W1)
MOVM T4,T2 ;SAVE A COPY OF ABS(W2)
MOVE P1,T4 ;
SETZ T5, ;ZERO T5
ASH P1,-30 ;SHIFT OFF MANTISSA
SUBI P1,1775 ;GET SHIFTING FACTOR
JUMPGE P1,GW ;IF P1 IS .GE. 0
;THEN GO TO GW
SETZ T4, ;OTHERWISE, SET W=0
JRST GW2 ;GO TO GW2
GW: AND T4,MSK3(P1) ;GET W
JUMPGE T2,GW2 ;IF W2 IS NEGATIVE
DMOVN T4,T4 ;NEGATE W
GW2: GFAD T0,T4 ;FORM W1 + W
EXTEND T0,[GSNGL T0] ;
FSC T0,4 ;*16
FIX T0,T0 ;IW1
GFSB T2,T4 ;W1 = W2 - W
JUMPLE T2,W2POS ;IF W2 .GT. 0
GFSB T2,SXTNTH ;W2 = W2-.0625
ADDI T0,1 ;IW1 = IW1+1
W2POS: MOVE T5,T0 ;SAVE A COPY OF IW1
JUMPGE T5,NPOS
ADDI T5,17
NPOS: ASH T5,-4 ;M1 = IW1/16
JUMPL T0,INEG ;IF IW1 .GE. 0
ADDI T5,1 ;M1 = M1+1
INEG: MOVE T4,T5 ;SAVE A COPY OF M1
ASH T4,4 ;P1 = 16*M1
SUB T4,T0 ;-IW1
ASH T4,1 ;
DMOVE T0,T2 ;SAVE A COPY OF W2
GFMP T2,Q7 ;Z = Q7*W2
GFAD T2,Q6 ;+Q6
GFMP T2,T0 ;*W2
GFAD T2,Q5 ;+Q5
GFMP T2,T0 ;*W2
GFAD T2,Q4 ;+Q4
GFMP T2,T0 ;*W2
GFAD T2,Q3 ;+Q3
GFMP T2,T0 ;*W2
GFAD T2,Q2 ;+Q2
GFMP T2,T0 ;*W2
GFAD T2,Q1 ;+Q1
GFMP T0,T2 ;*W2
GFMP T0,A1(T4) ;*A1(P1+1)
GFAD T0,A1(T4) ;+A1(P1+1)
EXTEND T0,[GFSC (T5)] ;ADD M1 TO THE EXP OF Z
JFCL EXCPT
RET3: POP P,P4 ;RESTORE ACCUMULATORS
POP P,P3
RET2: SKIPE IFLAG ;IF Y IS ODD NEGATIVE INTEGER
DMOVN T0,T0 ;NEGATE RESULT
POP P,P2
POP P,P1
POP P,T5
POP P,T4
RET1: POP P,T3
POP P,T2
GOODBY (1) ;RETURN
EXCPT: JUMPG T5,OVFL ;IF T5 >0 GO TO OVFL
JRST UNFL ;OTHERWISE, GO TO UNFL
SXTNTH: DOUBLE 177540000000,000000000000 ;.0625D0
BIGW: 201277740000 ;UPPER BOUND FOR WW
SMALLW: 576440000000 ;LOWER BOUND FOR W
RP1: DOUBLE 177552525252,252525251443 ;.833333333333332114D-1
RP2: DOUBLE 177263146314,314740401603 ;.125000000005037992D-01
RP3: DOUBLE 177044444441,161170665342 ;.223214212859242590D-2
RP4: DOUBLE 176570743756,332257605031 ;.434457756721631196D-3
C: DOUBLE 200156125073,051270137606 ;1.442695040888963407D0
Q7: DOUBLE 176076473357,034454617751 ;.149288526805956082D-4
Q6: DOUBLE 176450275727,001604616375 ;.154002904409897646D-3
Q5: DOUBLE 176753541760,306574663707 ;.133335413135857847D-02
Q4: DOUBLE 177247312533,160461663013 ;.961812905951724170D-02
Q3: DOUBLE 177470654106,270060115477 ;.555041086640855953D-1
Q2: DOUBLE 177675376757,374054231615 ;.240226506959095371D0
Q1: DOUBLE 200054271027,367643475706 ;.693147180559945296D0
A1: DOUBLE 200140000000,000000000000 ;A1(I), I=1,17 =
A12: DOUBLE 200075222575,025111033141 ;2**((1-I)/16). THIS
A13: DOUBLE 200072540306,347672220712 ;TABLE IS SEARCHED
A14: DOUBLE 200070146336,354125123411 ;TO DETERMINE P.
A15: DOUBLE 200065642374,312655165530 ;
A16: DOUBLE 200063422214,025077022007 ;
A17: DOUBLE 200061263452,021252033327 ;
A18: DOUBLE 200057204243,237260060666 ;
A19: DOUBLE 200055202363,063763571444 ;
A110: DOUBLE 200053254076,352205205126 ;
A111: DOUBLE 200051377326,251542504707 ;
A112: DOUBLE 200047572462,140443204215 ;
A113: DOUBLE 200046033760,121433342514 ;
A114: DOUBLE 200044341723,163526107032 ;
A115: DOUBLE 200042712701,343725057267 ;
A116: DOUBLE 200041325303,147630441731 ;
A117: DOUBLE 200040000000,000000000000 ;
A2: DOUBLE 170461734720,307535546011
A22: DOUBLE 607304062611,120221564632
A23: DOUBLE 170276106540,343712762662
A24: DOUBLE 607625040217,333155037217
A25: DOUBLE 170362230025,031075315002
A26: DOUBLE 170466404421,360475356413
A27: DOUBLE 607330176555,216025626545
A28: DOUBLE 607323056225,270604125171
MASK1: 000077777777 ;MASK FOR MANTISSA
MASK2: 200000000000 ;MASK FOR EXPONENT
REDMSK: 600000000000
RDMSK2: 700000000000
RDMSK3: 740000000000
RDMSK4: 760000000000
RDMSK5: 770000000000
RDMSK6: 774000000000
RDMSK7: 776000000000
RDMSK8: 777000000000
MASK: 777400000000 ;MASKS FOR FINDING W1
MSK1: 777600000000 ;
MSK2: 777700000000 ;
MSK3: 777740000000 ;
MSK4: 777760000000 ;
MSK5: 777770000000 ;
MSK6: 777774000000 ;
MSK7: 777776000000 ;
MSK8: 777777000000 ;
MSK9: 777777400000 ;
MSK10: 777777600000 ;
MSK11: 777777700000
MSK12: 777777740000
MSK13: 777777760000
MSK14: 777777770000
MSK15: 777777774000
MSK16: 777777776000
MSK17: 777777777000
MSK18: 777777777400
MSK19: 777777777600
MSK20: 777777777700
MSK21: 777777777740
MSK22: 777777777760
MSK23: 777777777770
MSK24: 777777777774
MSK25: 777777777776
MSK26: 777777777777
SEGMENT DATA
PTEMP: 0
IFLAG: 0
PRGEND
TITLE GLOG10 LOG BASE 10 FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. APRIL 8, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GLOG10
EXTERN GLG10.
GLOG10=GLG10.
PRGEND
TITLE GLOG NATURAL LOG FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. APRIL 8, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GLOG
EXTERN GLOG.
GLOG=GLOG.
PRGEND
TITLE GLOG. NATURAL LOG FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. APRIL 8, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DLOG10(X) AND DLOG(X) ARE CALCULATED AS FOLLOWS
; FOR X = 0 AN ERROR MESSAGE IS ISSUED AND -MACHINE INFINITY
; IS RETURNED AS THE RESULT.
; FOR X < 0 AN ERROR MESSAGE IS ISSUED, X IS SET T0 -X AND
; CALCULATION CONTINUES.
; FOR X > 0, X = F*2**F(M), 1/2 < F < 1
; DEFINE G AND N SO THAT F = G*2(-N), 1/SQRT(2) <= G < SQRT(2).
; NOW
; DLOG(X) = (K*M-N) * DLOG(2) + DLOG(G)
; AND
; DLOG10(X) = DLOG10(E) * DLOG(X) = DLOG(X)/DLOG(10)
;
; DLOG(G) IS EVALUATED BY DEFINING S = (G-1)/(G+1) AND Z = 2*S
; AND THEN CALCULATING DLOG(G) = DLOG((1+Z/2)/(1-Z/2)) USING
; A MINIMAX RATIONAL APPROXIMATION.
;THE RANGE OF DEFINITION FOR DLOG/DLOG10 IS GIVEN ABOVE, AND ERROR MESSAGES
; WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE
;REQUIRED (CALLED) ROUTINES: NONE
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; MOVEI L,ARG
; PUSHJ P,GLOG
; OR
; PUSHJ P,GLOG10
;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0.
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1.
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (GLG10.,GLOG10) ;ENTRY TO LOG BASE 10 ROUTINE.
SETZM FLAG ;CLEAR FLAG FOR DLOG10 ENTRY
JRST ALG ;GO TO ALGORITHM
HELLO (GLOG,.) ;ENTRY TO NATURAL LOG ROUTINE
SETOM FLAG ;SET FLAG FOR DLOG ENTRY
ALG: DMOVE T0,@(L) ;GET ARG
JUMPG T0,STRT ;IF ARG > ZERO GO TO STRT
JUMPN T0,ARGN ;OTHERWISE, IF ARG .NE. 0 GO TO ARGN
$LCALL AZM
;LERR (LIB,%,<DLOG or DLOG10: zero arg; result = -infinity>)
HRLZI T0,400000 ;SET RESULT TO
HRRZI T1,000001 ;LARGE NEGATIVE NUMBER
GOODBY (1)
ARGN: $LCALL NAA
;LERR (LIB,%,<DLOG or DLOG10: negative arg; result = LOG(ABS(arg))>)
DMOVN T0,T0 ;ARG = -ARG
STRT: PUSH P,T2 ;SAVE ACCUMULATORS
PUSH P,T3
PUSH P,T4
PUSH P,T5
HLRZ T2,T0 ;LEFT OF T0 TO RIGHT OF T2
LSH T2,-6 ;ISOLATE EXPONENT
SUBI T2,2000 ;SUBTRACT 1024 FROM THE EXP
DMOVE T4,T0 ;GET COPY OF ARG
AND T4,MASK1 ;EXTRACT MANTISSA
IOR T4,MASK2 ;SET EXP TO 0
CAMLE T4,HI ;IS HIGH PART OF F > SQRT(.5)
JRST FGT ;YES, GO TO FGT
CAME T4,HI ;NO, IS F = HIGH OF SQRT(.5)
JRST FLT ;NO, GO TO FLT
CAMLE T5,LOW ;YES, IS LOW PART > LOW PART OF SQRT(.5)
JRST FGT ;YES, GO TO FGT
FLT: SUBI T2,1 ;N = N-1
EXTEND T2,[GFLTR T2]
MOVEM T2,N ; SAVE N
GFAD T4,MHALF ;ZNUM = F-.5
DMOVE T2,T4 ;GET COPY OF ZNUM
EXTEND T4,[GFSC -1] ;ZDEN = ZNUM * .5
GFAD T4,HALF ; + .5
JRST EVALRZ
FGT: EXTEND T2,[GFLTR T2]
MOVEM T2,N ; SAVE N
DMOVE T2,T4
GFAD T2,B3 ;ZNUM = F - 1.0
EXTEND T4,[GFSC -1] ;ZDEN = F*.5
GFAD T4,HALF ; + .5
EVALRZ: GFDV T2,T4 ;Z = ZNUM/ZDEN
DMOVE T4,T2 ;
GFMP T4,T4 ;W = Z*Z
DMOVE T0,T4 ;SAVE COPY OF W
GFAD T4,B2 ; FORM B(W). B(W) = W + B2
GFMP T4,T0 ; * W
GFAD T4,B1 ; + B1
GFMP T4,T0 ; * W
GFAD T4,B0 ; + B0
DMOVEM T0,SAVEW ; SAVE A COPY OF W
GFMP T0,A2 ;FORM A(W). A(W)= A2*W
GFAD T0,A1 ; + A1
GFMP T0,SAVEW ; * W
GFAD T0,A0 ; + A0
GFDV T0,T4 ;R(Z) = A(W)/B(W)
GFMP T0,SAVEW ; * W
GFMP T0,T2 ; *Z
GFAD T0,T2 ; + Z
MOVE T2,N ;RETRIEVE N
MOVEI T3,0 ;ZERO OUT SECOND WORD
GFMP T2,C1 ;FORM N*C1
MOVE T4,N ;RETRIEVE N
MOVEI T5,0 ;ZERO OUT SECOND WORD
GFMP T4,C2 ;FORM N*C2
GFAD T0,T4 ;RESULT = N*C2 + R(Z)
GFAD T0,T2 ; + N*C1
SKIPN FLAG
GFMP T0,C3 ;IF DLOG10 ROUTINE, RESULT = RESULT*C3
POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4
POP P,T3
POP P,T2
GOODBY (1)
C1: DOUBLE 200054300000,000000000000 ;.693359375D0
C2: DOUBLE 601310277575,034757152744 ;-2.12194440054690583D-4
C3: DOUBLE 177767455730,251156241615 ;.434294481903251828D0
A0: DOUBLE 577037740007,152304514557 ;-.641249434237455811D2
A1: DOUBLE 200540611121,000552775450 ;.163839435630215342D2
A2: DOUBLE 577715357522,145224132710 ;-.789561128874912573D0
B0: DOUBLE 576517720013,037446761043 ;-.769499321084948798D3
B1: DOUBLE 201147002037,320522317573 ;.312032220919245328D3
B2: DOUBLE 577134251775,244603076112 ;-.356679777390346462D2
B3: DOUBLE 577640000000,000000000000 ;-.1D1
HALF: DOUBLE 200040000000,000000000000 ;.5D0
MHALF: DOUBLE 577740000000,000000000000 ;-.5D0
HI: 200055202363
LOW: 063763571444
MASK1: 000077777777
MASK2: 200000000000
SEGMENT DATA
N: 0
FLAG: 0
SAVEW: DOUBLE 000000000000,000000000000
PRGEND
TITLE GMOD DOUBLE PRECISION REMAINDER FUNCTION
SUBTTL MARY PAYNE /MHP/CKS 25-Jan-80
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GMOD
EXTERN GMOD.
GMOD=GMOD.
PRGEND
TITLE GMOD. DOUBLE PRECISION REMAINDER FUNCTION
SUBTTL MARY PAYNE /MHP/CKS 25-Jan-80
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
SEGMENT CODE
SALL
T6=6 ;EXTRA TEMP REGISTERS
T7=7
T8=10
T9=11
HELLO (GMOD,.)
DMOVEM T2,SAVE2 ;SAVE TEMP REGISTERS
DMOVEM T4,SAVE4
DMOVEM T6,SAVE6
DMOVEM T8,SAVE8
DMOVE T0,@0(L) ;GET ABS(A) IN T0-T1
CAIGE T0,0
DMOVN T0,T0
REPEAT: DMOVE T2,@1(L) ;GET ABS(B) IN T2-T3
CAIGE T2,0
DMOVN T2,T2
CAMG T0,T2 ;CHECK IF A .LE. B
JRST FAST ;YES, GO BE FAST
DMOVEM T2,B ;SAVE ABS(B)
DMOVE T8,T0 ;SAVE ABS(A) IN T8-T9
;BREAK B INTO HIGH, MIDDLE, AND LOW PARTS
;B IN T2/ XXXXHHHHHHHH MMMMMMMMLLLL
MOVE T6,T2 ;GET T6/ XXXXHHHHHHHH 000000000000
SETZ T7, ;
GFSB T2,T6 ; T2/ XXXXMMMMMMMM LLLL00000000
MOVE T4,T2 ; T4/ XXXXMMMMMMMM 000000000000
SETZ T5, ;
GFSB T2,T4 ; T2/ XXXXLLLL0000 000000000000
GFDV T8,B ;GET A/B
JFCL OVFL ;CATCH OVERFLOW
CAML T8,[204440000000] ;IS QUOTIENT OVER 2**35?
JRST [AND T9,[777700000000] ;YES, TRUNCATE IT TO 35 BITS
JRST BIG] ; TO KEEP PRODUCTS EXACT
EXTEND T9,[GFIX T8] ;UNDER 2**35, TRUNCATE TO NEAREST INTEGER
EXTEND T8,[GFLTR T9] ;DFLOAT THE INTEGER
BIG: GFMP T6,T8 ;GET A-B*[A/B] IN T0-T1
GFSB T0,T6
GFMP T4,T8
GFSB T0,T4
GFMP T2,T8
GFSB T0,T2
OVCONT: CAIGE T0,0 ;IF RESULT NEGATIVE, [A/B] WAS 1 TOO HIGH
GFAD T0,B ;CORRECT THE RESULT
CAMGE T0,B ;IF RESULT .LT. B, OK
JRST MRET
CAMG T0,B ;ELSE [A/B] WAS TOO LOW
CAML T1,B+1 ; AND MUST SUBTRACT SOME MORE MULTIPLES OF B
JRST REPEAT ;GO GET T0 MOD B IN T0
MRET: SKIPGE @0(L) ;GIVE RESULT SIGN OF A
DMOVN T0,T0
MRET1: DMOVE T8,SAVE8 ;RESTORE REGISTERS
DMOVE T6,SAVE6
DMOVE T4,SAVE4
DMOVE T2,SAVE2
POPJ P, ;DONE
FAST: CAMN T0,T2 ;HERE IF A HIGH .LE. B HIGH
CAMGE T1,T3 ;HIGH WORDS EQUAL, COMPARE LOW WORDS
JRST MRET ;A .LT. B, GO RETURN DMOD = A
GFSB T0,T2 ;ELSE RETURN DMOD = A - B
JRST MRET
OVFL: JUMPE T6,RTZ ;IF B IS ZERO, DMOD IS 0
EXTEND T0,[GFSC -1777] ;SCALE A TO AVOID OVERFLOW IN DIVIDE
DMOVE T8,T0 ;GET A
GFDV T8,B ;GET A/B
AND T9,[777700000000] ;TRUNCATE TO 35 BITS SO PRODUCTS ARE EXACT
GFMP T6,T8 ;GET A-B*[A/B]
GFSB T0,T6
GFMP T4,T8
GFSB T0,T4
GFMP T2,T8
GFSB T0,T2
EXTEND T0,[GFSC 1777] ;UNDO SCALING
JRST OVCONT ;CONTINUE
RTZ: SETZB T0,T1 ;A .EQ. B, RETURN DMOD = 0
JRST MRET1
SEGMENT DATA
B: BLOCK 2 ;ABS(B)
SAVE2: BLOCK 2 ;TEMP REGISTERS
SAVE4: BLOCK 2
SAVE6: BLOCK 2
SAVE8: BLOCK 2
PRGEND
TITLE GCOS SINE AND COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MAY 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GCOS
EXTERN GCOS.
GCOS=GCOS.
PRGEND
TITLE GSIN SINE AND COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MAY 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GSIN
EXTERN GSIN.
GSIN=GSIN.
PRGEND
TITLE GSIN. SINE AND COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MAY 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;COS(X) AND SIN(X) ARE CALCULATED AS FOLLOWS
; NOTE THAT COS(X) = COS(-X) AND SIN(X) = -SIN(-X)
; LET ABS(X) = N*PI + F WHERE ABS(F) <= PI/2
; WHEN ABS(F) < EPS (SEE CONSTANTS, BELOW), SIN(F) = F
; OTHERWISE, THE ARGUMENT REDUCTION IS:
; LET X1 = ABS(X)
; N = THE NEAREST INTEGER TO ABS(X), FOR SIN, OR
; TO ABS(X) + PI/2, FOR COS
; THEN XN = [N/PI], THE GREATEST INTEGER IN N/PI.
; THEN THE REDUCED ARGUMENT F = (((X1 - XN*C1) -XN*C2) -XN*C3)
; WHERE C1+C2+C3 = PI TO EXTRA PRECISION AND ARE GIVEN BELOW.
; LET G = F**2
; THEN R(G) = (G*XNUM/XDEN+RP1)*G
; WHERE XNUM = ((RP5*G+RP4)*G+RP3)*G+RP2
; XDEN = ((G*Q2)*G+Q1)*G+Q0
; AND RP5,RP4,RP3,RP2,RP1,Q2,Q1, AND Q0 ARE GIVEN BELOW
; AND SIN(F) = F + F*R(G). THE R(I) APPEAR BELOW
; FINALLY,
; SIN(X) = SIGN(X)*SIN(F)*(-1)**N, AND
; COS(X) = SIN(X + PI/2)
;THE RANGE OF DEFINITION FOR SIN IS ABS(X) <= 1686629713. ABS(X)+PI/2, IN
; COS(X), MUST BE LESS THAN 1686629713. SIN(X) = COS(X) = 0.0 AND AN
; ERROR MESSAGE WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE
;REQUIRED (CALLED) ROUTINES: NONE
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; MOVEI L,ARG
; PUSHJ P,GSIN
; OR
; PUSHJ P,GCOS
;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0.
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1.
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (GCOS,.) ;ENTRY TO COSINE ROUTINE
DMOVE T0,@(L) ;OBTAIN ARGUMENT
PUSH P,T2 ;SAVE AN ACCUMULATOR
HRLZI T2,201400 ;SET SGN TO 1.
JUMPGE T0,XPOS ;IF X IS NEGATIVE
DMOVN T0,T0 ;Y = -X
XPOS: DMOVEM T0,XDEN ;SAVE A COPY OF ABS(X)
SETOM FLAG ;SET FLAG FOR COSINE ENTRY
GFAD T0,PID2 ;Y = Y+PID2
JRST ALG ;PROCEED TO MAIN ALGORITHM
HELLO (GSIN,.) ;ENTRY TO SIN ROUTINE
DMOVE T0,@(L) ;GET ARG
PUSH P,T2 ;SAVE AN ACCUMULATOR
HRLZI T2,201400 ;SET SIGN TO 1.0
SETZM FLAG ;CLEAR FLAG FOR SINE ROUTINE
JUMPGE T0,ALG ;IF ARG IS NEGATIVE
DMOVN T0,T0 ;THEN, SET Y=-X AND
HRLZI T2,576400 ;SET VALUE OF SIGN TO -1.
ALG: MOVEM T2,SGN ;MOVE SIGN TO MEMORY
CAMGE T0,YMAX1 ;IF HI OF Y<HI OF YMAX
JRST YOK ;PROCEED TO YOK
CAME T0,YMAX1 ;IF HI OF Y > HI OF YMAX
JRST ERR1 ;GO TO ERR1
CAMGE T1,YMAX2 ;HI OF Y=HI OF YMAX, IS LO OF Y < LO OF YMAX?
JRST YOK ;YES, PROCEED TO YOK
ERR1: $LCALL ATZ
;LERR (LIB,%,<DSIN or DCOS: ABS(arg) too large; result = zero>)
MOVEI T0,0 ;SET RESULT
MOVEI T1,0 ;TO ZERO
POP P,T2 ;RESTORE ACCUMULATOR
GOODBY (1) ;
YOK: PUSH P,T3 ;SAVE ACCUMULATORS
PUSH P,T4 ;
PUSH P,T5
DMOVE T2,T0 ;SAVE A COPY OF ABS(X)
SKIPE FLAG ;IF THIS IS THE COS ROUTINE
DMOVE T2,XDEN ;RESTORE ABS(X) TO T2
GFMP T0,ODPI ;RN = Y/PI
JFCL
EXTEND T0,[DGFIXR T0] ;FIX RN
MOVE T5,T1 ;SAVE N
EXTEND T0,[DGFLTR T0] ;XN=FLOAT(N)
TRNE T5,1 ;IS N ODD?
MOVNS SGN ;YES,NEGATE SIGN
SKIPE FLAG ;IF THE COSINE IS WANTED
GFAD T0,PT5 ;THEN XN=XN-.5
DMOVE T4,T0 ;SAVE A COPY OF XN
GFMP T0,C1 ;FORM XN*C1
GFSB T2,T0 ;F=ABS(X)-(XN*C1)
DMOVE T0,T4 ;SAVE A COPY OF XN
GFMP T0,C3 ;FORM XN*C3
GFMP T4,C2 ;FORM XN*C2
GFSB T2,T4 ;-XN*C2
GFSB T2,T0 ;F=F-XN*C3
DMOVE T0,T2 ;SAVE A COPY OF F
JUMPGE T2,FPOS ;IF F IS NEGATIVE
DMOVN T2,T2 ;F = -F
FPOS: CAML T2,EPS ;IF F IS .GE. EPS
JRST GEEPS ;GO TO GEEPS
SKIPGE SGN ;IF SGN IS NEGATIVE
DMOVN T0,T0 ;F = -F
POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4
POP P,T3
POP P,T2
GOODBY (1) ;RETURN
GEEPS: GFMP T2,T2 ;G = F*F
DMOVE T4,T2 ;SAVE A COPY OF G
GFAD T2,Q2 ;XDEN = G+Q2
GFMP T2,T4 ;*G
GFAD T2,Q1 ;+Q1
GFMP T2,T4 ;*G
GFAD T2,Q0 ;+Q0
DMOVEM T2,XDEN ;MOVE XDEN TO MEMORY
DMOVE T2,T4 ;GET A COPY OF G
GFMP T2,RP5 ;XNUM = G*RP5
GFAD T2,RP4 ;+RP4
GFMP T2,T4 ;*G
GFAD T2,RP3 ;+RP3
GFMP T2,T4 ;*G
GFAD T2,RP2 ;+RP2
GFMP T2,T4 ;*G
GFDV T2,XDEN ;R(G) = XNUM/XDEN
GFAD T2,RP1 ;+RP1
GFMP T2,T4 ;*G
GFMP T2,T0 ;*F
GFAD T0,T2 ;+F
SKIPGE SGN ;IF SGN IS NEGATIVE
DMOVN T0,T0 ;THEN NEGATE T0
RET: POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4 ;
POP P,T3
POP P,T2
GOODBY (1)
EPS: 174340000000 ;2**(-30)
PID2: DOUBLE 200162207732,242102643022 ;PI/2
ODPI: DOUBLE 177750574603,156234420251 ;1/PI
C1: DOUBLE 200262207732,240000000000 ;HIGH 30 BITS OF PI
C2: DOUBLE 174442055060,200000000000 ;NEXT 28 BITS OF PI
C3: DOUBLE 171064611431,212134015604 ;C1+C2+C3=PI TO 120 BITS
Q2: DOUBLE 201161273135,127076131616 ;0.394924723520450141 D+3
Q1: DOUBLE 202142232235,112027153730 ;0.702492288221842518D+5
Q0: DOUBLE 202751252025,266402115312 ;0.541748285645351853D+7
RP5: DOUBLE 600516152672,335644427111 ;-0.121560740596710190D-1
RP4: DOUBLE 200342202301,360464200740 ;0.428183075897778265D+01
RP3: DOUBLE 576602640645,001056761116 ;-0.489487151969463797D+03
RP2: DOUBLE 202054054660,302527505074 ;0.451456904704461990D+05
RP1: DOUBLE 600125252525,125252525242 ;-.166666666666666667D0
PT5: DOUBLE 577740000000,000000000000 ;-.5
YMAX1: 203762207732 ;YMAX =
YMAX2: 242102643021 ; (PI*2**29)
SEGMENT DATA
FLAG: 0
SGN: DOUBLE 0,0
XDEN: DOUBLE 0,0
PRGEND
TITLE GSINH HYPERBOLIC SINE FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 7, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GSINH
EXTERN GSINH.
GSINH=GSINH.
PRGEND
TITLE GSINH. HYPERBOLIC SINE FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 7, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;GSINH(X) IS CALCULATED AS FOLLOWS
; LET V BE APPROXIMATELY 2 SO THAT LNV AND ABS(X)+LN V
; CAN BE EXACTLY REPRESENTED WHEN X IS EXACTLY REPRESENTABLE.
; THEN, LETTING Y = ABS(X), AND NOTING THAT -SINH(-X)=SINH(X),
; FOR
; 0 <= Y < EPS, GSINH = Y*SIGN(X)
; EPS <= Y <= 1, GSINH = X-X*(Z*R(Z)), WHERE Z = Y**2 AND
; R(Z) IS GIVEN BELOW.
; 1 < Y <= 709.089565, GSINH = SIGN(X)*(EXP(Y)-EXP(-Y))/2,
; 709.089565 <= Y < 1024 * LN(2)
; GSINH = SIGN(X)*((V/2)*EXP(Y - LN V)),
; Y >= 1024 * LN(2), SINH = +MACHINE INFINITY * SIGN(X)
; LET Z = Y**2. THEN
; R(Z) = (RP0 + Z*(RP1 + Z*(RP2 + Z*RP3)))/(Q0 + Z*(Q1 + Z*(Q2 + Z)))
; WHERE RP0, RP1, RP2, Q0, Q1, AND Q2 ARE GIVEN BELOW.
;THE RANGE OF DEFINITION FOR GSINH IS ABS(X) < 1024 * LN(2) AND ERROR MESSAGES
; WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE. GSINH WILL BE SET
; TO + MACHINE INFINITY * SIGN(X).
;REQUIRED (CALLED) ROUTINES: DEXP
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; MOVEI L,ARG
; PUSHJ P,GSINH
;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (GSINH,.) ;ENTRY TO GSINH ROUTINE
DMOVE T0,@(L) ;OBTAIN X
PUSH P,T5 ;SAVE REGISTER T5
HRLZI T5,200140 ;SET FLAG TO 1
JUMPGE T0,XPOS ;IF X IS NEGATIVE
MOVN T5,T5 ;NEGATE FLAG
DMOVN T0,T0 ;Y = -X
XPOS: CAMLE T0,ONE ;IF Y IS .GT. 1.0
JRST DCSH ;GO TO DCSH
LE: CAMG T0,TWOM30 ;IF Y IS .LE. 2**-30
JRST RET2 ;GO TO RET1
JUMPGE T5,NXT ;IF X IS NEGATIVE
DMOVN T0,T0 ;Y = X
NXT: PUSH P,T2 ;SAVE ACCUMULATORS
PUSH P,T3
PUSH P,T4
DMOVE T2,T0 ;SAVE A COPY OF X
DMOVNM T0,TEMP ;MOVE A COPY OF -X TO MEMORY
GFMP T2,T2 ;Z = X*X
DMOVE T4,T2 ;SAVE A COPY OF Z
GFAD T4,Q2 ;XDEN = Z + Q2
GFMP T4,T2 ;*Z
GFAD T4,Q1 ;+Q1
GFMP T4,T2 ;*Z
GFAD T4,Q0 ;+Q0
DMOVEM T2,Z ;MOVE A COPY OF Z TO MEMORY
GFMP T2,RP3 ;XNUM = Z*RP3
GFAD T2,RP2 ;+RP2
GFMP T2,Z ;*Z
GFAD T2,RP1 ;+RP1
GFMP T2,Z ;*Z
GFAD T2,RP0 ; +RP0
GFDV T2,T4 ;R(Z) = XNUM/XDEN
GFMP T2,Z ;*Z
GFMP T2,TEMP ;*(-X)
GFAD T0,T2 ;+X
POP P,T4 ;RESTORE ACCUMULATORS
POP P,T3
POP P,T2
POP P,T5
GOODBY (1) ;RETURN
DCSH: CAMGE T0,TWENT2 ;IF ABS(X) < 22, GO TO
JRST ALG2 ; FULL CALCULATION
CAMG T0,BIGX ;IF Y IS .LE. BIGX
JRST EASY ;THEN GO TO EASY
CAMGE T0,YMAXHI ;IF HI OF Y < YMAXHI
JRST GETW ; GO TO GETW
CAME T0,YMAXHI ;IF HI OF Y > YMAXHI
JRST MSTK ;GO TO MSTK
CAMGE T1,YMAXLO ;HI OF Y = YMAXHI
JRST GETW ;IF LO OF Y .LE. YMAXLO, GO TO GETW
MSTK: HRLOI T0,377777 ;SET RESULT TO
SETO T1, ;MACHINE INFINITY
JUMPG T5,DSNH ;IF X IS NEGATIVE
DMOVN T0,T0 ;RESULT = -RESULT
DSNH: $LCALL ROV
;LERR (LIB,%,<DSINH: result overflow>)
POP P,T5 ;RESTORE ACCUMULATOR
GOODBY (1) ;RETURN
GETW: GFAD T0,LNV ;W = Y-LNV
DMOVEM T0,TEMP ;MOVE W TO A TEMP REGISTER
FUNCT GEXP.,<TEMP> ;Z = EXP(W)
DMOVEM T0,TEMP ;SAVE A COPY OF Z
GFMP T0,CON1 ;RESULT = Z*CON1
GFAD T0,TEMP ;+Z
JFCL MSTK ;OVERFLOW POSSIBLE
RET2: JUMPGE T5,RET3 ;IF SGNFLG IS NEGATIVE
DMOVN T0,T0 ;RESULT = -RESULT
RET3: POP P,T5 ;RESTORE ACCUMULATOR
GOODBY (1) ;RETURN
EASY: DMOVEM T0,TEMP ;MOVE Y TO TEMP
FUNCT GEXP.,<TEMP> ;GET ITS EXP
EXTEND T0,[GFSC -1] ;DIV BY 2
JRST RET2 ;GET RIGHT SIGN FOR RESULT.
ALG2: PUSH P,T2 ;SAVE MORE ACCUMULATORS
PUSH P,T3
PUSH P,T4
DMOVEM T0,TEMP ;MOVE Y TO A TEMPORARY REGISTER
FUNCT GEXP.,<TEMP> ;Z = DEXP(Y)
DMOVN T2,T0 ;SAVE A COPY OF Z
MOVEM T5,SGNFLG ;MOVE FLAG TO MEMORY
HRLZI T4,200140 ;SET T4 TO 1.0
SETZ T5, ;CLEAR SECOND WORD
GFDV T4,T2 ;1/Z
GFAD T0,T4 ;+Z
EXTEND T0,[GFSC -1] ;*1/2
POP P,T4 ;RESTORE ACCUMULATORS
POP P,T3
POP P,T2
RET1: SKIPGE SGNFLG ;IF SGNFLG IS NEGATIVE
DMOVN T0,T0 ;RESULT = -RESULT
RET: POP P,T5 ;RESTORE ACCUMULATOR
GOODBY (1) ;RETURN
RP0: DOUBLE 202352744232,262463203065 ;.35181283430177117881D+6
RP1: DOUBLE 201655127025,264501221757 ;.11563521196851768270D+5
RP2: DOUBLE 201050741013,034133711232 ;.16375798202630751372D+3
RP3: DOUBLE 200062423475,303374403264 ;.78966127417357099479D+0
Q0: DOUBLE 575137624613,372031435521 ;-.21108770058106271242D+7
Q1: DOUBLE 202043241271,035545730675 ;.36162723109421836460D+5
Q2: DOUBLE 576635220743,361550001577 ;-.27773523119650701667D+3
LNV: DOUBLE 577723506750,010134300000 ;-.693147180559947174D0
ONE: 200140000000 ;1.0
BIGX: 201254242673 ;709.089565
YMAXHI: 201254271027 ;709.782712893383998706
YMAXLO: 367643475715 ;
TWOM30: 174340000000 ;2**-30
CON1: 172041452000 ;.186417721E-14
TWENT2: 200554000000 ;
SEGMENT DATA
SGNFLG: 0
TEMP: DOUBLE 0,0
Z: DOUBLE 0,0
PRGEND
TITLE GSQRT SQUARE ROOT FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MARCH 19, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GSQRT
EXTERN GSQRT.
GSQRT=GSQRT.
PRGEND
TITLE GSQRT. SQUARE ROOT FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MARCH 19, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;GSQRT(X) IS CALCULATED AS FOLLOWS
;THE ROUTINE CALCULATES THE SQUARE ROOT OF THE ABSOLUTE VALUE OF A
;DOUBLE PRECISION ARGUMENT BY DOING A LINEAR SINGLE PRECISION
;APPROXIMATION ON THE HIGH ORDER WORD, FOLLOWED BY TWO SINGLE PRECISION
;NEWTON ITERATIONS AND TWO DOUBLE PRECISION NEWTON ITERATIONS. AN ERROR
;MESSAGE RESULTS FOR NEGATIVE ARGUMENTS.
;REQUIRED (CALLED) ROUTINES: NONE
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; MOVEI L,ARG
; PUSHJ P,GSQRT
;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1
SEARCH MTHPRM
SEGMENT CODE
HELLO (GSQRT,.) ;ENTRY TO GSQRT ROUTINE
DMOVE T0,@(L) ;GET DP ARGUMENT
JUMPG T0,GSQRTP ;ARGUMENT IS GREATER THAN 0
JUMPE T0,GSQRT4 ;ARGUMENT IS ZERO
$LCALL NAA
;LERR (LIB,%,<GSQRT: negative arg; result = GSQRT(DABS(arg))>)
DMOVN T0,T0 ;ARG = -ARG
GSQRTP: PUSH P,T2 ;SAVE ACCUMULATORS
PUSH P,T3
PUSH P,T4
PUSH P,T5
DMOVE T4,T0 ;COPY ARG
LSH T4,-1 ;COMPUTE LINEAR APPROXIMATION
TLZE T4,40
JRST GSQRT2 ;YES, GO TO GSQRT2
ADD T4,[XWD 26,760700]
GFMP T4,[EXP 300145400000,0]
JRST GSQRT3
GSQRT2: ADD T4,[XWD 26,760700]
GFMP T4,[EXP 300165000000,0]
GSQRT3: DMOVE T2,T0 ;COPY ORIGINAL ARG
GFDV T2,T4 ;DO NEWTON ITERATIONS
GFAD T2,T4
EXTEND T2,[GFSC -1]
DMOVE T4,T0
GFDV T4,T2
GFAD T4,T2
EXTEND T4,[GFSC -1]
DMOVE T2,T0
GFDV T2,T4
GFAD T2,T4
EXTEND T2,[GFSC -1]
GFDV T0,T2
GFAD T0,T2
EXTEND T0,[GFSC -1]
POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4
POP P,T3
POP P,T2
GSQRT4: GOODBY (1) ;RETURN
PRGEND
TITLE GCOTAN TANGENT AND COTANGENT FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. APRIL 16, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GCOTAN
EXTERN GCOTN.
GCOTAN=GCOTN.
PRGEND
TITLE GTAN TANGENT AND COTANGENT FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC APRIL 16, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GTAN
EXTERN GTAN.
GTAN=GTAN.
PRGEND
TITLE GTAN. TANGENT AND COTANGENT FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC APRIL 16, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;TAN(X) AND COTAN(X) ARE CALCULATED AS FOLLOWS;
; THE IDENTITIES
; TAN(PI/2.0-G)=1.0/TAN(G)
; TAN(N*PI+H)=TAN(H), -PI/2.0 < H <= PI/2.0
; TAN(-X)=-TAN(X)
; COTAN(X)=1.0/TAN(X)
; COTAN(-X)=-COTAN(X)
; ARE USED TO REDUCE THE ORIGINAL PROBLEM TO A PROBLEM WITH
; AN ARGUMENT GREATER THAN -PI/2.0 AND LESS THAN OR EQUAL TO
; PI/2.0, PI=3.14159265358979324.
; THEN DEFINE N AND F SO THAT
; X=N*PI/4.0+F, 0.0 <= F <= PI/4.0, WITH PI=3.14159265358979324.
; THEN
; TAN(F) = F, IF F<2**(-30)
; = R(F), OTHERWISE
; WHERE
; R(F)=(((XP3*G+XP2)*G+XP1)*G+XP0)/((((Q4*G+Q3)*G+Q2)*G+Q1)*G+Q0)
; AND
; G = F*F
; XPi AND Qi ARE GIVEN BELOW
; THE APPROXIMATION IS DERIVED FROM ONE GIVEN IN CODY AND WAITE,
; "SOFTWARE MANUAL FOR THE ELEMENTARY FUNCTIONS"
; THE RESULT IS THEN RECIPROCATED, IF NECESSARY, AND GIVEN
; THE APPROPRIATE SIGN.
;THE RANGE OF DEFINITION FOR GTAN IS 0 < ABS(X) < ((2**29) * (PI/2)) = 843314856.D0
; AND FOR GCOTAN(X), (2**(-1023))*(1+(2**(-58))) < ABS(X) < ((2**29)*(PI/2)) IS NECESSARY.
; OTHERWISE, ERROR MESSAGES WILL RESULT.
;REQUIRED (CALLED) ROUTINES: NONE
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; MOVEI L,ARG
; PUSHJ P,GTAN
; OR
; PUSHJ P,GCOTAN
;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0,
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1.
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (GCOTN.,GCOTAN) ;ENTRY TO GCOTAN ROUTINE
PUSH P,T4 ;SAVE ACCUMULATORS
PUSH P,T5
SETZM FLAG ;CLEAR FLAG FOR GCOTAN
DMOVE T0,@(L) ;GET ARG
DMOVE T4,T0 ;Y = X
JUMPGE T0,XPOS ;IF X IS NEGATIVE
DMOVN T0,T0 ;SET Y = -X
XPOS: CAMGE T0,EPS1 ;IF THE HI ORDER PART OF Y IS TOO SMALL
JRST ERR1 ;THEN GO TO ERR1
CAME T0,EPS1 ;IF HI OF Y > HI OF EPS
JRST YOK ;THEN GO TO YOK
JUMPN T1,YOK ;HI OF Y = HI OF EPS, IS LO OF Y OK?
ERR1: $LCALL ROV
;LERR (LIB,%,DCOTAN: result overflow)
HRLOI T0,377777 ;SET RESULT TO MACHINE
HRLOI T1,377777 ;INFINITY
JUMPGE T4,RET ;IF ARG IS NEGATIVE
DMOVN T0,T0 ;RESULT = - RESULT
POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4
GOODBY (1) ;RETURN
HELLO (GTAN,.) ;ENTRY TO TAN ROUTINE
PUSH P,T4 ;SAVE ACCUMULATORS
PUSH P,T5
SETOM FLAG ;SET FLAG FOR GTAN ENTRY
DMOVE T0,@(L) ;GET COPY OF ARG
DMOVE T4,T0 ;Y = X
JUMPGE T0,YOK ;IF X IS NEGATIVE
DMOVN T0,T0 ;Y = -X
YOK: CAMGE T0,YMAX1 ;IF HI OF Y < HI OF YMAX
JRST ALG ;PROCEED TO ALG
CAME T0,YMAX1 ;IF HI OF Y > HI OF YMAX
JRST ERR2 ;GO TO ERR2
CAMG T1,YMAX2 ;HI OF Y=HI OF YMAX, IS LO OF Y .LE. LO OF YMAX?
JRST ALG ;YES, PROCEED TO ALG
ERR2: $LCALL ATZ
;LERR (LIB,%,<DTAN or DCOTAN: ABS(arg) too large; result = zero>)
MOVEI T0,0 ;SET RESULT TO ZERO
MOVEI T1,0 ;SET SECOND WORD TO ZERO
POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4
GOODBY (1) ;RETURN
ALG: PUSH P,T2 ;SAVE ACCUMULATORS
PUSH P,T3
CAML T0,PIO4HI ;COMPARE |ARG| WITH PI/4
JRST REDUCE ; REDUCTION IS NECESSARY
DMOVE T2,T0 ;|F| TO T2
DMOVE T0,T4 ;AND SIGNED F TO T0.
SETZM N ;STORE 0 FOR N
JRST FPOS ;BYPASS REDUCTION
REDUCE: GFMP T0,TWODPI ;Y*(2/PI)
EXTEND T0,[GFIXR T0] ;FIX IT
MOVEM T0,N ;MOVE N TO MEMORY
EXTEND T0,[GFLTR T0] ;FLOAT IT
JUMPLE T4,NEXT ;IF X IS POSITIVE
DMOVN T0,T0 ;NEGATE XN
NEXT: DMOVE T2,T0 ;GET A COPY OF -XN
GFMP T2,C1 ;FORM -XN*C1
GFAD T4,T2 ;F=X - XN*C1
DMOVE T2,T0 ;GET A COPY OF -XN
GFMP T2,C3 ;FORM -XN*C3
GFMP T0,C2 ;FORM -XN*C2
GFAD T4,T0 ;F=F - XN*C2
GFAD T2,T4 ;F=F - XN*C3
DMOVE T0,T2 ;MOVE COPY OF F INTO T0
JUMPGE T2,FPOS ;IF F IS NEGATIVE
DMOVN T2,T2 ;F = -F
FPOS: CAML T2,HIEPS ;IS F < EPS?
JRST NO ;NO, GO TO NO
HRLZI T4,200140 ;YES , F< EPS,
MOVEI T5,0 ;XDEN = 1.0
JRST FLGCHK ;GO TO CHECK FLAG
NO: GFMP T2,T2 ;NO, F .GE. EPS, SET G=F*F
DMOVE T4,T2 ;SAVE A COPY OF G
GFMP T2,XP3 ;XNUM = XP3*G +
GFAD T2,XP2 ; P2 *
GFMP T2,T4 ; G +
GFAD T2,XP1 ; P1 *
GFMP T2,T4 ; G *
GFMP T2,T0 ; F +
GFAD T0,T2 ; F
DMOVE T2,T4 ;SAVE A COPY OF G
GFMP T4,Q4 ;XDEN = Q4 * G +
GFAD T4,Q3 ; Q3 *
GFMP T4,T2 ; G +
GFAD T4,Q2 ; Q2 *
GFMP T4,T2 ; G +
GFAD T4,Q1 ; Q1 *
GFMP T4,T2 ; G +
GFAD T4,ONE ; 1.0
FLGCHK: MOVE T2,N ;GET COPY OF N
SKIPE FLAG ;IF FLAG IS NOT ZERO
JRST DA ;THEN GO TO DA
TRNN T2,1 ;OTHERWISE, IS N ODD?
JRST DOVN ;NO, GO GET RESULT
DMOVN T0,T0 ;XNUM = -XNUM
JRST NOVD ;GO TO NOVD
DA: TRNN T2,1 ;FLAG MUST BE ONE, IS N ODD?
JRST NOVD ;NO, GO TO NOVD
DMOVN T0,T0 ;XNUM = -XNUM
DOVN: GFDV T4,T0 ;RESULT = XDEN/XNUM
DMOVE T0,T4
JRST RET ;GO TO RET
NOVD: GFDV T0,T4 ;RESULT = XNUM/XDEN
RET: POP P,T3 ;RESTORE ACCUMULATORS
POP P,T2
POP P,T5
POP P,T4
GOODBY (1)
XP1: DOUBLE 600135665120,325457340031 ;-.133383500064219607D0
XP2: DOUBLE 177070072025,061672533663 ;.342488782358905900D-2
XP3: DOUBLE 601632425106,212723433423 ;-.178617073422544267D-4
Q1: DOUBLE 600004205175,300102305265 ;-.466716833397552942D0
Q2: DOUBLE 177364436365,013742053124 ;.256638322894401129D-1
Q3: DOUBLE 601227102333,067553367667 ;-.311815319070100273D-3
Q4: DOUBLE 175441335647,203547435105 ;.498194339937865123D-6
C1: DOUBLE 200162207732,240000000000 ;HIGH 30 BITS OF PI/2
C2: DOUBLE 174342055060,200000000000 ;C1+C2+C3=PI/2 TO EXTRA PREC
C3: DOUBLE 170764611431,212134015604
TWODPI: DOUBLE 200050574603,156234420251
ONE: DOUBLE 200140000000,000000000000
PIO4HI: 200062207733 ;HIGH ORDER PART OF PI/4
EPS1: 000240000000 ;HIGH ORDER PART OF EPS1
YMAX1: 203662207732 ;HIGH ORDER PART OF YMAX
YMAX2: 242102643021 ;LOW ORDER PART OF YMAX
HIEPS: 174340000000 ; 2**(-31)
SEGMENT DATA
N: 0
FLAG: 0
PRGEND
TITLE GTANH HYPERBOLIC TANGENT FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GTANH
EXTERN GTANH.
GTANH=GTANH.
PRGEND
TITLE GTANH. HYPERBOLIC TANGENT FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;TANH(X) IS CALCULATED AS FOLLOWS
; TANH(X) = -TANH(-X). LET F = ABS(X)
; IF F > 21.1409890, TANH = 1.0*SIGN(X)
; IF F > LN(3)/2 AND F <= 21.1409890, TANH = RESULT 1 =
; SIGN(X)*(1 - 2/(EXP(2*F) + 1)))
; IF F < 2**(-29), TANH = F*SIGN(X)
; OTHERWISE, LET G = F**2, AND OBTAIN RESULT 2 AS
; TANH = SIGN(X)*(F + F*R(G)), WHERE
; R(G) = G*(RP0+G*(RP1+RP2*G))/(Q0+G*(Q1+G*(Q2+G))).
; RP0, RP1, RP2, Q0, Q1, AND Q2 APPEAR BELOW.
;THE RANGE OF DEFINITION FOR TANH IS THE REPRESENTABLE REAL NUMBERS
;REQUIRED (CALLED) ROUTINES: DEXP
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; MOVEI L,ARG
; PUSHJ P,TANH
;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (GTANH,.) ;ENTRY TO GTANH ROUTINE
DMOVE T0,@(L) ;OBTAIN X
PUSH P,T5 ;SAVE AN ACCUMULATOR
MOVE T5,T0 ;SAVE A COPY OF X
JUMPGE T0,XPOS ;IF X IS NEGATIVE
DMOVN T0,T0 ;F = -X
XPOS: CAMG T0,LN2TC ;IF F IS .LE. LN2TC
JRST CALCF ;GO TO CALCF
HRLZI T0,200140 ;SET RESULT TO 1.0
SETZ T1, ;ZERO SECOND WORD
JUMPGE T5,RET1 ;IF X IS NEGATIVE
DMOVN T0,T0 ;RESULT = -RESULT
RET1: POP P,T5 ;RESTORE ACCUMULATORS
GOODBY (1) ;RETURN
CALCF: PUSH P,T2 ;SAVE MORE ACCUMULATORS
PUSH P,T3
PUSH P,T4
CAMG T0,LN3D2 ;IF F IS .LE. LN3D2
JRST ALG1 ;GO TO ALG1
EXTEND T0,[GFSC 1] ;F = F+F
DMOVEM T0,TEMP ;SAVE F IN A TEMP REGISTER
FUNCT GEXP.,<TEMP> ;EXP(2*F)
GFAD T0,ONE ;+1.0
HRLZI T2,577540 ;SET T2 TO -2.0
SETZ T3, ;ZERO SECOND WORD
GFDV T2,T0 ;-2/(EXP(2*F)+1)
DMOVE T0,T2
GFAD T0,ONE ;1. - 2/(EXP(2*F)+1)
JUMPGE T5,RET ;IF X IS NEGATIVE
DMOVN T0,T0 ;RESULT = -RESULT
POP P,T4 ;RESTORE ACCUMULATORS
POP P,T3
POP P,T2
POP P,T5
GOODBY (1) ;RETURN
ALG1: CAML T0,CON1 ;IF F IS .GE. 2**(-29)
JRST ALG2 ;GO TO ALG2
JUMPGE T5,RET ;IF X IS NEGATIVE
DMOVN T0,T0 ;RESULT = -RESULT
POP P,T4 ;RESTORE ACCUMULATORS
POP P,T3
POP P,T2
POP P,T5
GOODBY (1) ;RETURN
ALG2: JUMPGE T5,POSARG ;IF ARGUMENT < 0
DMOVN T0,T0 ; NEGATE |ARGUMENT|
POSARG: DMOVEM T0,TEMP ;SAVE SIGNED ARGUMENT.
GFMP T0,T0 ;G = F*F
DMOVE T2,T0 ;SAVE A COPY OF G
GFAD T2,Q2 ;XDEN = G+Q2
GFMP T2,T0 ;*G
GFAD T2,Q1 ;+Q1
GFMP T2,T0 ;*G
GFAD T2,Q0 ; + Q0
DMOVE T4,T0 ;SAVE A COPY OF G
GFMP T4,RP2 ;XNUM = G*RP2
GFAD T4,RP1 ; + RP1
GFMP T4,T0 ; * G
GFAD T4,RP0 ; + RP0
GFMP T0,T4 ; *G
GFDV T0,T2 ;R(G) = XNUM/XDEN
GFMP T0,TEMP ;RESULT = F*R
GFAD T0,TEMP ; + F
RET: POP P,T4 ;RESTORE ACCUMULATORS
POP P,T3
POP P,T2
POP P,T5
GOODBY (1) ;RETURN
RP0: DOUBLE 576415451321,262036074743 ;-.161341190239962281D+4
RP1: DOUBLE 577016306122,362132146345 ;-.992259296722360833D+2
RP2: DOUBLE 577702217271,210105420140 ;-.964374927772254698D0
Q0: DOUBLE 201545640742,272351322265 ;.484023570719886887D+4
Q1: DOUBLE 201442716132,150037036260 ;.22337720718962312926D+
Q2: DOUBLE 200770276517,017277377240 ;.112744743805349493D+3
LN2TC: 200552220277 ;21.1409890E0
CON1: 174340000000 ;2**(-30)
LN3D2: 200043117523 ;.549306144E0
ONE: DOUBLE 200140000000,000000000000 ;1.0
SEGMENT DATA
TEMP: DOUBLE 0,0
PRGEND
TITLE GINT DOUBLE PRECISION TRUNCATION TO INTEGER
SUBTTL CHRIS SMITH/CKS 29-Jan-80
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GINT
EXTERN GINT.
GINT=GINT.
PRGEND
TITLE GABS DOUBLE PRECISION ABSOLUTE VALUE
SUBTTL CHRIS SMITH/CKS 29-Jan-80
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GABS
EXTERN DABS.
GABS=DABS.
PRGEND
TITLE GABS. DOUBLE PRECISION ABSOLUTE VALUE
SUBTTL CHRIS SMITH/CKS 29-Jan-80
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GABS.
EXTERN DABS.
GABS.=DABS.
PRGEND
TITLE GMAX1 DOUBLE PRECISION MAXIMUM OF A SERIES OF ARGUMENTS
SUBTTL CHRIS SMITH/CKS 29-Jan-80
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GMAX1
EXTERN DMAX1.
GMAX1=DMAX1.
PRGEND
TITLE GMAX1. DOUBLE PRECISION MAXIMUM OF A SERIES OF ARGUMENTS
SUBTTL CHRIS SMITH/CKS 29-Jan-80
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GMAX1.
EXTERN DMAX1.
GMAX1.=DMAX1.
PRGEND
TITLE GMIN1 DOUBLE PRECISION MINIMUM OF A SERIES OF ARGUMENTS
SUBTTL CHRIS SMITH/CKS 29-Jan-80
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GMIN1
EXTERN DMIN1.
GMIN1=DMIN1.
PRGEND
TITLE GMIN1. DOUBLE PRECISION MINIMUM OF A SERIES OF ARGUMENTS
SUBTTL CHRIS SMITH/CKS 29-Jan-80
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GMIN1.
EXTERN DMIN1.
GMIN1.=DMIN1.
PRGEND
TITLE GSIGN DOUBLE PRECISION TRANSFER OF SIGN FUNCTION
SUBTTL CHRIS SMITH/CKS 29-Jan-80
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GSIGN
EXTERN DSIGN.
GSIGN=DSIGN.
PRGEND
TITLE GSIGN. DOUBLE PRECISION TRANSFER OF SIGN FUNCTION
SUBTTL CHRIS SMITH/CKS 29-Jan-80
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY GSIGN.
EXTERN DSIGN.
GSIGN.=DSIGN.
PRGEND
TITLE DTOG DOUBLE PRECISION TO GFLOAT DOUBLE PRECISION CONVERSION
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;ARG BLOCK OFFSETS
SRC==0 ;SOURCE ARRAY ADDRESS
DST==1 ;DESTINATION ARRAY ADDRESS
CNT==2 ;NUMBER OF ITEMS TO CONVERT
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (DTOG)
;CONVERT ONE NUMBER
DMOVE T0,@(L) ;GET THE D NUMBER
PUSHJ P,DTOGS ;CONVERT IT
GOODBY
HELLO (DTOGA)
;CONVERT AN ARRAY OF NUMBERS
PUSH P,T3 ;SAVE SOME AC'S
PUSH P,T4
PUSH P,T5
XMOVEI T3,@SRC(L) ;GET D PNTR
XMOVEI T4,@DST(L) ;GET G PNTR
MOVE T5,@CNT(L) ;GET COUNT
DTOGLP: DMOVE T0,(T3) ;GET A D NUMBER
PUSHJ P,DTOGS ;CONVERT IT
DMOVEM T0,(T4) ;SAVE IT
ADDI T3,2 ;INCR PNTRS
ADDI T4,2
SOJG T5,DTOGLP ;LOOP FOR ARRAY
POP P,T5
POP P,T4
POP P,T3
GOODBY
DTOGS: JUMPL T0,NEGDTG ;DO SEPARATELY IF NEGATIVE
JUMPE T0,DTGZER ;AND NOTHING IF ZERO
PUSHJ P,EXPDTG ;NOW DO THE MOST STUFF
POPJ P,
NEGDTG: DMOVN T0,T0 ;MAKE IT POSITIVE
PUSHJ P,EXPDTG ;DO MOST STUFF
DMOVN T0,T0 ;MAKE IT NEGATIVE AGAIN
POPJ P,
DTGZER: DMOVE T0,[EXP 0,0] ;LOAD ALL ZEROES
POPJ P,
EXPDTG: PUSH P,T2 ;SAVE AN AC
CAMN T0,[377777,,-1] ;'OVERFLOW' AMOUNT?
CAME T1,[377777,,-1]
JRST DTGOK ;NO
JRST DTGDON ;YES. LEAVE AS IS
DTGOK: LDB T2,[POINT 8,T0,8] ;GET THE EXPONENT
ADDI T2,1600 ;CONVERT TO G-TYPE EXP
TLZ T0,777000 ;CLEAR THE EXPONENT
DADD T0,[EXP 0,4] ;ROUND UP
ASHC T0,-3 ;SHIFT FRACTION
TLNN T0,100 ;DID WE OVERFLOW?
JRST EXPOK ;NO
ASHC T0,-1 ;YES. SHIFT 1 MORE
ADDI T2,1 ;AND INCR THE EXPONENT
EXPOK: DPB T2,[POINT 12,T0,11] ;DROP IN THE EXP
DTGDON: POP P,T2 ;RESTORE THE AC
POPJ P,
PRGEND
TITLE GTOD GFLOAT DOUBLE PRECISION TO DOUBLE PRECISION CONVERSION
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;ARG BLOCK OFFSETS
SRC==0 ;SOURCE ARRAY ADDRESS
DST==1 ;DESTINATION ARRAY ADDRESS
CNT==2 ;NUMBER OF ITEMS TO CONVERT
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (GTOD)
;CONVERT ONE NUMBER
DMOVE T0,@(L) ;GET THE G NUMBER
PUSHJ P,GTODS ;CONVERT IT
GOODBY
HELLO (GTODA)
;CONVERT AN ARRAY OF NUMBERS
PUSH P,T3 ;SAVE SOME AC'S
PUSH P,T4
PUSH P,T5
XMOVEI T3,@SRC(L) ;GET D PNTR
XMOVEI T4,@DST(L) ;GET G PNTR
MOVE T5,@CNT(L) ;GET COUNT
GTODLP: DMOVE T0,(T3) ;GET A G NUMBER
PUSHJ P,GTODS ;CONVERT IT
DMOVEM T0,(T4) ;SAVE IT
ADDI T3,2 ;INCR PNTRS
ADDI T4,2
SOJG T5,GTODLP ;LOOP FOR ARRAY
POP P,T5
POP P,T4
POP P,T3
GOODBY
GTODS: JUMPL T0,NEGGTD ;DO SEPARATELY IF NEGATIVE
JUMPE T0,GTDZER ;AND NOTHING IF ZERO
PUSHJ P,EXPGTD ;NOW DO THE MOST STUFF
POPJ P,
NEGGTD: DMOVN T0,T0 ;MAKE IT POSITIVE
PUSHJ P,EXPGTD ;DO MOST STUFF
DMOVN T0,T0 ;MAKE IT NEGATIVE AGAIN
POPJ P,
GTDZER: DMOVE T0,[EXP 0,0] ;LOAD ALL ZEROES
POPJ P,
EXPGTD: PUSH P,T2 ;SAVE AN AC
CAMN T0,[377777,,-1] ;'OVERFLOW' AMOUNT?
CAME T1,[377777,,-1]
JRST GTDOK ;NO
JRST GTDDON ;YES. LEAVE AS IS
GTDOK: LDB T2,[POINT 11,T0,11] ;GET THE EXPONENT
SUBI T2,1600 ;CONVERT TO D-TYPE EXP
JUMPL T2,GTDLOW ;EXPONENT TO SMALL
CAIL T2,400 ;TEST IF TOO BIG
JRST GTDHGH ;YUP. IT'S TO BIG
TLZ T0,777700 ;CLEAR THE EXPONENT
ASHC T0,3 ;SHIFT THE FRACTION
DPB T2,[POINT 9,T0,8] ;DROP IN THE EXP
JRST GTDDON
GTDLOW: $LCALL RUN
;LERR (LIB,%,<GTOD: result underflow>)
DMOVE T0,[EXP 0,0] ;LOAD ZEROES
JRST GTDDON ;AND LEAVE
GTDHGH: $LCALL ROV
;LERR (LIB,%,<GTOD: result overflow>)
DMOVE T0,[EXP 377777777777,377777777777] ;LOAD AN OVERFLOW
GTDDON: POP P,T2 ;RESTORE THE AC
POPJ P,
PRGEND
TITLE GSN.0
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GSN.0
GSN.0: GSNGL 0
PRGEND
TITLE GSN.2
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GSN.2
GSN.2: GSNGL 2
PRGEND
TITLE GSN.4
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GSN.4
GSN.4: GSNGL 4
PRGEND
TITLE GSN.6
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GSN.6
GSN.6: GSNGL 6
PRGEND
TITLE GSN.10
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GSN.10
GSN.10: GSNGL 10
PRGEND
TITLE GSN.12
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GSN.12
GSN.12: GSNGL 12
PRGEND
TITLE GSN.14
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GSN.14
GSN.14: GSNGL 14
PRGEND
TITLE GDB.0
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GDB.0
GDB.0: GDBLE 0
PRGEND
TITLE GDB.2
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GDB.2
GDB.2: GDBLE 2
PRGEND
TITLE GDB.4
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GDB.4
GDB.4: GDBLE 4
PRGEND
TITLE GDB.6
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GDB.6
GDB.6: GDBLE 6
PRGEND
TITLE GDB.10
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GDB.10
GDB.10: GDBLE 10
PRGEND
TITLE GDB.12
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GDB.12
GDB.12: GDBLE 12
PRGEND
TITLE GDB.14
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GDB.14
GDB.14: GDBLE 14
PRGEND
TITLE GFX.0
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GFX.0
GFX.0: GFIX 0
PRGEND
TITLE GFX.2
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GFX.2
GFX.2: GFIX 2
PRGEND
TITLE GFX.4
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GFX.4
GFX.4: GFIX 4
PRGEND
TITLE GFX.6
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GFX.6
GFX.6: GFIX 6
PRGEND
TITLE GFX.10
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GFX.10
GFX.10: GFIX 10
PRGEND
TITLE GFX.12
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GFX.12
GFX.12: GFIX 12
PRGEND
TITLE GFX.14
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GFX.14
GFX.14: GFIX 14
PRGEND
TITLE GFL.0
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GFL.0
GFL.0: GFLTR 0
PRGEND
TITLE GFL.2
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GFL.2
GFL.2: GFLTR 2
PRGEND
TITLE GFL.4
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GFL.4
GFL.4: GFLTR 4
PRGEND
TITLE GFL.6
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GFL.6
GFL.6: GFLTR 6
PRGEND
TITLE GFL.10
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GFL.10
GFL.10: GFLTR 10
PRGEND
TITLE GFL.12
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GFL.12
GFL.12: GFLTR 12
PRGEND
TITLE GFL.14
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GFL.14
GFL.14: GFLTR 14
PRGEND
TITLE GDIM G-floating positive difference
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
; OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.
;
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY GDIM
EXTERN GDIM.
GDIM=GDIM.
PRGEND
TITLE GDIM. G-Floating Positive Difference.
SUBTTL C. McCutcheon 6/26/81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
; Fortran Intrinsic Function for Positive Difference.
; Passed and returned arguments are G-floating double precision.
; Passed: A1,A2
; Returns: (from page 15-23 of standard)
; A1-A2 if A1 .GT. A2
; 0 if A1 .LE. A2
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (GDIM,.)
PUSH P,T2
PUSH P,T3
DMOVE T0,@0(L) ;A1
DMOVE T2,@1(L) ;A2
CAMN T0,T2 ;If A1 .LT. A2,
CAML T1,T3
CAMGE T0,T2
JRST RET0 ; return 0.
GFSB T0,@1(L) ;Subtract A1-A2
JFCL EXCEP ;Can underflow and overflow
RET: POP P,T3
POP P,T2
GOODBYE ;Return
RET0: SETZB T0,T1 ;Return zero
JRST RET
EXCEP: JUMPE T0,UNDER ;Underflow?
$LCALL ROV,RET ;No, result overflow
UNDER: $LCALL RUN,RET ;Result underflow
PRGEND
TITLE GINT. G-Floating truncation.
SUBTTL C. McCutcheon - 6/25/81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
; Fortran Instrinsic Function to return a G-floating double
; precision truncation of the double precision number passed.
; If Magnitude(A) .LT. 1, then 0 is returned,
; otherwise return the largest integer that does not exceed the
; magnitude of A, and whose sign is the same as that of A.
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (GINT,.) ;Enter routine
DMOVE T0,@0(L) ;Get argument to truncate.
JUMPN T0,NZERO ;If zero passed, return.
GOODBYE ;Return
NZERO: PUSH P,T2 ;Save ac.
CAIGE T0,0 ;If original number .LT. zero,
DMOVN T0,T0 ; negate it.
HLRZ T2,T0 ;Get exponent
LSH T2,-6 ;Put rightmost
CAIG T2,^O2000 ;[3215] If exponent .LE. 2000 then
JRST ZERO ; return zero (number .LT. 1)
CAIL T2,^O2073 ;If exponent .GE. 2073 then
JRST DONE ; return the number passed.
; Now shift out the fractional part of the number.
ASHC T0,-2073(T2) ;Shift into integer position.
MOVN T2,T2 ;Negate exponent.
ASHC T0,2073(T2) ;Shift back to where found.
SKIPGE @0(L) ;If original number .LT. zero,
DMOVN T0,T0 ; negate it again.
BYE: POP P,T2 ;Pop saved ac's
GOODBYE ;Return
DONE: DMOVE T0,@0(L) ;No fractional part, return
JRST BYE ; the number passed.
ZERO: SETZB T0,T1 ;Return 0
JRST BYE
PRGEND
TITLE GNINT. G-Floating nearest whole number.
SUBTTL C. McCutcheon - 6/25/81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
; Fortran Instrinsic Function to return a G-floating double
; precision nearest whole number of the double precision number
; passed.
; Number returned is: Integer(A+.5) if A .GE. 0
; Integer(A-.5) if A .LT. 0.
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (GNINT,.) ;Enter routine
DMOVE T0,@0(L) ;Get argument to truncate.
JUMPN T0,NZERO ;If zero passed,
GOODBYE ; return.
NZERO: PUSH P,T2 ;Save ac's
CAIGE T0,0 ;If original number .LT. zero,
DMOVN T0,T0 ; negate it.
; Now Shift out the fractional part of the number.
HLRZ T2,T0 ;Get exponent
LSH T2,-6 ;Put rightmost
CAIL T2,^O2073 ;If exponent .GE. 2073 then
JRST DONE ; return the number passed.
CAIGE T2,2000 ;number much .LT. 1.
JRST ZERO ; return 0.
GFAD T0,[200040,,0
0,,0] ;Add .5 before truncation.
HLRZ T2,T0 ;Get exponent
LSH T2,-6 ;Put rightmost
ASHC T0,-2073(T2) ;Shift into integer position.
MOVN T2,T2 ;Negate exponent.
ASHC T0,2073(T2) ;Shift back to where found.
SKIPGE @0(L) ;If original number .LT. zero,
DMOVN T0,T0 ; negate it again.
BYE: POP P,T2 ;Pop saved ac's
GOODBYE ;Return
DONE: DMOVE T0,@0(L) ;No fractional part, return
JRST BYE ; the number passed.
ZERO: SETZB T0,T1 ;Set to 0.
JRST BYE
PRGEND
TITLE GPROD. G-Floating product for single prec. factors.
SUBTTL C. McCutcheon - 6/25/81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
; Fortran Instrinsic Function to return a double precision
; product for two real arguments passed it.
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (GPROD,.) ;Enter routine
EXTEND T0,[GDBLE @1(L)] ;Convert 2nd arg to G-floating.
DMOVEM T0,MULT ;Save 2nd argument
EXTEND T0,[GDBLE @0(L)] ;Convert 1st arg to G-floating.
GFMP T0,MULT ;Multiply and leave result in AC0.
JFCL EXCEP ;Result can under/overflow
RET: POPJ P, ;Return
EXCEP: JUMPE T0,UNDER ;Underflow?
$LCALL ROV,RET
UNDER: $LCALL RUN,RET
SEGMENT DATA
MULT: BLOCK 2 ;Multiplier
PRGEND
TITLE IGNIN. Integer nearest whole number for G-Gloating.
SUBTTL C. McCutcheon - 6/25/81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
; Fortran Instrinsic Function to return the nearest
; integer to the G-floating double precision number passed.
; Number returned is: Integer(A+.5) if A .GE. 0
; Integer(A-.5) if A .LT. 0.
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (IGNIN.,IGNINT) ;Enter routine
DMOVE T0,@0(L) ;Get argument to round.
JUMPN T0,NZERO ;If number passed = 0,
GOODBYE ;Return
NZERO: CAIGE T0,0 ;If original number .LT. zero,
DMOVN T0,T0 ; negate it.
TLNN T0,200000 ;Is |number| .LT. 1/2 ?
JRST ZERO ;Yes, go return zero
PUSH P,T2 ;Save an ac
GFAD T0,[200040,,0
0,,0] ;Add 0.5G0 before truncation.
; Now shift out the fractional part of the number.
HLRZ T2,T0 ;Get exponent
LSH T2,-6 ;Put rightmost
CAILE T2,2043 ;If exponent .GE. 2043 then
JRST DONE ; return largest integer
TLZ T0,777700 ;Eliminate exponent.
ASHC T0,-2030(T2) ;Shift into integer position,
; integer result in T0.
SKIPGE @0(L) ;If original number .LT. zero,
MOVN T0,T0 ; negate it again.
BYE: POP P,T2 ;Pop saved ac's
GOODBYE ;Return
; Number is .LT. 1/2, return zero quickly (don't use the normal flow
; because GFAD 0.5G0 will round)
ZERO: SETZ T0, ;Return . . .
GOODBYE ; . . . 0
; Number too large to represent as integer. Return largest
; integer.
DONE: $LCALL ROV
;LERR (LIB,%,<IDNINT: Result overflow>)
HRLOI T0,377777 ;Largest positive integer.
SKIPGE @0(L) ;If original number .LT. zero,
HRLZI T0,400000 ; largest negative integer.
JRST BYE ;
END