Trailing-Edge
-
PDP-10 Archives
-
FORTRAN-10_Alpha_31-jul-86
-
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 ,2(4020)
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1986
;ALL RIGHTS RESERVED.
;
;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.
***** Begin Version 1A *****
3233 CKS 23-Mar-83
Rewrote GMOD to remove restrictions on the magnitudes of
its arguments.
3243 BCM 29-Mar-83
Get and clear PC flags with new GETFLG and RESFLG macros. Fixes
always clearing underflow in MTHTRP. Only applies to JOV branches.
***** End Revision History *****
***** Begin Version 2 *****
4002 JLC 3-Aug-83
Move GEXP. to after GTAN., which calls GEXP.
Make 1.0 a special case for GEXP3.
4011 JLC 4-May-84
Add code for MOD, AMOD, DMOD, and GMOD for 2nd arg=0.
4012 MRB 6-Jun-84
Change name of routines GTODA and DTOGA to GTODA. and DTOGA.
to allow flagger. These symbols will be defined in FORCOM.
4145 MRB 17-Aug-84
Change GABS., GMAX1., GMIN1. & GSIGN. to <nnn.+0> to get
macro to generate the correct polish stuff.
4020 DCE 17-Dec-85
The GFLOATING exponentiation routine uses an incorrect method of
getting the exponent of the magnitude of certain intermediate
GFLOATING numbers. A MOVM on the first word of the double word
value is used followed by stripping out the exponent field.
This does not work for all negative GFLOATING numbers where the
mantissa part in the first word is all zero (the first
non-zero bit of the mantissa is in the second word). Due to the
fact that these errors occur in an intermediate calculation,
it is impossible to specify the class of input numbers which
will result in this problem.
\
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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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 GEXP2. DOUBLE ** INTEGER EXPONENTIATION
SUBTTL MARY PAYNE/MHP/CKS
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1983, 1986
;ALL RIGHTS RESERVED.
;
;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 T3,T0 ;[3243] Copy result
DMOVE T0,ONE ;Get reciprocal of
GFDV T0,T3 ;[3243] 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 over/underflows.
GETFLG T2 ;[3243] get exception flags into t2
TLNE T2,(PC%FUF) ;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 OVMSGF ;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: GETFLG T2 ;[3243] Get exception flags into t2
TLNE T2,(PC%FUF) ;If underflow flag set
JRST UNDER ; underflow on product
SKIPL @1(L) ;Else, overflow on result if
JRST OVMSGF ; exponent > 0. Get message
DMOVE T3,T0 ;[3243] Copy result
DMOVE T0,ONE ;For exponent < 0, get
GFDV T0,T3 ;[3243] reciprocal of wrapped overflow
JOV RETF ;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
;[3243] entry point necessary because other branches join OVMSG
OVMSGF: RESFLG T2 ;[3243] restore the flags, t2-t3 are dbl wd PC
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,
RETF: RESFLG T2 ;[3243] clear the PC flags
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, 1986
;ALL RIGHTS RESERVED.
;
;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
CAMN T0,[200140,,0] ;IS IT EXACTLY 1.0?
JUMPE T1,POPRET ;YES. JUST RETURN EXACTLY 1.0
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: SKIPGE P1,T2 ;[4020] GET |W|
DMOVN P1,T2 ;[4020] THE HARD WAY (IF NEGATIVE)
SETZ P2, ;ZERO P2
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: CAIG P1,MSKTOP ;[4002] BEYOND TABLE?
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
SKIPGE T0,T4 ;[4020] SAVE A COPY OF ABS(W)
DMOVN T0,T4 ;[4020] THE HARD WAY (IF NEGATIVE)
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: CAIG P1,MSKTOP ;[4002] BEYOND TABLE?
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)
SKIPGE T4,T2 ;[4020] SAVE A COPY OF ABS(W2)
DMOVN T4,T2 ;[4020] THE HARD WAY (IF NEGATIVE)
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: CAIG P1,MSKTOP ;[4002] BEYOND TABLE?
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
MSKTOP==MSK26-MSK3 ;MAX INDEX TO USE FOR MSK3
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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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. G-FLOATING REMAINDER
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1983, 1986
;ALL RIGHTS RESERVED.
;
;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
; This routine was rewritten during edit 3233.
; G-floating MOD routine with no limit on the magnitudes
; of the arguments. This routine calculates
;
; |ARG1| - [|ARG1/ARG2|]*|ARG2|
;
; and returns the result with the sign of ARG1. Special code is
; provided for avoiding intermediate exceptions. The final result
; will not overflow. It may underflow, in which case a result of 0
; is returned and an error message is written.
;
; For A and B positive floating point numbers, let A = 2^n*F and
; B = 2^m*G. We would like to compute R = A - B*int[A/B]. We introduce
; the following sequences of numbers to assist in computing R. For
; p and d positive integers, let J = 2^p*G, I = 2^p*F and define
; I(0) = I - q*J, where q = 0 if I < J and 1 otherwise. For i >= 0 let
;
; L(i+1) = 2^d*I(i)
;
; I(i+1) = L(i+1) - J*int[L(i+1)/J].
;
; Using an induction argument it is not difficult to show that
;
; I(k) = 2^(kd)*I(0) - J*int[2^(kd)*I(0)/J].
;
; Now let 0 = < n - m = kd + r, where 0 =< r < d, then
;
; R = A - B*int[A/B]
; = 2^n*F - 2^m*G*int[2^(n-m)*F/G]
; ==> 2^p*R = 2^n*[I(0) + q*J] - 2^m*J*int[2^(n-m)*(I(0) + q*J)/J]
; = 2^n*I(0) - 2^m*J*int[2^(n-m)*I(0)]
; ==> 2^(p-m)*R = 2^(kd+r)*I(0) - J*int[2^(kd+r)*I(0)/J]
; = 2^r*I(k) - J*int[2^r*I(k)/J]
; = L - J*int[L/J],
;
; where L = 2^r*I(k).
;
; Based on the above, we have the following algorithm for computing R:
;
; Step 1: m <--- the exponent value of B
; c <--- the exponent value of A - m
; if c < 0, end with R = A
; Step 2: I <--- 2^p*(fraction field of A)
; J <--- 2^p*(fraction field of B)
; If I >= J, I <--- I - J
; go to step 5
; Step 3: L <--- 2^d*I
; Step 4: L <--- L - J*int[L/J]
; Step 5: c <--- c - d
; Step 6: if c > 0 go to step 3
; Step 7: if c = -d, exit with R = L
; Step 8: L <--- 2^(d+c)*I = 2^r*I
; Step 9: L <--- L - J*int[L/J]
; Step 10: R <--- 2^(m-p)*L
SEARCH MTHPRM
SEGMENT CODE
SALL
T6=6 ;Additional AC names
T7=7
T8=8
DEFINE DMOVM (AC,E) < ;Double absolute value
DMOVE AC,E
TLNE AC,(1B0)
DMOVN AC,AC
> ;DMOVM
DEFINE DCAML (AC,E) < ;Double integer compare
CAMN AC, E
CAMGE AC+1, E+1
CAMLE AC, E
> ;DCAML
HELLO (GMOD,.)
DMOVEM T2, SAV23 ;Save registers 2, 3,
DMOVEM T4, SAV45 ; 4, 5
DMOVEM T6, SAV67 ; 6, 7
MOVEM T8, SAV8 ; and 8
DMOVM T0, @0(L) ; T0 = |A|
DMOVM T4, @1(L) ; T4 = |B|
JUMPE T4,GMRETZ ;IF ARG2=0, GO RETURN ZERO
;
; Step 1
;
MOVE T6, T4 ; T6 = |B|hi
AND T6, [777700000000]
; T6 = Exponent field of B (including bias)
MOVE T7, T0 ; T7 = |A|hi
SUB T7, T6 ; High 12 bits of T7 = c
JUMPL T7, TSTSGN ; Done if c < 0
;
; Step 2
;
STEP2: ASHC T6, -30 ; Get c in the low bits of T7
; and m+2000 in low bits of T6
TLZ T0, 777700 ; T0 = I
TLZ T4, 777700 ; T4 = J
DCAML T0, T4 ; Compare I to J
DSUB T0, T4 ; If I >= J, I <--- I - J
JRST STEP5
;
; Steps 3 through 6
;
STEP3: SETZB T2, T3 ; T0/T3 = L = 2^35*I -- d = 35
DDIV T0, T4 ; T0 = int(L/J), T2 = L - J*int(L/J)
DMOVE T0, T2 ; T0 = L - J*int(L/J)
STEP5: SUBI T7, 106 ; c <--- c - d
JUMPG T7, STEP3 ; If c > 0 go to Step 3
;
; Steps 7 and 8 9: At this point c = r - d or r = c + d;
;
SETZB T2, T3 ; T0/T3 = 2^35*I
; Shift T0/T3 right T7 places
CAMG T7, [-43] ; More than 1 word to shift
JRST [EXCH T1, T2 ; Yes, do first 35 bits with word operations
EXCH T0, T1
MOVN T8, T7 ; Get negative shift count
ASHC T2, 43(T7) ; Shift T2-T3
ASH T2, -43(T8) ; Reposition T2
ASHC T1, 43(T7) ; Shift T1-T2
JRST STEP8]
MOVN T8, T7 ; Get negative shift count
ASHC T1, (T7) ; Shift T1-T2
ASH T1, (T8) ; Reposition T1
ASHC T0, (T7) ; Shift T0-T1
STEP8: DDIV T0, T4 ; T2 = 2^(p-m)*R
;
; Step 9 - Obtain R in floating point format and check for underflow
;
DMOVE T0, T2 ; Copy fraction into T0
EXTEND T0, [GFSC (T6)] ; Insert biased exponent and normalize
JFCL UNDER ; Can underflow
TSTSGN: SKIPGE @0(L) ;Remainder in T0. If A < 0
DMOVN T0,T0 ; negate it
RESTOR: DMOVE T2,SAV23 ;Restore registers 2, 3
DMOVE T4,SAV45 ; 4, 5
DMOVE T6,SAV67 ; 6, 7
MOVE T8,SAV8 ; and 8
POPJ P, ; Return
;
; If processing continues here, the remainder has underflowed.
;
UNDER: $LCALL RUN ;Result underflow
SETZB T0, T1 ;Store 0 for result
JRST RESTOR ;Restore registers and return
;
;IF ARG2 IS ZERO, RETURN 0 WITH A WARNING MESSAGE
;
GMRETZ: SETZB T0,T1 ;RETURN ZERO
$LCALL MZZ ;WITH MESSAGE
JRST RESTOR
SEGMENT DATA
SAV23: BLOCK 2
SAV45: BLOCK 2
SAV67: BLOCK 2
SAV8: BLOCK 1
PRGEND
TITLE GCOS SINE AND COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MAY 1, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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 GEXP EXPONENTIAL FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL.INC APRIL 3, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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 GINT DOUBLE PRECISION TRUNCATION TO INTEGER
SUBTTL CHRIS SMITH/CKS 29-Jan-80
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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.+0> ;[4015]
PRGEND
TITLE GMAX1 DOUBLE PRECISION MAXIMUM OF A SERIES OF ARGUMENTS
SUBTTL CHRIS SMITH/CKS 29-Jan-80
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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.+0> ;[4015]
PRGEND
TITLE GMIN1 DOUBLE PRECISION MINIMUM OF A SERIES OF ARGUMENTS
SUBTTL CHRIS SMITH/CKS 29-Jan-80
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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.+0> ;[4015]
PRGEND
TITLE GSIGN DOUBLE PRECISION TRANSFER OF SIGN FUNCTION
SUBTTL CHRIS SMITH/CKS 29-Jan-80
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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.+0> ;[4015]
PRGEND
TITLE DTOGA
SUBTTL M. R. BOUCHER/MRB 11-JUN-84
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1973, 1986
;ALL RIGHTS RESERVED.
;
;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 DTOGA
EXTERN DTOGA.
DTOGA=DTOGA.
PRGEND
TITLE DTOG DOUBLE PRECISION TO GFLOAT DOUBLE PRECISION CONVERSION
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1986
;ALL RIGHTS RESERVED.
;
;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.) ;[4012]
;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 GTODA
SUBTTL M. R. BOUCHER/MRB 11-JUN-84
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1973, 1986
;ALL RIGHTS RESERVED.
;
;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 GTODA
EXTERN GTODA.
GTODA=GTODA.
PRGEND
TITLE GTOD GFLOAT DOUBLE PRECISION TO DOUBLE PRECISION CONVERSION
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1986
;ALL RIGHTS RESERVED.
;
;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.) ;[4012]
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
;
;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, 1986
;ALL RIGHTS RESERVED.
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, 1986
;ALL RIGHTS RESERVED.
;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, 1986
;ALL RIGHTS RESERVED.
;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, 1986
;ALL RIGHTS RESERVED.
;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, 1986
;ALL RIGHTS RESERVED.
;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, 1986
;ALL RIGHTS RESERVED.
;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