Google
 

Trailing-Edge - PDP-10 Archives - BB-4157E-BM - fortran-ots-debugger/forcdx.mac
There is 1 other file named forcdx.mac in the archive. Click here to see a list.
	SEARCH	FORPRM
	TV	FORCDX	COMPLEX DOUBLE PRECISION ROUTINES ,6(2031)

;COPYRIGHT (C) 1979, 1981  BY  DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.

;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.

***** End Revision History *****

\

	PRGEND
TITLE	CDABS	COMPLEX ABSOLUTE VALUE FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
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) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
	TWOSEG	400000
	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
RET:	POP	P,T3
	POP	P,T2
	GOODBY	(1)		;RETURN

OVFL:	CAIN	T0,0		;OVERFLOW? - NOTE ERROR.  OTHERWISE,
	  JRST	RET		;UNDERFLOW, RETURN
	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)

	RELOC			;DATA
TEMP:	DOUBLE	0,0		;TEMPORARY STORAGE USED FOR SQRT CALL
	RELOC
	PRGEND
TITLE	CDEXP	COMPLEX EXPONENTIAL FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
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) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
	FSRCH
	TWOSEG	400000
	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:	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>)
	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>)
	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>)

RCHK:	JUMPN	T0,RET		;IF R = 0.0
	  LERR	(LIB,%,<entry CDEXP; underflow has occurred; DREAL(result)=0.0>)

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

	RELOC			;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
	RELOC
	PRGEND
TITLE	CDLOG	COMPLEX NATURAL LOG FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL  Mary Payne	8-Sep-80

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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  FORPRM
	TWOSEG  400000
	SALL

	HELLO	(CDLOG)		;ENTRY TO CDLOG ROUTINE
IFN FTGFL,<
	LDB	T1,[POINT 4,0(L),12] ;GET FIRST ARG TYPE
	CAIN	T1,TP%DPX	;G FLOATING?
	  JRST  CGLOG##		;YES
> ;END FTGFL
	PUSH	P,T2		;SAVE ACCUMULATORS
	PUSH	P,T3		
	PUSH	P,T4
	XMOVEI  T2,@(L)		;GET ADDRESS OF ARG
	DMOVE	T0,(T2)		;OBTAIN REAL PART OF ARGUMENT
	DMOVE	T2,2(T2)	;OBTAIN IMAG PART OF ARGUMENT

	JUMPE	T2,YZRO        ;Y = 0 IS A SPECIAL CASE

	JUMPE	T0,XZRO		;X = 0 IS A SPECIAL CASE

	PUSH	P,T5		;SAVE ANOTHER REGISTER
	DMOVEM  T0,TT0		;SAVE COPIES OF X AND
	DMOVEM  T2,TT2		; Y FOR FUTURE REFERENCE
	JUMPGE	T0,XPLUS	;IF X < 0
	DMOVN	T0,T0		;  NEGATE IT
XPLUS:	JUMPGE	T2,YPLUS	;IF Y < 0
	DMOVN	T2,T2		;  NEGATE IT
YPLUS:	DMOVEM	T0,TT4		;SAVE MAGNITUDES OF
	DMOVEM	T2,TT6		;  X AND Y
	AND	T0,MASK		;ISOLATE EXPONENT FIELDS
	AND	T2,MASK		;  OF X AND Y
	SUB	T2,T0		;EXPONENT OF Y - EXPONENT OF X
	CAML	T2,P32		;IF DIFFERENCE > 32
	JRST	BIG		;  |Y/X| IS LARGE
	CAMGE	T2,M32		;IF DIFFERENCE < -32
	JRST	SMALL		;  |Y/X| IS SMALL

	FUNCT	DATN2.,<TT2,TT0> ;NO EXCEPTIONS CAN OCCUR
	XMOVEI  T4,@1(L)	;GET POINTER AND STORE
	DMOVEM  T0,2(T4)	; IMAGINARY PART OF RESULT

	DMOVE	T4,TT4		;RECOVER |X|
	DMOVE	T2,TT6		;RECOVER |Y|
	CAML	T4,T2		;IF |X|HI .GE. |Y|HI,
	  JRST  NOXCH		;NO EXCHANGE NECESSARY
	EXCH	T4,T2		;|X|HI .LT. |Y|HI.  LARGER = A
	EXCH	T5,T3		;  TO T4. SMALLER = B TO T2
