Google
 

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