Trailing-Edge
-
PDP-10 Archives
-
bb-4157h-bm_fortran20_v10_16mt9
-
fortran-ots-debugger/mthcgx.mac
There are 9 other files named mthcgx.mac in the archive. Click here to see a list.
SEARCH MTHPRM
TV MTHCGX COMPLEX GFLOAT DOUBLE PRECISION ROUTINES ,1(4014)
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;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 *****
1404 EGM 6-Apr-81 --------
Separate CDX and CGX routines.
1405 DAW 6-Apr-81
Extended addressing support.
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 *****
3237 RVM 30-Mar-83
CGSQRT was not returning the primary square root (ie, the
result was not in the right half-plane). In addition,
there was an instruction which confused register T3 with
P3, and the underflow code didn't zero the result.
3255 20-19901 MRB 19-Nov-84
The CGLOG routine is returning an incorrect value for the case
(0,0) returns (+inifinity,0) should be (-infinity,0).
4014 MRB 18-Jul-84
The double precision sine and cosine functions destroy the
input variables. The result field gets written over the
input variable in some cases.
***** End Revision History *****
\
PRGEND
TITLE CGABS COMPLEX ABSOLUTE VALUE FUNCTION
; (DOUBLE PRECISION GFLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;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 CGABS
EXTERN CGABS.
CGABS=CGABS.
PRGEND
TITLE CGABS. COMPLEX ABSOLUTE VALUE FUNCTION
; (DOUBLE PRECISION GFLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;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.
;CGABS(Z) IS CALCULATED AS FOLLOWS
; LET Z = X + I*Y
; V = MAX(GABS(X),GABS(Y))
; W = MIN(GABS(X),GABS(Y))
; THEN
; CGABS = V*GSQRT(1.0 + (W/V)**2)
;THE RANGE OF DEFINITION FOR CGABS IS THE REPRESENTABLE REAL NUMBERS.
; AN OVERFLOW CAN OCCUR, IN WHICH CASE CGABS = + MACHINE INFINITY.
;REQUIRED (CALLED) ROUTINES: GSQRT
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; XMOVEI L,ARG
; PUSHJ P,CGABS
;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 (CGABS,.) ;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| IS GREATER, NO EXCHANGE NECESSARY
EXCH T2,T0 ;EXCHANGE ABS(X) AND ABS(Y)
EXCH T3,T1 ;
LT: JUMPE T2,RET ;Z = 0, HENCE CGABS = O
GFDV T0,T2 ;DIVIDE SMALLER BY LARGER; NO OVERFLOW
JFCL ANS ; RATIO NEGLIGIBLE IF UNDERFLOW.
CAMG T0,TWOM30 ;IF RATIO .LE. 2**(-30)
JRST ANS ; RATIO IS NEGLIBLE.
GFMP T0,T0 ;**2
GFAD T0,ONE ;+1.0
DMOVEM T0,TEMP
FUNCT GSQRT.,<TEMP> ;SQUARE ROOT IS IN AC T0
GFMP T0,T2 ;*V
JFCL OVFL ;NO UNDERFLOW, GET MESSAGE ON OVERFLOW
RET: POP P,T3
POP P,T2
GOODBY (1) ;RETURN
OVFL: $LCALL ATI
;LERR (LIB,%,<CGABS: CGABS(arg) too large; result=+infinity>)
RET1: POP P,T3
POP P,T2
GOODBY (1) ;RETURN
ANS: DMOVE T0,T2 ;ANSWER = ABS(LARGER) TO T0
POP P,T3
POP P,T2
GOODBY (1) ;RETURN
ONE: DOUBLE 200140000000,000000000000
TWOM30: 174340000000 ;2**(-30)
SEGMENT DATA
TEMP: DOUBLE 0,0 ;TEMPORARY STORAGE USED FOR SQRT CALL
PRGEND
TITLE CGEXP COMPLEX EXPONENTIAL FUNCTION
; (DOUBLE PRECISION GFLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;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 CGEXP
EXTERN CGEXP.
CGEXP=CGEXP.
PRGEND
TITLE CGEXP. COMPLEX EXPONENTIAL FUNCTION
; (DOUBLE PRECISION GFLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;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.
;CGEXP(Z), WHERE Z = X + I*Y, IS CALCULATED AS FOLLOWS
; CGEXP(Z) = GEXP(X)*(GCOS(Y)+I*GSIN(Y))
;THE RANGE OF DEFINITION FOR CGEXP IS AS FOLLOWS
;FOR Z = X+I*Y, IF GABS(Y) .GT. 6746518850.429 THE RESULT IS SET TO
; (0.0,0.0) AND AN ERROR MESSAGE IS RETURNED.
;FOR X.LT. -710.475860073943942 THE RESULT IS SET TO (0.0,0.0) AND AN ERROR
; MESSAGE IS RETURNED.
;FOR X .GT. 709.089565;
; IF Y = 0.0, THE RESULT IS SET TO (+INFINITY,0.0) AND AN
; ERROR MESSAGE IS RETURNED.
; IF 709.089565 < X < 1418.179131425648102, 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. 709.089565, THE GABS(DREAL(RESULT)) IS SET TO
; +INFINITY AND GABS(DIMAG(RESULT)) IS SET TO +INFINITY AND
; AN ERROR MESSAGE IS RETURNED.
;REQUIRED (CALLED) ROUTINES: GEXP,GSIN,GCOS
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; XMOVEI L,ARG
; PUSHJ P,CGEXP
;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 (CGEXP.,CGEXP) ;ENTRY TO CGEXP ROUTINE
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,YSAVQ ;SAVE Y
JUMPGE T2,QPOS
DMOVN T2,T2
QPOS: CAMGE T2,QYMAXH ;IF HI OF Y < QYMAXH
JRST QYOK ; GO TO QYOK
CAME T2,QYMAXH ;IF HI OF Y > QYMAXH
JRST ERQ1 ; GO TO ERQ1
CAMG T3,QYMAXL ;IF LO OF Y .LE. QYMAXL
JRST QYOK ; GO TO QYOK
ERQ1: $LCALL AIZ
;LERR (LIB,%,<CGEXP: 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)
QYOK: PUSH P,T5 ;SAVE ANOTHER REGISTER
DMOVEM T0,QXSAVE ;SAVE X
CAMLE T0,QNEGH ;IF HI OF X > QNEGH
JRST QXOK ; GO TO QXOK
CAME T0,QNEGH ;IF HI OF X < QNEGH
JRST QER ; GO TO QER
CAML T1,QNEGL ;IF LO OF X .GE. QNEGL
JRST QXOK ; GO TO QXOK
QER: SETZ T0, ;DREAL PART OF SOLUTION IS ZERO
SETZ T1,
SETZ T2, ;IMAG PART OF SOLUTION IS ZERO
SETZ T3,
JRST QYCHK ;CHECK Y
QXOK: FUNCT GCOS.,<YSAVQ> ;QCOSY=GCOS(Y)
DMOVEM T0,QCOSY
FUNCT GSIN.,<YSAVQ> ;GSINY=GSIN(Y)
DMOVEM T0,QSINY
DMOVE T4,QXSAVE ;T4=X
CAMG T4,QMAXH ;IF X IS NOT TOO LARGE
JRST QALG ;GO TO QALG
SKIPN YSAVQ ;ELSE, IF Y=0
JRST QMSTK ;GO TO QMSTK
CAMGE T4,QTBIGH ;IF HI OF S < QTBIGH
JRST QSOK ; GO TO QSOK
CAME T4,QTBIGH ;IF HI OF S > QTBIGH
JRST QSTIM ; GO TO QSTIM
CAMLE T5,QTBIGL ;IF LO OF S > QTBIGL
JRST QSTIM ; GO TO QSTIM
QSOK: EXTEND T4,[GFSC -1] ;ELSE, S=X/2
DMOVEM T4,QEXPX
FUNCT GEXP.,<QEXPX> ;T=EXP(S)
DMOVE T2,T0
GFMP T2,QSINY ;V=T*QSINY
JUMPGE T2,QVPOS ;IF V IS NEGATIVE
DMOVN T2,T2 ;NEGATE IT
QVPOS: HRLOI T4,377777 ;G3=XMAX
SETO T5,
GFDV T4,T0 ;Q=XMAX/T
DMOVEM T4,QQH ;SAVE Q
CAMLE T2,T4 ;IF V .GT. Q
JRST QSTIM ;THEN GO TO QSTIM
CAME T2,QQH ;IF V .LT. Q
JRST QGETI ;THEN GO TO QGETI
CAMG T3,QQL ;IF V .LE. Q
JRST QGETI ;THEN GO TO QGETI
QSTIM: HRLOI T2,377777 ;ELSE, SET IMAG SOLUTION TO XMAX
SETO T3,
$LCALL RTI
;LERR (LIB,%,<CGEXP: DREAL(arg) too large; DIMAG(result)=+infinity>)
JRST QD2
QGETI: GFMP T2,T0 ;IRES = V*T
QD2: SKIPGE QSINY ;IF SINY IS LESS THAN 0.0
DMOVN T2,T2 ;THEN NEGATE IRES
MOVE T4,QXSAVE
CAMGE T4,QTBIGH ;IF HI OF S < QTBIGH
JRST QOKS ; GO TO QOKS
CAME T4,QTBIGH ;IF HI OF S > QTBIGH
JRST QMSTK ; GO TO QMSTK
CAMLE T5,QTBIGL ;IF LO OF S > QTBIGL
JRST QMSTK ; GO TO QMSTK
QOKS: DMOVE T4,T0
GFMP T0,QCOSY ;V = T*QCOSY
JUMPGE T0,QVOK ;IF V IS NEGATIVE
DMOVN T0,T0 ;NEGATE IT
QVOK: CAMLE T0,QQH ;IF V .GT. Q
JRST QMSTK ;THEN GO TO QMSTK
CAME T0,QQH ;IF V .LT. Q
JRST QQLAB ;THEN GO TO QQLAB
CAMG T1,QQL ;IF V IS .LE. Q
JRST QQLAB ;THEN GO TO QQLAB
QMSTK: HRLOI T0,377777 ;RRES=XMAX
SETO T1,
$LCALL RTR
;LERR (LIB,%,<CGEXP: DREAL(arg) too large; DREAL(result)=+infinity>)
JRST DQD2
QQLAB: GFMP T0,T4 ;RRES = V*T
DQD2: SKIPGE QCOSY ;IF QCOSY .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
QALG: FUNCT GEXP.,<QXSAVE> ;QEXPX = DEXP(X)
DMOVE T2,T0
GFMP T2,QSINY ;IRES = QEXPX*QSINY
JFCL
GFMP T0,QCOSY ;RRES = QEXPX*QCOSY
JFCL
JUMPN T2,QRCHK ;IF IRES .NE. 0.0
;THEN GO CHECK RRES
QYCHK: SKIPN YSAVQ ;IF Y .EQ. 0
JRST QRCHK ;THEN GO CHECK RRES
$LCALL IPU
;LERR (LIB,%,<CGEXP: underflow has occurred; DIMAG(result)=0.0>)
QRCHK: JUMPN T0,QRET ;IF R = 0.0
$LCALL RPU
;LERR (LIB,%,<CGEXP: underflow has occurred; DREAL(result)=0.0>)
QRET: 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
QYMAXH: 203762207732 ;HIGH ORDER PART OF YMAX
QYMAXL: 242102643021 ;LOW ORDER PART OF YMAX
QTBIGH: 201354242673 ;1418.179131425648102
QTBIGL: 161647554056
QNEGH: 576523460613 ;-710.475860073943942
QNEGL: 202140360224
QMAXH: 201254242673 ;709.089565
SEGMENT DATA
QQH: 0
QQL: 0
QSINY: DOUBLE 0,0
QCOSY: DOUBLE 0,0
QXSAVE: DOUBLE 0,0
QEXPX: DOUBLE 0,0
YSAVQ: DOUBLE 0,0
PRGEND
TITLE CGLOG COMPLEX NATURAL LOG FUNCTION
; (DOUBLE PRECISION GFLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;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 CGLOG
EXTERN CGLOG.
CGLOG=CGLOG.
PRGEND
TITLE CGLOG COMPLEX NATURAL LOG FUNCTION
; (DOUBLE PRECISION FLOATING POINT)
SUBTTL Mary Payne 8-Sep-80
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;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 (CGLOG.,CGLOG) ;ENTRY TO CGLOG ROUTINE
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,P30 ;IF DIFFERENCE > 30
JRST BIG ; |Y/X| IS LARGE
CAMGE T2,M30 ;IF DIFFERENCE < -30
JRST SMALL ; |Y/X| IS SMALL
FUNCT GATN2.,<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 TO T4
EXCH T5,T3 ; SMALLER TO T2
NOXCH: MOVE T1,T4 ;HI OF LARGER TO T1
LSH T1,-30 ;ISOLATE ITS BIASED EXPONENT
SUBI T1,2000 ; AND UNBIAS IT
JUMPLE T1,LE ;IF UNBIASED EXPONENT POSITIVE,
SUBI T1,1 ; DECREMENT IT
LE: MOVN T1,T1 ;NEGATE SCALE INDEX
EXTEND T4,[GFSC (T1)] ;SCALE LARGER TO [1/2,2)
EXTEND T2,[GFSC (T1)] ; AND SMALLER BY SAME FACTOR
JFCL ;NO OVERFLOW, UNDERFLOW WON'T MATTER
LSH T1,1 ;MULTIPLY NEGATED SCALE INDEX BY 2
DMOVEM T4,BIGGER ;SAVE SCALED LARGER OF |X| AND |Y|
GFMP T4,T4 ;SQUARE LARGER. NO EXCEPTIONS
GFMP T2,T2 ;SQUARE SMALLER. NO OVERFLOW AND
JFCL ; UNDERFLOW WON'T MATTER
GFAD 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
DMOVE T4,BIGGER ;RECOVER LARGER
DMOVE T0,T4 ; AND COPY INTO T0
GFSB T0,ONE ;GET A - 1
GFAD T4,ONE ;AND A + 1
GFMP T4,T0 ;AND THEIR PRODUCT
GFAD T2,T4 ;T2 NOW HAS A**2 + B**2 - 1
DMOVE T4,T2 ;COPY IT INTO T4
GFAD T2,TWO ;GET DENOMINATOR OF Z IN T2
EXTEND T4,[GFSC 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,-30 ;ISOLATE ITS BIASED EXPONENT
SUBI T2,2000 ;UNBIAS THE EXPONENT
MOVN T2,T2 ; AND NEGATE IT
EXTEND T4,[GFSC (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
EXTEND T4,[GFSC 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 ;STORE FINAL INDEX FOR FUTURE USE
DMOVE T2,T4 ;GET COPY OF SCALED (A**2 + B**2)
GFSB T4,ONE ;SUBTRACT ONE AND MULTIPLY BY TWO
EXTEND T4,[GFSC 1] ; TO GET NUMERATOR OF Z IN T4
GFAD T2,ONE ;GET DENOMINATOR OF Z IN T2
;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: GFDV T4,T2 ;Z = ZNUM / ZDEN
DMOVEM T4,SAVEZ ;SAVE A COPY OF Z
GFMP T4,T4 ;W = Z**2
DMOVE T0,T4 ; AND MAKE COPY
GFAD T4,B2 ;FORM B(W) = W + B2
GFMP T4,T0 ; * W
GFAD T4,B1 ; + B1
GFMP T4,T0 ; * W
GFAD T4,B0 ; + B0
DMOVE T2,T0 ;MAKE ANOTHER COPY OF 2
GFMP T0,A2 ;FORM A(W) = W * A2
GFAD T0,A1 ; + A1
GFMP T0,T2 ; * W
GFAD T0,A0 ; + A0
GFDV T0,T4 ;R(Z) = A(W) / B(W)
GFMP T0,SAVEZ ; * Z
GFMP 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.
EXTEND T2,[GFLTR NN] ;RECOVER SCALE INDEX
JUMPN T2,REST ;IF NN = 0,
GFAD 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: DMOVE T4,T2 ;GET COPY IN T4
GFMP T4,C2 ;GET NN * LN(2)LO
GFMP T2,C1 ;GET NN * LN(2)HI
GFAD T0,T4 ;R(Z) + NN * LN(2)LO
GFAD T0,SAVEZ ; + Z
GFAD T0,T2 ; + NN * LN(2)HI
DONE: EXTEND T0,[GFSC -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 GLOG.,<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|
GFDV T4,TT2 ;GET |X/Y|. NO OVERFLOW
JFCL UND1 ; UNDERFLOW IS POSSIBLE
JRST SMERG ;MERGE WITH SMALL
SMALL: FUNCT GLOG.,<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
GFDV T4,TT4 ;GET Y/X. NO OVERFLOW
JFCL UND2 ; UNDERFLOW IS POSSIBLE
SMERG: GFAD T2,T4 ;SMALL CORRECTION FOR IMAG
GFMP T4,T4 ;GET (Y/X)**2. NO OVERFLOW
JFCL UND1 ; INDERFLOW IS POSSIBLE
EXTEND T4,[GFSC -1] ;(1/2)*(Y/X)**2. NO OVERFLOW
JFCL UND1 ; UNDERFLOW IS POSSIBLE
GFAD 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
$LCALL IPU
;LERR (LIB,%,<CDLOG: Imaginary part underflow>)
UND1: JUMPN T0,STOR ;NO MESSAGE IF LN NOT ZERO
$LCALL RPU
;LERR (LIB,%,<CDLOG: Real part underflow>)
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 GLOG.,<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 GLOG.,<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 200140000000,000000000000 ;1.0
TWO: EXP 200240000000,000000000000 ;2.0
MASK: EXP 377700000000 ;EXPONENT MASK
P30: EXP 003600000000 ;HIGH 12 BITS = 36 OCTAL
M30: EXP 774200000000 ;HIGH 12 BITS = -36 OCTAL
CPI: EXP 200262207732,242102643022 ;PI
DPIO2: EXP 200162207732,242102643022 ;PI / 2
RT2O2: EXP 200055202363 ;SQRT(2) / 2
RT2: EXP 200155202363 ;SQRT(2)
A0: EXP 577037740007,152304514557 ;-.641249434237455811D2
A1: EXP 200540611121,000552775450 ;.163839435630215342D2
A2: EXP 577715357522,145224132710 ;-.789561128874912573D0
B0: EXP 576517720013,037446761043 ;-.769499321084948798D3
B1: EXP 201147002037,320522317572 ;.312032220919245328D3
B2: EXP 577134251775,244603076112 ;-.356679777390346462D2
B3: EXP 577640000000,000000000000 ;-1
C1: EXP 200054300000,000000000000 ;LN(2)HI
C2: EXP 601310277575,034757152745 ;LN(2)LO
;DATA
SEGMENT DATA
TT0: BLOCK 2 ;FOR GATN2 AND GLOG CALLS
TT2: BLOCK 2 ;FOR GATN2 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 CGSIN COMPLEX SINE FUNCTION
; (DOUBLE PRECISION GFLOATING)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;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 CGSIN
EXTERN CGSIN.
CGSIN=CGSIN.
PRGEND
TITLE CGCOS COMPLEX COSINE FUNCTION
; (DOUBLE PRECISION GFLOATING)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;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 CGCOS
EXTERN CGCOS.
CGCOS=CGCOS.
PRGEND
TITLE CGSIN. COMPLEX SINE FUNCTION
; (DOUBLE PRECISION GFLOATING)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;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.
;CGSIN(Z) AND CGCOS(Z), WHERE Z = X + I*Y, ARE CALCULATED AS FOLLOWS
; CGSIN(Z) =GSIN(X)*GCOSH(Y) + I*GCOS(X)*GSINH(Y)
; CGCOS(Z) =GCOS(X)*GCOSH(Y) - I*GSIN(X)*GSINH(Y)
; THE FUNCTIONS GSINH AND GCOSH ARE CODED IN LINE.
;THE RANGE OF DEFINITION FOR GDSIN AND GDCOS IS AS FOLLOWS
;FOR
; Z = X+I*Y. IF GABS(X) > ((2**29)*PI - PI/2) THE RESULT IS SET TO 0.0 AND
; AN ERROR MESSAGE IS RETURNED.
; THE REAL PART OF THE RESULT CANNOT UNDERFLOW BECAUSE
; COSH(Y) IS NEVER LESS THAN 1.
; THE IMAGINARY PART OF THE RESULT MAY UNDERFLOW; IF THIS HAPPENS THE
; IMAGINARY PART IS SET TO 0 AND A MESSAGE IS PROVIDED.
;FOR GABS(Y) > 88.029692 CALCULATIONS PROCEED AS FOLLOWS:
; FOR THE REAL PART OF THE RESULT
; LET T = GABS(GSIN(X)), THEN
; FOR T = 0.0 THE REAL PART OF THE RESULT IS SET TO 0.0
; FOR GLOG(T) + GABS(Y) - GLOG(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 = GABS(GCOS(X)), THEN
; FOR T = 0.0 THE IMAGINARY PART OF THE RESULT IS SET TO 0.0
; FOR GLOG(T) + GABS(Y) - GLOG(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: GSIN,GCOS,GEXP,GLOG
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; XMOVEI L,ARG
; PUSHJ P,CGSIN
; OR
; PUSHJ P,CGCOS
;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 (CGCOS.,CGCOS) ;ENTRY TO CGCOS ROUTINE
PUSH P,T5
MOVEI T5,1 ;SET FLAG TO 1 FOR CGCOS
JRST GSTRT ;GO TO CGSIN ROUTINE
HELLO (CGSIN.,CGSIN) ;ENTRY TO CGSIN ROUTINE
PUSH P,T5
GBEG: SETZI T5, ;SET FLAG TO 0 FOR CGSIN
GSTRT: 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 T3,T4
JUMPGE T3,QXPOS ;IF REAL PART OF ARG < 0
DMOVN T3,T3 ;REAL PART=ABS(REAL PART)
QXPOS: CAMLE T3,XMAXQH ;IF ABS(X) IS TOO LARGE
JRST QERR1 ;GO TO QERR1
CAME T3,XMAXQH ;IF ABS(X) IS LESS THAN XMAXQH
JRST QXOK ;GO TO QXOK
CAMG T4,XMAXQL ;IF ABS(X) IS NOT TOO LARGE
JRST QXOK ;GO TO QXOK
QERR1: $LCALL ARZ
;LERR (LIB,%,<CGSIN or CGCOS: GABS(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
QXOK: DMOVE T2,2(T2) ;Y = IMAG(Z)
DMOVEM T2,YSAVQ ;SAVE Y
DMOVEM T0,ARGSVQ
FUNCT GSIN.,<ARGSVQ> ;CALCULATE GSIN(X)
DMOVEM T0,SXQ ;SXQ = GSIN(X)
FUNCT GCOS.,<ARGSVQ> ;CALCULATE GCOS(X)
DMOVEM T0,CXQ ;CXQ = GCOS(X)
JUMPE T5,NQSWT ;IF THIS IS THE CGSIN ROUTINE
;THEN GO TO NQSWT
DMOVN T2,SXQ ;OTHERWISE, PUT -GSIN(X)
EXCH T2,CXQ ;IN CXQ, AND COS(X)
EXCH T3,CXQL
DMOVEM T2,SXQ ;IN SXQ.
NQSWT: DMOVE T2,YSAVQ ;GET A COPY OF Y
JUMPG T2,QYPOS ;IF Y IS NEGATIVE
DMOVN T2,T2 ;NEGATE IT
QYPOS: CAMGE T2,XQHI ;IF HI OF AY < XQHI
JRST YSMALQ ; GO TO YSMALQ
CAME T2,XQHI ;IF HI OF AY > XQHI
JRST OVFLQ ; GO TO OVFLQ
CAMLE T3,XQLO ;IF LO OF AY > XQLO
JRST OVFLQ ; GO TO OVFLQ
YSMALQ: CAMG T2,TWQM29 ;IF ABS(Y) .LE. 2**(-29)
JRST QEASY ; GO TO QEASY CODE.
CAMGE T2,ONEQ ;IS ABS(Y) .GE. 1?
JRST LSTH7 ; BRANCH IF NOT
DMOVEM T2,ARGSVQ ;GET COPY OF ABS(Y)
FUNCT GEXP.,<ARGSVQ> ;GET EXP(ABS(Y))
DMOVE T4,ONEQ ;GET EXP(-ABS(Y)) BY
GFDV T4,T0 ; RECIPROCATION.
GFSB T0,T4 ;T0 GETS
EXTEND T0,[GFSC -1] ; SINH(ABS(Y))
JRST CCCQ ;TO REST OF CALCULATION
LSTH7: DMOVNM T2,ARGSVQ ; SO AS TO GET
FUNCT GEXP.,<ARGSVQ> ; EXP(-ABS(Y))
DMOVE T4,T0 ;SAVE 1/T IN T4
GFMP T2,T2 ;SS = AY**2
DMOVEM T2,TEMPQ ;SAVE A COPY OF SS
GFAD T2,QQ2 ;XDEN = QQ2 +SS
GFMP T2,TEMPQ ;*SS
GFAD T2,QQ1 ;+QQ1
GFMP T2,TEMPQ ;*SS
GFAD T2,QQ0 ;+QQ0
DMOVE T0,TEMPQ ;SAVE A COPY OF SS
GFMP T0,RPQ3 ;XNUM = RPQ3*SS
GFAD T0,RPQQ2 ;+RPQQ2
GFMP T0,TEMPQ ;*SS
GFAD T0,RPQQ1 ;+RPQQ1
GFMP T0,TEMPQ ;*SS
GFAD T0,RPQQ0 ;+RPQQ0
GFDV T0,T2 ;XNUM/XDEN
DMOVN T2,ARGSVQ ;GET AY BACK
GFMP T0,TEMPQ ;*SS
GFMP T0,ARGSVQ ;*(-AY)
GFAD T0,T2 ;+AY
CCCQ: DMOVE T2,T0
GFAD T0,T4 ;CC = SS+(1/T)
DMOVE T4,YSAVQ ;RESTORE Y
JUMPGE T4,YOKQ ;IF Y IS LESS THAN 0.0
DMOVN T2,T2 ;THEN NEGATE SS
YOKQ: GFMP T0,SXQ ;GET REAL PART OF RESULT
GFMP T2,CXQ ;GET IMAG PART OF RESULT
JFCL QIMUND ;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
OVFLQ: DMOVE T0,SXQ ;T=SXQ
JUMPE T0,GOTRQ ;IF T IS EQUAL TO 0.
JUMPG T0,SQXPOS ;IF SXQ IS NEGATIVE
DMOVN T0,T0 ;NEGATE SXQ
SQXPOS: DMOVEM T0,ARGSVQ ;OTHERWISE
FUNCT GLOG.,<ARGSVQ> ;CALCULATE LOG(T)
GFAD T0,T2 ;T = LOG(T)+AY
GFAD T0,LN2Q ;T = T-LOG(2.0)
CAMGE T0,XQHI ;IF HI OF T < XQHI
JRST QONTIN ; GO TO QONTIN
CAME T0,XQHI ;IF HI OF T > XQHI
JRST QERR2 ; GO TO QERR2
CAMG T1,XQLO ;IF LO OF T .LE. XQLO
JRST QONTIN ; GO TO QONTIN
QERR2: $LCALL AIR
;LERR (LIB,%,<CGSIN or CGCOS: GABS(DIMAG(arg)) too large; DREAL(result) = infinity>)
HRLOI T0,377777 ;REAL PART OF RESULT SET TO INFINITY
SETO T1, ;SET SECOND WORD
JRST SKPQQ2 ;SKIP NEXT 2 INSTRUCTIONS
QONTIN: DMOVEM T0,ARGSVQ
FUNCT GEXP.,<ARGSVQ> ;RRES = EXP(T)
SKPQQ2: SKIPGE SXQ ;IF SXQ IS LESS THAN 0.
DMOVN T0,T0 ;THEN NEGATE RRES
GOTRQ: SKIPN T4,CXQ ;IF CXQ = 0,
JRST IMAG7 ;THEN PREPARE TO RETURN
MOVE T5,CXQL ;OTHERWISE GET LOWER WORD
DMOVEM T0,SXQ ;SAVE RRES
JUMPGE T4,TPOSQ ;IF T IS NEGATIVE
DMOVN T4,T4 ;THEN GET ABS(T)
TPOSQ: DMOVEM T4,ARGSVQ
FUNCT GLOG.,<ARGSVQ> ;CALCULATE T=LOG(T)
GFAD T2,T0 ;T = T+AY
GFAD T2,LN2Q ;T = T-LOG(2)
CAMGE T2,XQHI ;IF HI OF T < XQHI
JRST QONTN2 ; GO TO QONTN2
CAME T2,XQHI ;IF HI OF T > XQHI
JRST QERR3 ;GO TO QERR3
CAMG T3,XQLO ;IF LO OF T .LE. XQLO
JRST QONTN2 ; GO TO QONTN2
QERR3: $LCALL AII
;LERR (LIB,%,<CGSIN or CGCOS: GABS(DIMAG(arg)) too large; DIMAG(result) = infinity>)
HRLOI T2,377777 ;SET IRES TO INFINITY
SETO T3,
JRST SKPQ3 ;SKIP NEXT 3 INSTRUCTIONS
QONTN2: DMOVEM T2,ARGSVQ
FUNCT GEXP.,<ARGSVQ> ;IRES = EXP(T)
DMOVE T2,T0
SKPQ3: DMOVE T0,SXQ
MOVE T4,YSAVQ
XOR T4,CXQ ;T4 HAS THE SIGN OF CXQ*Y
JUMPGE T4,RETQ ;IF T4 IS LESS THAN 0.
DMOVN T2,T2 ;THEN NEGATE IRES
JRST RETQ ;RETURN
IMAG7: SETZ T2, ;SET IMAGINARY PART
SETZ T3, ; OF RESULT TO ZERO.
JRST RETQ ;RETURN
QEASY: DMOVE T2,YSAVQ ;GET SIGNED VALUE OF Y
GFMP T2,CXQ ;SINH(Y)*COS(X) = Y*COS(X)
JFCL QIMUND ; UNDERFLOW POSSIBLE
DMOVE T0,SXQ ;COSH(Y)*SIN(X) = SIN(X).
RETQ: XMOVEI T4,@1(L) ;[4014]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)
QIMUND: $LCALL IPU,RETQ
;LERR (LIB,%,<CDSIN or CDCOS: imaginary part underflow>,RETQ)
TWQM29: 174440000000 ;2**(-29)
XQHI: 201254242673 ;709.089565712824051
XQLO: 161647554056 ;
XMAXQH: 203762207732 ;HIGH ORDER 2**29 * PI
XMAXQL: 236772245275 ;LOW ORDER 2**29 * PI
ONEQ: DOUBLE 200140000000,000000000000
LN2Q: DOUBLE 577723506750,010134302063 ;-(NATURAL LOG OF 2.0)
RPQQ0: DOUBLE 202352744232,262463203065 ;.35181283430177117881D+6
RPQQ1: DOUBLE 201655127025,264501221757 ;.11563521196851768270D+5
RPQQ2: DOUBLE 201050741013,034133711232 ;.16375798202630751372D+3
RPQ3: DOUBLE 200062423475,303374403264 ;.78966127417357099479D+0
QQ0: DOUBLE 575137624613,372031435521 ;-.21108770058106271242D+7
QQ1: DOUBLE 202043241271,035545730675 ;.36162723109421836460D+5
QQ2: DOUBLE 576635220743,361550001577 ;-.27773523119650701667D+3
SEGMENT DATA
CXQ: 0
CXQL: 0
SXQ: DOUBLE 0,0
YSAVQ: DOUBLE 0,0
ARGSVQ: DOUBLE 0,0
TEMPQ: DOUBLE 0,0
PRGEND
TITLE CGSQRT COMPLEX SQUARE ROOT FUNCTION
; (DOUBLE PRECISION GFLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;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 CGSQRT
EXTERN CGSQT.
CGSQRT=CGSQT.
PRGEND
TITLE CGSQT. COMPLEX SQUARE ROOT FUNCTION
; (DOUBLE PRECISION GFLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;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 CGSQRT(Z) = U + I*V IS DEFINED AS FOLLOWS
; IF X >= 0.0, THEN
; U = GSQRT((GABS(X)+CGABS(Z))/2.0)
; V = Y/(2.0*U)
; IF X < 0.0 AND Y >= 0.0, THEN
; U = Y/(2.0*V)
; V = GSQRT((GABS(X)+CGABS(Z))/2.0)
; IF X AND Y ARE BOTH < 0.0, THEN
; U = Y/(2.0*V)
; V = GSQRT((GABS(X)+CGABS(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: GSQRT
;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,CGSQRT
;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 (CGSQT.,CGSQRT) ;ENTRY TO COMPLEX SQUARE ROOT ROUTINE.
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: GFDV P3,T4 ;CALC S/L.
JFCL ;NO UNDERFLOW ERROR MESSAGE ALLOWED.
GFMP P3,P3 ;CALC (S/L)**2.
JFCL ;NO UNDERFLOW ERROR MESSAGE ALLOWED.
GFAD P3,ONE ;HAVE (1+(S/L)**2) IN AC P3.
DMOVEM P3,TEMP
FUNCT GSQRT.,<TEMP> ;CALC. THE GSQRT OF
;(1+(S/L)**2)=1+EPS.
JUMPE P1,YGTX ;GO TO YGTX IF DABS(Y) > DABS(X).
XGTY: GFAD T0,ONE ;CALC. 2 + EPS.
EXTEND T0,[GFSC -1] ;CALC. 1+EPS/2.
DMOVEM T0,TEMP ;PUT IN TEMP FOR
FUNCT GSQRT.,<TEMP> ;CALL TO GSQRT
DMOVE P3,T0 ;SAVE GSQRT(1+EPS/2) IN AC P3.
DMOVEM T4,TEMP
FUNCT GSQRT.,<TEMP> ;CALC.
GFMP T0,P3 ;CALC. GSQRT(DABS(X)*(1+EPS/2)).
JRST OUT1 ;GO TO REST OF CALC.
YGTX: CAIGE T4,200140 ;IF DABS(Y)>1, GO TO POSSIBLE OVERFLOW CASE.
JRST YXSMAL ;ELSE, GO TO POSSIBLE UNDERFLOW.
EXTEND T0,[GFSC -3] ;CALC. (1+EPS)/8.
GFMP T4,T0 ;CALC. DABS(Y)*(1+EPS)/8
DMOVE T0,XSAVE ;DABS(X) TO AC T0
JUMPGE T0,NXT
DMOVN T0,T0
NXT: EXTEND T0,[GFSC -3] ;CALC. DABS(X)/8
JFCL ;SUPPRESS UNDERFLOW ERROR MESSAGE.
GFAD T0,T4 ;CALC.(DABS(X)/8)+(DABS(Y)*(1+EPS)/8).
DMOVEM T0,TEMP
FUNCT GSQRT.,<TEMP> ;CALC.
EXTEND T0,[GFSC 1] ;GET DSQRT(|X|/2 + |Y|*(1+EPS)/2)
OUT1: JUMPGE T2,POSY
DMOVN T2,T2
POSY: GFDV T2,T0 ;DIVIDE DABS(Y)/2 BY THE
JFCL UNDFL
EXTEND T2,[GFSC -1] ;SQRT TERM.
JFCL UNDFL
JRST SIGNS ;GO TO REST OF CALC.
YXSMAL: GFMP 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: GFAD P3,T0 ;GET |X| + CDABS(Z)
EXTEND P3,[GFSC 1] ;[3237] MULTIPLY ARG BY 2
DMOVEM P3,TEMP ;SAVE FOR DSQRT
FUNCT GSQRT.,<TEMP> ;CALC. GSQRT OF 2 * (|X| + CDABS(Z))
DMOVE P3,T0 ;GET DSQRT((|X| + CDABS(Z))*2)
EXTEND T0,[GFSC -1] ;GET DSQRT((|X| + CDABS(Z))/2)
GFDV 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 ;[3237] 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 GSQRT
DMOVN T0,T0
POSX: DMOVEM T0,TEMP
FUNCT GSQRT.,<TEMP> ;T0,T1 GET GSQRT(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: EXTEND T2,[GFSC -1] ;GET |Y|/2
DMOVEM T2,TEMP ; AND STORE IT FOR GSQRT
FUNCT GSQRT.,<TEMP> ;T0,T1 GETS GSQRT(DABS(Y/2))
DMOVE T2,T0 ;T2,T3 ALSO
JUMPG T4,DONE ;DONE IF Y > 0
DMOVN T2,T2 ;[3237] 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 ;[3237] Go to UFCONT to do fixup when through
; LERR (LIB,%,<CDSQRT: real part underflow>,UFCONT)
$LCALL IPU ;[3237] Go to UFCONT to do fixup when through
; LERR (LIB,%,<CDSQRT: imaginary part underflow>)
UFCONT: SETZB T2,T3 ;[3237] STORE 0 FOR RESULT OF UNDERFLOW
JRST SIGNS ;[3237]
ONE: DOUBLE 200140000000,000000000000 ;1.0D0
SEGMENT DATA
TEMP: DOUBLE 0,0 ;TEMPORARY STORAGE FOR GSQRT ARGS
XSAVE: DOUBLE 0,0
YSAVE: DOUBLE 0,0
END