NOXCH:  MOVE	T1,T4		;HI OF LARGER TO T1
	LSH	T1,-33		;ISOLATE ITS BIASED EXPONENT
	SUBI	T1,200		;  AND UNBIAS IT
	JUMPLE  T1,LE		;IF UNBIASED EXPONENT POSITIVE,
	  SUBI  T1,1		;  DECREMENT IT
LE:	MOVN	T1,T1		;NEGATE SCALE INDEX
	FSC	T4,(T1)		;SCALE A TO [1/2,2) AND B BY
	FSC	T2,(T1)		;  SAME FACTOR. NO EXCEPTIONS
	LSH	T1,1		;MULTIPLY NEGATED SCALE INDEX BY 2
	DMOVEM  T4,BIGGER	;SAVE SCALED LARGER OF |X| AND |Y|
	DFMP	T4,T4		;SQUARE LARGER.  NO EXCEPTIONS
	DFMP	T2,T2		;SQUARE SMALLER. NO EXCEPTIONS
	DFAD	T4,T2		;GET SCALED SQUARE OF MODULUS
	CAMLE	T4,RT2		;IF T4 .GT. SQRT(2)HI
	  JRST  MORE		;  GO TO MORE SCALING
	CAMGE	T4,RT2O2	;IF T4 .LT. SQRT(1/2)HI
	  JRST  MORE		;  GO TO MORE SCALING
	JUMPN	T1,JOIN		;IF SCALE INDEX NOT 0, GO TO JOIN
	MOVEM	T1,NN		;  ELSE STORE 0 FOR FUTURE USE

;At this point the scale index is zero, which means that (A**2 + B**2)
;is unscaled. Also (A**2 + B**2) is between SQRT(1/2) and SQRT(2),
;which implies that a damaging loss of significance can occur in the
;calculation of (A**2 + B**2 - 1), which is one factor of Z, used
;in the rational approximation. Growth of error due to loss of
;significance can be avoided by calculating instead
;
;	     (A - 1) * (A + 1) + B**2
;
;which is done in the following lines of code:

	DMOVE	T4,BIGGER	;RECOVER LARGER
	DMOVE	T0,T4		;  AND COPY INTO T0
	DFSB	T0,ONE		;GET A - 1
	DFAD	T4,ONE		;AND A + 1
	DFMP	T4,T0		;AND THEIR PRODUCT
	DFAD	T2,T4		;T2 NOW HAS A**2 + B**2 - 1
	DMOVE	T4,T2		;COPY IT INTO T4
	DFAD	T2,TWO		;GET DENOMINATOR OF Z IN T2
	FSC	T4,1		;  AND NUMERATOR IN T4
	JRST	MERGE		;MERGE FOR RATIONAL APPROX

;The following code applies to scaled values of (A**2 + B**2) lying
;outside the interval (SQRT(1/2),SQRT(2)); it carries out more
;scaling to get the scaled (A**2 + B**2) inside this interval.

MORE:	MOVE	T2,T4		;HI OF T4 TO T2 FOR RESCALING
	LSH	T2,-33		;ISOLATE ITS BIASED EXPONENT
	SUBI	T2,200		;UNBIAS THE EXPONENT
	MOVN	T2,T2		;  AND NEGATE IT
	FSC	T4,(T2)		;T4 NOW IN [1/2,1)
	ADD	T1,T2		;TOTAL NEGATED SCALE INDEX TO T1
	CAML	T4,RT2O2	;IF T4 .GT. SQRT(1/2)
	  JRST  JOIN		;  SCALING IS DONE
	FSC	T4,1		;SCALE T4 TO [1,SQRT(2)]
	ADDI	T1,1		;  AND INCREMENT NEGATED INDEX

