Trailing-Edge
-
PDP-10 Archives
-
FORTRAN-10_Alpha_31-jul-86
-
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 ,1(4100)
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1986
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SUBTTL REVISION HISTORY
COMMENT \
***** Begin Revision History *****
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.
***** Begin Version 1A *****
3240 RVM 30-Mar-83
CDSQRT used an old standard for register names, this edit
updated the routine to be more modern. Also, memory
locations allocated but not used were deleted.
3246 The double precision sine and cosine functions destroy the
input variables. The result field gets written over the
input variable in some cases.
3247 10-33795 MRB 24-May-83
Fix in CDLOG.
3255 20-19901 MRB 23-Jan-84
The CDLOG routine is returning an incorrect value for the case
(0,0) returns (+inifinity,0) should be (-infinity,0).
***** End Revision History *****
***** Begin Version 10 *****
4100 MRB 16-Apr-84
Change name of CDLOG routine to CDLOG.
(include a symbol named CDLOG equivalenced to CDLOG.)
\
PRGEND
TITLE CDABS COMPLEX ABSOLUTE VALUE FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1986
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY 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, 1986
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;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, 1986
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY 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, 1986
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;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 LOG FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL Mike Boucher/MRB 16-Apr-84
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1984, 1986
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CDLOG
EXTERN CDLOG.
CDLOG=CDLOG.
PRGEND
TITLE CDLOG COMPLEX NATURAL LOG FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL Mary Payne 8-Sep-80
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1986
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;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,TT0 ;[3247]RECOVER X
DFDV T4,TT2 ;[3247]GET X/Y. NO OVERFLOW
JFCL UND1 ; UNDERFLOW IS POSSIBLE
DMOVN T4,T4 ;[3247]GET -X/Y
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,TT2 ;[3247]RECOVER Y
DFDV T4,TT0 ;[3247]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)>);[3255]
HRLOI T0,777777 ;[3255]REAL ANSWER IS NEGATIVE 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, 1986
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY 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, 1986
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY 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, 1986
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;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,@1(L) ;[3246]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, 1986
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY 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, 1986
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
; 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 >= 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 = 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, P1, P3, AND P4 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.
;[3240] Use the new names of the preserved registers, P1-P4,
;[3240] rather than the old names G1-G4.
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,P1 ;SAVE REGISTER P1.
PUSH P,P3 ;SAVE REGISTER P3
PUSH P,P4 ;SAVE REGISTER P4
DMOVEM T0,XSAVE ;SAVE X IN XSAVE
DMOVEM T2,YSAVE ;SAVE Y IN YSAVE
DMOVE P3,T2 ;DABS(Y) TO AC P3.
JUMPGE T2,YPOS ;IF Y IS NEGATIVE
DMOVN P3,T2 ;NEGATE IT
YPOS: HRRZI P1,1 ;SET P1 TO 1
CAML T4,P3 ;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,P3 ;PUT MAX(ABS(X),ABS(Y)) IN T4
EXCH T5,P4 ;AND MIN(ABS(X),ABS(Y)) IN P3,P4
SETZ P1, ;SET P1 TO 0
DWN: DFDV P3,T4 ;CALC S/L.
JFCL ;NO UNDERFLOW ERROR MESSAGE ALLOWED.
DFMP P3,P3 ;CALC (S/L)**2.
JFCL ;NO UNDERFLOW ERROR MESSAGE ALLOWED.
DFAD P3,ONE ;HAVE (1+(S/L)**2) IN AC P3.
DMOVEM P3,TEMP
FUNCT DSQRT.,<TEMP> ;CALC. THE DSQRT OF
;(1+(S/L)**2)=1+EPS.
JUMPE P1,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 P3,T0 ;SAVE DSQRT(1+EPS/2) IN AC P3.
DMOVEM T4,TEMP
FUNCT DSQRT.,<TEMP> ;CALC.
DFMP T0,P3 ;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 P3,XSAVE ;DABS(X) TO AC P3.
JUMPGE P3,XOK ;IF X IS NEGATIVE
DMOVN P3,P3 ;NEGATE IT
XOK: DFAD P3,T0 ;GET |X| + CDABS(Z)
FSC P3,1 ;MULTIPLY ARG BY 2
DMOVEM P3,TEMP ;SAVE FOR DSQRT
FUNCT DSQRT.,<TEMP> ;CALC. DSQRT OF 2 * (|X| + CDABS(Z))
DMOVE P3,T0 ;GET DSQRT((|X| + CDABS(Z))*2)
FSC T0,-1 ;GET DSQRT((|X| + CDABS(Z))/2)
DFDV T4,P3 ;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 RESULT
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,P4 ;RESTORE AC P4
POP P,P3 ;RESTORE AC P3
POP P,P1 ;RESTORE AC P1
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 IMAGINARY 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
SEGMENT DATA
TEMP: DOUBLE 0,0 ;TEMPORARY STORAGE FOR DSQRT ARGS
XSAVE: DOUBLE 0,0
YSAVE: DOUBLE 0,0
END