Trailing-Edge
-
PDP-10 Archives
-
BB-4157E-BM
-
fortran-ots-debugger/fordbl.mac
There are 3 other files named fordbl.mac in the archive. Click here to see a list.
SEARCH FORPRM
TV FORDBL DOUBLE PRECISION ROUTINES ,6(2031)
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 *****
1101 SWG 16-Aug-79
Cleanup FORDAR for V6 - take out F40 and KA conditionals
remove following macros: FLADD, DFA, DAM, DFS, DSM
FLMUL, DFM, DMM, FLDIV, DFD, DDM, all of which are KA
only library routines. Replace DFN in IDF macro
with DMOVN. remove KI conditionals because no
longer necessary
1351 EDS 16-Mar-81 Q10-04786
Fix TWOSEG and RELOC problems.
1405 DAW 6-Apr-81
Extended addressing support.
1464 DAW 12-May-81
New error message format.
1673 CDM 9-Sept-81
Added new double precision routines.
(DDIM, DINT, DNINT, DPROD, IDNINT)
1713 BL 15-Sep-81
Added 'DFL.3'.
1721 CDM 15-Sep-81
Edited DDIM to eliminate DDIM(-INF,+INF) overflow message.
1724 BL 17-Sep-81
Clean up error messages.
1746 CDM 25-Sep-81
Simple change to DDIM.
***** End Revision History *****
\
PRGEND
TITLE DACOS ARC SINE AND ARC COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 19, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DACOS
EXTERN DACOS.
DACOS=DACOS.
PRGEND
TITLE DASIN ARC SINE AND ARC COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 19, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DASIN
EXTERN DASIN.
DASIN=DASIN.
PRGEND
TITLE DASIN. ARC SINE AND ARC COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 19, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 + 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:
; XMOVEI L,ARG
; PUSHJ P,DASIN
; OR
; PUSHJ P,DACOS
;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 FORPRM
TWOSEG 400000
SALL
HELLO (DACOS,.) ;ENTRY TO DACOS ROUTINE
HRRZI T0,2 ;SET FLAG TO TWO
MOVEM T0,FLAG ;MOVE FLAG TO MEMORY
JRST GETX ;GO TO GETX
HELLO (DASIN,.) ;ENTRY TO DASIN 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: CAMGE T0,ONE ;IF Y IS .LT. ONE
JRST ALG ;GO TO ALG
CAMN T0,ONE ;IF Y IS .GT. ONE, GO TO ERR
JUMPE T1,ALG
ERR0: LERR (LIB,%,DASIN or DACOS; ABS(arg) > one; 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,201400 ;SET T2,T3 TO
SETZ T3, ;ONE
DFSB T2,T0 ;G = 1-Y
FSC T2,-1 ;* 1/2
DMOVEM T2,TEMP ;SAVE A COPY OF G
FUNCT DSQRT.,<TEMP> ;Y = DSQRT(G)
FSC T0,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,TNM11H ;IF Y < TNM11H
JRST GOTRES ;GO TO GOTRES
DMOVE T2,T0 ;SAVE A COPY OF Y
DFMP T2,T2 ;G = Y*Y
GOTY: DMOVE T4,T2 ;SAVE A COPY OF G
DFAD T4,Q4 ;XDEN = G + Q4
DFMP T4,T2 ;*G
DFAD T4,Q3 ;+Q3
DFMP T4,T2 ;*G
DFAD T4,Q2 ;+Q2
DFMP T4,T2 ;*G
DFAD T4,Q1 ;+Q1
DFMP T4,T2 ;*G
DFAD T4,Q0 ;+Q0
DMOVEM T2,TEMP ;SAVE A COPY OF G
DFMP T2,RP5 ;XNUM = G*RP5
DFAD T2,RP4 ;+RP4
DFMP T2,TEMP ;*G
DFAD T2,RP3 ;+RP3
DFMP T2,TEMP ;*G
DFAD T2,RP2 ;+RP2
DFMP T2,TEMP ;*G
DFAD T2,RP1 ;+RP1
DFMP T2,TEMP ;*G
DFDV T2,T4 ;RESULT = XNUM/XDEN
DFMP T2,T0 ;*Y
DFAD T0,T2 ;+Y
GOTRES: MOVE T5,I ;GET A COPY OF I
SKIPE FLAG ;IF FLAG IS .NE. 0
JRST NE ;GO TO NE
DFAD T0,A(T5) ;RESULT = RESULT + A(I)
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
DFAD T0,A(T5) ;+A(I)
POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4
POP P,T3
POP P,T2
GOODBY (1) ;RETURN
NEGX: DFAD T0,B(T5) ;RESULT = RESULT + B(I)
POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4
POP P,T3
POP P,T2
GOODBY (1) ;RETURN
RP1: DOUBLE 572112065225,361372635617 ;-.27368494524164255994D+2
RP2: DOUBLE 206711524715,202552212153 ;.57208227877891731407D+2
RP3: DOUBLE 571302372325,226222540767 ;-.39688862997504877339D+2
RP4: DOUBLE 204504702731,073034447126 ;.10152522233806463645D+2
RP5: DOUBLE 577233210222,206360633365 ;-.69674573447350646411D+0
Q0: DOUBLE 567267447760,165074063632 ;-.16421096714498560795D+3
Q1: DOUBLE 211641111704,007534102703 ;.41714430248260412556D+3
Q2: DOUBLE 566202106100,352253670303 ;-.38186303361750149284D+3
Q3: DOUBLE 210455717445,226216157457 ;.15095270841030604719D+3
Q4: DOUBLE 572202642744,101474740606 ;-.23823859153670238830D+2
TNM11H: 134537657770 ;HIGH ORDER PART OF 1.D-11
A: DOUBLE 0,0
A1: DOUBLE 201622077325,021026430215 ;PI/2
B: DOUBLE 202622077325,021026430215 ;PI
B1: DOUBLE 201622077325,021026430215 ;PI/2
ONE: 201400000000 ;1.0
HALF: 200400000000 ;1/2
RELOC ;DATA
TEMP: DOUBLE 0,0
I: 0
FLAG: 0
RELOC
PRGEND
TITLE DATAN2 ARC TAN FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MAY 9, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DATAN2
EXTERN DATN2.
DATAN2=DATN2.
PRGEND
TITLE DATAN ARC TAN FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MAY 9, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DATAN
EXTERN DATAN.
DATAN=DATAN.
PRGEND
TITLE DATAN. ARC TAN FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
;COPYRIGHT (C) 1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
TWOSEG 400000
SALL
T6=6 ;MORE AC NAMES
T7=7
T8=10
HELLO (DATAN,.) ;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,9 ;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-2000+24(T2) ;GET OFFSET INTO XHI TABLES
DMOVE T3,T0 ;GET A COPY OF X
DFSB T0,XHI(T2) ;GET X-XHI
DFMP T3,XHI(T2) ; X*XHI
DFAD T3,ONE ; 1 + X*XHI
DFDV 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
DFMP T3,T3 ;GET Z**2
DMOVE T5,P06 ;GET P(Z**2)
DFMP T5,T3
DFAD T5,P05
DFMP T5,T3
DFAD T5,P04
DFMP T5,T3
DFAD T5,P03
DFMP T5,T3
DFAD T5,P02
DFMP T5,T3
DFAD T5,P01
DFMP T3,T5
DFMP T3,T0 ; * Z
DFAD T0,T3 ; + Z = ATAN(Z)
SMALLX: DFAD T0,ATANLO(T2) ; + ATAN(XHI) LOW
DFAD 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
DFDV 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 (DATN2.,DATAN2) ;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
DFDV 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,9 ;GET INDEX INTO OFFSET TABLES
ASHC T2,3
HLRZ T2,OFFSET-2000+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
DFSB T7,T5 ;GET LOW 35 BITS OF V = VLO
DFMP T5,XHI(T2) ;GET V*XHI = VHI * XHI
DFMP T7,XHI(T2) ; + VLO * XHI
DFSB T0,T5 ;GET (U - VHI*XHI)
DFSB T0,T7 ; - VLO*XHI
DMOVE T5,USAVE ;GET U BACK
DFMP T5,XHI(T2) ;GET U * XHI
DFAD T3,T5 ;GET V + U*XHI
DFDV 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
DFDV 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
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
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 175657740000,0 ; .1054534912109375000 X1
EXP 176407740000,0 ; .1288757324218750000 X2
EXP 176477740000,0 ; .1562194824218750000
EXP 176617640000,0 ; .1952209472656250000
EXP 176777400000,0 ; .2497558593750000000
EXP 177477540000,0 ; .3121948242187500000
EXP 177617200000,0 ; .3898925781250000000
EXP 177776300000,0 ; .4984130859375000000
EXP 200515740000,0 ; .6522216796875000000
EXP 200674040000,0 ; .8673095703125000000
EXP 201453440000,0 ; 1.170166015625000000
EXP 201645100000,0 ; 1.645019531250000000
EXP 202511040000,0 ; 2.570800781250000000
EXP 203527000000,0 ; 5.359375000000000000 X14
ATANHI: EXP 000000000000,000000000000 ; .0000000000000000000 ATAN(0)
EXP 175656261520,251717717115 ; .1050651824695016827 ATAN(X1)
EXP 176406373154,305103547431 ; .1281692623534198424
EXP 176475276500,226075230141 ; .1549669515088296698
EXP 176612661305,353341503017 ; .1927961221331547523
EXP 176765175625,123431412160 ; .2447488705187413410
EXP 177465675077,207755453362 ; .3026068193134274483
EXP 177574536624,346575570716 ; .3717628303965550497
EXP 177731362665,337061624571 ; .4623772720674932799
EXP 200447716236,341667321523 ; .5779354489017924548
EXP 200555632634,072362412013 ; .7144577222199199908
EXP 200672140433,153152445005 ; .8636495725719767213
EXP 201406227204,252435004161 ; 1.024591515455200379
EXP 201463117110,112456461235 ; 1.199822549393342831
EXP 201542714675,142761043345 ; 1.386328658261967261 ATAN(X14)
PI2: EXP 201622077325,021026430215 ; 1.570796326794896619 PI/2
MPI: EXP 575155700452,356751347563 ;-3.141592653589793238 -PI
EXP 575173246105,164207746145 ;-3.036527471120291556 -PI + ATAN(X1)
EXP 575176220221,273215536144 ;-3.013423391236373396 -PI + ATAN(X2)
EXP 575201554376,370255221171 ;-2.986625702080963569
EXP 575206433527,115527433724 ;-2.948796531456638486
EXP 575215150344,104133030272 ;-2.896843783071051897
EXP 575224470162,337747115121 ;-2.838985834276365790
EXP 575235354335,213631126655 ;-2.769829823193238188
EXP 575251036741,252657532242 ;-2.679215381522299958
EXP 575267664122,247327234107 ;-2.563657204688000784
EXP 575311247221,375446052165 ;-2.427134931369873248
EXP 575334330561,311604060764 ;-2.277943081017816517
EXP 575361014155,104167751653 ;-2.117001138134592860
EXP 576016720236,050401400603 ;-1.941770104196450407
EXP 576076516023,100703762713 ;-1.755263995327825977 -PI + ATAN(X14)
EXP 576155700452,356751347563 ;-1.570796326794896619 -PI/2
ATANLO: EXP 000000000000,000000000000 ; .0000000000000000000
EXP 701136265552,021570267746 ;-.1105497383317563340E-19
EXP 075632563733,177450755317 ; .5435918950106886920E-20
EXP 077603757522,264676361123 ; .2053885961335437547E-19
EXP 077723405617,212207745036 ; .2474984160195000159E-19
EXP 700044202740,142243000710 ;-.2518569148307390217E-19
EXP 100702445617,347031410060 ; .4770635578227287401E-19
EXP 100414524241,267545144713 ; .2844597940229199763E-19
EXP 100625025353,262726507370 ; .4288548085099372680E-19
EXP 676255576133,135473517361 ;-.7162797697818810113E-19
EXP 676323642342,201773750626 ;-.6356616556133364840E-19
EXP 077413353720,052424757433 ; .1415925447541254247E-19
EXP 100744743210,367241225634 ; .5134543068800780929E-19
EXP 101423305250,056513174673 ; .5831512827058072309E-19
EXP 102617464410,310763545272 ; .1692382723892392993E-18
EXP 101611431424,134015603464 ; .8333742918520878328E-19
EXP 675166346353,243762174314 ;-.1666748583704175666E-18
EXP 102634261642,105051607712 ; .1746358738541957411E-18
EXP 103601511025,077765605721 ; .3266520381981663157E-18
EXP 676115710653,365124065054 ;-.9192589013278796940E-19
EXP 675060707135,225203171017 ;-.1961351253927427867E-18
EXP 675072766647,260206474405 ;-.1918605498534914687E-18
EXP 101716137637,073361174657 ; .9787193190895619423E-19
EXP 674134635612,010745612637 ;-.3550693134652264558E-18
EXP 700335536464,205477162116 ;-.1536916027087339636E-19
EXP 103746522614,251310022042 ; .4122184681426969927E-18
EXP 103760133656,162370070313 ; .4202802795595514454E-18
EXP 101457607713,122451564536 ; .6432483060209586270E-19
EXP 103567657506,360715220731 ; .3183514413117920163E-18
EXP 676000222177,166457565523 ;-.1083597300998368435E-18
EXP 074603276433,074574160521 ; .2563414018821732668E-20
EXP 676166346353,243762174314 ;-.8333742918520878328E-19
;COEFFICIENTS OF APPROXIMATION POLYNOMIAL ATAN(X) = X*P(X**2)
P06: EXP 175461762413,330264551705 ; .0747006049800000000
P05: EXP 602213603465,225464265500 ;-.0908796288218500000
P04: EXP 175707070366,207135266005 ; .1111109168530032000
P03: EXP 601333333333,310141050127 ;-.1428571421988482590
P02: EXP 176631463146,146202001672 ; .1999999999989370798
P01: EXP 600252525252,252525266207 ;-.3333333333333326895
ONE: EXP 201400000000,000000000000 ; 1.000000000000000000
EPS: EXP 142561156421 ;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 175623327343 ;TAN(PI/32) (HIGH WORD, ROUNDED DOWN)
MAXX: EXP 204504715427 ;1/TAN(PI/32) (HIGH WORD, ROUNDED UP)
RELOC
SGNFLG: BLOCK 1 ;SIGN TO BE ATTACHED TO RESULT
USAVE: BLOCK 2 ;TEMP STORAGE FOR U
RELOC
PRGEND
TITLE DCOSH HYPERBOLIC COSINE FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 7, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DCOSH
EXTERN DCOSH.
DCOSH=DCOSH.
PRGEND
TITLE DCOSH. HYPERBOLIC COSINE FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 7, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 <= 88.029678, DCOSH = (EXP(Y)+EXP(-Y))/2,
; 88.029678 <= Y < 128 * LN(2)
; DCOSH = (V/2)*EXP(Y - LN V),
; Y >= 128 * LN(2), DCOSH = +MACHINE INFINITY
;THE RANGE OF DEFINITION FOR DCOSH IS ABS(X) < 128 * LN(2) AND ERROR MESSAGES
; WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE. DCOSH WILL BE SET
; TO + MACHINE INFINITY.
;REQUIRED (CALLED) ROUTINES: DEXP
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; XMOVEI L,ARG
; PUSHJ P,DCOSH
;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 FORPRM
TWOSEG 400000
SALL
HELLO (DCOSH,.) ;ENTRY FOR DCOSH 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 ERRD ;GO TO ERRD
CAMGE T1,YMAXLO ;HI OF Y = YMAXHI
JRST GETW ; IF LO OF Y .LE. YMAXLO, GO TO GETW
ERRD: HRLOI T0,377777 ;SET RESULT TO
SETO T1, ;MACHINE INFINITY
LERR (LIB,%,DCOSH: result overflow)
GOODBY (1) ;RETURN
GETW: DFAD T0,LNV ;W = Y-LNV
DMOVEM T0,TEMP ;MOVE W TO A TEMP REGISTER
FUNCT DEXP.,<TEMP> ;Z = EXP(W)
DMOVEM T0,TEMP ;SAVE A COPY OF Z
FMP T0,CON1 ;RESULT = Z*CON1
SETZ T1, ;ZERO SECOND WORD
DFAD T0,TEMP ;+Z
JFCL ERRD
GOODBY (1) ;RETURN
EASY: DMOVEM T0,TEMP ;MOVE Y TO TEMP
FUNCT DEXP.,<TEMP> ;GET ITS EXPONENT
FSC T0,-1 ;DIV BY 2
GOODBY (1) ;RETURN
ALG2: CAMG T0,TM32 ;IF ABS(X) < 2**(-32)
JRST TINY ; THE RESULT = 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 DEXP.,<TEMP> ;Z = DEXP(Y)
DMOVE T2,T0 ;SAVE A COPY OF Z
HRLZI T4,201400 ;SET T4 TO 1.0
SETZ T5, ;CLEAR SECOND WORD
DFDV T4,T2 ;1/Z
DFAD T0,T4 ;+Z
FSC T0,-1 ;*1/2
POP P,T5
POP P,T4 ;RESTORE ACCUMULATORS
POP P,T3
POP P,T2
GOODBY (1) ;RETURN
TINY: HRLZI T0,201400 ;RESULT IS 1.
SETZ T1,
GOODBY (1) ;RETURN
ONE: 201400000000 ;1.0
BIGX: 207540074635 ;88.029691
YMAXHI: 207542710277 ;88.722839111672999604
YMAXLO: 276434757153 ;
LNV: DOUBLE 577235067500,101343000000 ;-.693147180559947174D0
CON1: 120414520000 ;.186417721E-14
TWENT2: 205540000000 ;22
TM32: 141400000000 ;2**(-32)
RELOC ;DATA
TEMP: DOUBLE 0,0
RELOC
PRGEND
TITLE DEXP EXPONENTIAL FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL.INC APRIL 3, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DEXP
EXTERN DEXP.
DEXP=DEXP.
PRGEND
TITLE DEXP. EXPONENTIAL FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. APRIL 3, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 <= -89.415986292232944914, EXP = 0
; IF X > 88.029691931113054295, 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:
; XMOVEI L,ARG
; PUSHJ P,DEXP
;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 FORPRM
TWOSEG 400000
SALL
HELLO (DEXP,.) ;ENTRY TO DEXP 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 ERRD ; GO TO ERR
CAMG T1,BIGXLO ;IF LO OF X .LE. BIGXLO,
JRST EXP1 ; GO TO EXP1
ERRD: LERR (LIB,%,DEXP: result overflow)
HRLOI T0,377777 ;DEXP = +MACHINE INFINITY
HRLOI T1,377777
GOODBY (1) ;RETURN
OUT2: LERR (LIB,%,DEXP: result underflow)
MOVEI T0,0 ;EXP = 0
MOVEI T1,0
GOODBY (1) ;RETURN
EXP1: JUMPE T0,[MOVSI T0,(1.0) ;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 |ARGUMENT| IN T2
CAML T2,LN2OV2 ;IS IT .LT. (1/2)*LN(2)?
JRST REDUCE ;IF SO, IT MUST BE REDUCED.
SETZ T4, ;MAKE SCALE INDEX 0.
MOVEM T4,SAVEN ; AND SAVE IT.
DMOVE T4,T0 ;GET COPY OF ARGUMENT.
JRST MERGE ;MERGE WITH MAIN FLOW.
REDUCE: DMOVE T2,T0 ;GET COPY OF ARG
DFMP T2,RNDLN2 ;ARG/LN2
FIXR T2,T2 ;NEAREST INTEGER = N
FLTR T2,T2 ;FLOAT N
MOVEI T3,0 ;
FIX T4,T0 ;[X]
FLTR T4,T4 ;X1 = FLOAT[X]
MOVEI T5,0
DFSB T0,T4 ;X2 = X - X1
MOVEM T2,SAVEN ;SAVEN
DFMP T2,C1 ;N*C1
DFAD T4,T2 ;X1 + (N*C1)
DFAD T4,T0 ;+X2
MOVE T2,SAVEN ;RETRIEVE N
MOVEI T3,0 ;ZERO SECOND WORD
DFMP T2,C2 ;FORM N*C2
DFAD T4,T2 ;(N*C2)+X1+(N*C1)+X2
DMOVE T0,T4 ;SAVE G
MOVM T2,T4 ;GET ABS(G)
MERGE: CAML T2,TWOM32 ;IF REDUCED ARG IS>= 2**-32
JRST APPRX ; GO TO APPRX
DFAD T0,ONE ;R(G) = 1. + G
FSC T0,-1 ;*.5
JRST BRNCH ;GO TO BRNCH
APPRX: DFMP T4,T4 ;Z = G*G
DMOVE T2,T4 ;SAVE Z
DFMP T4,XP2 ;Z*XP2
DFAD T4,XP1 ;+XP1
DFMP T4,T2 ;* Z
DFAD T4,XP0 ;+ XP0
DFMP T0,T4 ;* G
DMOVE T4,T2 ;SAVE Z
DFMP T4,XQ3 ;XQ3*Z
DFAD T4,XQ2 ;+XQ2
DFMP T4,T2 ;*Z
DFAD T4,XQ1 ;+ XQ1
DFMP T4,T2 ;* Z
DFAD T4,XQ0 ; + XQ0
DFSB T4,T0 ; XQ - G*XP
DFDV T0,T4 ;(G*XP)/(XQ-G*XP)
DFAD T0,XQ0 ; + .5
BRNCH: MOVE T4,SAVEN ;RETRIEVE N
FIX T4,T4
ADDI T4,1 ;N = N+1
FSC T0,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: 570232254036 ;-89.415986292232944914
SMLXLO: 301750634730 ;
BIGXHI: 207540074636 ;88.029691931113054295
BIGXLO: 077042573256 ;
LN2OV2: 177542710300 ;(1/2) * LN(2)
RNDLN2: DOUBLE 201561250731,112701376057 ;1.44269504088896341 = 1/LN2
ONE: DOUBLE 201400000000,000000000000 ;1.0D0
TWOM32: DOUBLE 141400000000,000000000000 ;2**-32
C1: DOUBLE 577235000000,000000000000 ;-0.693359375D0
C2: DOUBLE 164675002027,030206250331 ;2.12194440054690583D-4
XP0: DOUBLE 177400000000,000000000000 ;0.250D0
XP1: DOUBLE 171760351374,212411245446 ;0.757531801594227767D-2
XP2: DOUBLE 162410550412,271511036101 ;0.315551927656846464D-4
XQ0: DOUBLE 200400000000,000000000000 ;0.5D0
XQ1: DOUBLE 174721345024,167754776545 ;0.568173026985512218D-1
XQ2: DOUBLE 166512741427,012152337316 ;0.631218943743985036D-3
XQ3: DOUBLE 154623154303,071202125117 ;0.751040283998700461D-6
RELOC ;DATA
SAVEN: 0
RELOC
PRGEND
TITLE DEXP2. DOUBLE ** INTEGER EXPONENENTIATION
SUBTTL CHRIS SMITH/CKS 28-Jan-80
;COPYRIGHT (C) 1980,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
TWOSEG 400000
SALL
HELLO (DEXP2,.)
PUSH P,T2 ;SAVE TEMP REGISTERS
PUSH P,T3
PUSH P,T4
DMOVE T0,[EXP 201400000000,0] ;SET RESULT TO 1
MOVM T4,@1(L) ;GET ABS(EXPONENT)
JUMPE T4,D2RET ;EXPONENT ZERO, GO RETURN 1
DMOVE T2,@0(L) ;GET BASE
JUMPN T2,D2POS ;BASE NOT ZERO, GO GET BASE**ABS(EXP)
JRST D2ZERO ;BASE ZERO, GO HANDLE
D2LP: DFMP T2,T2 ;SQUARE BASE
JFCL D2OVUN ;CATCH OVERFLOW
D2POS: TRNE T4,1 ;CHECK LOW BIT OF EXPONENT
DFMP T0,T2 ;MULTIPLY ANSWER BY BASE
JFCL D2OVUN ;CATCH OVERFLOW
LSH T4,-1 ;DIVIDE EXPONENT BY 2
JUMPN T4,D2LP ;HANDLE ALL BITS OF EXPONENT
SKIPL @1(L) ;NEGATIVE EXPONENT?
JRST D2RET ;NO, WE HAVE RESULT
DMOVE T2,[EXP 201400000000,0] ;YES, GET RECIPROCAL
DFDV T2,T0
DMOVE T0,T2
D2RET: POP P,T4 ;RESTORE TEMP REGISTERS
POP P,T3
POP P,T2
POPJ P, ;DONE
D2ZERO: SKIPL T4,@1(L) ;BASE 0, RESULT IS 0 IF EXPONENT POSITIVE
JRST D2RTZ ;POSITIVE EXPONENT, GO RETURN 0
JRST D2OVFL ;ELSE RESULT OVERFLOWS
D2OVUN: SKIPG T4,@1(L) ;NEGATIVE EXPONENT?
JRST D2UNFL ;YES, RESULT UNDERFLOWS
D2OVFL: LERR (LIB,%,DEXP2: result overflow)
HRLOI T0,377777 ;OVERFLOW, GUESS POSITIVE RESULT
HRLOI T1,377777
SKIPL @0(L) ;CHECK SIGN OF BASE
JRST D2RET ;BASE NONNEGATIVE, GO RETURN +INFINITY
TRNE T4,1 ;ODD EXPONENT?
DMOVN T0,T0 ;NEGATIVE**ODD, RETURN -INFINITY
JRST D2RET
D2UNFL: LERR (LIB,%,DEXP2: result underflow)
D2RTZ: SETZB T0,T1 ;RETURN 0
JRST D2RET
PRGEND
TITLE DEXP3. POWER FUNCTION
; (DOUBLE PRECISION)
SUBTTL IMSL, INC. JUNE 4, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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.
;DEXP3 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.)
; -128.9375 < FLOAT(INT((Y*LOG2(X))*16))/16 < 127.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 DEXP3 IS GIVEN ABOVE, AND ERROR MESSAGES
; WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE
;REQUIRED (CALLED) ROUTINES: NONE
;REGISTERS T2, T3, T4, T5, G1, G2, G3, AND G4 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; XMOVEI 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 FORPRM
TWOSEG 400000
SALL
HELLO (DEXP3,.) ;ENTRY TO DEXP3. ROUTINE
PUSH P,T2 ;SAVE ACCUMULATORS
PUSH P,T3
DMOVE T0,@(L) ;GET THE BASE
DMOVE T2,@1(L) ;GET THE EXPONENT
SETZM IFLAG ;SET INTEGER EXP FLAG TO 0
JUMPG T0,XOK ;IF X IS NOT GT 0
JUMPE T0,X0 ;AND IF X IS NOT = 0
PUSH P,T4 ;SAVE ACCUMULATORS
PUSH P,T5
PUSH P,G1
PUSH P,G2
DMOVE T4,T2 ;GET A COPY OF Y
JUMPGE T4,YP ;IF Y IS NEGATIVE
DMOVN T4,T4 ;NEGATE IT
YP: MOVE G2,T4 ;GET HIGH ORDER OF Y
SETZ G1, ;SET G1 TO ZERO
LSHC G1,11 ;GET EXPONENT OF Y
SUBI G1,170 ;GET SHIFTING FACTOR
LSHC T4,(G1) ;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,ERR0 ;IF LO OF Y .NE. 0, GO TO ERR
JUMPE T4,YINT ;IF HI OF Y = 0, GO TO YINT
ERR0: LERR (LIB,%,<DEXP3: negative ** non-integer; ABS(base) used instead of base>)
JRST CONT ;GO TO CONT
X0: JUMPG T2,RET1 ;IF Y IS NOT GT 0
LERR (LIB,%,<DEXP3: base is 0 and exponent is le 0; result = infinity>)
HRLOI T0,377777 ;RESULT = INFINITY
HRLOI T1,377777
POP P,T3 ;RESTORE ACCUMULATORS
POP P,T2
GOODBY (1) ;RETURN
XOK: PUSH P,T4 ;SAVE ACCUMULATORS
PUSH P,T5
PUSH P,G1
PUSH P,G2
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,G3
PUSH P,G4
MOVE T4,T0 ;OBTAIN THE EXPONENT
ASH T4,-33 ;SHIFT MANTISSA OFF
SUBI T4,200 ;SUBTRACT 128 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 G1,T0 ;SAVE A COPY OF F
DFAD G1,A1(T5) ;F+A1(P+1)
DFSB T0,A1(T5) ;F-A1(P+1)
SUBI T5,2 ;FORM P + 1
ASH T5,-1 ;(P+1)/2
DFSB T0,A2(T5) ;Z=(F-A1(P+1))-A2((P+1)/2)
DFDV T0,G1 ;/(F+A1(P+1))
FSC T0,1 ;
DMOVE G1,T0 ;SAVE A COPY OF Z
DFMP G1,G1 ;FORM Z**2
DMOVE G3,G1 ;SAVE A COPY OF Z**2
DFMP G3,RP4 ;R(Z)=RP4*Z**2
DFAD G3,RP3 ;+RP3
DFMP G3,G1 ;*Z**2
DFAD G3,RP2 ;+RP2
DFMP G3,G1 ;*Z**2
DFAD G3,RP1 ;+RP1
DFMP G3,G1 ;*Z**2
DFMP G3,T0 ;*Z
DFAD G3,T0 ;+Z
DMOVE G1,G3 ;SAVE A COPY OF R(Z)
DFMP G3,C ;U2 = R(Z) * C
DFAD G3,G1 ;+ R(Z)
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, ;
FSC T4,-4 ;
DMOVE T0,T2 ;SAVE A COPY OF Y
JUMPGE T2,YPOS ;IF Y IS NEGATIVE
DMOVN T0,T0 ;NEGATE IT
YPOS: MOVE G1,T0 ;GET COPY OF HIGH ORDER OF Y
ASH G1,-33 ;SHIFT OFF MANTISSA
CAIGE G1,227 ;IF THE EXPONENT IS LESS THAN 227
JRST SMALLY ;GO TO SMALLY
CAIE G1,227 ;IF THE EXPONENT IS > 227
JRST LARGEY ;GO TO LARGEY
SETZ T1, ;ZERO SECOND WORD
JRST SGNCHK ;GO TO SGNCHK
SMALLY: SUBI G1,175 ;GET SHIFTING FACTOR
SETZ T1, ;SET Y1 TO 0
JUMPGE G1,GETY1 ;IF G1 IS .GE. 0
SETZ T0, ;THEN GO TO GETY1
JRST GETY2 ;
GETY1: AND T0,MASK(G1) ;GET Y1
JRST SGNCHK ;GO TO CHECK SIGN
LARGEY: CAIL G1,272 ;IF EXPONENT IS .GE. 272
JRST SGNCHK ;GO TO SGNCHK
SUBI G1,230 ;GET SHIFTING FACTOR
AND T1,REDMSK(G1) ;GET LOW ORDER PART OF Y1
SGNCHK: JUMPGE T2,GETY2 ;IF Y IS NEGATIVE
DMOVN T0,T0 ;NEGATE Y1
GETY2: DMOVE G1,T2 ;SAVE A COPY OF Y
DFSB T2,T0 ;Y2 = Y-Y1
DFMP T2,T4 ;FORM Y2*U1
JFCL
DFMP G1,G3 ;FORM Y*U2
JFCL
DFAD T2,G1 ;W = Y*U2+Y2*U1
DFMP T4,T0 ;FORM U1*Y1
JFCL
DMOVE T0,T2 ;RECONSTRUCT W
DFAD T0,T4
JFCL
CAMG T0,BIGW ;IF W IS NOT TOO BIG
JRST WOK ;GO TO WOK
LERR (LIB,%,DEXP3: 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
LERR (LIB,%,DEXP3: result underflow)
SETZ T0, ; RESULT = 0
SETZ T1,
POP P,G4 ;RESTORE ACCUMULATORS
POP P,G3
POP P,G2
POP P,G1
POP P,T5
POP P,T4
POP P,T3
POP P,T2
GOODBY (1) ;RETURN
WOK2: SETZ G2, ;ZERO G2
MOVM G1,T2 ;GET HIGH ORDER OF W
MOVE G3,G1 ;
ASH G3,-33 ;SHIFT OFF MANTISSA
SUBI G3,175 ;GET SHIFTING FACTOR
JUMPGE G3,GETW1 ;IF G3 .GE. 0
;THEN GO TO GETW1
SETZ G1, ;OTHERWISE, SET W1 = 0
JRST GETW2 ;GO TO GETW2
GETW1: AND G1,MASK(G3) ;GETW1
JUMPG T2,GETW2 ;IF W IS NEGATIVE
DMOVN G1,G1 ;NEGATE W1
GETW2: DFSB T2,G1 ;W2 = W-W1
DFAD T4,G1 ;W = W1+U1*Y1
MOVM T0,T4 ;SAVE A COPY OF ABS(W)
MOVE G1,T0 ;
SETZ T1, ;ZERO T1
ASH G1,-33 ;SHIFT OFF MANTISSA
SUBI G1,175 ;GET SHIFTING FACTOR
JUMPGE G1,GTW1 ;IF G1 IS .GE. 0
;THEN GO TO GTW1
SETZ T0, ;OTHERWISE, SET W1 TO 0
JRST GTW2 ;GO TO GET W2
GTW1: AND T0,MASK(G1) ;GET W1
JUMPGE T4,GTW2 ;IF W IS NEGATIVE
DMOVN T0,T0 ;NEGATE W1
GTW2: DFSB T4,T0 ;FORM W-W1
DFAD T2,T4 ;W2 = W2+(W-W1)
MOVM T4,T2 ;SAVE A COPY OF ABS(W2)
MOVE G1,T4 ;
SETZ T5, ;ZERO T5
ASH G1,-33 ;SHIFT OFF MANTISSA
SUBI G1,175 ;GET SHIFTING FACTOR
JUMPGE G1,GW ;IF G1 IS .GE. 0
;THEN GO TO GW
SETZ T4, ;OTHERWISE, SET W=0
JRST GW2 ;GO TO GW2
GW: AND T4,MASK(G1) ;GET W
JUMPGE T2,GW2 ;IF W2 IS NEGATIVE
DMOVN T4,T4 ;NEGATE W
GW2: DFAD T0,T4 ;FORM W1 + W
FSC T0,4 ;*16
FIX T0,T0 ;IW1
DFSB T2,T4 ;W1 = W2 - W
JUMPLE T2,W2POS ;IF W2 .GT. 0
DFSB 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
DFMP T2,Q7 ;Z = Q7*W2
DFAD T2,Q6 ;+Q6
DFMP T2,T0 ;*W2
DFAD T2,Q5 ;+Q5
DFMP T2,T0 ;*W2
DFAD T2,Q4 ;+Q4
DFMP T2,T0 ;*W2
DFAD T2,Q3 ;+Q3
DFMP T2,T0 ;*W2
DFAD T2,Q2 ;+Q2
DFMP T2,T0 ;*W2
DFAD T2,Q1 ;+Q1
DFMP T0,T2 ;*W2
DFMP T0,A1(T4) ;*A1(P1+1)
DFAD T0,A1(T4) ;+A1(P1+1)
FSC T0,(T5) ;ADD M1 TO THE EXP OF Z
RET3: POP P,G4 ;RESTORE ACCUMULATORS
POP P,G3
RET2: SKIPE IFLAG ;IF Y IS ODD NEGATIVE INTEGER
DMOVN T0,T0 ;NEGATE RESULT
POP P,G2
POP P,G1
POP P,T5
POP P,T4
RET1: POP P,T3
POP P,T2
GOODBY (1) ;RETURN
SXTNTH: DOUBLE 175400000000,000000000000 ;.0625D0
BIGW: 127.0E0 ;UPPER BOUND FOR WW
SMALLW: -128.9375E0 ;LOWER BOUND FOR W
RP1: DOUBLE 175525252525,125252514430 ;.833333333333332114D-1
RP2: DOUBLE 172631463146,147404016032 ;.125000000005037992D-01
RP3: DOUBLE 170444444413,211706653423 ;.223214212859242590D-2
RP4: DOUBLE 165707437566,322576050305 ;.434457756721631196D-3
C: DOUBLE 177705243545,053405770274 ;.442695040888963407D0
Q7: DOUBLE 160764733570,344546177513 ;.149288526805956082D-4
Q6: DOUBLE 164502757270,016046163745 ;.154002904409897646D-3
Q5: DOUBLE 167535417606,065746637065 ;.133335413135857847D-02
Q4: DOUBLE 172473125333,204616630130 ;.961812905951724170D-02
Q3: DOUBLE 174706541065,300601154766 ;.555041086640855953D-1
Q2: DOUBLE 176753767577,340542316147 ;.240226506959095371D0
Q1: DOUBLE 200542710277,276434757056 ;.693147180559945296D0
A1: DOUBLE 201400000000,000000000000 ;A1(I), I=1,17 =
A12: DOUBLE 200752225750,251110331413 ;2**((1-I)/16). THIS
A13: DOUBLE 200725403067,076722207114 ;TABLE IS SEARCHED
A14: DOUBLE 200701463367,141251234104 ;TO DETERMINE P.
A15: DOUBLE 200656423746,126551655275 ;
A16: DOUBLE 200634222140,250770220071 ;
A17: DOUBLE 200612634520,212520333267 ;
A18: DOUBLE 200572042434,372600606660 ;
A19: DOUBLE 200552023631,237635714441 ;
A110: DOUBLE 200532540767,122052051262 ;
A111: DOUBLE 200513773265,115425047073 ;
A112: DOUBLE 200475724623,004432042153 ;
A113: DOUBLE 200460337602,214333425134 ;
A114: DOUBLE 200443417233,235261070316 ;
A115: DOUBLE 200427127017,037250572672 ;
A116: DOUBLE 200413253033,076304417305 ;
A117: DOUBLE 200400000000,000000000000 ;
A2: DOUBLE 077756350307,256663004531
A22: DOUBLE 101406261124,022156463124
A23: DOUBLE 702043260343,345371331424
A24: DOUBLE 676250402176,331550372175
A25: DOUBLE 676111401243,043654640102
A26: DOUBLE 101640442174,047535641316
A27: DOUBLE 676017655543,202562654434
A28: DOUBLE 100613445334,141025236276
MASK1: 000777777777 ;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
RELOC ;DATA
PTEMP: 0
IFLAG: 0
RELOC
PRGEND
TITLE DLOG10 LOG BASE 10 FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. APRIL 8, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DLOG10
EXTERN DLG10.
DLOG10=DLG10.
PRGEND
TITLE DLOG NATURAL LOG FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. APRIL 8, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DLOG
EXTERN DLOG.
DLOG=DLOG.
PRGEND
TITLE DLOG. NATURAL LOG FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. APRIL 8, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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:
; XMOVEI L,ARG
; PUSHJ P,DLOG
; OR
; PUSHJ P,DLOG10
;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 FORPRM
TWOSEG 400000
SALL
HELLO (DLG10.,DLOG10) ;ENTRY TO LOG BASE 10 ROUTINE.
SETZM FLAG ;CLEAR FLAG FOR DLOG10 ENTRY
JRST ALG ;GO TO ALGORITHM
HELLO (DLOG,.) ;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
LERR (LIB,%,DLOG or DLOG10: zero arg; result = -infinity)
HRLZI T0,400000 ;SET RESULT TO
HRRZI T1,000001 ;LARGE NEGATIVE NUMBER
GOODBY (1)
ARGN: 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
MOVE T2,T0 ;GET COPY OF HIGH ORDER PART OF ARG
ASH T2,-33 ;SHIFT OFF THE MANTISSA
SUBI T2,200 ;SUBTRACT 128 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
FLTR T2,T2
MOVEM T2,N ; SAVE N
DFAD T4,MHALF ;ZNUM = F-.5
DMOVE T2,T4 ;GET COPY OF ZNUM
FSC T4,-1 ;ZDEN = ZNUM * .5
DFAD T4,HALF ; + .5
JRST EVALRZ
FGT: FLTR T2,T2
MOVEM T2,N ; SAVE N
DMOVE T2,T4
DFAD T2,B3 ;ZNUM = F - 1.0
FSC T4,-1 ;ZDEN = F*.5
DFAD T4,HALF ; + .5
EVALRZ: DFDV T2,T4 ;Z = ZNUM/ZDEN
DMOVE T4,T2 ;
DFMP T4,T4 ;W = Z*Z
DMOVE T0,T4 ;SAVE COPY OF W
DFAD T4,B2 ; FORM B(W). B(W) = W + B2
DFMP T4,T0 ; * W
DFAD T4,B1 ; + B1
DFMP T4,T0 ; * W
DFAD T4,B0 ; + B0
DMOVEM T0,SAVEW ; SAVE A COPY OF W
DFMP T0,A2 ;FORM A(W). A(W)= A2*W
DFAD T0,A1 ; + A1
DFMP T0,SAVEW ; * W
DFAD T0,A0 ; + A0
DFDV T0,T4 ;R(Z) = A(W)/B(W)
DFMP T0,SAVEW ; * W
DFMP T0,T2 ; *Z
JFCL
DFAD T0,T2 ; + Z
MOVE T2,N ;RETRIEVE N
MOVEI T3,0 ;ZERO OUT SECOND WORD
DFMP T2,C1 ;FORM N*C1
MOVE T4,N ;RETRIEVE N
MOVEI T5,0 ;ZERO OUT SECOND WORD
DFMP T4,C2 ;FORM N*C2
DFAD T0,T4 ;RESULT = N*C2 + R(Z)
DFAD T0,T2 ; + N*C1
SKIPN FLAG
DFMP 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 200543000000,000000000000 ;.693359375D0
C2: DOUBLE 613102775750,347571527443 ;-2.12194440054690583D-4
C3: DOUBLE 177674557305,111562416145 ;.434294481903251828D0
A0: DOUBLE 570377400073,123045145570 ;-.641249434237455811D2
A1: DOUBLE 205406111210,005527754477 ;.163839435630215342D2
A2: DOUBLE 577153575223,052241327104 ;-.789561128874912573D0
B0: DOUBLE 565177200130,374467610432 ;-.769499321084948798D3
B1: DOUBLE 211470020376,205223175724 ;.312032220919245328D3
B2: DOUBLE 571342517755,046030761117 ;-.356679777390346462D2
B3: DOUBLE 576400000000,000000000000 ;-.1D1
HALF: DOUBLE 200400000000,000000000000 ;.5D0
MHALF: DOUBLE 577400000000,000000000000 ;-.5D0
HI: 200552023631
LOW: 237635714441
MASK1: 000777777777
MASK2: 200000000000
RELOC ;DATA
N: 0
FLAG: 0
SAVEW: DOUBLE 000000000000,000000000000
RELOC
PRGEND
TITLE DMOD DOUBLE PRECISION
SUBTTL MARY PAYNE /MHP/CKS 25-Jan-80
;COPYRIGHT (C) 1980,1981 BY DIGITAL EQUIPMENT CORPORATION
;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.
NOSYM
ENTRY DMOD
EXTERN DMOD.
DMOD=DMOD.
PRGEND
TITLE DMOD. DOUBLE PRECISION
SUBTTL MARY PAYNE /MHP/CKS 25-Jan-80
;COPYRIGHT (C) 1980,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
TWOSEG 400000
SALL
T6=6 ;EXTRA TEMP REGISTERS
T7=7
T8=10
T9=11
HELLO (DMOD,.)
DMOVEM T2,SAVE2 ;SAVE TEMP REGISTERS
DMOVEM T4,SAVE4
DMOVEM T6,SAVE6
DMOVEM T8,SAVE8
DMOVE T0,@0(L) ;GET ABS(A) IN T0-T1
CAIGE T0,0
DMOVN T0,T0
REPEAT: DMOVE T2,@1(L) ;GET ABS(B) IN T2-T3
CAIGE T2,0
DMOVN T2,T2
CAMG T0,T2 ;CHECK IF A .LE. B
JRST FAST ;YES, GO BE FAST
DMOVEM T2,B ;SAVE ABS(B)
DMOVE T8,T0 ;SAVE ABS(A) IN T8-T9
;BREAK B INTO HIGH, MIDDLE, AND LOW PARTS
;B IN T2/ XXXHHHHHHHHH MMMMMMMMMLLL
MOVE T6,T2 ;GET T6/ XXXHHHHHHHHH 000000000000
SETZ T7, ;
DFSB T2,T6 ; T2/ XXXMMMMMMMMM LLL000000000
MOVE T4,T2 ; T4/ XXXMMMMMMMMM 000000000000
SETZ T5, ;
DFSB T2,T4 ; T2/ XXXLLL000000 000000000000
DFDV T8,B ;GET A/B
JFCL OVFL ;CATCH OVERFLOW
CAML T8,[244400000000] ;IS QUOTIENT OVER 2**35?
JRST [AND T9,[777000000000] ;YES, TRUNCATE IT TO 35 BITS
JRST BIG] ; TO KEEP PRODUCTS EXACT
FIX T9,T8 ;UNDER 2**35, TRUNCATE TO NEAREST INTEGER
MOVSI T8,276000 ;DFLOAT THE INTEGER
DFAD T8,[EXP 0,0]
BIG: DFMP T6,T8 ;GET A-B*[A/B] IN T0-T1
DFSB T0,T6
DFMP T4,T8
DFSB T0,T4
DFMP T2,T8
DFSB T0,T2
OVCONT: CAIGE T0,0 ;IF RESULT NEGATIVE, [A/B] WAS 1 TOO HIGH
DFAD T0,B ;CORRECT THE RESULT
CAMGE T0,B ;IF RESULT .LT. B, OK
JRST MRET
CAMG T0,B ;ELSE [A/B] WAS TOO LOW
CAML T1,B+1 ; AND MUST SUBTRACT SOME MORE MULTIPLES OF B
JRST REPEAT ;GO GET T0 MOD B IN T0
MRET: SKIPGE @0(L) ;GIVE RESULT SIGN OF A
DMOVN T0,T0
MRET1: DMOVE T8,SAVE8 ;RESTORE REGISTERS
DMOVE T6,SAVE6
DMOVE T4,SAVE4
DMOVE T2,SAVE2
POPJ P, ;DONE
FAST: CAMN T0,T2 ;HERE IF A HIGH .LE. B HIGH
CAMGE T1,T3 ;HIGH WORDS EQUAL, COMPARE LOW WORDS
JRST MRET ;A .LT. B, GO RETURN DMOD = A
DFSB T0,T2 ;ELSE RETURN DMOD = A - B
JRST MRET
OVFL: JUMPE T6,RTZ ;IF B IS ZERO, DMOD IS 0
FSC T0,-177 ;SCALE A TO AVOID OVERFLOW IN DIVIDE
DMOVE T8,T0 ;GET A
DFDV T8,B ;GET A/B
AND T9,[777000000000] ;TRUNCATE TO 35 BITS SO PRODUCTS ARE EXACT
DFMP T6,T8 ;GET A-B*[A/B]
DFSB T0,T6
DFMP T4,T8
DFSB T0,T4
DFMP T2,T8
DFSB T0,T2
FSC T0,177 ;UNDO SCALING
JRST OVCONT ;CONTINUE
RTZ: SETZB T0,T1 ;A .EQ. B, RETURN DMOD = 0
JRST MRET1
RELOC ;DATA
B: BLOCK 2 ;ABS(B)
SAVE2: BLOCK 2 ;TEMP REGISTERS
SAVE4: BLOCK 2
SAVE6: BLOCK 2
SAVE8: BLOCK 2
RELOC
PRGEND
TITLE DCOS SINE AND COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MAY 1, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DCOS
EXTERN DCOS.
DCOS=DCOS.
PRGEND
TITLE DSIN SINE AND COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MAY 1, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DSIN
EXTERN DSIN.
DSIN=DSIN.
PRGEND
TITLE DSIN. SINE AND COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MAY 1, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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)], THE GREATEST INTEGER IN ABS(X)
; X2 = ABS(X) - X1
; N = THE NEAREST INTEGER TO ABS(X), FOR SIN, OR
; TO ABS(X) + PI/2, FOR COS
; THEN THE REDUCED ARGUMENT F = ((X1 - N*C1) + X2) - N*C2
; WHERE C1 = 3.1104 (OCTAL) AND C2 IS 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) <= 6746518852. ABS(X)+PI/2, IN
; COS(X), MUST BE LESS THAN 6746518852. 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:
; XMOVEI L,ARG
; PUSHJ P,DSIN
; OR
; PUSHJ P,DCOS
;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 FORPRM
TWOSEG 400000
SALL
HELLO (DCOS,.) ;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
DFAD T0,PID2 ;Y = Y+PID2
JRST ALG ;PROCEED TO MAIN ALGORITHM
HELLO (DSIN,.) ;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: 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
DFMP T0,ODPI ;RN = Y/PI
JFCL
DMOVE T4,T0 ;SAVE A COPY OF RN
CAMG T0,TWO26 ;IF RN IS .LE. 2**26
JRST LESS ;THEN GO TO LESS
ASH T4,-33 ;SHIFT OFF MANTISSA
SUBI T4,275 ;OBTAIN SHIFTING FACTOR
ASH T1,(T4) ;SHIFT SECOND WORD APPROPRIATELY
MOVE T5,T1 ;SAVE BITS OF SECOND WORD
ASH T1,-1 ;SHIFT OFF REMAINDER
MOVN T4,T4 ;PREPARE TO SHIFT BACK
ADDI T4,1 ;ADD 1 TO SHIFTING FACTOR
ASH T1,(T4) ;SHIFT BACK SECOND WORD
TRNN T5,1 ;IF RIGHTMOST OF SAVED BITS IS OFF
JRST NORND ;THEN GO TO NORND
DFAD T0,ONE ;OTHERWISE ADD ONE TO XN
ASH T5,-1 ;SHIFT REMAINDER OFF OF N
ADDI T5,1 ;RND N UP
JRST GOTN ;GO TO GOTN
NORND: ASH T5,-1 ;SHIFT OFF FRACTIONAL PART
JRST GOTN ;GO TO GOTN
LESS: FIXR T0,T0 ;FIX RN
MOVE T5,T0 ;SAVE N
FLTR T0,T0 ;XN=FLOAT(N)
MOVEI T1,0 ;ZERO SECOND WORD
GOTN: TRNE T5,1 ;IS N ODD?
MOVNS SGN ;YES,NEGATE SIGN
SKIPE FLAG ;IF THE COSINE IS WANTED
DFAD T0,PT5 ;THEN XN=XN-.5
DMOVE T4,T0 ;SAVE A COPY OF XN
DFMP T0,C1 ;FORM XN*C1
DFSB T2,T0 ;GET |X| - XN * C1
DMOVE T0,T4 ;GET A COPY OF XN
DFMP T4,C2 ;GET XN * C2
DFMP T0,C3 ;GET XN * C3
DFSB T2,T4 ;F = (|X| - XN*C1) - XN*C2
DFSB T2,T0 ; - 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: DFMP T2,T2 ;G = F*F
DMOVE T4,T2 ;SAVE A COPY OF G
DFAD T2,Q2 ;XDEN = G+Q2
DFMP T2,T4 ;*G
DFAD T2,Q1 ;+Q1
DFMP T2,T4 ;*G
DFAD T2,Q0 ;+Q0
DMOVEM T2,XDEN ;MOVE XDEN TO MEMORY
DMOVE T2,T4 ;GET A COPY OF G
DFMP T2,RP5 ;XNUM = G*RP5
DFAD T2,RP4 ;+RP4
DFMP T2,T4 ;*G
DFAD T2,RP3 ;+RP3
DFMP T2,T4 ;*G
DFAD T2,RP2 ;+RP2
DFMP T2,T4 ;*G
DFDV T2,XDEN ;R(G) = XNUM/XDEN
DFAD T2,RP1 ;+RP1
DFMP T2,T4 ;*G
DFMP T2,T0 ;*F
DFAD 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)
PID2: DOUBLE 201622077325,021026430215 ;PI/2
ODPI: DOUBLE 177505746033,162344202512 ;1/PI
EPS: DOUBLE 142400000000,000000000000 ;2**(-31)
C1: DOUBLE 202622077325,020000000000 ;HIGH 34 BITS OF PI
C2: DOUBLE 140413214106,220000000000 ;NEXT 31 BITS OF PI
C3: DOUBLE 101423063050,270033407151 ;C1+C2+C3=PI TO XTRA PREC
Q2: DOUBLE 211612731352,270761316154 ;0.394924723520450141D+3
Q1: DOUBLE 221422322352,120271537275 ;0.702492288221842518D+5
Q0: DOUBLE 227512520255,264021153121 ;0.541748285645351853D+7
RP5: DOUBLE 605161526726,356444271113 ;-0.121560740596710190D-1
RP4: DOUBLE 203422023017,204642007402 ;0.428183075897778265D+01
RP3: DOUBLE 566026406450,010567611157 ;-0.489487151969463797D+03
RP2: DOUBLE 220540546606,025275050734 ;0.451456904704461990D+05
RP1: DOUBLE 601252525252,252525252421 ;-.166666666666666667D0
TWO26: DOUBLE 233400000000,000000000000 ;2**26
ONE: DOUBLE 201400000000,000000000000 ;1.0
PT5: DOUBLE 577400000000,000000000000 ;-.5
YMAX1: 241622077325 ;YMAX = (PI*2**31)
YMAX2: 021026430215
RELOC ;DATA
FLAG: 0
SGN: DOUBLE 0,0
XDEN: DOUBLE 0,0
RELOC
PRGEND
TITLE DSINH HYPERBOLIC SINE FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 7, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DSINH
EXTERN DSINH.
DSINH=DSINH.
PRGEND
TITLE DSINH. HYPERBOLIC SINE FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 7, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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.
;DSINH(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, DSINH = Y*SIGN(X)
; EPS <= Y <= 1, DSINH = X-X*(Z*R(Z)), WHERE Z = Y**2 AND
; R(Z) IS GIVEN BELOW.
; 1 < Y <= 88.029678, DSINH = SIGN(X)*(EXP(Y)-EXP(-Y))/2,
; 88.029678 <= Y < 128 * LN(2)
; DSINH = SIGN(X)*((V/2)*EXP(Y - LN V)),
; Y >= 128 * 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 DSINH IS ABS(X) < 128 * LN(2) AND ERROR MESSAGES
; WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE. DSINH 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:
; XMOVEI L,ARG
; PUSHJ P,DSINH
;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 FORPRM
TWOSEG 400000
SALL
HELLO (DSINH,.) ;ENTRY TO DSINH ROUTINE
DMOVE T0,@(L) ;OBTAIN X
PUSH P,T5 ;SAVE REGISTER T5
HRLZI T5,201400 ;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,TWOM31 ;IF Y IS .LE. 2**-31
JRST RET2 ;GO RETURN
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
DFMP T2,T2 ;Z = X*X
DMOVE T4,T2 ;SAVE A COPY OF Z
DFAD T4,Q2 ;XDEN = Z + Q2
DFMP T4,T2 ;*Z
DFAD T4,Q1 ;+Q1
DFMP T4,T2 ;*Z
DFAD T4,Q0 ;+Q0
DMOVEM T2,Z ;MOVE A COPY OF Z TO MEMORY
DFMP T2,RP3 ;XNUM = Z*RP3
DFAD T2,RP2 ;+RP2
DFMP T2,Z ;*Z
DFAD T2,RP1 ;+RP1
DFMP T2,Z ;*Z
DFAD T2,RP0 ; +RP0
DFDV T2,T4 ;R(Z) = XNUM/XDEN
DFMP T2,Z ;*Z
DFMP T2,TEMP ;*(-X)
DFAD 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 ERRD ;GO TO ERRD
CAMGE T1,YMAXLO ;HI OF Y = YMAXHI
JRST GETW ;IF LO OF Y .LE. YMAXLO, GO TO GETW
ERRD: HRLOI T0,377777 ;SET RESULT TO
SETO T1, ;MACHINE INFINITY
JUMPG T5,DSNH ;IF X IS NEGATIVE
DMOVN T0,T0 ;RESULT = -RESULT
DSNH: LERR (LIB,%,DSINH: result overflow)
POP P,T5 ;RESTORE ACCUMULATOR
GOODBY (1) ;RETURN
GETW: DFAD T0,LNV ;W = Y-LNV
DMOVEM T0,TEMP ;MOVE W TO A TEMP REGISTER
FUNCT DEXP.,<TEMP> ;Z = EXP(W)
DMOVEM T0,TEMP ;SAVE A COPY OF Z
FMP T0,CON1 ;RESULT = Z*CON1
SETZ T1, ;ZERO SECOND WORD
DFAD T0,TEMP ;+Z
JFCL ERRD ;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 DEXP.,<TEMP> ;GET ITS EXP
FSC T0,-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 DEXP.,<TEMP> ;Z = DEXP(Y)
DMOVN T2,T0 ;SAVE A COPY OF Z
MOVEM T5,SGNFLG ;MOVE FLAG TO MEMORY
HRLZI T4,201400 ;SET T4 TO 1.0
SETZ T5, ;CLEAR SECOND WORD
DFDV T4,T2 ;1/Z
DFAD T0,T4 ;+Z
FSC T0,-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 223527442325,224632030652 ;.35181283430177117881D+6
RP1: DOUBLE 216551270255,245012217565 ;.11563521196851768270D+5
RP2: DOUBLE 210507410130,341337112321 ;.16375798202630751372D+3
RP3: DOUBLE 200624234756,033744032643 ;.78966127417357099479D+0
Q0: DOUBLE 551376246137,720314355210 ;-.21108770058106271242D+7
Q1: DOUBLE 220432412710,355457306744 ;.36162723109421836460D+5
Q2: DOUBLE 566352207437,615500015770 ;-.27773523119650701667D+3
ONE: 201400000000 ;1.0
BIGX: 207540074635 ;88.029691
YMAXHI: 207542710277 ;88.722839111672999604
YMAXLO: 276434757153 ;
TWOM31: 142400000000 ;2**-31
LNV: DOUBLE 577235067500,101343000000 ;-.693147180559947174D0
CON1: 120414520000 ;.186417721E-14
TWENT2: 205540000000 ;22
RELOC ;DATA
SGNFLG: 0
TEMP: DOUBLE 0,0
Z: DOUBLE 0,0
RELOC
PRGEND
TITLE DSQRT SQUARE ROOT FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MARCH 19, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DSQRT
EXTERN DSQRT.
DSQRT=DSQRT.
PRGEND
TITLE DSQRT. SQUARE ROOT FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MARCH 19, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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.
;DSQRT(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:
; XMOVEI L,ARG
; PUSHJ P,DSQRT
;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 FORPRM
TWOSEG 400000
SALL
HELLO (DSQRT,.) ;ENTRY TO DSQRT ROUTINE
DMOVE T0,@(L) ;GET DP ARGUMENT
JUMPG T0,DSQRTP ;ARGUMENT IS GREATER THAN 0
JUMPE T0,DSQRT4 ;ARGUMENT IS ZERO
LERR (LIB,%,DSQRT: negative arg; result = SQRT(ABS(arg)))
DMOVN T0,T0 ;ARG = -ARG
DSQRTP: PUSH P,T2 ;SAVE ACCUMULATORS
PUSH P,T3
PUSH P,T4
PUSH P,T5
MOVE T2,T0 ;GET SPARE COPY OF HIGH ORDER BITS
MOVE T3,T0 ;GET SPARE COPY OF HIGH ORDER BITS
MOVE T4,T0 ;GET SPARE COPY OF HIGH ORDER BITS
LSH T3,-1
TLZE T3,400 ;WAS EXPONENT ODD?
JRST DSQRT2 ;YES, GO TO DSQRT2
ADD T3,[XWD 267,607000] ;NO,COMPUTE LINEAR APPROXIMATION
FMPRI T3,301454 ;RESCALE
JRST DSQRT3
DSQRT2: ADD T3,[XWD 267,607000] ;ODD EXPONENT, COMPUTE LINEAR APPROXIMATION
FMPRI T3,301650 ;RESCALE
DSQRT3: FDV T4,T3 ;ORIGINAL ARG / INITIAL GUESS
FAD T4,T3 ;AVERAGE THEM
FSC T4,-1 ;
FDV T2,T4 ;ORIGINAL ARG / SECOND GUESS
FAD T2,T4 ;AVERAGE THEM
FSC T2,-1
MOVEI T3,0 ;LOW HALF OF 1ST APPROX. IS 0
DMOVE T4,T0 ;GET ARG BACK
DFDV T4,T2 ;ORIGINAL ARG / THIRD GUESS
DFAD T2,T4 ;AVERAGE THEM
FSC T2,-1 ;
DFDV T0,T2 ;ORIGINAL ARG / FOURTH GUESS
DFAD T0,T2 ;AVERAGE THEM
FSC T0,-1
POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4
POP P,T3
POP P,T2
DSQRT4: GOODBY (1) ;RETURN
PRGEND
TITLE DCOTAN TANGENT AND COTANGENT FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. APRIL 16, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DCOTAN
EXTERN DCOTN.
DCOTAN=DCOTN.
PRGEND
TITLE DTAN TANGENT AND COTANGENT FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC APRIL 16, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DTAN
EXTERN DTAN.
DTAN=DTAN.
PRGEND
TITLE DTAN. TANGENT AND COTANGENT FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC APRIL 16, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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**(-31)
; = 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, Qi ARE GIVEN BELOW
; THE RESULT IS THEN RECIPROCATED, IF NECESSARY, AND GIVEN
; THE APPROPRIATE SIGN.
; THE APPROXIMATION USED IS A MODIFICATION OF HART 4287.
;THE RANGE OF DEFINITION FOR DTAN IS 0 < ABS(X) < ((2**31) * (PI/2)) = 3373259426.D0
; AND FOR DCOTAN(X), (2**(-127))*(1+(2**(-61))) < ABS(X) < ((2**31)*(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:
; XMOVEI L,ARG
; PUSHJ P,DTAN
; OR
; PUSHJ P,DCOTAN
;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 FORPRM
TWOSEG 400000
SALL
HELLO (DCOTN.,DCOTAN) ;ENTRY TO DCOTAN ROUTINE
PUSH P,T4 ;SAVE ACCUMULATORS
PUSH P,T5
SETZM FLAG ;CLEAR FLAG FOR DCOTAN
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: LERR (LIB,%,DCOTAN: ABS(arg) too large; result = zero)
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 (DTAN,.) ;ENTRY TO TAN ROUTINE
PUSH P,T4 ;SAVE ACCUMULATORS
PUSH P,T5
SETOM FLAG ;SET FLAG FOR DTAN 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: LERR (LIB,%,DTAN or DCOTAN: result underflow)
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 T0
DMOVE T0,T4 ;AND SIGNED F TO T0.
SETZM N ;STORE 0 FOR N.
JRST FPOS ;BYPASS REDUCTION
REDUCE: DFMP T0,TWODPI ;Y*(2/PI)
MOVE T2,T0 ;SAVE A COPY OF Y*(2/PI)
CAMG T2,TWO26 ;IF Y*(2/PI) IS .LE. 2**26
JRST LES ;THEN GO TO LES
ASH T2,-33 ;SHIFT OFF MANTISSA
SUBI T2,275 ;
ASH T1,(T2) ;SHIFT SECOND WORD APPROPRIATELY
MOVE T3,T1 ;SAVE BITS OF SECOND WORD
ASH T1,-1 ;SHIFT REMAINDER OFF
MOVN T2,T2 ;PREPARE TO SHIFT BACK
ADDI T2,1 ;
ASH T1,(T2) ;SHIFT BACK
TRNN T3,1 ;IF RIGHTMOST OF SAVED BITS IS OFF
JRST NORND ;THEN GO TO NORND
DFAD T0,ONE ;OTHERWISE, ADD 1 TO XN
ASH T3,-1 ;PREPARE TO LOOK AT NEXT BIT
ADDI T3,1 ;ADD 1 TO N
JRST RND ;GO TO RND
NORND: ASH T3,-1 ;SHIFT OFF FRACTIONAL PART
RND: MOVEM T3,N ;MOVE N TO MEMORY
JRST GOTIT ;GO TO GOTIT
LES: FIXR T0,T0 ;FIX IT
MOVEM T0,N ;MOVE N TO MEMORY
FLTR T0,T0 ;FLOAT IT
MOVEI T1,0 ;ZERO SECOND WORD
GOTIT: JUMPLE T4,NEXT ;IF X IS POSITIVE
DMOVN T0,T0 ;NEGATE XN
NEXT: DMOVE T2,T0 ;GET A COPY OF -XN
DFMP T2,C1 ;FORM -XN*C1
DFAD T4,T2 ;GET |X| - N*C1
DMOVE T2,T0 ;GET A COPY OF -XN
DFMP T0,C2 ;GET -XN*C2
DFMP T2,C3 ;GET -XN*C3
DFAD T2,T4 ;F = (|X| - XN*C1) - XN*C2
DFAD T2,T0 ; - 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,201400 ;YES , F< EPS,
MOVEI T5,0 ;XDEN = 1.0
JRST FLGCHK ;GO TO CHECK FLAG
NO: DFMP T2,T2 ;NO, F .GE. EPS, SET G=F*F
DMOVE T4,T2 ;SAVE A COPY OF G
DFMP T2,XP4 ;XNUM = XP4 * G +
DFAD T2,XP3 ; P3 *
DFMP T2,T4 ; G +
DFAD T2,XP2 ; P2 *
DFMP T2,T4 ; G +
DFAD T2,XP1 ; P1 *
DFMP T2,T4 ; G *
DFMP T2,T0 ; F +
DFAD T0,T2 ; F
DMOVE T2,T4 ;SAVE A COPY OF G
DFMP T4,Q4 ;XDEN = Q4 * G +
DFAD T4,Q3 ; Q3 *
DFMP T4,T2 ; G +
DFAD T4,Q2 ; Q2 *
DFMP T4,T2 ; G +
DFAD T4,Q1 ; Q1 *
DFMP T4,T2 ; G +
DFAD 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: DFDV T4,T0 ;RESULT = XDEN/XNUM
DMOVE T0,T4
JRST RET ;GO TO RET
NOVD: DFDV T0,T4 ;RESULT = XNUM/XDEN
RET: POP P,T3 ;RESTORE ACCUMULATORS
POP P,T2
POP P,T5
POP P,T4
GOODBY (1)
XP1: DOUBLE 601346652066,115432245623 ;-.1372889460941120802
XP2: DOUBLE 171401224404,126140152543 ; .3925934686364577602E-02
XP3: DOUBLE 616034314474,026356056117 ;-.2882482747560198194E-04
XP4: DOUBLE 147766720604,360442362056 ; .2927308283322907641E-07
Q1: DOUBLE 600036052305,321342375437 ;-.4706222794274454135
Q2: DOUBLE 173702007252,226306364361 ; .2746669449551304872E-01
Q3: DOUBLE 612131325464,137044675153 ;-.4030063705745304384E-03
Q4: DOUBLE 155540343710,045637644377 ; .1312960309685759549E-05
C1: DOUBLE 201622077325,020000000000 ;FIRST 34 BITS OF PI/2
C2: DOUBLE 137413214106,220000000000 ;NEXT 31 BITS OF PI/2
C3: DOUBLE 100423063050,270033407151 ;NEXT 62 BITS OF PI/2.
EPS1: 002400000000 ;HIGH ORDER PART OF EPS1
YMAX1: 240622077325 ;HIGH ORDER PART OF YMAX
YMAX2: 021026430215 ;LOW ORDER PART OF YMAX
HIEPS: 142400000000 ; 2**(-31)
PIO4HI: 200622077325 ;HIGH ORDER PART OF PI/4
TWODPI: DOUBLE 200505746033,162344202512
ONE: DOUBLE 201400000000,000000000000
TWO26: 233400000000 ;2**26
RELOC ;DATA
N: 0
FLAG: 0
RELOC
PRGEND
TITLE DTANH HYPERBOLIC TANGENT FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DTANH
EXTERN DTANH.
DTANH=DTANH.
PRGEND
TITLE DTANH. HYPERBOLIC TANGENT FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 > 22.1807100, TANH = 1.0*SIGN(X)
; IF F > LN(3)/2 AND F <= 22.1807100, TANH = RESULT 1 =
; SIGN(X)*(1 - 2/(EXP(2*F) + 1)))
; IF F < 2**(-31)*SQRT(3), 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:
; XMOVEI 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 FORPRM
TWOSEG 400000
SALL
HELLO (DTANH,.) ;ENTRY TO DTANH 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,LN2T32 ;IF F IS .LE. LN2T32
JRST CALCF ;GO TO CALCF
HRLZI T0,201400 ;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
FSC T0,1 ;F = F+F
DMOVEM T0,TEMP ;SAVE F IN A TEMP REGISTER
FUNCT DEXP.,<TEMP> ;EXP(2*F)
DFAD T0,ONE ;+1.0
HRLZI T2,575400 ;SET T2 TO -2.0
SETZ T3, ;ZERO SECOND WORD
DFDV T2,T0 ;-2/(EXP(2*F)+1)
DMOVE T0,T2
DFAD 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. .806549008E-9
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
DFMP T0,T0 ;G = F*F
DMOVE T2,T0 ;SAVE A COPY OF G
DFAD T2,Q2 ;XDEN = G+Q2
DFMP T2,T0 ;*G
DFAD T2,Q1 ;+Q1
DFMP T2,T0 ;*G
DFAD T2,Q0 ; + Q0
DMOVE T4,T0 ;SAVE A COPY OF G
DFMP T4,RP2 ;XNUM = G*RP2
DFAD T4,RP1 ; + RP1
DFMP T4,T0 ; * G
DFAD T4,RP0 ; + RP0
DFMP T0,T4 ; *G
DFDV T0,T2 ;R(G) = XNUM/XDEN
DFMP T0,TEMP ;RESULT = F*R
DFAD T0,TEMP ; + F
RET: POP P,T4 ;RESTORE ACCUMULATORS
POP P,T3
POP P,T2
POP P,T5
GOODBY (1) ;RETURN
RP0: DOUBLE 564154513215,220360747434 ;-.161341190239962281D+4
RP1: DOUBLE 570163061227,221321463447 ;-.992259296722360833D+2
RP2: DOUBLE 577022172714,101054201376 ;-.964374927772254698D0
Q0: DOUBLE 215456407425,323513222652 ;.484023570719886887D+4
Q1: DOUBLE 214427161323,100370362574 ;.22337720718962312926D+4
Q2: DOUBLE 207702765170,172773772374 ;.112744743805349493D+3
LN2T32: 205542710300 ;22.1807100E0
CON1: 141673317272 ;2**-32 * SQRT(3)
LN3D2: 200431175236 ;.549306144E0
HALF: DOUBLE 200400000000,000000000000 ;1/2
ONE: DOUBLE 201400000000,000000000000 ;1.0
RELOC ;DATA
TEMP: DOUBLE 0,0
RELOC
PRGEND
COMMENT | Not in version 6
TITLE DINT DOUBLE PRECISION TRUNCATION TO INTEGER
SUBTTL CHRIS SMITH/CKS 29-Jan-80
;COPYRIGHT (C) 1980,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DINT
EXTERN DINT.
DINT=DINT.
PRGEND
TITLE DINT. DOUBLE PRECISION TRUNCATION TO INTEGER
SUBTTL CHRIS SMITH/CKS 29-Jan-80
;COPYRIGHT (C) 1980,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
TWOSEG 400000
SALL
HELLO (DINT,.)
PUSH P,T2 ;SAVE TEMP REGISTERS
PUSH P,T3
PUSH P,T4
DMOVE T0,@(L) ;GET ABS(ARG) IN T0-T1
JUMPGE T0,DINT1
DMOVN T0,T0
DINT1: CAMGE T0,[201400000000] ;IS NUMBER .LT. 1?
JRST [SETZB T0,T1 ;YES, DINT IS 0
JRST R1]
CAML T0,[277400000000] ;IS NUMBER ALL INTEGER?
JRST R ;YES, FAST RETURN
LDB T4,[POINT 8,T0,8] ;GET EXPONENT
TRC T4,-1 ;GET ONE'S COMPLEMENT FOR RIGHT SHIFT
MOVSI T2,777000 ;GET INITIAL MASK TO GET RID OF FRACTION
SETZ T3,
ASHC T2,201(T4) ;MOVE MASK OVER INTEGER PART
AND T0,T2 ;CLEAR FRACTION BITS
AND T1,T3
R: SKIPGE @(L) ;ATTACH ORIGINAL SIGN
DMOVN T0,T0
R1: POP P,T4 ;RESTORE TEMP REGISTERS
POP P,T3
POP P,T2
POPJ P, ;RETURN
PRGEND
End of code not in version 6 |
TITLE DABS DOUBLE PRECISION ABSOLUTE VALUE
SUBTTL H. P. WEISS/HPW 11-DEC-73
;COPYRIGHT (C) 1973,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DABS
EXTERN DABS.
DABS=DABS.
PRGEND
TITLE DABS. DOUBLE PRECISION ABSOLUTE VALUE
SUBTTL /KK/DMN/TWE
;COPYRIGHT (C) 1980,1981 BY DIGITAL EQUIPMENT CORPORATION
;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.
;DOUBLE PRECISION ABSOLUTE VALUE FUNCTION
;THIS ROUTINE RETURNS THE ABSOLUTE VALUE OF A DOUBLE PRECISION
;ARGUMENT
SEARCH FORPRM
TWOSEG 400000
SALL
HELLO (DABS,.) ;[235] ENTRY TO DABS ROUTINE
DMOVE T0,@(L) ;GET ARGUMENT
JUMPGE T0,DABSX ;IF POSITIVE, RETURN
DMOVN T0,T0 ;ELSE NEGATE IT
DABSX: GOODBY (1) ;RETURN
PRGEND
TITLE DMAX1 DOUBLE PRECISION MAXIMUM OF A SERIES OF ARGUMENTS
SUBTTL H. P. WEISS/HPW 11-DEC-73
;COPYRIGHT (C) 1973,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DMAX1
EXTERN DMAX1.
DMAX1=DMAX1.
PRGEND
TITLE DMAX1. DOUBLE PRECISION MAXIMUM OF A SERIES OF ARGUMENTS
SUBTTL /TWE/CKS 29-Jan-80
;COPYRIGHT (C) 1980,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
TWOSEG 400000
SALL
HELLO (DMAX1,.) ;[235] ENTRY TO DMAX1 ROUTINE
PUSH P,T2 ;SAVE AC T2
PUSH P,T3 ;Save AC t3
HLLZ T3,-1(L) ;Get the F10 arg count
DMOVE T0,@(L) ;GET FIRST ARGUMENT
JRST DMAX.2 ;ADDRESS OF NEXT, START CHECKING
DMAX.1: XMOVEI T2,@0(L) ;GET ADDRESS OF NEXT ARGUMENT
CAMLE T0,(T2) ;IS HIGH ORDER .GT. THIS ONE?
JRST DMAX.2 ;YES, GET ADDRESS OF NEXT IN LIST
CAME T0,(T2) ;ARE HIGH ORDER WORDS EQUAL?
JRST DMAX.3 ;NO, REPLACEMENT NECESSARY
CAMGE T1,1(T2) ;SKIP IF PRESENT ARG IS LARGER
DMAX.3: DMOVE T0,(T2) ;PICK UP NEW ARG
DMAX.2: ADDI L,1 ;Point to next arg
AOBJN T3,DMAX.1 ;Loop for all args.
POP P,T3 ;End of arg list
POP P,T2 ;, restore acs
GOODBY (2) ;Return answer in 0 and 1
PRGEND
TITLE DMIN1 DOUBLE PRECISION MINIMUM OF A SERIES OF ARGUMENTS
SUBTTL H. P. WEISS/HPW 11-DEC-73
;COPYRIGHT (C) 1973,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DMIN1
EXTERN DMIN1.
DMIN1=DMIN1.
PRGEND
TITLE DMIN1. DOUBLE PRECISION MINIMUM OF A SERIES OF ARGUMENTS
SUBTTL /TWE/CKS 29-Jan-80
;COPYRIGHT (C) 1980,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
TWOSEG 400000
SALL
HELLO (DMIN1,.) ;[235] ENTRY TO DMIN1 ROUTINE
PUSH P,T2 ;SAVE AC T2
PUSH P,T3 ;And T3
HLLZ T3,-1(L) ;Get the F10 arg count
DMOVE T0,@(L) ;GET FIRST ARGUMENT
JRST DMIN.2 ;ADDRESS OF NEXT, START CHECKING
DMIN.1: XMOVEI T2,@0(L) ;GET ADDRESS OF NEXT ARGUMENT
CAMGE T0,(T2) ;IS HIGH ORDER .LT. THIS ONE?
JRST DMIN.2 ;YES, GET ADDRESS OF NEXT IN LIST
CAME T0,(T2) ;ARE HIGH ORDER WORDS EQUAL?
JRST DMIN.3 ;NO, REPLACEMENT NECESSARY
CAMLE T1,1(T2) ;SKIP IF PRESENT ARG IS LARGER
DMIN.3: DMOVE T0,(T2) ;PICK UP NEW ARG
DMIN.2: ADDI L,1 ;Point to next arg
AOBJN T3,DMIN.1 ;End of arg list?
POP P,T3 ;Yes, restore T3
POP P,T2 ; and T2
GOODBY (2) ;Return answer in 0 and 1
PRGEND
TITLE DSIGN DOUBLE PRECISION SIGN ROUTINE
SUBTTL H. P. WEISS/HPW 11-DEC-73
;COPYRIGHT (C) 1973,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DSIGN
EXTERN DSIGN.
DSIGN=DSIGN.
PRGEND
TITLE DSIGN. DOUBLE PRECISION SIGN ROUTINE
SUBTTL /KK/DMN/CKS 29-Jan-80
;COPYRIGHT (C) 1980,1981 BY DIGITAL EQUIPMENT CORPORATION
;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.
;FROM LIB40 V.006.
;DOUBLE PRECISION TRANSFER OF SIGN
;THIS ROUTINE RETURNS ABSF(ARG1)*SIGN(ARG2)
;THE ROUTINE MAKES USE OF THE FOLLOWING TABLE:
;ARG1 ARG2 RESULT CHANGE OF SIGN?
;+ + + NO
;+ - - YES
;- + + YES
;- - - NO
SEARCH FORPRM
TWOSEG 400000
SALL
HELLO (DSIGN,.) ;[235] ENTRY TO DSIGN ROUTINE
DMOVE T0,@(L) ;GET FIRST ARGUMENT
SKIPGE @1(L) ;THEN
JUMPL T0,OUT ;CHOOSE
SKIPL @1(L) ;THE
JUMPGE T0,OUT ;CORRECT
DMOVN T0,T0 ;SIGN.
OUT: GOODBY (2) ;EXIT
PRGEND
TITLE DFLOAT INTEGER TO DOUBLE PRECISION CONVERSION
SUBTTL H. P. WEISS/HPW 11-DEC-73
;COPYRIGHT (C) 1973,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DFLOAT
EXTERN DFLOT.
DFLOAT=DFLOT.
PRGEND
TITLE DFLOT. INTEGER TO DOUBLE PRECISION CONVERSION
SUBTTL ED YOURDON/KK/VAA/TWE/DMN
;COPYRIGHT (C) 1973,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
EXTERN DFL.0
NOSYM
TWOSEG 400000
SALL
HELLO (DFLOT.,DFLOAT) ;[235] ENTRY TO DFLOAT ROUTINE
MOVE T0,@(L) ;PICK UP THE ARGUMENT
PJRST DFL.0 ;USE DFL.0 ROUTINE
PRGEND
TITLE IDINT DOUBLE PRECISION TO INTEGER CONVERSION
SUBTTL H. P. WEISS/HPW 11-DEC-73
;COPYRIGHT (C) 1973,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY IDINT
EXTERN IDINT.
IDINT=IDINT.
PRGEND
TITLE IDINT. DOUBLE PRECISION TO INTEGER CONVERSION
SUBTTL TOM EGGERS/DMN
;COPYRIGHT (C) 1973,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
EXTERN IDF.0
NOSYM
TWOSEG 400000
SALL
HELLO (IDINT,.) ;[235] ENTRY TO IDINT
DFIX1: DMOVE T0,@(L) ;GET THE ARGUMENT
PJRST IDF.0 ;USE IDF.0 ROUTINE
PRGEND
TITLE DBLE REAL TO DOUBLE PRECISION CONVERSION
SUBTTL H. P. WEISS/HPW 11-DEC-73
;COPYRIGHT (C) 1973,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY DBLE
EXTERN DBLE.
DBLE=DBLE.
PRGEND
TITLE DBLE. REAL TO DOUBLE PRECISION CONVERSION
SUBTTL ED YOURDON/KK
;COPYRIGHT (C) 1973,1981 BY DIGITAL EQUIPMENT CORPORATION
;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.
;FROM LIB40 V.005.
;SINGLE PRECISION TO DOUBLE PRECISION CONVERSION
;THIS ROUTINE CONVERTS A SINGLE PRECISION ARGUMENT TO A
;DOUBLE PRECISION ANSWER WITH THE LOW ORDER WORD SET TO 0
SEARCH FORPRM
NOSYM
TWOSEG 400000
SALL
HELLO (DBLE,.) ;[235] ENTRY TO DOUBLE ROUTINE
MOVE T0,@(L) ;GET THE ARGUMENT IN HIGH ORDER
SETZ T1, ;SET LOW ORDER PORTION TO ZERO
GOODBY (1) ;EXIT
PRGEND
TITLE SNGL DOUBLE PRECISION TO REAL CONVERSION
SUBTTL H. P. WEISS/HPW 11-DEC-73
;COPYRIGHT (C) 1973,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM
NOSYM
ENTRY SNGL
EXTERN SNGL.
SNGL=SNGL.
PRGEND
TITLE SNGL. DOUBLE PRECISION TO REAL CONVERSION
SUBTTL /DMN/TWE
;COPYRIGHT (C) 1973,1981 BY DIGITAL EQUIPMENT CORPORATION
;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.
;DOUBLE PRECISION TO SINGLE PRECISION CONVERSION FUNCTION
;THIS ROUTINE OBTAINS THE MOST SIGNIFICANT PART OF A DOUBLE
;PRECISION ARGUMENT.
SEARCH FORPRM
EXTERN SNG.0
NOSYM
TWOSEG 400000
SALL
HELLO (SNGL,.) ;[235] ENTRY TO SNGL ROUTINE
DMOVE T0,@(L) ;GET THE ARGUMENT IN AC 0
PJRST SNG.0 ;USE AC0 SINGLE ROUTINE
PRGEND
UNIVERSAL FORDAR GENERAL DOUBLE PRECISION ARITHMETICS
SUBTTL TOM EGGERS/DMN/DRT/JNG/SWG 14-Aug-79
;COPYRIGHT (C) 1972,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 MACROS FOR FORDAR UNIVERSAL FILE
IF1,< ;ONLY ONCE
; DOUBLE PRECISION FLOAT FUNCTION "DFLOAT"
DEFINE DFL (X) <
XALL
ENTRY DFL.'X ;ENTRY POINT TO DFL.'X
SIXBIT /DFL.'X/
DFL.'X: MOVEI X+1,0 ;CLEAR LOW ORDER WORD
ASHC X,-8 ;MAKE ROOM FOR EXPONENT IN HI WORD
TLC X,243000 ;SET EXP TO 27+8 DECIMAL
DFAD X,[EXP 0,0] ;NORMALIZE
POPJ P, ;RETURN X=THE DOUBLE PRECISION RESULT
>; END DFL
; DOUBLE PRECISION FIX FUNCTION "IDINT"
; DOUBLE TO INTEGER
DEFINE IDF (X) <
XALL
ENTRY IDF.'X
SIXBIT /IDF.'X/
IDF.'X: PUSH P,L ;SAVE THE SCRATCH REG
HLRE L,X ;GET THE EXPONENT
ASH L,-9 ;RIGHT 8 BITS
JUMPGE X,IDF.XT ;JUMP IF POS.
DMOVN X,X ;NEGATE
TRC L,-1 ;COMPLEMENT THE EXPONENT
IDF.XT: TLZ X,777000 ;CLEAR THE EXPONENT
ASHC X,-201-^D26(L) ;CHANGE FRACTION TO INTEGER
TLNE L,400000 ;SKIP IF POS.
MOVN X,X ;NEGATE
POP P,L ;RESTORE THE SCRATCH REG
POPJ P, ;RETURN X=FIXED NUMBER
>; END IDF
; DOUBLE PRECISION TO SINGLE FUNCTION
DEFINE SNG (X)<
XALL
ENTRY SNG.'X
SIXBIT /SNG.'X/
SNG.'X: JUMPL X,SNG3 ;NEGATIVE ARGUMENT?
TLNE X+1,(1B1) ;POSITIVE. ROUND REQUIRED?
TRON X,1 ;YES, TRY TO ROUND BY SETTING LSB
POPJ P, ;WE WON, FINISHED
MOVE X+1,X ;COPY HIGH PART OF ARG
AND X,[777000,,1] ;MAKE UNNORMALIZED LSB, SAME EXPONENT
FAD X,X+1 ;ROUND & RENORMALIZE
POPJ P,
;HERE IF ARG IS NEGATIVE
SNG3: DMOVN X,X ;MAKE POSITIVE
TLNE X+1,(1B1) ;NEED ROUNDING?
TRON X,1 ;YES, TRY TO DO IT BY SETTING LSB
JRST SNG4 ;DONE
MOVN X+1,X ;MAKE RE-NEGATED COPY OF HIGH PART
ORCA X,[777,,-1] ;GET UNNORM NEG LSB WITH SAME EXPONENT
FADR X,X+1 ;ROUND & NORMALIZE
POPJ P,
SNG4: MOVN X,X ;RE-NEGATE
POPJ P, ;EXIT
>; END SNG
>; END IF1
PRGEND
TITLE DFL.0
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM,FORDAR
NOSYM
TWOSEG 400000
DFL 0
PRGEND
TITLE DFL.2
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
DFL 2
PRGEND
TITLE DFL.3
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
DFL 3
PRGEND
TITLE DFL.4
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
DFL 4
PRGEND
TITLE DFL.6
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
DFL 6
PRGEND
TITLE DFL.10
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
DFL 10
PRGEND
TITLE DFL.12
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
DFL 12
PRGEND
TITLE DFL.14
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
DFL 14
PRGEND
TITLE IDF.0
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM,FORDAR
NOSYM
TWOSEG 400000
IDF 0
PRGEND
TITLE IDF.2
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
IDF 2
PRGEND
TITLE IDF.4
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
IDF 4
PRGEND
TITLE IDF.6
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
IDF 6
PRGEND
TITLE IDF.10
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
IDF 10
PRGEND
TITLE IDF.12
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
IDF 12
PRGEND
TITLE IDF.14
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
IDF 14
PRGEND
TITLE SNG.0
;COPYRIGHT (C) 1979,1981 BY DIGITAL EQUIPMENT CORPORATION
;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 FORPRM,FORDAR
NOSYM
TWOSEG 400000
SNG 0
PRGEND
TITLE SNG.2
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
SNG 2
PRGEND
TITLE SNG.4
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
SNG 4
PRGEND
TITLE SNG.6
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
SNG 6
PRGEND
TITLE SNG.10
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
SNG 10
PRGEND
TITLE SNG.12
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
SNG 12
PRGEND
TITLE SNG.14
SEARCH FORPRM,FORDAR
NOSYM
TWOSEG 400000
SNG 14
PRGEND
TITLE DDIM Double precision 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) 1981 BY DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.
SEARCH FORPRM
TWOSEG 400000
NOSYM
ENTRY DDIM
EXTERN DDIM.
DDIM=DDIM.
PRGEND
TITLE DDIM. Double Precision Positive Difference.
SUBTTL C. McCutcheon 6/26/81
;COPYRIGHT (C) 1981 BY DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.
;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 double precision.
; Passed: A1,A2
; Returns:
; A1-A2 if A1 .GT. A2
; 0 if A1 .LE. A2
SEARCH FORPRM
TWOSEG 400000
SALL
HELLO (DDIM,.)
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.
DFSB T0,T2 ;A1-A2
RET: POP P,T3
POP P,T2
GOODBYE ;Return
RET0: SETZB T0,T1 ;Return zero
JRST RET
PRGEND
TITLE DINT Double precision truncation
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
; OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.
;
;COPYRIGHT (C) 1981 BY DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.
SEARCH FORPRM
TWOSEG 400000
NOSYM
ENTRY DINT
EXTERN DINT.
DINT=DINT.
PRGEND
TITLE DINT. Double Precision Truncation.
SUBTTL C. McCutcheon - 6/25/81
;COPYRIGHT (C) 1981 BY DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.
;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
; 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 FORPRM
TWOSEG 400000
SALL
HELLO (DINT,.) ;Enter routine
DMOVE T0,@0(L) ;Get argument to truncate.
JUMPN T0,NZERO ;If zero
GOODBYE ;Return.
NZERO: PUSH P,T2 ;Save ac's
SKIPGE @0(L) ;If original number .LT. zero,
DMOVN T0,T0 ; negate it.
HLRZ T2,T0 ;Get exponent
LSH T2,-^D9 ;Put rightmost
CAIG T2,200 ;If number .LT. 1.0D0,
JRST ZERO ; return zero.
CAIL T2,276 ;If exponent .GE. 276 then there is
JRST DONE ; no fraction, return the number passed.
; Now shift out the fractional part of the number.
ASHC T0,-276(T2) ;Shift into integer position.
MOVN T2,T2 ;Negate exponent.
ASHC T0,276(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
; No fractional part, return number passed to function
DONE: DMOVE T0,@0(L)
JRST BYE
ZERO: SETZB T0,T1 ;Return a 0.
JRST BYE
PRGEND
TITLE DNINT Double precision nearest whole number
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
; OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.
;
;COPYRIGHT (C) 1981 BY DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.
SEARCH FORPRM
TWOSEG 400000
NOSYM
ENTRY DNINT
EXTERN DNINT.
DNINT=DNINT.
PRGEND
TITLE DNINT. Double Precision nearest whole number.
SUBTTL C. McCutcheon - 6/25/81
;COPYRIGHT (C) 1981 BY DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.
;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 double precision
; nearest whole number to the double precision number passed.
; Number returned is: Integer(A+.5) if A .GE. 0
; Integer(A-.5) if A .LT. 0.
SEARCH FORPRM
TWOSEG 400000
HELLO (DNINT,.) ;Enter routine
DMOVE T0,@0(L) ;Get argument to truncate.
JUMPN T0,NZERO ;Is number passed = 0?
GOODBYE ;Yes.
NZERO: PUSH P,T2 ;Save an ac
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.
HLRZ T2,T0 ;Get exponent
LSH T2,-^D9 ;Put rightmost
CAIL T2,276 ;If no fraction possible,
JRST DONE ; return the number passed.
; Now shift out the fractional part of the number.
DFAD T0,[200400,,0
0,,0] ;Add .5D0 before truncation.
HLRZ T2,T0 ;Get new exponent.
LSH T2,-^D9 ;Put rightmost.
ASHC T0,-276(T2) ;Shift into integer position.
MOVN T2,T2 ;Negate exponent.
ASHC T0,276(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.
GOODBYE ;Return.
; No fractional part found, return number passed to function.
DONE: DMOVE T0,@0(L)
JRST BYE
ZERO: SETZB T0,T1 ;Return 0.
JRST BYE
PRGEND
TITLE DPROD Double precision product for single precision
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
; OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.
;
;COPYRIGHT (C) 1981 BY DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.
SEARCH FORPRM
TWOSEG 400000
NOSYM
ENTRY DPROD
EXTERN DPROD.
DPROD=DPROD.
PRGEND
TITLE DPROD. Double Precision product for single precision.
SUBTTL C. McCutcheon - 6/25/81
;COPYRIGHT (C) 1981 BY DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.
;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 FORPRM
TWOSEG 400000
SALL
HELLO (DPROD,.) ;Enter routine
MOVE T1,@1(L) ;Get 2nd argument
MOVEM T1,MULT ;Save 2nd argument
MOVE T0,@0(L) ;Get 1st argument
SETZB T1,MULT+1 ;Set both second words to 0
DFMP T0,MULT ;Multiply and leave result in AC0.
GOODBYE ;Return
RELOC
MULT: BLOCK 2 ;Multiplier
RELOC
PRGEND
TITLE IDNINT Integer nearest whole number for double precision
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
; OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.
;
;COPYRIGHT (C) 1981 BY DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.
SEARCH FORPRM
TWOSEG 400000
NOSYM
ENTRY IDNINT
EXTERN IDNIN.
IDNINT=IDNIN.
PRGEND
TITLE IDNIN. Integer nearest whole number for double prec.
SUBTTL C. McCutcheon - 6/25/81
;COPYRIGHT (C) 1981 BY DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.
;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 to return the nearest
; integer to the double precision number passed.
; Number returned is: Integer(A+.5) if A .GE. 0
; Integer(A-.5) if A .LT. 0.
SEARCH FORPRM
TWOSEG 400000
SALL
HELLO (IDNIN.,IDNINT) ;Enter routine
DMOVE T0,@0(L) ;Get argument to round.
JUMPN T0,NZERO ;If number passed is 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
DFAD T0,[200400,,0 ;Add 0.5D0 before truncation
0,,0] ; to round the number
; Now shift out the fractional part of the number.
HLRZ T2,T0 ;Get exponent.
LSH T2,-^D9 ;Put rightmost.
CAILE T2,243 ;If not representable, (overflow imminent)
JRST DONE ; return largest integer
TLZ T0,777000 ;Eliminate exponent.
ASHC T0,-233(T2) ;Shift into integer position, leaving
; 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
GOODBYE ;Return
; Number is .LT. 1/2, return zero quickly (don't use the normal flow
; because DFAD 0.5D0 will round)
ZERO: SETZ T0, ;Return . . .
GOODBYE ; . . . 0
; Number too large to represent as integer. Return largest
; integer.
DONE: LERR (LIB,%,<IDNINT: result overflow>)
HRLOI T0,377777 ;Largest positive number.
SKIPGE @0(L) ;If original number .LT. zero,
HRLZI T0,400000 ; largest negative number.
JRST BYE
END