;At this point the scaled (A**2 + B**2) is in the interval
;[SQRT(1/2),SQRT(2)]. T1 contains the negative of the index
;necessary to compensate later for the scaling.  In the following
;lines of code, the index, together with its proper sign is
;temporarily stored.  Then a rational approximation is used to
;evaluate the natural log of the scaled (A**2 + B**2).

JOIN:	MOVNM	T1,NN		;RESTORE SIGN OF INDEX AND SAVE
	DMOVE	T2,T4		;GET COPY OF SCALED (A**2 + B**2)
	FSC	T4,1		;  AND MULTIPLY BY 2
	DFSB	T4,TWO		;SUBTRACT TWO TO GET NUMERATOR
	DFAD	T2,ONE		;ADD 1 TO GET DENOMINATOR OF Z

;The following code is taken from the DLOG routine.  The constants
;A0, A1, A2, B0, B1, and B2 are those given in the DLOG routine.
;Note that all tests restricting the scaled (A**2 + B**2) to
;the interval (SQRT(1/2),SQRT(2)) have been carried out only
;on the "high" words of these limiting values.  This is
;adequate for the validity of the approximation to be used,
;since only low order bits, and not the order of magnitude,
;of the error bound of the approximation are affected.

MERGE:  DFDV	T4,T2		;Z = ZNUM / ZDEN. NO EXCEPTIONS
	DMOVEM  T4,SAVEZ	;SAVE A COPY OF Z
	DFMP	T4,T4		;W = Z**2
	DMOVE	T0,T4		;  AND MAKE COPY
	DFAD	T4,B2		;FORM B(W) = W + B2
	DFMP	T4,T0		;  * W
	DFAD	T4,B1		;  + B1
	DFMP	T4,T0		;  * W
	DFAD	T4,B0		;  + B0
	DMOVE	T2,T0		;MAKE ANOTHER COPY OF 2
	DFMP	T0,A2		;FORM A(W) = W * A2
	DFAD	T0,A1		;  + A1
	DFMP	T0,T2		;  * W
	DFAD	T0,A0		;  + A0
	DFDV	T0,T4		;R(Z) = A(W) / B(W)
	DFMP	T0,SAVEZ	;  * Z
	DFMP	T0,T2		;  * W
	  JFCL			;NO OVERFLOW, UNDERFLOW WON'T HURT

;As indicated in the DLOG routine, Z still needs to be added into R(Z). 
;To increase accuracy, this addition is deferred until after the
;addition, below, of NN * LN(2))LO to the part of R(Z) so far
;obtained.

	FLTR	T2,NN		;RECOVER SCALE INDEX
	JUMPN	T2,REST		;IF NN = 0,
	DFAD	T0,SAVEZ	;  COMPLETE R(Z) BY ADDING IN Z
	JRST	DONE		;GO STORE RESULT AND RESTORE REGISTERS

;The following lines of code complete the calculation of the real
;part for NN not zero.  The real part is the sum of NN*(LN(2)HI),
;Z, R(Z), and NN*(LN(2)LO).  Of these terms the first is guaranteed
;to have the largest magnitude.  THe second has magnitude larger
;than the third by a factor of at least 2**7, and the last has
;magnitude smaller than the first by at least a factor of 2**(-10).
;The maximum value of |Z| is about 2/7, and this is the worst case
;from the point of view of accuracy, because it is only slightly
;smaller than (LN(2)HI).  The calculation of Z is accurate to a
;small multiple of its LSB, and hence the calculation of the sum
;is limited to a small multiple of ITS LSB.  Note that if |Z| is
;appreciably less than its bound of 2/7, the error bound of the
;sum decreases towards its LSB/2.  THe overall accuracy will be
;maximized by summing the terms in the reverse order in which
;they are listed above.

REST:	SETZ	T3,		;GET DFLOAT(NN) 
	DMOVE	T4,T2		;GET COPY IN T4
	DFMP	T4,C2		;GET NN * LN(2)LO
	DFMP	T2,C1		;GET NN * LN(2)HI
	DFAD	T0,T4		;R(Z) + NN * LN(2)LO
	DFAD	T0,SAVEZ	;  + Z
	DFAD	T0,T2		;  + NN * LN(2)HI
	
