Trailing-Edge
-
PDP-10 Archives
-
BB-4157F-BM_1983
-
fortran/ots-debugger/mthcdx.mac
There are 9 other files named mthcdx.mac in the archive. Click here to see a list.
SEARCH MTHPRM
TV MTHCDX COMPLEX DOUBLE PRECISION ROUTINES ,7(3200)
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SUBTTL REVISION HISTORY
COMMENT \
***** Begin Revision History *****
1351 EDS 18-Mar-81 Q10-04786
Fix TWOSEG and RELOC problems.
1404 EGM 6-Apr-81 --------
Separate CDX and CGX routines, add GFL checks to CDX routines.
1405 DAW 6-Apr-81
Support extended addressing.
1464 DAW 12-May-81
Error messages.
***** Begin Version 6A *****
2055 TGS 10-May-82 10-32471
Fix typo in CDLOG so STOR POPJs instead of POPing P and
plowing into YZRO.
***** Begin Mathlib *****
3200 JLC 10-May-82
Mathlib integration. Change all LERRs to $LCALLs. Change
TWOSEG/RELOC to SEGMENT macros.
***** End Revision History *****
\
PRGEND
TITLE CDABS COMPLEX ABSOLUTE VALUE FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CDABS
EXTERN CDABS.
CDABS=CDABS.
PRGEND
TITLE CDABS. COMPLEX ABSOLUTE VALUE FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;CDABS(Z) IS CALCULATED AS FOLLOWS
; LET Z = X + I*Y
; V = MAX(ABS(X),ABS(Y))
; W = MIN(ABS(X),ABS(Y))
; THEN
; CDABS = V*SQRT(1.0 + (W/V)**2)
;THE RANGE OF DEFINITION FOR CDABS IS THE REPRESENTABLE REAL NUMBERS.
; AN OVERFLOW CAN OCCUR, IN WHICH CASE CDABS = + MACHINE INFINITY.
;REQUIRED (CALLED) ROUTINES: DSQRT
;REGISTERS T2, T3, T2, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; XMOVEI L,ARG
; PUSHJ P,CDABS
;THE ROUTINE IS CALLED WITH A DOUBLE PRECISION VECTOR.
;THE REAL PART OF THE ARGUMENT IS EXPECTED TO BE STORED IN THE
;FIRST DOUBLE PRECISION WORD. THE IMAGINARY PART OF THE ARGUMENT
;IS EXPECTED TO BE STORED IN THE SECOND DOUBLE PRECISION WORD.
;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN T0
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN T1
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (CDABS,.) ;ENTRY TO MODULUS ROUTINE
PUSH P,T2
PUSH P,T3
XMOVEI T2,@(L)
DMOVE T0,(T2) ;ABS(X)
JUMPGE T0,XPOS ;IF X IS NEGATIVE
DMOVN T0,T0 ;X = -X
XPOS: DMOVE T2,2(T2) ;ABS(Y)
JUMPGE T2,YPOS ;IF Y IS NEGATIVE
DMOVN T2,T2 ;Y = -Y
;OBTAIN MAX(ABS(X),ABS(Y))IN T2 AND MIN IN T0.
;ONLY THE HIGH WORDS ARE COMPARED, WHICH CAN RESULT IN THE MAX AND MIN
;BEING INTERCHANGED IF |X| AND |Y| ARE VERY NEARLY EQUAL. THIS IS OK.
;IT CAN RESULT IN SMALLER/LARGER BEING SLIGHTLY GREATER THAN 1.
YPOS: CAMG T0,T2 ;COMPARE |X|HI WITH |Y|HI
JRST LT ;|X| GREATER, NO EXCHANGE NECESSARY
EXCH T2,T0 ;EXCHANGE ABS(X) AND ABS(Y)
EXCH T3,T1 ;
LT: JUMPE T2,RET ;Z = 0, HENCE CDABS = O
DFDV T0,T2 ;U/V
JFCL ANS ;NO OVERFLOW, RATIO NEGLIGIBLE IF UNDERFLOW.
CAMG T0,TWOM32 ;IF RATIO .LE. 2**(-32)
JRST ANS ; RATIO IS NEGLIGIBLE.
DFMP T0,T0 ;**2
DFAD T0,ONE ;+1.0
DMOVEM T0,TEMP
FUNCT DSQRT.,<TEMP> ;SQUARE ROOT IS IN AC T0
DFMP T0,T2 ;*V
JFCL OVFL ;NO UNDERFLOW, GET MESSAGE ON OVERFLOW
RET: POP P,T3
POP P,T2
GOODBY (1) ;RETURN
OVFL: $LCALL ROV,RET
; LERR (LIB,%,<CDABS: result overflow>,RET)
ANS: DMOVE T0,T2 ;ANSWER = ABS(LARGER) TO T0
POP P,T3
POP P,T2
GOODBY (1) ;RETURN
ONE: DOUBLE 201400000000,000000000000
TWOM32: 141400000000 ;2**(-32)
SEGMENT DATA
TEMP: DOUBLE 0,0 ;TEMPORARY STORAGE USED FOR SQRT CALL
PRGEND
TITLE CDEXP COMPLEX EXPONENTIAL FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CDEXP
EXTERN CDEXP.
CDEXP=CDEXP.
PRGEND
TITLE CDEXP. COMPLEX EXPONENTIAL FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;CDEXP(Z), WHERE Z = X + I*Y, IS CALCULATED AS FOLLOWS
; CDEXP FIRST CHECKS THE TYPE OF THE ARGUMENT SO THAT THE
; APPROPRIATE RESULT (DOUBLE PRECISION OR EXTENDED
; PRECISION) IS RETURNED.
; CDEXP(Z) = DEXP(X)*(DCOS(Y)+I*DSIN(Y))
;THE RANGE OF DEFINITION FOR CDEXP IS AS FOLLOWS
;FOR Z = X+I*Y, IF DABS(Y) .GT. 6746518850.429 THE RESULT IS SET TO
; (0.0,0.0) AND AN ERROR MESSAGE IS RETURNED.
;FOR X.LT.-89.415986 THE RESULT IS SET TO (0.0,0.0) AND AN ERROR
; MESSAGE IS RETURNED.
;FOR X .GT. 88.029692;
; IF Y = 0.0, THE RESULT IS SET TO (+INFINITY,0.0) AND AN
; ERROR MESSAGE IS RETURNED.
; IF 88.029 < X < 176.058, AND A COMPONENT OF THE RESULT
; IS OUT OF RANGE, THEN AN ERROR MESSAGE IS RETURNED AND
; ONLY THAT COMPONENT IS SET TO +INFINITY.
; IF X/2. IS .GT. 88.029692, THE DABS(DREAL(RESULT)) IS SET TO
; +INFINITY AND DABS(DIMAG(RESULT)) IS SET TO +INFINITY AND
; AN ERROR MESSAGE IS RETURNED.
;REQUIRED (CALLED) ROUTINES: DEXP,DSIN,DCOS
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; XMOVEI L,ARG
; PUSHJ P,CDEXP
;THE ROUTINE IS CALLED WITH TWO DOUBLE PRECISION VECTORS.
;THE REAL PART OF THE ARGUMENT IS EXPECTED TO BE STORED IN THE
;FIRST DOUBLE PRECISION WORD OF THE FIRST VECTOR.
;THE IMAGINARY PART OF THE ARGUMENT IS EXPECTED TO BE STORED IN THE
;SECOND DOUBLE PRECISION WORD OF THE FIRST VECTOR.
;THE REAL PART OF THE SOLUTION IS RETURNED IN THE FIRST DOUBLE
;PRECISION WORD OF THE SECOND VECTOR.
;THE IMAGINARY PART OF THE SOLUTION IS RETURNED IN THE SECOND DOUBLE
;PRECISION WORD OF THE SECOND VECTOR.
SEARCH MTHPRM
FSRCH
SEGMENT CODE
SALL
HELLO (CDEXP.,CDEXP) ;ENTRY TO CDEXP ROUTINE
IFN FTGFL,<
LDB T1,[POINTR 0(L),ARGTYP] ;GET FIRST ARG TYPE
CAIN T1,TP%DPX ;G FLOATING?
JRST CGEXP## ;YES
> ;END FTGFL
PUSH P,T2 ;SAVE REGISTER T2
PUSH P,T3 ;SAVE REGISTER T3
PUSH P,T4 ;SAVE REGISTER T4
XMOVEI T2,@(L) ;GET ADDRESS OF Z
DMOVE T0,(T2) ;X=DREAL(Z)
DMOVE T2,2(T2) ;Y=IMAG(Z)
DMOVEM T2,YSAVE ;SAVE Y
JUMPGE T2,YPOS
DMOVN T2,T2
YPOS: CAMGE T2,YMAXH ;IF HI OF Y < YMAXH
JRST YOK ; GO TO YOK
CAME T2,YMAXH ;IF HI OF Y > YMAXH
JRST ERR1 ; GO TO ERR1
CAMG T3,YMAXL ;IF LO OF Y .LE. YMAXL
JRST YOK ; GO TO YOK
ERR1: $LCALL AIZ
; LERR (LIB,%,<entry CDEXP; DIMAG(arg) too large in absolute value; result=(0.0,0.0)>)
SETZ T0, ;REAL PART OF SOLUTION IS ZERO
SETZ T1,
SETZ T2, ;IMAG PART OF SOLUTION IS ZERO
SETZ T3,
XMOVEI T4,@1(L) ;GET ADDRESS OF 2ND ARGUMENT
DMOVEM T0,(T4) ;SAVE REAL PART OF RESULT
DMOVEM T2,2(T4) ;SAVE IMAG PART OF RESULT
POP P,T4 ;RESTORE ACCUMULATORS
POP P,T3
POP P,T2
GOODBY (1)
YOK: PUSH P,T5 ;SAVE ANOTHER REGISTER
DMOVEM T0,XSAVE ;SAVE X
CAMLE T0,NEGH ;IF HI OF X > NEGH
JRST XOK ; GO TO XOK
CAME T0,NEGH ;IF HI OF X < NEGH
JRST ZER ; GO TO ZER
CAML T1,NEGL ;IF LO OF X .GE. NEGL
JRST XOK ; GO TO XOK
ZER: SETZ T0, ;DREAL PART OF SOLUTION IS ZERO
SETZ T1,
SETZ T2, ;IMAG PART OF SOLUTION IS ZERO
SETZ T3,
JRST YCHK ;CHECK Y
XOK: FUNCT DCOS.,<YSAVE> ;DCOSY=DCOS(Y)
DMOVEM T0,DCOSY
FUNCT DSIN.,<YSAVE> ;DSINY=DSIN(Y)
DMOVEM T0,DSINY
DMOVE T4,XSAVE ;T4=X
CAMG T4,XMAXH ;IF X IS NOT TOO LARGE
JRST ALG ;GO TO ALG
SKIPN YSAVE ;ELSE, IF Y=0
JRST ERR0 ;GO TO ERR
CAMGE T4,TBIGH ;IF HI OF S < TBIGH
JRST SOK ; GO TO SOK
CAME T4,TBIGH ;IF HI OF S > TBIGH
JRST STIM ; GO TO STIM
CAMLE T5,TBIGL ;IF LO OF S > TBIGL
JRST STIM ; GO TO STIM
SOK: FSC T4,-1 ;ELSE, S=X/2
DMOVEM T4,DEXPX
FUNCT DEXP.,<DEXPX> ;T=EXP(S)
DMOVE T2,T0
DFMP T2,DSINY ;V=T*DSINY
JUMPGE T2,VPOS ;IF V IS NEGATIVE
DMOVN T2,T2 ;NEGATE IT
VPOS: HRLOI T4,377777 ;G3=XMAX
SETO T5,
DFDV T4,T0 ;Q=XMAX/T
DMOVEM T4,QH ;SAVE Q
CAMLE T2,T4 ;IF V .GT. Q
JRST STIM ;THEN GO TO STIM
CAME T2,QH ;IF V .LT. Q
JRST GETI ;THEN GO TO GETI
CAMG T3,QL ;IF V .LE. Q
JRST GETI ;THEN GO TO GETI
STIM: HRLOI T2,377777 ;ELSE, SET IMAG SOLUTION TO XMAX
SETO T3,
; LERR (LIB,%,<entry CDEXP; DREAL(arg) too large; DIMAG(result)=+infinity>)
$LCALL RTI
JRST D2
GETI: DFMP T2,T0 ;IRES = V*T
D2: SKIPGE DSINY ;IF SINY IS LESS THAN 0.0
DMOVN T2,T2 ;THEN NEGATE IRES
MOVE T4,XSAVE
CAMGE T4,TBIGH ;IF HI OF S < TBIGH
JRST OKS ; GO TO OKS
CAME T4,TBIGH ;IF HI OF S > TBIGH
JRST ERR0 ; GO TO ERR
CAMLE T5,TBIGL ;IF LO OF S > TBIGL
JRST ERR0 ; GO TO ERR
OKS: DMOVE T4,T0
DFMP T0,DCOSY ;V = T*DCOSY
JUMPGE T0,VOK ;IF V IS NEGATIVE
DMOVN T0,T0 ;NEGATE IT
VOK: CAMLE T0,QH ;IF V .GT. Q
JRST ERR0 ;THEN GO TO ERR
CAME T0,QH ;IF V .LT. Q
JRST LAB ;THEN GO TO LAB
CAMG T1,QL ;IF V IS .LE. Q
JRST LAB ;THEN GO TO LAB
ERR0: HRLOI T0,377777 ;RRES=XMAX
SETO T1,
; LERR (LIB,%,<entry CDEXP; DREAL(arg) too large; DREAL(result)=+infinity>)
$LCALL RTR
JRST DD2
LAB: DFMP T0,T4 ;RRES = V*T
DD2: SKIPGE DCOSY ;IF DCOSY .LT. 0.0
DMOVN T0,T0 ;THEN NEGATE RRES
XMOVEI T4,@1(L) ;GET ADDRESS OF 2ND ARGUMENT
DMOVEM T0,(T4) ;SAVE REAL PART OF RESULT
DMOVEM T2,2(T4) ;SAVE IMAG PART OF RESULT
POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4
POP P,T3
POP P,T2
GOODBY (1) ;RETURN
ALG: FUNCT DEXP.,<XSAVE> ;DEXPX = DEXP(X)
DMOVE T2,T0
DFMP T2,DSINY ;IRES = DEXPX*DSINY
JFCL
DFMP T0,DCOSY ;RRES = DEXPX*DCOSY
JFCL
JUMPN T2,RCHK ;IF IRES .NE. 0.0
;THEN GO CHECK RRES
YCHK: SKIPN YSAVE ;IF Y .EQ. 0
JRST RCHK ;THEN GO CHECK RRES
; LERR (LIB,%,<entry CDEXP; underflow has occurred; DIMAG(result)=0.0>)
$LCALL IPU
RCHK: JUMPN T0,RET ;IF R = 0.0
; LERR (LIB,%,<entry CDEXP; underflow has occurred; DREAL(result)=0.0>)
$LCALL RPU
RET: XMOVEI T4,@1(L) ;GET ADDRESS OF 2ND ARGUMENT
DMOVEM T0,(T4) ;SAVE REAL PART OF RESULT
DMOVEM T2,2(T4) ;SAVE IMAG PART OF RESULT
POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4
POP P,T3
POP P,T2
GOODBY (1) ;RETURN
YMAXH: 241622077325 ;HIGH ORDER PART OF YMAX
YMAXL: 011556457065 ;LOW ORDER PART OF YMAX
TBIGH: 210540074636 ;176.059383862226109
TBIGL: 077042573256
NEGH: 570232254036 ;-89.415986292232944914
NEGL: 301750634730
XMAXH: 207540074636 ;88.029691
SEGMENT DATA
QH: 0
QL: 0
DSINY: DOUBLE 0,0
DCOSY: DOUBLE 0,0
XSAVE: DOUBLE 0,0
DEXPX: DOUBLE 0,0
YSAVE: DOUBLE 0,0
PRGEND
TITLE CDLOG COMPLEX NATURAL LOG FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL Mary Payne 8-Sep-80
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;CDLOG(Z) IS CALCULATED AS FOLLOWS
; CDLOG FIRST CHECKS THE TYPE OF THE ARGUMENT SO THAT THE
; APPROPRIATE RESULT (D OR G) IS RETURNED.
; LET Z = X + I*Y
; IF Y IS NOT ZERO, THEN
; CDLOG(Z) = U + I*V
; U = (1/2) * (DLOG (X**2 + Y**2))
; V = DATAN2(Y,X)
; IF X IS NOT ZERO AND Y IS ZERO, THEN
; IF X IS POSITIVE, CDLOG(Z) = DLOG(X)
; IF X IS NEGATIVE, CDLOG(Z) = DLOG(ABS(X)) + I*PI
; IF X AND Y ARE ZERO, THEN
; CDLOG(Z) = +MACHINE INFINITY
; AND AN ERROR MESSAGE IS RETURNED
;REQUIRED (CALLED) ROUTINES: DLOG, DATAN, DATAN2
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; XMOVEI L,ARG
; PUSHJ P,CDLOG
;THE ROUTINE IS CALLED WITH TWO DOUBLE PRECISION VECTORS.
;THE REAL PART OF THE ARGUMENT IS EXPECTED TO BE STORED IN THE
;FIRST DOUBLE PRECISION WORD OF THE FIRST VECTOR.
;THE IMAGINARY PART OF THE ARGUMENT IS EXPECTED TO BE STORED IN THE
;SECOND DOUBLE PRECISION WORD OF THE FIRST VECTOR.
;THE REAL PART OF THE SOLUTION IS RETURNED IN THE FIRST DOUBLE
;PRECISION WORD OF THE SECOND VECTOR.
;THE IMAGINARY PART OF THE SOLUTION IS RETURNED IN THE SECOND DOUBLE
;PRECISION WORD OF THE SECOND VECTOR.
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (CDLOG) ;ENTRY TO CDLOG ROUTINE
IFN FTGFL,<
LDB T1,[POINT 4,0(L),12] ;GET FIRST ARG TYPE
CAIN T1,TP%DPX ;G FLOATING?
JRST CGLOG## ;YES
> ;END FTGFL
PUSH P,T2 ;SAVE ACCUMULATORS
PUSH P,T3
PUSH P,T4
XMOVEI T2,@(L) ;GET ADDRESS OF ARG
DMOVE T0,(T2) ;OBTAIN REAL PART OF ARGUMENT
DMOVE T2,2(T2) ;OBTAIN IMAG PART OF ARGUMENT
JUMPE T2,YZRO ;Y = 0 IS A SPECIAL CASE
JUMPE T0,XZRO ;X = 0 IS A SPECIAL CASE
PUSH P,T5 ;SAVE ANOTHER REGISTER
DMOVEM T0,TT0 ;SAVE COPIES OF X AND
DMOVEM T2,TT2 ; Y FOR FUTURE REFERENCE
JUMPGE T0,XPLUS ;IF X < 0
DMOVN T0,T0 ; NEGATE IT
XPLUS: JUMPGE T2,YPLUS ;IF Y < 0
DMOVN T2,T2 ; NEGATE IT
YPLUS: DMOVEM T0,TT4 ;SAVE MAGNITUDES OF
DMOVEM T2,TT6 ; X AND Y
AND T0,MASK ;ISOLATE EXPONENT FIELDS
AND T2,MASK ; OF X AND Y
SUB T2,T0 ;EXPONENT OF Y - EXPONENT OF X
CAML T2,P32 ;IF DIFFERENCE > 32
JRST BIG ; |Y/X| IS LARGE
CAMGE T2,M32 ;IF DIFFERENCE < -32
JRST SMALL ; |Y/X| IS SMALL
FUNCT DATN2.,<TT2,TT0> ;NO EXCEPTIONS CAN OCCUR
XMOVEI T4,@1(L) ;GET POINTER AND STORE
DMOVEM T0,2(T4) ; IMAGINARY PART OF RESULT
DMOVE T4,TT4 ;RECOVER |X|
DMOVE T2,TT6 ;RECOVER |Y|
CAML T4,T2 ;IF |X|HI .GE. |Y|HI,
JRST NOXCH ;NO EXCHANGE NECESSARY
EXCH T4,T2 ;|X|HI .LT. |Y|HI. LARGER = A
EXCH T5,T3 ; TO T4. SMALLER = B TO T2
NOXCH: MOVE T1,T4 ;HI OF LARGER TO T1
LSH T1,-33 ;ISOLATE ITS BIASED EXPONENT
SUBI T1,200 ; AND UNBIAS IT
JUMPLE T1,LE ;IF UNBIASED EXPONENT POSITIVE,
SUBI T1,1 ; DECREMENT IT
LE: MOVN T1,T1 ;NEGATE SCALE INDEX
FSC T4,(T1) ;SCALE A TO [1/2,2) AND B BY
FSC T2,(T1) ; SAME FACTOR. NO EXCEPTIONS
LSH T1,1 ;MULTIPLY NEGATED SCALE INDEX BY 2
DMOVEM T4,BIGGER ;SAVE SCALED LARGER OF |X| AND |Y|
DFMP T4,T4 ;SQUARE LARGER. NO EXCEPTIONS
DFMP T2,T2 ;SQUARE SMALLER. NO EXCEPTIONS
DFAD T4,T2 ;GET SCALED SQUARE OF MODULUS
CAMLE T4,RT2 ;IF T4 .GT. SQRT(2)HI
JRST MORE ; GO TO MORE SCALING
CAMGE T4,RT2O2 ;IF T4 .LT. SQRT(1/2)HI
JRST MORE ; GO TO MORE SCALING
JUMPN T1,JOIN ;IF SCALE INDEX NOT 0, GO TO JOIN
MOVEM T1,NN ; ELSE STORE 0 FOR FUTURE USE
;At this point the scale index is zero, which means that (A**2 + B**2)
;is unscaled. Also (A**2 + B**2) is between SQRT(1/2) and SQRT(2),
;which implies that a damaging loss of significance can occur in the
;calculation of (A**2 + B**2 - 1), which is one factor of Z, used
;in the rational approximation. Growth of error due to loss of
;significance can be avoided by calculating instead
;
; (A - 1) * (A + 1) + B**2
;
;which is done in the following lines of code:
DMOVE T4,BIGGER ;RECOVER LARGER
DMOVE T0,T4 ; AND COPY INTO T0
DFSB T0,ONE ;GET A - 1
DFAD T4,ONE ;AND A + 1
DFMP T4,T0 ;AND THEIR PRODUCT
DFAD T2,T4 ;T2 NOW HAS A**2 + B**2 - 1
DMOVE T4,T2 ;COPY IT INTO T4
DFAD T2,TWO ;GET DENOMINATOR OF Z IN T2
FSC T4,1 ; AND NUMERATOR IN T4
JRST MERGE ;MERGE FOR RATIONAL APPROX
;The following code applies to scaled values of (A**2 + B**2) lying
;outside the interval (SQRT(1/2),SQRT(2)); it carries out more
;scaling to get the scaled (A**2 + B**2) inside this interval.
MORE: MOVE T2,T4 ;HI OF T4 TO T2 FOR RESCALING
LSH T2,-33 ;ISOLATE ITS BIASED EXPONENT
SUBI T2,200 ;UNBIAS THE EXPONENT
MOVN T2,T2 ; AND NEGATE IT
FSC T4,(T2) ;T4 NOW IN [1/2,1)
ADD T1,T2 ;TOTAL NEGATED SCALE INDEX TO T1
CAML T4,RT2O2 ;IF T4 .GT. SQRT(1/2)
JRST JOIN ; SCALING IS DONE
FSC T4,1 ;SCALE T4 TO [1,SQRT(2)]
ADDI T1,1 ; AND INCREMENT NEGATED INDEX
;At this point the scaled (A**2 + B**2) is in the interval
;[SQRT(1/2),SQRT(2)]. T1 contains the negative of the index
;necessary to compensate later for the scaling. In the following
;lines of code, the index, together with its proper sign is
;temporarily stored. Then a rational approximation is used to
;evaluate the natural log of the scaled (A**2 + B**2).
JOIN: MOVNM T1,NN ;RESTORE SIGN OF INDEX AND SAVE
DMOVE T2,T4 ;GET COPY OF SCALED (A**2 + B**2)
FSC T4,1 ; AND MULTIPLY BY 2
DFSB T4,TWO ;SUBTRACT TWO TO GET NUMERATOR
DFAD T2,ONE ;ADD 1 TO GET DENOMINATOR OF Z
;The following code is taken from the DLOG routine. The constants
;A0, A1, A2, B0, B1, and B2 are those given in the DLOG routine.
;Note that all tests restricting the scaled (A**2 + B**2) to
;the interval (SQRT(1/2),SQRT(2)) have been carried out only
;on the "high" words of these limiting values. This is
;adequate for the validity of the approximation to be used,
;since only low order bits, and not the order of magnitude,
;of the error bound of the approximation are affected.
MERGE: DFDV T4,T2 ;Z = ZNUM / ZDEN. NO EXCEPTIONS
DMOVEM T4,SAVEZ ;SAVE A COPY OF Z
DFMP T4,T4 ;W = Z**2
DMOVE T0,T4 ; AND MAKE COPY
DFAD T4,B2 ;FORM B(W) = W + B2
DFMP T4,T0 ; * W
DFAD T4,B1 ; + B1
DFMP T4,T0 ; * W
DFAD T4,B0 ; + B0
DMOVE T2,T0 ;MAKE ANOTHER COPY OF 2
DFMP T0,A2 ;FORM A(W) = W * A2
DFAD T0,A1 ; + A1
DFMP T0,T2 ; * W
DFAD T0,A0 ; + A0
DFDV T0,T4 ;R(Z) = A(W) / B(W)
DFMP T0,SAVEZ ; * Z
DFMP T0,T2 ; * W
JFCL ;NO OVERFLOW, UNDERFLOW WON'T HURT
;As indicated in the DLOG routine, Z still needs to be added into R(Z).
;To increase accuracy, this addition is deferred until after the
;addition, below, of NN * LN(2))LO to the part of R(Z) so far
;obtained.
FLTR T2,NN ;RECOVER SCALE INDEX
JUMPN T2,REST ;IF NN = 0,
DFAD T0,SAVEZ ; COMPLETE R(Z) BY ADDING IN Z
JRST DONE ;GO STORE RESULT AND RESTORE REGISTERS
;The following lines of code complete the calculation of the real
;part for NN not zero. The real part is the sum of NN*(LN(2)HI),
;Z, R(Z), and NN*(LN(2)LO). Of these terms the first is guaranteed
;to have the largest magnitude. THe second has magnitude larger
;than the third by a factor of at least 2**7, and the last has
;magnitude smaller than the first by at least a factor of 2**(-10).
;The maximum value of |Z| is about 2/7, and this is the worst case
;from the point of view of accuracy, because it is only slightly
;smaller than (LN(2)HI). The calculation of Z is accurate to a
;small multiple of its LSB, and hence the calculation of the sum
;is limited to a small multiple of ITS LSB. Note that if |Z| is
;appreciably less than its bound of 2/7, the error bound of the
;sum decreases towards its LSB/2. THe overall accuracy will be
;maximized by summing the terms in the reverse order in which
;they are listed above.
REST: SETZ T3, ;GET DFLOAT(NN)
DMOVE T4,T2 ;GET COPY IN T4
DFMP T4,C2 ;GET NN * LN(2)LO
DFMP T2,C1 ;GET NN * LN(2)HI
DFAD T0,T4 ;R(Z) + NN * LN(2)LO
DFAD T0,SAVEZ ; + Z
DFAD T0,T2 ; + NN * LN(2)HI
DONE: FSC T0,-1 ;DIVIDE BY 2 TO GET LOG OF MODULUS
DMOVEM T0,@1(L) ;STORE REAL PART OF RESUUT
POP P,T5
POP P,T4
POP P,T3
POP P,T2
POPJ P, ;DONE
BIG: FUNCT DLOG.,<TT6> ;GET LN(|Y|)
DMOVE T2,DPIO2 ;IMAGINARY PART NEAR +/- PI/2
SKIPGE TT2 ; + FOR Y POSITIVE
DMOVN T2,T2 ; - FOR Y NEGATIVE
DMOVE T4,TT4 ;RECOVER |X|
DFDV T4,TT2 ;GET |X/Y|. NO OVERFLOW
JFCL UND1 ; UNDERFLOW IS POSSIBLE
JRST SMERG ;MERGE WITH SMALL
SMALL: FUNCT DLOG.,<TT4> ;GET LN(|X|)
SETZ T2, ;IMAGINARY PART IS NEAR
SETZ T3, ; 0 OR +/- PI
SKIPGE TT0 ;0 FOR X POSITIVE
DMOVE T2,CPI ; +/-PI FOR X NEGATIVE
SKIPGE TT2 ; + PI FOR Y POSITIVE
DMOVN T2,T2 ; - PI FOR Y NEGATIVE
DMOVE T4,TT6 ;RECOVER Y
DFDV T4,TT4 ;GET Y/X. NO OVERFLOW
JFCL UND2 ; UNDERFLOW IS POSSIBLE
SMERG: DFAD T2,T4 ;SMALL CORRECTION FOR IMAG
DFMP T4,T4 ;GET (Y/X)**2. NO OVERFLOW
JFCL UND1 ; INDERFLOW IS POSSIBLE
FSC T4,-1 ;(1/2)*(Y/X)**2. NO OVERFLOW
JFCL UND1 ; UNDERFLOW IS POSSIBLE
DFAD T0,T4 ;REAL = LN(|X|) + (1/2)*(Y/X)**2
JRST STOR ;NO EXCEPTIONS. GO STORE RESULTS
UND2: JUMPN T2,UND1 ;NO MESSAGE IF |IMAG| NEAR PI OR PI/2
; LERR (LIB,%,<CDLOG: Imaginary part underflow>)
$LCALL IPU
UND1: JUMPN T0,STOR ;NO MESSAGE IF LN NOT ZERO
; LERR (LIB,%,<CDLOG: Real part underflow>)
$LCALL RPU
STOR: XMOVEI T4,@1(L) ;GET ADRESS FOR RESULTS
DMOVEM T0,(T4) ;STORE REAL PART
DMOVEM T2,2(T4) ;STORE IMAGINARY PART
POP P,T5 ;RESTORE REGISTERS T5
POP P,T4
POP P,T3
POP P,T2 ;THROUGH T2
POPJ P, ;RETURN
YZRO: JUMPE T0,XYZRO ;X = Y = 0 IS A SPECIAL CASE
JUMPGE T0,RPOS ;Y IS 0, X ISN'T. IF X < 0
DMOVN T0,T0 ;X=ABS(X)
DMOVE T2,CPI ;V=PI
RPOS: DMOVEM T0,TT0 ;MOVE X TO MEMORY
FUNCT DLOG.,<TT0> ;COMPUTE U=DLOG(X); NO EXCEPTIONS
JRST DONE1
XZRO: DMOVE T0,T2 ;X IS 0, Y ISN'T. COPY Y TO T0
DMOVE T2,DPIO2 ;IMAG PART IS +/- PI/2
JUMPGE T0,YPOS ;IF Y < 0,
DMOVN T0,T0 ; MAKE IT POSITIVE, AND
DMOVN T2,T2 ; MAKE IMAG PART NEGATIVE
YPOS: DMOVEM T0,TT0 ;PREPARE TO GET LOG
FUNCT DLOG.,<TT0> ;GET LOG(|Y|), NO EXCEPTIONS
DONE1: XMOVEI T4,@1(L) ;GET THE ADDRESS OF THE SECOND ARGUMENT
DMOVEM T0,(T4) ;SAVE THE REAL PART OF THE ANSWER
DMOVEM T2,2(T4) ;SAVE THE IMAG PART OF THE ANSWER
POP P,T4
POP P,T3 ;RESTORE ACCUMULATORS
POP P,T2
POPJ P, ;RETURN
XYZRO: $LCALL ZIZ
;LERR (LIB,%,<CDLOG: arg is (0.0,0.0); result=(+infinity,0.0)>)
HRLOI T0,377777 ;REAL ANSWER IS POSITIVE INFINITY
SETO T1, ;
XMOVEI T4,@1(L) ;GET THE ADDRESS OF THE SECOND ARGUMENT
DMOVEM T0,(T4) ;SAVE THE REAL PART OF THE ANSWER
DMOVEM T2,2(T4) ;SAVE THE IMAG PART OF THE ANSWER
POP P,T4
POP P,T3 ;RESTORE ACCUMULATORS
POP P,T2
POPJ P, ;RETURN
;CONSTANTS
ONE: EXP 201400000000,000000000000 ;1.0
TWO: EXP 202400000000,000000000000 ;2.0
MASK: EXP 377000000000 ;EXPONENT MASK
P32: EXP 040000000000 ;HIGH 9 BITS = 40 OCTAL
M32: EXP 740000000000 ;HIGH 9 BITS = -40 OCTAL
CPI: EXP 202622077325,021026430215 ;PI
DPIO2: EXP 201622077325,021026430215 ;PI / 2
RT2O2: EXP 200552023632 ;SQRT(2) / 2
RT2: EXP 201552023632 ;SQRT(2)
A0: EXP 570377400073,123045145570 ;-.641249434237455811D2
A1: EXP 205406111210,005527754477 ;.163839435630215342D2
A2: EXP 577153575223,052241327104 ;-.789561128874912573D0
B0: EXP 565177200130,374467610432 ;-.769499321084948798D3
B1: EXP 211470020376,205223175724 ;.312032220919245328D3
B2: EXP 571342517755,046030761117 ;-.356679777390346462D2
B3: EXP 576400000000,000000000000 ;-.1D1
C1: EXP 200543000000,000000000000 ;LN(2)HI
C2: EXP 613102775750,347571527450 ;LN(2)LO
;DATA
SEGMENT DATA
TT0: BLOCK 2 ;FOR DATN2 AND DLOG CALLS
TT2: BLOCK 2 ;FOR ATAN2 CALL
TT4: BLOCK 2 ;FOR ABS(X)
TT6: BLOCK 2 ;FOR ABS(Y)
NN: BLOCK 2 ;FOR SCALE INDEX
BIGGER: BLOCK 2
SAVEZ: BLOCK 2
PRGEND
TITLE CDSIN COMPLEX SINE FUNCTION
; (DOUBLE PRECISION)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CDSIN
EXTERN CDSIN.
CDSIN=CDSIN.
PRGEND
TITLE CDCOS COMPLEX COSINE FUNCTION
; (DOUBLE PRECISION)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CDCOS
EXTERN CDCOS.
CDCOS=CDCOS.
PRGEND
TITLE CDSIN. COMPLEX SINE AND COSINE ROUTINES
; (DOUBLE PRECISION)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;CDSIN(Z) AND CDCOS(Z), WHERE Z = X + I*Y, ARE CALCULATED AS FOLLOWS
; CDSIN FIRST CHECKS THE TYPE OF THE ARGUMENT SO THAT THE
; APPROPRIATE RESULT (DOUBLE PRECISION OR EXTENDED
; PRECISION) WILL BE RETURNED.
; CDSIN(Z) =DSIN(X)*DCOSH(Y) + I*DCOS(X)*DSINH(Y)
; CDCOS(Z) =DCOS(X)*DCOSH(Y) - I*DSIN(X)*DSINH(Y)
; THE FUNCTIONS DSINH AND DCOSH ARE CODED IN LINE.
; THE REAL PART OF CDSIN AND THE IMAGINARY PART OF CDCOS CANNOT UNDERFLOW
; BECAUSE COSH(Y) IS NEVER LESS THAN 1.
; THE OTHER PART OF THE RESULT MAY UNDERFLOW; IF THIS HAPPENS IT IS SET TO 0
; AND A MESSAGE IS PROVIDED.
;THE RANGE OF DEFINITION FOR CDSIN AND CDCOS IS AS FOLLOWS
;FOR
; Z = X+I*Y. IF DABS(X) > ((2**31)*PI - PI/2) THE RESULT IS SET TO 0.0 AND
; AN ERROR MESSAGE IS RETURNED.
;FOR DABS(Y) > 88.029692 CALCULATIONS PROCEED AS FOLLOWS:
; FOR THE REAL PART OF THE RESULT
; LET T = ABS(DSIN(X)), THEN
; FOR T = 0.0 THE REAL PART OF THE RESULT IS SET TO 0.0
; FOR DLOG(T) + DABS(Y) - DLOG(2.0) > 88.722839 THE REAL PART OF THE
; RESULT IS SET TO PLUS OR MINUS MACHINE INFINITY AND
; AN ERROR MESSAGE IS RETURNED. (88.722839=88.029692+LN2)
; FOR THE IMAGINARY PART OF THE RESULT
; LET T = DABS(DCOS(X)), THEN
; FOR T = 0.0 THE IMAGINARY PART OF THE RESULT IS SET TO 0.0
; FOR DLOG(T) + DABS(Y) - DLOG(2.0) > 88.722839 THE IMAGINARY PART OF THE
; RESULT IS SET TO PLUS OR MINUS MACHINE INFINITY AND
; AN ERROR MESSAGE IS RETURNED. (88.722839=88.029692+LN2)
;REQUIRED (CALLED) ROUTINES: DSIN,DCOS,DEXP,DLOG
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; XMOVEI L,ARG
; PUSHJ P,CDSIN
; OR
; PUSHJ P,CDCOS
;THE ROUTINE IS CALLED WITH TWO DOUBLE PRECISION VECTORS.
;THE REAL PART OF THE ARGUMENT IS EXPECTED TO BE STORED IN THE
;FIRST DOUBLE PRECISION WORD OF THE FIRST VECTOR.
;THE IMAGINARY PART OF THE ARGUMENT IS EXPECTED TO BE STORED IN THE
;SECOND DOUBLE PRECISION WORD OF THE FIRST VECTOR.
;THE REAL PART OF THE SOLUTION IS RETURNED IN THE FIRST DOUBLE
;PRECISION WORD OF THE SECOND VECTOR.
;THE IMAGINARY PART OF THE SOLUTION IS RETURNED IN THE SECOND DOUBLE
;PRECISION WORD OF THE SECOND VECTOR.
SEARCH MTHPRM
FSRCH
SEGMENT CODE
SALL
HELLO (CDCOS.,CDCOS) ;ENTRY TO CDCOS ROUTINE
IFN FTGFL,<
LDB T1,[POINTR 0(L),ARGTYP] ;GET FIRST ARG TYPE
CAIN T1,TP%DPX ;G FLOATING?
JRST CGCOS## ;YES
> ;END FTGFL
PUSH P,T5
MOVEI T5,1 ;SET FLAG TO 1 FOR CDCOS
JRST STRT ;GO TO CDSIN ROUTINE
HELLO (CDSIN.,CDSIN) ;ENTRY TO CMPLX SINE ROUTINE
IFN FTGFL,<
LDB T1,[POINTR 0(L),ARGTYP] ;GET FIRST ARG TYPE
CAIN T1,TP%DPX ;G FLOATING?
JRST CGSIN## ;YES
> ;END FTGFL
PUSH P,T5
BEG: SETZI T5, ;SET FLAG TO 0 FOR CDSIN
STRT: PUSH P,T2
PUSH P,T3
PUSH P,T4
XMOVEI T2,@(L) ;GET ADDRESS OF ARG
DMOVE T0,(T2) ;X = REAL(Z)
DMOVE T3,T0 ;MOVE REAL PART OF ARG TO T4,T5
JUMPGE T3,XPOS ;IF REAL PART OF ARG < 0
DMOVN T3,T3 ;REAL PART=ABS(REAL PART)
XPOS: CAMLE T3,XMAXH ;IF ABS(X) IS TOO LARGE
JRST ERR1 ;GO TO ERR1
CAME T3,XMAXH ;IF ABS(X) IS LESS THAN XMAXH
JRST XOK ;GO TO XOK
CAMG T4,XMAXL ;IF ABS(X) IS NOT TOO LARGE
JRST XOK ;GO TO XOK
ERR1: $LCALL ARZ
;LERR (LIB,%,<CDSIN or CDCOS; DABS(DREAL(arg)) too large; result = (0.0,0.0)>)
SETZI T0, ;SET REAL(RESULT) TO 0.0
SETZ T1,
SETZI T2, ;SET IMAG(RESULT) TO 0.0
SETZ T3,
XMOVEI T4,@1(L) ;GET ADDRESS OF 2ND ARGUMENT
DMOVEM T0,(T4) ;SAVE REAL PART OF RESULT
DMOVEM T2,2(T4) ;SAVE IMAGINARY PART OF RESULT
POP P,T4 ;RESTORE ACCUMULATORS
POP P,T3
POP P,T2
POP P,T5
GOODBY (1) ;RETURN
XOK: DMOVE T2,2(T2) ;Y = IMAG(Z)
DMOVEM T2,YSAVE ;SAVE Y
DMOVEM T0,ARGSAV
FUNCT DSIN.,<ARGSAV> ;CALCULATE DSIN(X)
DMOVEM T0,SX ;SX = DSIN(X)
FUNCT DCOS.,<ARGSAV> ;CALCULATE DCOS(X)
DMOVEM T0,CX ;CX = DCOS(X)
JUMPE T5,NOSWT ;IF THIS IS THE CDSIN ROUTINE
;THEN GO TO NOSWT
DMOVN T2,SX ;OTHERWISE, PUT -DSIN(X)
EXCH T2,CX ;IN CX, AND COS(X)
EXCH T3,CXL
DMOVEM T2,SX ;IN SX.
NOSWT: DMOVE T2,YSAVE ;GET A COPY OF Y
NE: JUMPG T2,YPOS ;IF Y IS NEGATIVE
DMOVN T2,T2 ;NEGATE IT
YPOS: CAMGE T2,XHI ;IF HI OF AY < XHI
JRST YSMALL ; GO TO YSMALL
CAME T2,XHI ;IF HI OF AY > XHI
JRST OVFL ; GO TO OVFL
CAMLE T3,XLO ;IF LO OF AY > XLO
JRST OVFL ; GO TO OVFL
YSMALL: CAMG T2,TWOM31 ;IF ABS(Y) .LE. 2**(-31)
JRST EASY ; GO TO EASY CODE.
CAMGE T2,ONE ;IS ABS(Y) .GE. ONE?
JRST LSTH1 ; BRANCH IF NOT
DMOVEM T2,ARGSAV ;GET COPY OF ABS(Y)
FUNCT DEXP.,<ARGSAV> ;GET EXP(ABS(Y))
DMOVE T4,ONE ;GET EXP(-ABS(Y)) BY
DFDV T4,T0 ; RECIPROCATION.
DFSB T0,T4 ;T0 GETS
FSC T0,-1 ; SINH(ABS(Y))
JRST CCCC ;TO REST OF CALCULATION
LSTH1: DMOVNM T2,ARGSAV ;GET -ABS(Y)
FUNCT DEXP.,<ARGSAV> ; EXP(-ABS(Y))
DMOVE T4,T0
DFMP T2,T2 ;SS = AY**2
DMOVEM T2,TEMP ;SAVE A COPY OF SS
DFAD T2,Q2 ;XDEN = Q2 +SS
DFMP T2,TEMP ;*SS
DFAD T2,Q1 ;+Q1
DFMP T2,TEMP ;*SS
DFAD T2,Q0 ;+Q0
DMOVE T0,TEMP ;SAVE A COPY OF SS
DFMP T0,RP3 ;XNUM = RP3*SS
DFAD T0,RP2 ;+RP2
DFMP T0,TEMP ;*SS
DFAD T0,RP1 ;+RP1
DFMP T0,TEMP ;*SS
DFAD T0,RP0 ;+RP0
DFDV T0,T2 ;XNUM/XDEN
DMOVN T2,ARGSAV ;GET AY BACK
DFMP T0,TEMP ;*SS
DFMP T0,ARGSAV ;*AY
DFAD T0,T2 ;+AY
CCCC: DMOVE T2,T0
DFAD T0,T4 ;CC = SS+(1/T)
DMOVE T4,YSAVE ;RESTORE Y
JUMPGE T4,YOK ;IF Y IS LESS THAN 0.0
DMOVN T2,T2 ;THEN NEGATE SS
YOK: DFMP T0,SX ;GET REAL PART OF RESULT
DFMP T2,CX ;GET IMAG PART OF RESULT
JFCL IMUND ;UNDERFLOW POSSIBLE
XMOVEI T4,@1(L) ;GET ADDRESS OF 2ND ARGUMENT
DMOVEM T0,(T4) ;SAVE REAL PART OF RESULT
DMOVEM T2,2(T4) ;SAVE IMAGINARY PART OF RESULT
POP P,T4 ;RESTORE ACCUMULATORS
POP P,T3
POP P,T2
POP P,T5
GOODBY (1) ;RETURN
OVFL: DMOVE T0,SX ;T=SX
JUMPE T0,GOTR ;IF T IS EQUAL TO 0.
JUMPG T0,SXPOS ;IF SX IS NEGATIVE
DMOVN T0,T0 ;NEGATE SX
SXPOS: DMOVEM T0,ARGSAV ;OTHERWISE
FUNCT DLOG.,<ARGSAV> ;CALCULATE LOG(T)
DFAD T0,T2 ;T = LOG(T)+AY
DFAD T0,LN2 ;T = T-LOG(2.0)
CAMGE T0,XHI ;IF HI OF T < XHI
JRST CONTIN ; GO TO CONTIN
CAME T0,XHI ;IF HI OF T > XHI
JRST ERR2 ; GO TO ERR2
CAMG T1,XLO ;IF LO OF T .LE. XLO
JRST CONTIN ; GO TO CONTIN
ERR2: $LCALL AIR
;LERR (LIB,%,<CDSIN or CDCOS; DABS(DIMAG(arg)) too large; DREAL(result) = infinity>)
HRLOI T0,377777 ;REAL PART OF RESULT SET TO INFINITY
SETO T1, ;SET SECOND WORD
JRST SKP2 ;SKIP NEXT 2 INSTRUCTIONS
CONTIN: DMOVEM T0,ARGSAV
FUNCT DEXP.,<ARGSAV> ;RRES = EXP(T)
SKP2: SKIPGE SX ;IF SX IS LESS THAN 0.
DMOVN T0,T0 ;THEN NEGATE RRES
GOTR: SKIPN T4,CX ;IF CX = 0,
JRST IMAG0 ;THEN PREPARE TO RETURN
MOVE T5,CXL ;OTHERWISE GET LOWER WORD
DMOVEM T0,SX ;SAVE RRES
JUMPGE T4,TPOS ;IF T IS NEGATIVE
DMOVN T4,T4 ;THEN GET ABS(T)
TPOS: DMOVEM T4,ARGSAV
FUNCT DLOG.,<ARGSAV> ;CALCULATE T=LOG(T)
DFAD T2,T0 ;T = T+AY
DFAD T2,LN2 ;T = T-LOG(2)
CAMGE T2,XHI ;IF HI OF T < XHI
JRST CONTN2 ; GO TO CONTN2
CAME T2,XHI ;IF HI OF T > XHI
JRST ERR3 ;GO TO ERR3
CAMG T3,XLO ;IF LO OF T .LE. XLO
JRST CONTN2 ; GO TO CONTN2
ERR3: $LCALL AII
;LERR (LIB,%,CDSIN or CDCOS; DABS(DIMAG(arg)) too large; DIMAG(result) = infinity)
HRLOI T2,377777 ;SET IRES TO INFINITY
SETO T3,
JRST SKP3 ;SKIP NEXT 3 INSTRUCTIONS
CONTN2: DMOVEM T2,ARGSAV
FUNCT DEXP.,<ARGSAV> ;IRES = EXP(T)
DMOVE T2,T0
SKP3: DMOVE T0,SX
MOVE T4,YSAVE
XOR T4,CX ;T4 HAS THE SIGN OF CX*Y
JUMPGE T4,RET ;IF T4 IS LESS THAN 0.
DMOVN T2,T2 ;THEN NEGATE IRES
JRST RET ;RETURN
IMAG0: SETZ T2, ;SET IMAGINARY PART
SETZ T3, ; OF RESULT TO ZERO.
JRST RET ;RETURN
EASY: DMOVE T2,YSAVE ;GET SIGNED VALUE OF Y
DFMP T2,CX ;SINH(Y)*COS(X) = Y*COS(X)
JFCL IMUND ;UNDERFLOW POSSIBLE.
DMOVE T0,SX ;COSH(Y)*SIN(X) = SIN(X).
RET: XMOVEI T4,@(L) ;GET ADDRESS OF 2ND ARGUMENT
DMOVEM T0,(T4) ;SAVE REAL PART OF RESULT
DMOVEM T2,2(T4) ;SAVE IMAGINARY PART OF RESULT
POP P,T4
POP P,T3
POP P,T2
POP P,T5
GOODBY (1)
IMUND: $LCALL IPU,RET
; LERR (LIB,%,<CDSIN or CDCOS: imaginary part underflow>,RET)
TWOM31: 142400000000 ;2**(-31)
LN2: DOUBLE 577235067500,101343020625 ;-(NATURAL LOG OF 2.0)
XHI: 207540074636 ;88.029691931113054295
XLO: 077042573256 ;
XMAXH: 241622077325 ;HIGH ORDER OF 2**31 * PI
XMAXL: 012605434744 ;LOW ORDER OF 2**31 * PI
ONE: DOUBLE 201400000000,000000000000
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,320314355210 ;-.21108770058106271242D+7
Q1: DOUBLE 220432412710,355457306744 ;.36162723109421836460D+5
Q2: DOUBLE 566352207437,215500015770 ;-.27773523119650701667D+3
SEGMENT DATA
CX: 0
CXL: 0
SX: DOUBLE 0,0
YSAVE: DOUBLE 0,0
ARGSAV: DOUBLE 0,0
TEMP: DOUBLE 0,0
PRGEND
TITLE CDSQRT COMPLEX SQUARE ROOT FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CDSQRT
EXTERN CDSQT.
CDSQRT=CDSQT.
PRGEND
TITLE CDSQT. COMPLEX SQUARE ROOT FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
; LET Z = X + I*Y
; THEN THE ANSWER CDSQRT(Z) = U + I*V IS DEFINED AS FOLLOWS
; CDSQRT FIRST CHECKS THE TYPE OF THE ARGUMENT SO THAT THE
; APPROPRIATE RESULT (DOUBLE PRECISION OR EXTENDED
; PRECISION) IS RETURNED.
; IF X AND Y ARE BOTH >= 0.0, THEN
; U = DSQRT((DABS(X)+CDABS(Z))/2.0)
; V = Y/(2.0*U)
; IF X >= 0.0 AND Y < 0.0, THEN
; U = SQRT((DABS(X)+CDABS(Z))/2.0)
; V = Y/(2.0*U)
; IF X < 0.0 AND Y >= 0.0, THEN
; U = Y/(2.0*V)
; V = DSQRT((DABS(X)+CDABS(Z))/2.0)
; IF X AND Y ARE BOTH < 0.0, THEN
; U = Y/(2.0*V)
; V = -DSQRT((DABS(X)+CDABS(Z))/2.0)
; THE RESULT IS IN THE RIGHT HALF PLANE, THAT IS, THE POLAR ANGLE OF THE
; RESULT LIES IN THE INTERVAL [-PI/2,+PI/2].
;REQUIRED (CALLED) ROUTINES: DSQRT
;REGISTERS T2, T3, T4, T5, G1, G3, AND G4 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; XMOVEI L,ARG
; PUSHJ P,CDSQRT
;THE ROUTINE IS CALLED WITH TWO DOUBLE PRECISION VECTORS.
;THE REAL PART OF THE ARGUMENT IS EXPECTED TO BE STORED IN THE
;FIRST DOUBLE PRECISION WORD OF THE FIRST VECTOR.
;THE IMAGINARY PART OF THE ARGUMENT IS EXPECTED TO BE STORED IN THE
;SECOND DOUBLE PRECISION WORD OF THE FIRST VECTOR.
;THE REAL PART OF THE SOLUTION IS RETURNED IN THE FIRST DOUBLE
;PRECISION WORD OF THE SECOND VECTOR.
;THE IMAGINARY PART OF THE SOLUTION IS RETURNED IN THE SECOND DOUBLE
;PRECISION WORD OF THE SECOND VECTOR.
SEARCH MTHPRM
FSRCH
SEGMENT CODE
SALL
HELLO (CDSQT.,CDSQRT) ;ENTRY TO COMPLEX SQUARE ROOT ROUTINE.
IFN FTGFL,<
LDB T1,[POINTR 0(L),ARGTYP] ;GET FIRST ARG TYPE
CAIN T1,TP%DPX ;G FLOATING?
JRST CGSQRT## ;YES
> ;END FTGFL
PUSH P,T2 ;SAVE ACCUMULATORS
PUSH P,T3
XMOVEI T2,@0(L) ;PICK UP ADDRESS OF ARGUMENT
DMOVE T0,(T2) ;PICK UP X IN AC T0.
DMOVE T2,2(T2) ;PICK UP Y IN AC T2.
PUSH P,T4 ;SAVE AC T4.
PUSH P,T5
JUMPE T2,ZREAL ;JUMP IF Y=0
JUMPE T0,ZIMAG ;JUMP IF X=0
DMOVE T4,T0 ;DABS(X) TO AC T4.
JUMPGE T0,XPOS ;IF X IS NEGATIVE
DMOVN T4,T4 ;NEGATE
XPOS: PUSH P,G1 ;SAVE REGISTER G1.
PUSH P,G3 ;SAVE REGISTER G3
PUSH P,G4 ;SAVE REGISTER G4
DMOVEM T0,XSAVE ;SAVE X IN XSAVE
DMOVEM T2,YSAVE ;SAVE Y IN YSAVE
DMOVE G3,T2 ;DABS(Y) TO AC G3.
JUMPGE T2,YPOS ;IF Y IS NEGATIVE
DMOVN G3,T2 ;NEGATE IT
YPOS: HRRZI G1,1 ;SET G1 TO 1
CAML T4,G3 ;IF ABS(X) > ABS(Y)
JRST DWN ;THEN GO TO DWN
;IF HIGH WORDS OF |X| AND |Y| ARE EQUAL,
; DON'T BOTHER DECIDING WHETHER TO EXCHANGE
; THEM. IF |X| AND |Y| ARE NEARLY EQUAL,
; THEIR MAX AND MIN MAY BE SWITCHED BELOW.
; THIS MEANS THAT S/L MAY BE SLIGHTLY GREATER
; THAN 1, BUT THIS DOES NOT HURT.
XCHNG: EXCH T4,G3 ;PUT MAX(ABS(X),ABS(Y)) IN T4
EXCH T5,G4 ;AND MIN(ABS(X),ABS(Y)) IN G3,G4
SETZ G1, ;SET G1 TO 0
DWN: DFDV G3,T4 ;CALC S/L.
JFCL ;NO UNDERFLOW ERROR MESSAGE ALLOWED.
DFMP G3,G3 ;CALC (S/L)**2.
JFCL ;NO UNDERFLOW ERROR MESSAGE ALLOWED.
DFAD G3,ONE ;HAVE (1+(S/L)**2) IN AC G3.
DMOVEM G3,TEMP
FUNCT DSQRT.,<TEMP> ;CALC. THE DSQRT OF
;(1+(S/L)**2)=1+EPS.
JUMPE G1,YGTX ;GO TO YGTX IF DABS(Y) > DABS(X).
XGTY: DFAD T0,ONE ;CALC. 2 + EPS.
FSC T0,-1 ;CALC. 1+EPS/2.
DMOVEM T0,TEMP ;PUT IN TEMP FOR
FUNCT DSQRT.,<TEMP> ;CALL TO DSQRT
DMOVE G3,T0 ;SAVE DSQRT(1+EPS/2) IN AC G3.
DMOVEM T4,TEMP
FUNCT DSQRT.,<TEMP> ;CALC.
DFMP T0,G3 ;CALC. DSQRT(DABS(X)*(1+EPS/2)).
JRST OUT1 ;GO TO REST OF CALC.
YGTX: CAIGE T4,201400 ;IF DABS(Y)>1, GO TO POSSIBLE OVERFLOW CASE.
JRST YXSMAL ;ELSE, GO TO POSSIBLE UNDERFLOW.
FSC T0,-3 ;CALC. (1+EPS)/8
DFMP T4,T0 ;CALC. DABS(Y)*(1+EPS)/8
DMOVE T0,XSAVE ;DABS(X) TO AC T0
JUMPGE T0,NXT
DMOVN T0,T0
NXT: FSC T0,-3 ;CALC. DABS(X)/8
JFCL ;SUPPRESS UNDERFLOW ERROR MESSAGE.
DFAD T0,T4 ;CALC.(DABS(X)/8)+(DABS(Y)*(1+EPS)/8).
DMOVEM T0,TEMP
FUNCT DSQRT.,<TEMP> ;CALC.
FSC T0,1 ;DSQRT(DABS(X)/2 + (DABS(Y)*(1+EPS)/2)).
OUT1: JUMPGE T2,POSY
DMOVN T2,T2
POSY: DFDV T2,T0 ;DIVIDE DABS(Y)/2 BY THE
JFCL UNDFL
FSC T2,-1 ;SQRT TERM.
JFCL UNDFL
JRST SIGNS ;GO TO REST OF CALC.
YXSMAL: DFMP T0,T4 ;DABS(Y)*(1+EPS) = CDABS(Z).
DMOVE G3,XSAVE ;DABS(X) TO AC G3.
JUMPGE G3,XOK ;IF X IS NEGATIVE
DMOVN G3,G3 ;NEGATE IT
XOK: DFAD G3,T0 ;GET |X| + CDABS(Z)
FSC G3,1 ;MULTIPLY ARG BY 2
DMOVEM G3,TEMP ;SAVE FOR DSQRT
FUNCT DSQRT.,<TEMP> ;CALC. DSQRT OF 2 * (|X| + CDABS(Z))
DMOVE G3,T0 ;GET DSQRT((|X| + CDABS(Z))*2)
FSC T0,-1 ;GET DSQRT((|X| + CDABS(Z))/2)
DFDV T4,G3 ;GET |Y| / DSQRT (2 * (|X| + CDABS(Z)))
DMOVE T2,T4 ;PUT A TEMP ANSWER IN AC T2.
SIGNS: SKIPL XSAVE ;IF NEGATIVE X, SWITCH REAL & IMAG
JRST OK ; PARTS OF RESTLT
EXCH T2,T0 ;EXCHANGE
EXCH T1,T3 ;EXCHANGE
;SET UP THE REAL AND
;THE IMAGINARY ANSWERS WITH
OK: SKIPGE YSAVE ;THE APPROPRIATE
DMOVN T2,T2 ;SIGNS.
XMOVEI T4,@1(L) ;GET THE ADDRESS OF THE SECOND ARGUMENT
DMOVEM T0,(T4) ;SAVE REAL PART OF RESULT
DMOVEM T2,2(T4) ;SAVE IMAG PART OF RESULT
POP P,G4 ;RESTORE AC G4
POP P,G3 ;RESTORE AC G3
POP P,G1 ;RESTORE AC G1
POP P,T5 ;RESTORE AC T5
POP P,T4 ;RESTORE AC T4
POP P,T3 ;RESTORE AC T3
POP P,T2 ;RESTORE AC T2
GOODBY (1) ;RETURN
ZREAL: JUMPE T0,DONE ;Y = 0, HOW ABOUT X?
MOVE T4,T0 ;X NOT ZERO, SAVE IT.
JUMPGE T0,POSX ;GET DABS(X) FOR DSQRT
DMOVN T0,T0
POSX: DMOVEM T0,TEMP ;
FUNCT DSQRT.,<TEMP> ;T0,T1 GET DSQRT(DABS(X))
JUMPG T4,DONE ;T0,T1 AND T2,T3 OK FOR X .GE. 0
EXCH T0,T2 ;INTERCHANGE FOR X<0
EXCH T1,T3 ;
JRST DONE ;FINISHED
ZIMAG: MOVE T4,T2 ;X=0, SAVE Y
JUMPGE T2,PSY ;IF Y < 0
DMOVN T2,T2 ; NEGATE IT
PSY: FSC T2,-1 ;GET |Y|/2
DMOVEM T2,TEMP ; AND STORE IT FOR DSQRT
FUNCT DSQRT.,<TEMP> ;T0,T1 GETS DSQRT(DABS(Y/2))
DMOVE T2,T0 ;T2,T3 ALSO
JUMPG T4,DONE ;DONE IF Y > 0
DMOVN T2,T2 ;NEGATE IMAG PART IF Y < 0
DONE: XMOVEI T4,@1(L) ;GET THE ADDRESS OF THE SECOND ARGUMENT
DMOVEM T0,(T4) ;SAVE REAL PART OF RESULT
DMOVEM T2,2(T4) ;SAVE IMAG PART OF RESULT
POP P,T5 ;RESTORE ACCUMULATORS
POP P,T4
POP P,T3
POP P,T2 ;RESTORE AC T2.
GOODBY (1) ;RETURN
UNDFL: SKIPGE XSAVE ;DOES RESULT HAVE REAL & IMAG PARTS SWITCHED?
$LCALL RPU,UFCONT
; LERR (LIB,%,<CDSQRT: real part underflow>,UFCONT)
$LCALL IPU
; LERR (LIB,%,<CDSQRT: imaginary part underflow>)
UFCONT: SETZB T2,T3 ;STORE 0 FOR RESULT OF UNDERFLOW
JRST SIGNS
ONE: DOUBLE 201400000000,000000000000 ;1.0D0
SQ2: DOUBLE 201552023631,237635714441 ;SQRT(2).
SEGMENT DATA
TEMP: DOUBLE 0,0 ;TEMPORARY STORAGE FOR DSQRT ARGS
XSAVE: DOUBLE 0,0
YSAVE: DOUBLE 0,0
END