Trailing-Edge
-
PDP-10 Archives
-
fortv11
-
mthdbl.mac
There are 9 other files named mthdbl.mac in the archive. Click here to see a list.
SEARCH MTHPRM
TV MTHDBL DOUBLE PRECISION ROUTINES ,2(4020)
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SUBTTL REVISION HISTORY
COMMENT \
***** Begin Revision History *****
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.
***** Begin Mathlib *****
3055 CKS/JLC 17-Mar-82
Fix DEXP2 overflow/underflow message output.
3200 JLC 10-May-82
Mathlib integration. Change all LERRs to $LCALLs. Change
TWOSEG/RELOC to SEGMENT macros.
3216 JLC 3-Sep-82
Remove extraneous MTHDAR universal file and SEARCHes of it,
since it was moved to MTHPRM long ago.
3221 BCM 29-Oct-82
Moved SNGL, SNGL., and SNG.nn routines to the end of file
MTHSNG and out of MTHDBL. Resolves forward reference problem.
***** Begin Version 1A *****
3232 CKS 23-Mar-83
Rewrote DMOD to remove restrictions on the magnitudes of
its arguments.
3243 BCM 29-Mar-83
Get and clear PC flags with new GETFLG and RESFLG macros. Fixes
always clearing underflow in MTHTRP. Only applies to JOV branches.
3244 TGS 07-Mar-83
Retract edit 3232 to MTHDBL. DMOD was returning inaccurate results
in the 2d word.
3245 MRB 11-Apr-83
DCOTAN pops to many AC's off the stack if an result to large occurs.
3254 20-19812 MRB 13-Dec-83
The DEXP2. routine is failing to clear the result field on an
underflow. (This is edit 3374 to FOROT7.)
***** End Revision History *****
***** Begin Version 2 *****
4002 JLC 3-Aug-83
Move DEXP. to after DTAN., which calls DEXP.
Make 1.0 a special case for DEXP3.
4011 JLC 4-May-84
Add code for MOD, AMOD, DMOD, and GMOD for 2nd arg=0.
4020 DCE 17-Dec-85
The GFLOATING exponentiation routine uses an incorrect method of
getting the exponent of the magnitude of certain intermediate
GFLOATING numbers. A MOVM on the first word of the double word
value is used followed by stripping out the exponent field.
This does not work for all negative GFLOATING numbers where the
mantissa part in the first word is all zero (the first
non-zero bit of the mantissa is in the second word). Due to the
fact that these errors occur in an intermediate calculation,
it is impossible to specify the class of input numbers which
will result in this problem.
\
PRGEND
TITLE DACOS ARC SINE AND ARC COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 19, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY 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) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY 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) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DASIN(X) AND DACOS(X) ARE CALCULATED AS FOLLOWS
; LET R(G) = (G*(RP1+G*(RP2+G*(RP3+G*(RP4+G*RP5)))))/(Q0+G*(Q1
; +G*(Q2+G*(Q3+G*(Q4+G)))))
; (RP1,RP2,RP3,RP4,RP5,Q0,Q1,Q2,Q3,Q4, AND Q5 ARE GIVEN BELOW)
; LET S = Y + Y*R(G)
; FOR SITUATIONS 1. AND 3. BELOW,
; G = Y**2, WHERE Y = ABS(X)
; FOR SITUATIONS 2. AND 4., G = (1.0 - ABS(X))/2.0
; AND Y = -2.0*SQRT(G)
; LET W = ABS(X) AND TERM THE RESULT T.
; THEN, CONSIDER THE FOUR SITUATIONS:
; 1. DASIN, FOR W <= 0.5. T = S, AND T IS NEGATED
; FOR NEGATIVE X
; 2. DASIN, FOR W > 0.5. T = PI/2 + S, AND T IS NEGATED
; FOR NEGATIVE X
; 3. DACOS, FOR W <= 0.5. T = PI/2 - S IF X > 0
; = PI/2 + S IF X < 0
; 4. DACOS FOR W > 0.5. T = -S IF X > 0
; = PI + S IF X < 0.
;THE RANGE OF DEFINITION FOR DASIN/DACOS IS ABS(X) <= 1.0. AND ERROR MESSAGES
; WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE. DASIN/DACOS WILL BE SET
; TO + MACHINE INFINITY.
;REQUIRED (CALLED) ROUTINES: DSQRT
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; 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 MTHPRM
SEGMENT CODE
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: $LCALL AOI
;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
SEGMENT DATA
TEMP: DOUBLE 0,0
I: 0
FLAG: 0
PRGEND
TITLE DATAN2 ARC TAN FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MAY 9, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DATAN2
EXTERN DATN2.
DATAN2=DATN2.
PRGEND
TITLE DATAN ARC TAN FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL Chris Smith/CKS
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DATAN
EXTERN DATAN.
DATAN=DATAN.
PRGEND
TITLE DATAN. ARC TAN FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DATAN(X) is computed as follows:
;
;If X < 0, compute DATAN(|X|) below, then DATAN(X) = -DATAN(-X).
;
;If X > 0, use the identity
;
; DATAN(X) = DATAN(XHI) + DATAN(Z)
; Z = (X - XHI) / (1 + X*XHI)
;
;where XHI is chosen so that |Z| <= tan(pi/32).
;
;XHI is chosen to be exactly representable as a double precision number,
;and so Z can be calculated without loss of significance.
;
;DATAN(XHI) is found by table lookup. It is stored as ATANHI + ATANLO to
;provide guard bits for the final addition to DATAN(Z).
;
;DATAN(Z) is evaluated by means of a polynomial approximation from Hart et al.
;(formula 4904).
;
;If X < tan(pi/32), DATAN(X) = DATAN(Z).
;If X > 1/tan(pi/32), DATAN(X) = pi/2 - DATAN(1/X).
;
;If tan(pi/32) < X < 1/tan(pi/32), an appropriate XHI is obtained by indexing
;into a table. The table tells which XHI to use for various ranges of X.
;The index into the table is formed from the low 3 exponent bits and the high
;3 fraction bits of X.
SEARCH MTHPRM
SEGMENT CODE
SALL
T6=6 ;MORE AC NAMES
T7=7
T8=10
HELLO (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
$LCALL RUN,RET
; LERR (LIB,%,<DATAN2: result underflow>,,RET)
;IF SECOND ARG (X) POSITIVE, RESULT UNDERFLOWS
DMOVN T0,MPI ;ELSE RESULT IS PI WITH SIGN OF FIRST ARG
JRST YSIGN ;GO ATTACH SIGN AND RETURN
DIVCHK: SKIPE @0(L) ;CHECK FOR BOTH ARGS 0
JRST OVER ;ATAN2(NONZERO,0) IS SAME AS OVERFLOW
$LCALL BAZ
;LERR (LIB,%,<DATAN2: both arguments are zero, result=0.0>)
SETZB T0,T1 ;RETURN 0
JRST RET
OVER: DMOVE T0,PI2 ;OVERFLOW, RESULT IS PI/2 WITH SIGN OF FIRST ARG
YSIGN: SKIPGE @0(L) ;ATTACH SIGN OF FIRST ARGUMENT (Y)
DMOVN T0,T0
JRST RET ;RETURN
SUBTTL TABLES
;This table is indexed by the low 3 exponent bits and the high 3 fraction bits
;of X, where MINX < X < MAXX. It gives the offset into XHI, ATANHI, and ATANLO
;of a suitable XHI.
DEFINE OFFS (X) <
IRP X,<XWD 2*X,X>
>
RADIX 10
OFFSET: OFFS <1,1,1,2,2,3,3,4,4,4,5,5,5,6,6,7,7,7,8,8,8,9,9,9>
OFFS <10,10,10,10,11,11,11,12,12,12,12,12,13,13,13,13>
OFFS <13,13,14,14,14,14,14,14,14,14,14,14,14,14,14>
RADIX 8
PI2OFFS=2*^D15 ;OFFSET INTO ATANHI AND ATANLO OF PI/2
MPIOFFS=2*^D16 ;OFFSET OF -PI
MPI2OFFS=2*^D31 ;OFFSET OF -PI/2
XHI: EXP 000000000000,0 ; .0000000000000000000 not used
EXP 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)
SEGMENT DATA
SGNFLG: BLOCK 1 ;SIGN TO BE ATTACHED TO RESULT
USAVE: BLOCK 2 ;TEMP STORAGE FOR U
PRGEND
TITLE DCOSH HYPERBOLIC COSINE FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 7, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DCOSH
EXTERN DCOSH.
DCOSH=DCOSH.
PRGEND
TITLE DCOSH. HYPERBOLIC COSINE FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 7, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DCOSH(X) IS CALCULATED AS FOLLOWS
; LET V BE APPROXIMATELY 2 SO THAT LNV AND ABS(X)+LN V
; CAN BE EXACTLY REPRESENTED WHEN X IS EXACTLY REPRESENTABLE.
; THEN, LETTING Y = ABS(X), FOR
; Y <= 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 MTHPRM
SEGMENT CODE
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
$LCALL ROV
;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)
SEGMENT DATA
TEMP: DOUBLE 0,0
PRGEND
TITLE DEXP2. DOUBLE ** INTEGER EXPONENTIATION
SUBTTL MARY PAYNE/MHP/CKS
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1983, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
SEGMENT CODE
SALL
; Mary Payne, July 30, 1982
HELLO (DEXP2,.)
DMOVEM T2,SAVE2 ;Save T2-T4
MOVEM T4,SAVE4
DMOVE T0,ONE ;Floating 1 to T0-T1
MOVM T2,@1(L) ;|exponent| to T2
JUMPE T2,EXP0 ;Exponent = 0 is special
DMOVE T3,@0(L) ;Base to T3-T4.
JUMPN T3,STEP1 ;If base not 0 go to main flow
JRST BASE0 ;Else to special code
LOOP: DFMP T3,T3 ;Square current result
JOV OVER2 ;Over/underflow possible
STEP1: TRNE T2,1 ;If exponent is odd
DFMP T0,T3 ; update current result
JOV OVER ; Branch on over/underflow
LSH T2,-1 ;Discard low bit of exponent
JUMPN T2,LOOP ;Iterate if not 0
SKIPL @1(L) ;If exponent > 0
JRST RET ; return
DMOVE T2,T0 ;Copy result
DMOVE T0,ONE ;Get reciprocal of
DFDV T0,T2 ;result; Underflow impossible
JOV OVMSG ; On overflow get message
JRST RET ;Else return
EXP0: SKIPE @0(L) ;Exponent 0. If base not
JRST RET ; 0, result is 1. Return
$LCALL ZZZ ;zero**zero is indeterminate, result=zero
SETZB T0,T1 ;Store 0
JRST RET
BASE0: SKIPL @1(L) ;If exponent > 0
JRST ZERO ; result is 0
$LCALL ROV ;Else overflow
HRLOI T0,377777 ;Store +biggest
HRLOI T1,377777
JRST RET
ZERO: SETZB T0,T1 ;Result is 0
JRST RET ;Return
;
;The following block of code deals with over/underflow in the
;square operation at LOOP:. Note that the "exponent" cannot be
;0 -- LOOP: is entered only if T2 is not 0. Moreover, if T2 is
;not 1 subsequent operations will aggravate the over/underflow
;condition in such a way that both the result of the iteration
;and its reciprocal will have the same exception as currently
;indicated. If, however, T2 = 1, and the square overflowed, it
;is possible that its reciprocal will be in range. We therefore
;complete the current pass through the loop, and if the LSH of T2
;makes it zero, we join the handling at OVER: for overflow/underflow
;on the MUL of T0 by T3. Note that no exception can occur on the
;MUL of T0 by a wrapped over/underflow of the square, so that the
;exception flags will still be valid after this step.
;
OVER2: DFMP T0,T3 ;No over/underflow. Hence flags
; from square of T3 still valid
LSH T2,-1 ;Discard low bit of exponent
JUMPE T2,OVER ;If T2 = 0, T0 has wrapped final
; result or its reciprocal
; which may be in range
;Final product surely over/underflows.
GETFLG T2 ;[3243] get exception flags into t2
TLNE T2,(PC%FUF) ;[3243] If underflow flag set, reciprocal
JRST UNDER ; overflows. Go test sign of exponent
SKIPL @1(L) ;For overflow, if exponent > 0
JRST UNDMSG ; final result underflows.
JRST OVMSGF ;[3243] Else reciprocal gives overflow
;
;The rest of the code handles over/underflow on the product of
;T0 by T3 and calculation of the reciprocal, if this is done.
;
OVER: GETFLG T2 ;[3243] get exception flags into t2
TLNE T2,(PC%FUF) ;[3243] If underflow flag set
JRST UNDER ; underflow on product
SKIPL @1(L) ;Else, overflow on result if
JRST OVMSGF ; exponent > 0. Get message
DMOVE T3,T0 ;[3243] Copy result
DMOVE T0,ONE ;For exponent < 0, get
DFDV T0,T3 ;[3243] reciprocal of wrapped overflow
JOV RETF ;[3243] Underflow impossible; overflow
; compensates previous overflow
JRST UNDMSG ;Else, get underflow message
UNDER: SKIPL @1(L) ;Product underflowed. If exponent
JRST UNDMSG ; >/= 0, result underflows
;Else reciprocal overflows
;[3243] entry point necessary because other branches join OVMSG
OVMSGF: RESFLG T2 ;[3243] restore the flags, t2-t3 are dbl wd PC
OVMSG: $LCALL ROV ;Result overflow
JUMPL T0,NEGOV ;If result > 0
HRLOI T0,377777 ;Store +BIGGEST
HRLOI T1,377777
JRST RET ; and return
NEGOV: MOVSI T0,400000 ;If result < 0, store -BIGGEST
MOVEI T1,1
JRST RET ; and return
UNDMSG: $LCALL RUN ;Result underflow
SETZB T0,T1 ;[3254]
RETF: RESFLG T2 ;[3243] reset all flags, t2-t3 are dbl wd PC
RET: DMOVE T2,SAVE2 ;Restore T2-T4
MOVE T4,SAVE4
POPJ P,
ONE: EXP 1.0,0 ;D-floating 1.0
SEGMENT DATA
SAVE2: BLOCK 2 ;Temp for T2-T3
SAVE4: BLOCK 1 ;Temp for T4
PRGEND
TITLE DEXP3. POWER FUNCTION
; (DOUBLE PRECISION)
SUBTTL IMSL, INC. JUNE 4, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;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 MTHPRM
SEGMENT CODE
SALL
HELLO (DEXP3,.) ;ENTRY TO DEXP3. ROUTINE
PUSH P,T2 ;SAVE ACCUMULATORS
PUSH P,T3
DMOVE T0,@(L) ;GET THE BASE
CAMN T0,[1.0] ;IS IT EXACTLY 1.0?
JUMPE T1,POPRET ;YES. JUST RETURN EXACTLY 1.0
DMOVE T2,@1(L) ;GET THE EXPONENT
SETZM IFLAG ;SET INTEGER EXP FLAG TO 0
JUMPG T0,XOK ;IF X IS NOT GT 0
JUMPE T0,X0 ;AND IF X IS NOT = 0
PUSH P,T4 ;SAVE ACCUMULATORS
PUSH P,T5
PUSH P,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: $LCALL NNA
;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
JUMPE T2,ZERZER ;SEPARATE OUT 0**0 CASE
$LCALL ZNI
;LERR (LIB,%,<DEXP3: base is 0 and exponent is le 0; result = infinity>)
HRLOI T0,377777 ;RESULT = INFINITY
HRLOI T1,377777
POPRET: POP P,T3 ;RESTORE ACCUMULATORS
POP P,T2
GOODBY (1) ;RETURN
ZERZER: $LCALL ZZZ ;"0**0 undefined, result = 0"
JRST POPRET
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
OVFL: $LCALL ROV ;RESULT OVERFLOW
HRLOI T0,377777 ; RESULT = INFINITY
HRLOI T1,377777
JRST RET3 ;GO TO RET3
WOK: CAML T0,SMALLW ;IF W IS NOT TOO SMALL
JRST WOK2 ;THEN PROCEED
UNFL: $LCALL RUN ;RESULT UNDERFLOWN
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: SKIPGE G1,T2 ;[4020] GET |W|
DMOVN G1,T2 ;[4020] THE HARD WAY (IF NEGATIVE)
SETZ G2, ;ZERO G2
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: CAIG G1,MSKTOP ;[4002] IF OUTSIDE TABLE, -1 IS MASK
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
SKIPGE T0,T4 ;[4020] SAVE A COPY OF ABS(W)
DMOVN T0,T4 ;[4020] THE HARD WAY (IF NEGATIVE)
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: CAIG G1,MSKTOP ;[4002] IF OUTSIDE TABLE, -1 IS MASK
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)
SKIPGE T4,T2 ;[4020] SAVE A COPY OF ABS(W2)
DMOVN T4,T2 ;[4020] THE HARD WAY (IF NEGATIVE)
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: CAIG G1,MSKTOP ;[4002] IF OUTSIDE TABLE, -1 IS MASK
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
JFCL EXCPT
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
EXCPT: JUMPG T5,OVFL ;IF T5 > 0 GO TO OVFL
JRST UNFL ;OTHERWISE, GO TO UNFL
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
MSKTOP==MSK26-MASK ;TOP INDEX ALLOWABLE TO MASK TABLE
SEGMENT DATA
PTEMP: 0
IFLAG: 0
PRGEND
TITLE DLOG10 LOG BASE 10 FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. APRIL 8, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DLOG10
EXTERN DLG10.
DLOG10=DLG10.
PRGEND
TITLE DLOG NATURAL LOG FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. APRIL 8, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DLOG
EXTERN DLOG.
DLOG=DLOG.
PRGEND
TITLE DLOG. NATURAL LOG FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. APRIL 8, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DLOG10(X) AND DLOG(X) ARE CALCULATED AS FOLLOWS
; FOR X = 0 AN ERROR MESSAGE IS ISSUED AND -MACHINE INFINITY
; IS RETURNED AS THE RESULT.
; FOR X < 0 AN ERROR MESSAGE IS ISSUED, X IS SET T0 -X AND
; CALCULATION CONTINUES.
; FOR X > 0, X = F*2**F(M), 1/2 < F < 1
; DEFINE G AND N SO THAT F = G*2(-N), 1/SQRT(2) <= G < SQRT(2).
; NOW
; DLOG(X) = (K*M-N) * DLOG(2) + DLOG(G)
; AND
; DLOG10(X) = DLOG10(E) * DLOG(X) = DLOG(X)/DLOG(10)
;
; DLOG(G) IS EVALUATED BY DEFINING S = (G-1)/(G+1) AND Z = 2*S
; AND THEN CALCULATING DLOG(G) = DLOG((1+Z/2)/(1-Z/2)) USING
; A MINIMAX RATIONAL APPROXIMATION.
;THE RANGE OF DEFINITION FOR DLOG/DLOG10 IS GIVEN ABOVE, AND ERROR MESSAGES
; WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE
;REQUIRED (CALLED) ROUTINES: NONE
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; 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 MTHPRM
SEGMENT CODE
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
$LCALL AZM
;LERR (LIB,%,DLOG or DLOG10: zero arg; result = -infinity)
HRLZI T0,400000 ;SET RESULT TO
HRRZI T1,000001 ;LARGE NEGATIVE NUMBER
GOODBY (1)
ARGN: $LCALL NAA
; LERR (LIB,%,DLOG or DLOG10: negative arg;result = LOG(ABS(arg)))
DMOVN T0,T0 ;ARG = -ARG
STRT: PUSH P,T2 ;SAVE ACCUMULATORS
PUSH P,T3
PUSH P,T4
PUSH P,T5
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
SEGMENT DATA
N: 0
FLAG: 0
SAVEW: DOUBLE 000000000000,000000000000
PRGEND
TITLE DMOD DOUBLE PRECISION
SUBTTL MARY PAYNE /MHP/CKS 07-Apr-83
;[3244] Replace DMOD routine
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
NOSYM
ENTRY DMOD
EXTERN DMOD.
DMOD=DMOD.
PRGEND
TITLE DMOD. DOUBLE PRECISION
SUBTTL MARY PAYNE /MHP/CKS 07-Apr-83
;[3244] Replace DMOD. routine
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
SEGMENT CODE
SALL
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
SKIPN @1(L) ;IF ARG2 IS ZERO
JRST DMRETZ ;GO RETURN 0 WITH MSG
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
DMRETZ: SETZB T0,T1 ;RETURN ZERO IF ARG2=0
$LCALL MZZ
JRST MRET ;GO RESTORE REGISTERS
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
SEGMENT DATA
B: BLOCK 2 ;ABS(B)
SAVE2: BLOCK 2 ;TEMP REGISTERS
SAVE4: BLOCK 2
SAVE6: BLOCK 2
SAVE8: BLOCK 2
PRGEND
TITLE DCOS SINE AND COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MAY 1, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DCOS
EXTERN DCOS.
DCOS=DCOS.
PRGEND
TITLE DSIN SINE AND COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MAY 1, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DSIN
EXTERN DSIN.
DSIN=DSIN.
PRGEND
TITLE DSIN. SINE AND COSINE FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MAY 1, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;COS(X) AND SIN(X) ARE CALCULATED AS FOLLOWS
; NOTE THAT COS(X) = COS(-X) AND SIN(X) = -SIN(-X)
; LET ABS(X) = N*PI + F WHERE ABS(F) <= PI/2
; WHEN ABS(F) < EPS (SEE CONSTANTS, BELOW), SIN(F) = F
; OTHERWISE, THE ARGUMENT REDUCTION IS:
; LET X1 = [ABS(X)], 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 MTHPRM
SEGMENT CODE
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: $LCALL ATZ
;LERR (LIB,%,DSIN or DCOS: ABS(arg) too large; result = zero)
MOVEI T0,0 ;SET RESULT
MOVEI T1,0 ;TO ZERO
POP P,T2 ;RESTORE ACCUMULATOR
GOODBY (1) ;
YOK: PUSH P,T3 ;SAVE ACCUMULATORS
PUSH P,T4 ;
PUSH P,T5
DMOVE T2,T0 ;SAVE A COPY OF ABS(X)
SKIPE FLAG ;IF THIS IS THE COS ROUTINE
DMOVE T2,XDEN ;RESTORE ABS(X) TO T2
DFMP T0,ODPI ;RN = Y/PI
JFCL
DMOVE T4,T0 ;SAVE A COPY OF RN
CAMGE T0,TWO26 ;IF RN IS .LT. 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
SEGMENT DATA
FLAG: 0
SGN: DOUBLE 0,0
XDEN: DOUBLE 0,0
PRGEND
TITLE DSINH HYPERBOLIC SINE FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 7, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DSINH
EXTERN DSINH.
DSINH=DSINH.
PRGEND
TITLE DSINH. HYPERBOLIC SINE FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 7, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;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 MTHPRM
SEGMENT CODE
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: $LCALL ROV
;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
SEGMENT DATA
SGNFLG: 0
TEMP: DOUBLE 0,0
Z: DOUBLE 0,0
PRGEND
TITLE DSQRT SQUARE ROOT FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MARCH 19, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DSQRT
EXTERN DSQRT.
DSQRT=DSQRT.
PRGEND
TITLE DSQRT. SQUARE ROOT FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. MARCH 19, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;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 MTHPRM
SEGMENT CODE
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
$LCALL NAA
;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) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DCOTAN
EXTERN DCOTN.
DCOTAN=DCOTN.
PRGEND
TITLE DTAN TANGENT AND COTANGENT FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC APRIL 16, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DTAN
EXTERN DTAN.
DTAN=DTAN.
PRGEND
TITLE DTAN. TANGENT AND COTANGENT FUNCTIONS
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC APRIL 16, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;TAN(X) AND COTAN(X) ARE CALCULATED AS FOLLOWS;
; THE IDENTITIES
; TAN(PI/2.0-G)=1.0/TAN(G)
; TAN(N*PI+H)=TAN(H), -PI/2.0 < H <= PI/2.0
; TAN(-X)=-TAN(X)
; COTAN(X)=1.0/TAN(X)
; COTAN(-X)=-COTAN(X)
; ARE USED TO REDUCE THE ORIGINAL PROBLEM TO A PROBLEM WITH
; AN ARGUMENT GREATER THAN -PI/2.0 AND LESS THAN OR EQUAL TO
; PI/2.0, PI=3.14159265358979324.
; THEN DEFINE N AND F SO THAT
; X=N*PI/4.0+F, 0.0 <= F <= PI/4.0, WITH PI=3.14159265358979324.
; THEN
; TAN(F) = F, IF F<2**(-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 MTHPRM
SEGMENT CODE
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: $LCALL ROV
;LERR (LIB,%,DCOTAN: Result overflow)
HRLOI T0,377777 ;SET RESULT TO MACHINE
HRLOI T1,377777 ;INFINITY
JUMPGE T4,ERR3 ;[3245]IF ARG IS NEGATIVE
DMOVN T0,T0 ;RESULT = - RESULT
ERR3: POP P,T5 ;[3245]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: $LCALL RUN
;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 T2
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
SEGMENT DATA
N: 0
FLAG: 0
PRGEND
TITLE DTANH HYPERBOLIC TANGENT FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DTANH
EXTERN DTANH.
DTANH=DTANH.
PRGEND
TITLE DTANH. HYPERBOLIC TANGENT FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;TANH(X) IS CALCULATED AS FOLLOWS
; TANH(X) = -TANH(-X). LET F = ABS(X)
; IF F > 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 MTHPRM
SEGMENT CODE
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
SEGMENT DATA
TEMP: DOUBLE 0,0
PRGEND
TITLE DEXP EXPONENTIAL FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL.INC APRIL 3, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DEXP
EXTERN DEXP.
DEXP=DEXP.
PRGEND
TITLE DEXP. EXPONENTIAL FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. APRIL 3, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DEXP(X) IS CALCULATED AS FOLLOWS
; IF X <= -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 MTHPRM
SEGMENT CODE
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: $LCALL ROV
;LERR (LIB,%,DEXP: result overflow)
HRLOI T0,377777 ;DEXP = +MACHINE INFINITY
HRLOI T1,377777
GOODBY (1) ;RETURN
OUT2: $LCALL RUN
;LERR (LIB,%,DEXP: result underflow)
MOVEI T0,0 ;EXP = 0
MOVEI T1,0
GOODBY (1) ;RETURN
EXP1: JUMPE T0,[MOVSI T0,(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 NOT, 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
SEGMENT DATA
SAVEN: 0
PRGEND
TITLE DABS DOUBLE PRECISION ABSOLUTE VALUE
SUBTTL H. P. WEISS/HPW 11-DEC-73
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1973, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DABS
EXTERN DABS.
DABS=DABS.
PRGEND
TITLE DABS. DOUBLE PRECISION ABSOLUTE VALUE
SUBTTL /KK/DMN/TWE
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;DOUBLE PRECISION ABSOLUTE VALUE FUNCTION
;THIS ROUTINE RETURNS THE ABSOLUTE VALUE OF A DOUBLE PRECISION
;ARGUMENT
SEARCH MTHPRM
SEGMENT CODE
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) DIGITAL EQUIPMENT CORPORATION 1973, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DMAX1
EXTERN DMAX1.
DMAX1=DMAX1.
PRGEND
TITLE DMAX1. DOUBLE PRECISION MAXIMUM OF A SERIES OF ARGUMENTS
SUBTTL /TWE/CKS 29-Jan-80
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
SEGMENT CODE
SALL
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) DIGITAL EQUIPMENT CORPORATION 1973, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DMIN1
EXTERN DMIN1.
DMIN1=DMIN1.
PRGEND
TITLE DMIN1. DOUBLE PRECISION MINIMUM OF A SERIES OF ARGUMENTS
SUBTTL /TWE/CKS 29-Jan-80
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
SEGMENT CODE
SALL
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) DIGITAL EQUIPMENT CORPORATION 1973, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DSIGN
EXTERN DSIGN.
DSIGN=DSIGN.
PRGEND
TITLE DSIGN. DOUBLE PRECISION SIGN ROUTINE
SUBTTL /KK/DMN/CKS 29-Jan-80
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;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 MTHPRM
SEGMENT CODE
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) DIGITAL EQUIPMENT CORPORATION 1973, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DFLOAT
EXTERN DFLOT.
DFLOAT=DFLOT.
PRGEND
TITLE DFLOT. INTEGER TO DOUBLE PRECISION CONVERSION
SUBTTL ED YOURDON/KK/VAA/TWE/DMN
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1973, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
EXTERN DFL.0
NOSYM
SEGMENT CODE
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) DIGITAL EQUIPMENT CORPORATION 1973, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY IDINT
EXTERN IDINT.
IDINT=IDINT.
PRGEND
TITLE IDINT. DOUBLE PRECISION TO INTEGER CONVERSION
SUBTTL TOM EGGERS/DMN
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1973, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
EXTERN IDF.0
NOSYM
SEGMENT CODE
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) DIGITAL EQUIPMENT CORPORATION 1973, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY DBLE
EXTERN DBLE.
DBLE=DBLE.
PRGEND
TITLE DBLE. REAL TO DOUBLE PRECISION CONVERSION
SUBTTL ED YOURDON/KK
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1973, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;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 MTHPRM
NOSYM
SEGMENT CODE
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 DFL.0
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
SEGMENT CODE
DFL 0
PRGEND
TITLE DFL.2
SEARCH MTHPRM
NOSYM
SEGMENT CODE
DFL 2
PRGEND
TITLE DFL.3
SEARCH MTHPRM
NOSYM
SEGMENT CODE
DFL 3
PRGEND
TITLE DFL.4
SEARCH MTHPRM
NOSYM
SEGMENT CODE
DFL 4
PRGEND
TITLE DFL.6
SEARCH MTHPRM
NOSYM
SEGMENT CODE
DFL 6
PRGEND
TITLE DFL.10
SEARCH MTHPRM
NOSYM
SEGMENT CODE
DFL 10
PRGEND
TITLE DFL.12
SEARCH MTHPRM
NOSYM
SEGMENT CODE
DFL 12
PRGEND
TITLE DFL.14
SEARCH MTHPRM
NOSYM
SEGMENT CODE
DFL 14
PRGEND
TITLE IDF.0
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
SEGMENT CODE
IDF 0
PRGEND
TITLE IDF.2
SEARCH MTHPRM
NOSYM
SEGMENT CODE
IDF 2
PRGEND
TITLE IDF.4
SEARCH MTHPRM
NOSYM
SEGMENT CODE
IDF 4
PRGEND
TITLE IDF.6
SEARCH MTHPRM
NOSYM
SEGMENT CODE
IDF 6
PRGEND
TITLE IDF.10
SEARCH MTHPRM
NOSYM
SEGMENT CODE
IDF 10
PRGEND
TITLE IDF.12
SEARCH MTHPRM
NOSYM
SEGMENT CODE
IDF 12
PRGEND
TITLE IDF.14
SEARCH MTHPRM
NOSYM
SEGMENT CODE
IDF 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) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY DDIM
EXTERN DDIM.
DDIM=DDIM.
PRGEND
TITLE DDIM. Double Precision Positive Difference.
SUBTTL C. McCutcheon 6/26/81
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
; Fortran Intrinsic Function for Positive Difference.
; Passed and returned arguments are double precision.
; Passed: A1,A2
; Returns:
; A1-A2 if A1 .GT. A2
; 0 if A1 .LE. A2
SEARCH MTHPRM
SEGMENT CODE
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
JFCL EXCEP ;Can overflow and underflow
RET: POP P,T3
POP P,T2
GOODBYE ;Return
RET0: SETZB T0,T1 ;Return zero
JRST RET
EXCEP: JUMPE T0,UNDER
$LCALL ROV,RET ;Result overflow
UNDER: $LCALL ROV,RET ;Result underflow
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) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY DINT
EXTERN DINT.
DINT=DINT.
PRGEND
TITLE DINT. Double Precision Truncation.
SUBTTL C. McCutcheon - 6/25/81
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
; Fortran Instrinsic Function to return a Double Precision
; truncation of the Double Precision number passed.
; If Magnitude(A) .LT. 1, then 0 is returned,
; otherwise return the largest integer that does not exceed the
; magnitude of A, and whose sign is the same as that of A.
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (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) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY DNINT
EXTERN DNINT.
DNINT=DNINT.
PRGEND
TITLE DNINT. Double Precision nearest whole number.
SUBTTL C. McCutcheon - 6/25/81
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
; Fortran Instrinsic Function to return the 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 MTHPRM
SEGMENT CODE
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) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.
SEARCH MTHPRM
SEGMENT CODE
NOSYM
ENTRY DPROD
EXTERN DPROD.
DPROD=DPROD.
PRGEND
TITLE DPROD. Double Precision product for single precision.
SUBTTL C. McCutcheon - 6/25/81
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
; Fortran Instrinsic Function to return a double precision
; product for two real arguments passed it.
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (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.
JFCL EXCEP ;Over/underflow can occur
RET: POPJ P,
EXCEP: JUMPE T0,UNDER ;Underflow?
$LCALL ROV,RET
UNDER: $LCALL RUN,RET
SEGMENT DATA
MULT: BLOCK 2 ;Multiplier
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) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.
SEARCH MTHPRM
SEGMENT CODE
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) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
; Fortran Intrinsic Function 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 MTHPRM
SEGMENT CODE
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: $LCALL ROV
;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