DONE:	FSC	T0,-1		;DIVIDE BY 2 TO GET LOG OF MODULUS
	DMOVEM  T0,@1(L)	;STORE REAL PART OF RESUUT

	POP	P,T5
	POP	P,T4
	POP	P,T3
	POP	P,T2
	POPJ	P,		;DONE


BIG:	FUNCT	DLOG.,<TT6>	;GET LN(|Y|)
	DMOVE	T2,DPIO2	;IMAGINARY PART NEAR +/- PI/2
	SKIPGE	TT2		;  + FOR Y POSITIVE
	DMOVN	T2,T2		;  - FOR Y NEGATIVE
	DMOVE	T4,TT4		;RECOVER |X|
	DFDV	T4,TT2		;GET |X/Y|. NO OVERFLOW
	JFCL	UND1		;  UNDERFLOW IS POSSIBLE
	JRST	SMERG		;MERGE WITH SMALL

SMALL:	FUNCT	DLOG.,<TT4>	;GET LN(|X|)
	SETZ	T2,		;IMAGINARY PART IS NEAR
	SETZ	T3,		;  0 OR +/- PI
	SKIPGE	TT0		;0 FOR X POSITIVE
	DMOVE	T2,CPI		;  +/-PI FOR X NEGATIVE
	SKIPGE	TT2		;    + PI FOR Y POSITIVE
	DMOVN	T2,T2		;    - PI FOR Y NEGATIVE
	DMOVE	T4,TT6		;RECOVER Y
	DFDV	T4,TT4		;GET Y/X. NO OVERFLOW
	JFCL	UND2		;  UNDERFLOW IS POSSIBLE
SMERG:	DFAD	T2,T4		;SMALL CORRECTION FOR IMAG
	DFMP	T4,T4		;GET (Y/X)**2. NO OVERFLOW
	  JFCL	UND1		;  INDERFLOW IS POSSIBLE
	FSC	T4,-1		;(1/2)*(Y/X)**2. NO OVERFLOW
	JFCL	UND1		;  UNDERFLOW IS POSSIBLE
	DFAD	T0,T4		;REAL = LN(|X|) + (1/2)*(Y/X)**2
	JRST	STOR		;NO EXCEPTIONS. GO STORE RESULTS

UND2:	JUMPN	T2,UND1		;NO MESSAGE IF |IMAG| NEAR PI OR PI/2
	LERR	(LIB,%,<CDLOG: Imaginary part underflow>)

UND1:	JUMPN	T0,STOR		;NO MESSAGE IF LN NOT ZERO
	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

	POP	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:  LERR	(LIB,%,<CDLOG: arg is (0.0,0.0); result=(+infinity,0.0)>
)
	HRLOI	T0,377777	;REAL ANSWER IS POSITIVE INFINITY
	SETO	T1,		;
	XMOVEI  T4,@1(L)	;GET THE ADDRESS OF THE SECOND ARGUMENT
	DMOVEM  T0,(T4)		;SAVE THE REAL PART OF THE ANSWER 
	DMOVEM  T2,2(T4)	;SAVE THE IMAG PART OF THE ANSWER
	POP	P,T4
	POP	P,T3		;RESTORE ACCUMULATORS
	POP	P,T2	
	POPJ	P,		;RETURN

;CONSTANTS

