Google
 

Trailing-Edge - PDP-10 Archives - CFS_TSU04_19910205_1of1 - update/ftnosr/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, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

	SUBTTL	REVISION HISTORY

COMMENT \

***** Begin Revision History *****

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, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	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, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

;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, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	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, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;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, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	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, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;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, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	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, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	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, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;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, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	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, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

;  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