ONE:	EXP 201400000000,000000000000	;1.0
TWO:	EXP 202400000000,000000000000	;2.0
MASK:	EXP 377000000000		;EXPONENT MASK
P32:	EXP 040000000000		;HIGH 9 BITS = 40 OCTAL
M32:	EXP 740000000000		;HIGH 9 BITS = -40 OCTAL
CPI:	EXP 202622077325,021026430215	;PI
DPIO2:  EXP 201622077325,021026430215	;PI / 2
RT2O2:  EXP 200552023632		;SQRT(2) / 2
RT2:	EXP 201552023632		;SQRT(2)
A0:	EXP 570377400073,123045145570	;-.641249434237455811D2
A1:	EXP 205406111210,005527754477	;.163839435630215342D2
A2:	EXP 577153575223,052241327104	;-.789561128874912573D0
B0:	EXP 565177200130,374467610432	;-.769499321084948798D3
B1:	EXP 211470020376,205223175724	;.312032220919245328D3
B2:	EXP 571342517755,046030761117	;-.356679777390346462D2
B3:	EXP 576400000000,000000000000	;-.1D1
C1:	EXP 200543000000,000000000000	;LN(2)HI
C2:	EXP 613102775750,347571527450	;LN(2)LO

;DATA
	RELOC		
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
	RELOC
	PRGEND
TITLE	CDSIN	COMPLEX SINE FUNCTION
;		(DOUBLE PRECISION)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	CDSIN
EXTERN	CDSIN.
CDSIN=CDSIN.
PRGEND
TITLE	CDCOS	COMPLEX COSINE FUNCTION
;		(DOUBLE PRECISION)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
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) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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 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.

;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	FORPRM
	FSRCH
	TWOSEG	400000
	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:	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:	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:	LERR	(LIB,%,CDSIN or CDCOS; DABS(DIMAG(arg)) too large; DIMAG(result) = infinity)
	HRLOI	T2,377777		;SET IRES TO INFINITY
	SETO	T3,
	JRST	SKP3			;SKIP NEXT 3 INSTRUCTIONS
	
CONTN2:	DMOVEM	T2,ARGSAV
	FUNCT	DEXP.,<ARGSAV>		;IRES = EXP(T)
	DMOVE	T2,T0

SKP3:	DMOVE	T0,SX
	MOVE	T4,YSAVE
	XOR	T4,CX			;T4 HAS THE SIGN OF CX*Y
	JUMPGE	T4,RET			;IF T4 IS LESS THAN 0.
	  DMOVN	T2,T2			;THEN NEGATE IRES
	JRST    RET			;RETURN
IMAG0:	SETZ	T2,			;SET IMAGINARY PART 
	SETZ	T3,			; OF RESULT TO ZERO.
	JRST	RET			;RETURN

EASY:	DMOVE	T2,YSAVE		;GET SIGNED VALUE OF Y
	DFMP	T2,CX			;SINH(Y)*COS(X) = Y*COS(X)
	  JFCL	IMUND			;UNDERFLOW POSSIBLE.
	DMOVE	T0,SX			;COSH(Y)*SIN(X) = SIN(X).

RET:	XMOVEI	T4,@(L)			;GET ADDRESS OF 2ND ARGUMENT
	DMOVEM	T0,(T4)			;SAVE REAL PART OF RESULT
	DMOVEM	T2,2(T4)		;SAVE IMAGINARY PART OF RESULT
	POP	P,T4
	POP	P,T3
	POP	P,T2
	POP	P,T5
	GOODBY	(1)

IMUND:	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

	RELOC					;DATA
CX:	0
CXL:	0
SX:	DOUBLE	0,0
YSAVE:	DOUBLE	0,0
ARGSAV: DOUBLE	0,0
TEMP:	DOUBLE	0,0
	RELOC
	PRGEND
TITLE	CDSQRT	COMPLEX SQUARE ROOT FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
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) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

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

;  LET Z = X + I*Y
;  THEN THE ANSWER CDSQRT(Z) = U + I*V IS DEFINED AS FOLLOWS

;  CDSQRT FIRST CHECKS THE TYPE OF THE ARGUMENT SO THAT THE
;  APPROPRIATE RESULT (DOUBLE PRECISION OR EXTENDED
;  PRECISION) IS RETURNED.

;  IF X AND Y ARE BOTH >= 0.0, THEN
;    U = DSQRT((DABS(X)+CDABS(Z))/2.0)
;    V = Y/(2.0*U)

;  IF X >= 0.0 AND Y < 0.0, THEN
;    U = SQRT((DABS(X)+CDABS(Z))/2.0)
;    V = Y/(2.0*U)

;  IF X < 0.0 AND Y >= 0.0, THEN
;    U = Y/(2.0*V)
;    V = DSQRT((DABS(X)+CDABS(Z))/2.0)

;  IF X AND Y ARE BOTH < 0.0, THEN
;    U = Y/(2.0*V)
;    V = -DSQRT((DABS(X)+CDABS(Z))/2.0)

;  THE RESULT IS IN THE RIGHT HALF PLANE, THAT IS, THE POLAR ANGLE OF THE
;    RESULT LIES IN THE INTERVAL [-PI/2,+PI/2].

;REQUIRED (CALLED) ROUTINES:  DSQRT

;REGISTERS T2, T3, T4, T5, G1, G3, AND G4 WERE SAVED, USED, AND RESTORED

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  XMOVEI	L,ARG
;  PUSHJ	P,CDSQRT

;THE ROUTINE IS CALLED WITH TWO DOUBLE PRECISION VECTORS.
;THE REAL PART OF THE ARGUMENT IS EXPECTED TO BE STORED IN THE
;FIRST DOUBLE PRECISION WORD OF THE FIRST VECTOR.
;THE IMAGINARY PART OF THE ARGUMENT IS EXPECTED TO BE STORED IN THE
;SECOND DOUBLE PRECISION WORD OF THE FIRST VECTOR.
;THE REAL PART OF THE SOLUTION IS RETURNED IN THE FIRST DOUBLE
;PRECISION WORD OF THE SECOND VECTOR.
;THE IMAGINARY PART OF THE SOLUTION IS RETURNED IN THE SECOND DOUBLE
;PRECISION WORD OF THE SECOND VECTOR.

	SEARCH	FORPRM
	FSRCH
	TWOSEG	400000
	SALL

	HELLO	(CDSQT.,CDSQRT)	;ENTRY TO COMPLEX SQUARE ROOT ROUTINE.
IFN FTGFL,<
	LDB	T1,[POINTR 0(L),ARGTYP]	;GET FIRST ARG TYPE
	CAIN	T1,TP%DPX	;G FLOATING?
	  JRST	CGSQRT##	;YES
> ;END FTGFL
	PUSH	P,T2		;SAVE ACCUMULATORS
	PUSH	P,T3
	XMOVEI	T2,@0(L)	;PICK UP ADDRESS OF ARGUMENT
	DMOVE	T0,(T2)		;PICK UP X IN AC T0.
	DMOVE	T2,2(T2)	;PICK UP Y IN AC T2.
	PUSH	P,T4		;SAVE AC T4.
	PUSH	P,T5
	JUMPE	T2,ZREAL	;JUMP IF Y=0
	JUMPE	T0,ZIMAG	;JUMP IF X=0
	DMOVE	T4,T0		;DABS(X) TO AC T4.
	JUMPGE	T0,XPOS		;IF X IS NEGATIVE
	DMOVN	T4,T4		;NEGATE

XPOS:   PUSH	P,G1		;SAVE REGISTER G1.
	PUSH	P,G3		;SAVE REGISTER G3
	PUSH	P,G4		;SAVE REGISTER G4
	DMOVEM	T0,XSAVE	;SAVE X IN XSAVE
	DMOVEM	T2,YSAVE	;SAVE Y IN YSAVE
	DMOVE	G3,T2		;DABS(Y) TO AC G3.
	JUMPGE	T2,YPOS		;IF Y IS NEGATIVE
	  DMOVN	G3,T2		;NEGATE IT
YPOS:	HRRZI	G1,1		;SET G1 TO 1
	CAML	T4,G3		;IF ABS(X) > ABS(Y)
	  JRST	DWN		;THEN GO TO DWN
				;IF HIGH WORDS OF |X| AND |Y| ARE EQUAL,
				; DON'T BOTHER DECIDING WHETHER TO EXCHANGE
				; THEM.  IF |X| AND |Y| ARE NEARLY EQUAL,
				; THEIR MAX AND MIN MAY BE SWITCHED BELOW.
				; THIS MEANS THAT S/L MAY BE SLIGHTLY GREATER
				; THAN 1, BUT THIS DOES NOT HURT.
XCHNG:	EXCH	T4,G3		;PUT MAX(ABS(X),ABS(Y)) IN T4
	EXCH	T5,G4		;AND MIN(ABS(X),ABS(Y)) IN G3,G4
	SETZ	G1,		;SET G1 TO 0
DWN:	DFDV	G3,T4		;CALC S/L.
	JFCL			;NO UNDERFLOW ERROR MESSAGE ALLOWED.
	DFMP	G3,G3		;CALC (S/L)**2.
	JFCL			;NO UNDERFLOW ERROR MESSAGE ALLOWED.
	DFAD	G3,ONE		;HAVE (1+(S/L)**2) IN AC G3.
	DMOVEM	G3,TEMP
	FUNCT	DSQRT.,<TEMP>	;CALC. THE DSQRT OF
				;(1+(S/L)**2)=1+EPS.
	JUMPE	G1,YGTX		;GO TO YGTX IF DABS(Y) > DABS(X).

XGTY:	DFAD	T0,ONE		;CALC. 2 + EPS.
	FSC	T0,-1		;CALC. 1+EPS/2.
	DMOVEM	T0,TEMP		;PUT IN TEMP FOR
	FUNCT	DSQRT.,<TEMP>	;CALL TO DSQRT
	DMOVE	G3,T0		;SAVE DSQRT(1+EPS/2) IN AC G3.
	DMOVEM	T4,TEMP
	FUNCT	DSQRT.,<TEMP>	;CALC.
	DFMP	T0,G3		;CALC. DSQRT(DABS(X)*(1+EPS/2)).
	JRST	OUT1		;GO TO REST OF CALC.

YGTX:	CAIGE	T4,201400	;IF DABS(Y)>1, GO TO POSSIBLE OVERFLOW CASE.
	JRST	YXSMAL		;ELSE, GO TO POSSIBLE UNDERFLOW.
	FSC	T0,-3		;CALC. (1+EPS)/8
	DFMP	T4,T0		;CALC. DABS(Y)*(1+EPS)/8
	DMOVE	T0,XSAVE	;DABS(X) TO AC T0
	JUMPGE	T0,NXT
	  DMOVN	T0,T0
NXT:	FSC	T0,-3		;CALC. DABS(X)/8
	JFCL			;SUPPRESS UNDERFLOW ERROR MESSAGE.
	DFAD	T0,T4		;CALC.(DABS(X)/8)+(DABS(Y)*(1+EPS)/8).
	DMOVEM	T0,TEMP
	FUNCT	DSQRT.,<TEMP>	;CALC.
	FSC	T0,1		;DSQRT(DABS(X)/2 + (DABS(Y)*(1+EPS)/2)).

OUT1:	JUMPGE	T2,POSY
	  DMOVN	T2,T2
POSY:	DFDV	T2,T0		;DIVIDE DABS(Y)/2 BY THE
	  JFCL	UNDFL
	FSC	T2,-1		;SQRT TERM.
	  JFCL	UNDFL
	JRST	SIGNS		;GO TO REST OF CALC.

YXSMAL:	DFMP	T0,T4		;DABS(Y)*(1+EPS) = CDABS(Z).
	DMOVE	G3,XSAVE	;DABS(X) TO AC G3.
	JUMPGE	G3,XOK		;IF X IS NEGATIVE
	  DMOVN	G3,G3		;NEGATE IT
XOK:	DFAD	G3,T0		;GET |X| + CDABS(Z)
	FSC	G3,1		;MULTIPLY ARG BY 2
	DMOVEM	G3,TEMP		;SAVE FOR DSQRT
	FUNCT	DSQRT.,<TEMP>	;CALC. DSQRT OF 2 * (|X| + CDABS(Z))
	DMOVE	G3,T0		;GET DSQRT((|X| + CDABS(Z))*2)
	FSC	T0,-1		;GET DSQRT((|X| + CDABS(Z))/2)
	DFDV	T4,G3		;GET |Y| / DSQRT (2 * (|X| + CDABS(Z)))
	DMOVE	T2,T4		;PUT A TEMP ANSWER IN AC T2.

SIGNS:  SKIPL	XSAVE		;IF NEGATIVE X, SWITCH REAL & IMAG
	  JRST	OK		;  PARTS OF RESTLT
	EXCH	T2,T0		;EXCHANGE
	EXCH	T1,T3		;EXCHANGE
				;SET UP THE REAL AND
  	  			;THE IMAGINARY ANSWERS WITH
OK:	SKIPGE	YSAVE		;THE APPROPRIATE
	  DMOVN	T2,T2		;SIGNS.
	XMOVEI	T4,@1(L)	;GET THE ADDRESS OF THE SECOND ARGUMENT
	DMOVEM	T0,(T4)		;SAVE REAL PART OF RESULT
	DMOVEM	T2,2(T4)	;SAVE IMAG PART OF RESULT
	POP	P,G4		;RESTORE AC G4
	POP	P,G3		;RESTORE AC G3
	POP	P,G1		;RESTORE AC G1
	POP	P,T5		;RESTORE AC T5
	POP	P,T4		;RESTORE AC T4
	POP	P,T3		;RESTORE AC T3
	POP	P,T2		;RESTORE AC T2
	GOODBY	(1)		;RETURN

ZREAL:	JUMPE	T0,DONE		;Y = 0, HOW ABOUT X?
	MOVE	T4,T0		;X NOT ZERO, SAVE IT.
	JUMPGE	T0,POSX		;GET DABS(X) FOR DSQRT
	  DMOVN	T0,T0
POSX:	DMOVEM	T0,TEMP		;
	FUNCT	DSQRT.,<TEMP>	;T0,T1 GET DSQRT(DABS(X))
	JUMPG	T4,DONE		;T0,T1 AND T2,T3 OK FOR X .GE. 0
	  EXCH	T0,T2		;INTERCHANGE FOR X<0
	  EXCH	T1,T3		;
	JRST	DONE		;FINISHED

ZIMAG:	MOVE	T4,T2		;X=0, SAVE Y
	JUMPGE	T2,PSY		;IF Y < 0
	  DMOVN	T2,T2		;  NEGATE IT
PSY:	FSC	T2,-1		;GET |Y|/2
	DMOVEM	T2,TEMP		; AND STORE IT FOR DSQRT
	FUNCT	DSQRT.,<TEMP>	;T0,T1 GETS DSQRT(DABS(Y/2))
	DMOVE	T2,T0		;T2,T3  ALSO
	JUMPG	T4,DONE		;DONE IF Y > 0
	  DMOVN	T2,T2		;NEGATE IMAG PART IF Y < 0

DONE:	XMOVEI	T4,@1(L)	;GET THE ADDRESS OF THE SECOND ARGUMENT
	DMOVEM	T0,(T4)		;SAVE REAL PART OF RESULT
	DMOVEM	T2,2(T4)	;SAVE IMAG PART OF RESULT
	POP	P,T5		;RESTORE ACCUMULATORS
	POP	P,T4
	POP	P,T3
	POP	P,T2		;RESTORE AC T2.
	GOODBY	(1)		;RETURN

UNDFL:	SKIPGE	XSAVE		;DOES RESULT HAVE REAL & IMAG PARTS SWITCHED?
	  LERR	(LIB,%,<CDSQRT: real part underflow>,UFCONT)
	LERR	(LIB,%,<CDSQRT: imaginary part underflow>)
UFCONT:	SETZB	T2,T3		;STORE 0 FOR RESULT OF UNDERFLOW
	JRST	SIGNS



ONE:	DOUBLE	201400000000,000000000000	;1.0D0
SQ2:	DOUBLE	201552023631,237635714441	;SQRT(2).

	RELOC					;DATA
TEMP:	DOUBLE	0,0				;TEMPORARY STORAGE FOR DSQRT ARGS
XSAVE:	DOUBLE	0,0
YSAVE:	DOUBLE	0,0
	RELOC
	END