Google
 

Trailing-Edge - PDP-10 Archives - fortv11 - mthcpx.mac
There are 9 other files named mthcpx.mac in the archive. Click here to see a list.

	SEARCH	MTHPRM
	TV	MTHCPX	COMPLEX ROUTINES ,7(4002)

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 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	16-Mar-81	Q10-04786
	Fix TWOSEG and RELOC problems.

1405	DAW	6-Apr-81
	Support extended addressing.

1464	DAW	12-May-81
	Error messages.

1673	CDM	9-Sept-81
	Added complex routines.
	(CMPL.I, CMPL.C, CMPL.D)

***** Begin Mathlib *****

3024	CDM	17-Nov-81
	Changed REAL. to REAL.C, defining REAL. to have integer
	arguments as in ANSI-77 standard.

3200	JLC	10-May-82
	Mathlib integration. Change all LERRs to $LCALLs. Change
	TWOSEG/RELOC to SEGMENT macros.

3201	RVM	20-May-82
	Add the routine CMPL.G (convert two gfloating numbers to
	complex).

3217	JLC	8-Oct-82
	Modify SIXBIT names of the complex multiply and divide routines
	so that TRACE won't look for an argument block.

3220	PLB	12-Oct-82
	Remove IFIWs from invocations of FUNCT macro, so that we
	can put IFIW into the defn of FUNCT.

***** Begin Version 1A *****

3234	CKS	23-Mar-83
	The maximum value of the arg to SIN was given erroneously
	as 36394.429 instead of the correct value of 823549.66.  
	Fix this in the code and any of the comments which it
	appears.

3247	10-33795	MRB	24-May-83
	Fixed bug in the interface of CEXP2 with CEXP3, and to
	fix bug in block 2 (of CEXP3) at XISO.

3251	10-34015	MRB	30-Aug-83
	Added a special case in CEXP2. When the imaginary part of 
	the complex number is zero and the real part in an integer
	call EXP2. instead. (This is edit 3344 to FOROT7) 

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

***** Begin Version 2 *****

4002	JLC	8-Aug-83
	Install the beginnings of user-fixable library routine
	error results. Fix an overflow msg.

\

	PRGEND
TITLE	CABS	COMPLEX ABSOLUTE VALUE FUNCTION  
;		(SINGLE 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	CABS
EXTERN	CABS.
CABS=CABS.
PRGEND
TITLE	CABS.	COMPLEX ABSOLUTE VALUE FUNCTION
;		(SINGLE 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.

;CABS(Z) IS CALCULATED AS FOLLOWS

;  LET Z = X + I*Y
;  V = MAX(ABS(X),ABS(Y))
;  W = MIN(ABS(X),ABS(Y))
;  THEN
;	CABS = V*SQRT(1.0 + (W/V)**2)

;THE RANGE OF DEFINITION FOR CABS IS THE REPRESENTABLE REAL NUMBERS.
;  AN OVERFLOW CAN OCCUR, IN WHICH CASE CABS = + MACHINE INFINITY.

;REQUIRED (CALLED) ROUTINES:  SQRT

;REGISTER T2 WAS SAVED, USED, AND RESTORED

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  XMOVEI	L,ARG
;  PUSHJ	P,CABS

;THE ANSWER IS RETURNED IN ACCUMULATOR T0

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(CABS,.)	;ENTRY TO MODULUS ROUTINE
	PUSH	P,T2
	XMOVEI	T2,@(L)		
	MOVM	T0,(T2)		;ABS(X)
	MOVM	T2,1(T2)	;ABS(Y)

;				 OBTAIN MAX(ABS(X),ABS(Y))IN T2
;				 AND MIN IN T0
	CAML	T0,T2
	  EXCH	T2,T0		;THEN CALCULATE CABS
	JUMPE	T2,RET		;Z = 0, HENCE CABS = O
	FDVR	T0,T2		;U/V
	  JFCL	ANS		;NO OVERFLOW, RATIO NEGLIGIBLE IF UNDERFLOW
	CAMG	T0,TWOM14	;IF RATIO .LE. 2**-14,
	  JRST	ANS		;RATIO IS NEGLIGIBLE.
	FMPR	T0,T0		;**2
	FADRI	T0,201400	;+1.0
	MOVEM	T0,TEMP		
	FUNCT	SQRT.,<TEMP>	;[3220] SQUARE ROOT IS IN AC T0
     	FMPR	T0,T2		;*V 
	  JFCL	OVFL		;NO UNDERFLOW, GET MESSAGE ON OVERFLOW
RET:	POP	P,T2
	GOODBY	(1)		;RETURN

OVFL:
	$SSPR			;SAVE FIXED-UP RESULT
	$LCALL	ROV		;RESULT OVERFLOW
	$RSPR			;GET FIXED-UP RESULT
	JRST	RET		;RETURN

ANS:	MOVE	T0,T2		;ANSWER = ABS(LARGER) TO T0
	POP	P,T2		;RESTORE T2
	GOODBY	(1)		;RETURN

TWOM14:	163400000000		;2**-14

	SEGMENT	DATA
TEMP:	0			;TEMPORARY STORAGE USED FOR SQRT CALL

	PRGEND
TITLE	CEXP	COMPLEX EXPONENTIAL FUNCTION  
;		(SINGLE 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	CEXP
EXTERN	CEXP.
CEXP=CEXP.
PRGEND
TITLE	CEXP.	COMPLEX EXPONENTIAL FUNCTION
;		(SINGLE 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.

;CEXP(Z), WHERE Z = X + I*Y, IS CALCULATED AS FOLLOWS
;    CEXP(Z) = EXP(X)*(COS(Y)+I*SIN(Y))
;THE RANGE OF DEFINITION FOR CEXP IS AS FOLLOWS

;FOR Z = X+I*Y, IF ABS(Y) .GT. 823549.66 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 X/2. IS .GT. 88.029692, THE ABS(REAL(RESULT)) IS SET TO
;    +INFINITY AND ABS(AIMAG(RESULT)) IS SET TO +INFINITY AND
;    AN ERROR MESSAGE IS RETURNED.

;REQUIRED (CALLED) ROUTINES:  EXP,SIN,COS

;REGISTERS T2, T3, AND T4 WERE SAVED, USED, AND RESTORED

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  MOVEI	L,ARG
;  PUSHJ	P,CEXP

;THE REAL PART OF THE SOLUTION IS RETURNED IN ACCUMULATOR T0
;THE IMAG PART OF THE SOLUTION IS RETURNED IN ACCUMULATOR T1

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(CEXP,.)	;ENTRY TO NAME ROUTINE
	PUSH	P,T2		;SAVE REGISTER T2
	PUSH	P,T3		;SAVE REGISTER T3
	PUSH	P,T4		;SAVE REGISTER T4
	XMOVEI	T1,@(L)		;GET ADDRESS OF Z
	MOVE	T0,(T1)		;X=REAL(Z)
	MOVE	T1,1(T1)	;Y=IMAG(Z)
	MOVEM	T1,YSAVE	;SAVE Y
	MOVM	T1,T1		;T1=ABS(T1)
	CAMG	T1,YMAX		;IF ABS(Y) IS NOT TOO LARGE
	  JRST	YOK		;GO TO YOK
	$SCPX			;SAVE FIXED-UP RESULT
	$LCALL	AIZ
;LERR	(LIB,%,<CEXP: AIMAG(arg) too large in absolute value; result=(0,0)>)
	$RCPX			;GET FIXED-UP RESULT
	SETZ	T0,		;REAL PART OF SOLUTION IS ZERO
	SETZ	T1,		;IMAG PART OF SOLUTION IS ZERO
	JRST	RET		;RETURN

YOK:	MOVEM	T0,XSAVE	;SAVE X
	CAML	T0,NEGCON	;IF X IS NOT TOO SMALL
	  JRST	XOK		;GO TO XOK
	SETZ	T0,		;REAL PART OF SOLUTION IS ZERO
	SETZ	T1,		;IMAG PART OF SOLUTION IS ZERO
	JRST	YCHK		;CHECK Y

XOK:	FUNCT	COS.,<YSAVE>	;[3220] COSY=COS(Y)
	MOVEM	T0,COSY		
	FUNCT	SIN.,<YSAVE>	;[3220] SINY=SIN(Y)
	MOVEM	T0,SINY	
	MOVE	T2,XSAVE	;T2=X
	MOVE	T1,YSAVE	;T1=Y
	CAMG	T2,TBIG		;IF X IS NOT TOO LARGE
	  JRST	ALG		;GO TO ALG
	CAIN	T1,0		;ELSE, IF Y=0
	  JRST	ERR0		;GO TO ERR
	FSC	T2,-1		;ELSE, S=X/2
	CAMLE	T2,TBIG		;IF S IS TOO LARGE
	  JRST	STIM		;GO TO STIM
	MOVEM	T2,XSAVE
	FUNCT	EXP.,<XSAVE>	;T=EXP(S)
	MOVE	T1,T0		
	FMPR	T1,SINY		;V=T*SINY
	MOVM	T1,T1		;V=ABS(V)
	HRLOI	T4,377777	;T4=XMAX
	FDVR	T4,T0		;Q=XMAX/T
	CAIG	T1,T4		;IF V IS LESS THAN OR EQUAL TO Q
	  JRST	GETI		;GO TO GETI

STIM:	HRLOI	T1,377777	;ELSE, SET IMAG SOLUTION TO XMAX
	$LCALL	IPO
;LERR	(LIB,%,<CEXP: imaginary part overflow>)
	JRST	D2
	
GETI:	FMPR	T1,T0		;IRES = V*T
D2:	MOVE	T4,SINY			
	CAIGE	T4,		;IF SINY IS LESS THAN 0.0
	  MOVN	T1,T1		;THEN NEGATE IRES
	CAMLE	T2,TBIG		;IF S IS TOO LARGE
	  JRST 	ERR0		;GO TO ERR
	MOVE	T2,T0
	FMPR	T0,COSY		;V = T*COSY
	MOVM	T0,T0		;V = ABS(V)
	HRLOI	T3,377777	;T3=XMAX
	FDVR	T3,T2		;Q = XMAX/T
	CAIG	T0,T3		;IF V IS LESS THAN OR EQUAL TO Q
	  JRST	LAB		; THEN GO TO LAB

ERR0:	HRLOI	T0,377777	;RRES=XMAX
	$LCALL	RPO
;LERR	(LIB,%,<CEXP: real part overflow>)
	JRST	DD2

LAB:	FMPR	T0,T2		;RRES = V*T
DD2:	MOVE	T2,COSY
	CAIGE	T2,		;IF COSY.LT. 0.0
	  MOVN	T0,T0		;THEN NEGATE RRES
	JRST RET		;RETURN

ALG:	MOVEM	T2,XSAVE
    	FUNCT	EXP.,<XSAVE>	;EXPX = EXP(X)
	MOVEM	T0,EXPX
	MOVE	T1,EXPX
	FMPR	T1,SINY		;IRES = EXPX*SINY
	JFCL
	FMPR	T0,COSY		;RRES = EXPX*COSY
	JFCL
	CAIE	T1,		;IF IRES .NE. 0.0
	  JRST	RCHK		;THEN GO CHECK RRES

YCHK:	MOVE	T2,YSAVE		
	CAIN	T2,		;IF Y .NE. 0
	  JRST RCHK		;THEN GO CHECK RRES
	$LCALL	IPU
;LERR	(LIB,%,<CEXP: imaginary part underflow>)

RCHK:	CAIN	T0,		;IF R = 0.0
	$LCALL	RPU
;LERR	(LIB,%,<CEXP: real part underflow>)

RET:	POP	P,T4
	POP	P,T3
	POP	P,T2
    	GOODBY	(1)		;RETURN

YMAX:	823549.66		;[3234] MAX ARG TO SIN/COS
TBIG:   207540074636		;88.029692
NEGCON: 570232254037		;-89.415986

	SEGMENT	DATA
SINY:   0
COSY:   0
XSAVE:  0
EXPX:   0
YSAVE:  0

	PRGEND
TITLE	CEXP2	COMPLEX ** INTEGER

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 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	CEXP2
EXTERN	CEXP2.
CEXP2=CEXP2.
PRGEND
TITLE	CEXP3	COMPLEX ** COMPLEX

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 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	CEXP3
EXTERN	CEXP3.
CEXP3=CEXP3.
PRGEND
	TITLE	CEXP3.	Complex ** complex exponenentiation
	SUBTTL	Mary Payne /MHP			29-May-81


;	New Routine for CEXP3; Mary Payne, May 29, 1981.
;
;This routine calculates (X + iY)**(A + iB), where i is SQRT(-1).

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

;[3247]THE CEXP2 ENTRY:
;[3247]
;[3247]CEXP2 CALCULATES (X + IY)**N, WHERE N IS A (REAL) INTEGER.
;[3247]THE INITIALIZATION FOR THE CEXP2 ENTRY IS GIVEN BELOW. THIS
;[3247]CODE COMPLETES THE CALCULATION FOR THE  EXPONENT = 0 OR 1,
;[3247]AND ALSO COMPLETES THE CASE FOR WHICH THE COMPLEX BASE IS
;[3247]0 + I*0. IF X = 0, THE ROUTINE MERGES WITH CEXP3 AT XTINY;
;[3247]IN BLOCK 2. ALL OTHER CASES MERGE WITH THE CEXP3 ROUTINE AT
;[3247]JOIN3: IN BLOCK 2.

	HELLO	(CEXP2,.)	;[3247]CEXP2, COMPLEX**INTEGER
	DMOVEM	T4,SAVT4	;[3247]SAVE REGISTERS T4,T5
	DMOVE	T0,@(L)		;[3247]X + IY TO T0,T1
	MOVE	T4,@1(L)	;[3247]INTEGER EXPONENT TO T4
	JUMPE	T4,I0		;[3247]SEPARATE OUT EXPONENT = 0
	MOVE	T5,T4		;[3247]COPY EXPONENT TO T5
	SOJE	T5,I1		;[3247]SEPARATE OUT EXPONENT = 1
	  JFCL			;[3247]  SUPPRESS ERROR MESSAGE
	JUMPE	T1,CX1		;[3251]IMAGINARY PART = 0
	JUMPN	T0,CX2		;[3251]ARGS (<>0,<>0) = NORMAL CASE 
	JRST	NOREAL		;[3251]NO REAL PART = SPECIAL CASE
CX1:	JUMPL	T0,NOIMAG	;[3251]ARGS (-REAL,0) = SPECIAL CASE
	JUMPE	T0,ARG0		;[3251]ARGS (0,0) = SPECIAL CASE
;
;[3247]THE CX2: BLOCK OF CODE DOES THE PROPER INITIALIZATION
;[3247]FOR MERGING WITH CEXP3 AT JOIN3: OR XTINY: IN BLOCK 2.

CX2:	DMOVEM	T2,SAVT2	;[3247]SAVE REGISTERS T2,T3
	PUSHJ	P,DFL.4##	;[3247]D. P. FLOAT EXPONENT
	DMOVEM	T4,REXP		;[3247]SAVE REAL PART OF EXPONENT
	SETZB	T4,ISIGN	;[3247]INITIALIZE IMAGINARY SIGN FLAG
	SETZB	T5,RSIGN	;[3247]  REAL SIGN FLAG AND PHI
	MOVEM	T4,PHIFLG	;[3247]  FLAG TO ZERO
	DMOVEM	T4,IEXP		;[3247]Save 0 IMAG. PART OF EXPONENT

	MOVE	T4,T0		;[3247]COPY X TO T4 WITH 0 EXTEND
	MOVE	T2,T1		;[3247]  IN T5. COPY Y TO T2
	SETZ	T3,		;[3247]  AND ZERO EXTEND
	EXCH	T0,T1		;[3247]Y TO T0, X TO T1
	JUMPE	T1,XTINY	;[3247]IF X = 0, MERGE WITH CEXP3 AT XTINY
	JRST	JOIN3		;[3247] ELSE AT JOIN3, BOTH IN BLOCK 2.
;[3247]THE REST OF THE CODE DEALS WITH THE SPECIAL CASES OF
;[3247]EXPONENT = 0 OR 1, AND X = Y = 0.

I1:	DMOVE	T4,SAVT4	;[3247]EXPONENT IS 1. RESTORE REGISTERS
				;[3247]T4, T5, AND RETURN THE COMPLEX
				;[3247]ARGUMENT AS THE RESULT
	POPJ	P,		;[3247]RETURN

I0:	JUMPN	T0,ANSRVM	;[3247]EXPONENT IS 0
	JUMPE	T1,INDET	;[3247]IF X = Y = 0, RESULTS INDETERMINATE
ANSRVM:	MOVE	T0,SPONE	;[3247]IF X AND Y NOT BOTH 0,
	SETZ	T1,		;[3247]  (X+IY)**0 IS 1.0 + I**0
	DMOVE	T4,SAVT4	;[3247]RESTORE T4,T5 AND
	POPJ	P,		;[3247]RETURN

INDET:	SETZB	T0,T1		;[3247]FOR 0**0, RETURN 0 + I*0
	$SCPX			;SAVE FIXED-UP RESULT
;[3247]	LERR	(LIB,%,<CEXP2:(0,0)**0 IS INDETERMINATE, RESULT IS (0,0)>)
	$LCALL	ZZZ		;[3247]ZERO**ZERO IS INDETERMINATE
	$RCPX			;GET FIXED-UP RESULT
	DMOVE	T4,SAVT4	;[3247]RESTORE REGISTERS T4, T5
	POPJ	P,		;[3247]RETURN
;[3251]
;[3251]	>>>>> Special case (0,iY)**n <<<<<
;[3251]
NOREAL:	MOVEM	T1,EXPARG	;[3251]COPY OF THE IMAGINARY PART
	MOVEM	T4,EXPARG+1	;[3251]COPY OF THE EXPONENT
	FUNCT	EXP2.,<EXPARG,EXPARG+1>;[3251]DO THE EXPONENT THE EASY WAY
	TRNE	T4,2		;[3251]NUMBER NEED TO BE NEGATED?
	MOVN	T0,T0		;[3251]YES; NEGATE THE RESULT.
	SETZ	T1,		;[3251]ASSUME EXPONENT EVEN; 
				;[3251]RESULT = (IY^N,0)
	TRNE	T4,1		;[3251]IS EXPONENT ODD?
	EXCH	T0,T1		;[3251]YES; RESULT = (0,IY^N)
	DMOVE	T4,SAVT4	;[3251]RESTORE REGISTERS T4, T5
	POPJ	P,		;[3251]RETURN
;[3251]
;[3251]	>>>>> SPECIAL CASE (-X,0)**N <<<<<
;[3251]
NOIMAG:	MOVEM	T0,EXPARG	;[3251]COPY OF THE REAL PART
	MOVEM	T4,EXPARG+1	;[3251]COPY OF THE EXPONENT
	FUNCT	EXP2.,<EXPARG,EXPARG+1>;[3251]DO THE EXPONENT THE EASY WAY
	SETZ	T1,		;[3251]NO IMAGINARY PART
	DMOVE	T4,SAVT4	;[3251]RESTORE REGISTERS T4, T5
	POPJ	P,		;[3251]RETURN

ARG0:	JUMPL	T4,INFIN	;[3247]ARG = (0,0). JUMP IF EXPONENT < 0
	DMOVE	T4,SAVT4	;[3247]ELSE RESULT IS 0 + I*0. RESTORE
	POPJ	P,		;[3247]REGISTERS T4,T5 AND RETURN

INFIN:	HRLOI	T0,377777	;[3247](0 + I*0)**NEGATIVE IS IMPLIED DIV
	MOVE	T1,T0		;[3247]  BY 0. BOTH PARTS = 377777777777
;[3247]	LERR	(LIB,%,<CEXP2: (0,0)**NEGATIVE: RESULT = (BIGGEST,BIGGEST)>)
	$LCALL	ZNI		;[3247]ZERO**NEGATIVE IS INFINITY
	DMOVE	T4,SAVT4	;[3247]RESTORE REGISTERS T4,T5
	POPJ	P,		;[3247]Return

;[3247]
;[3247]	CEXP3 ROUTINE STARTS HERE:
;[3247]
	HELLO	(CEXP3,.)	;CEXP3, complex ** complex
;
;			*****************
;			*    BLOCK 1    *
;			*****************
;
;Separate out the special cases: X = Y = 0 and A = B = 0:


	MOVEM	T2,SAVT2	;Save Registers T2
	MOVEM	T3,SAVT3	;  and T3
	MOVEM	T4,SAVT4	;  and T4
	DMOVE	T1,@(L)		;X and Y to T1 and T2.
	DMOVE	T3,@1(L)	;A and B to T3 and T4.
	JUMPN	T1,NXTST	;[3247]If X not zero, continue
				;  with main flow below.
	JUMPN	T2,XIS0		;Special case: X = 0, Y not.
				;  Treated in BLOCK 2 after
				;  its main flow is complete.

;	>>>>  Special Case: X = Y = 0  <<<<
;
;In the following block of code X = Y = 0. No detailed
;calculations are needed. If A > 0, both the real and
;imaginary parts of the result are zero. If A = 0
;the complex result has modulus one, and indeterminate
;direction, so both its real and imaginary parts should
;be identified as indeterminate. If A < 0, the complex
;result has modulus infinity, and indeterminate direction.
;The real and imaginary parts should probably be identified
;as indeterminate. A is in T3.
;
	JUMPG	T3,ZERO		;If A > 0, result is (0,0).
	JUMPL	T3,BTHIND	;If A = 0, results indeterminate.
	$LCALL	ZZZ		;zero**zero is indeterminate, result = zero
	SETZB	T0,T1		;Return (0,0)
	JRST	RST234		;Restore registers 2, 3, and 4.
BTHIND:	HRLOI	T0,377777	;Else results are indeterminate
	HRLOI	T1,377777	;Provide something for storage
	JUMPE	T4,CPXZNI	;[4002]If no complex part, treat like
				;[4002] zero**neg real
	$LCALL	ZCI		;(0,0)**(negative,non-zero) is indeterminate
	JRST	RST234		;Restore registers 2, 3, and 4.

CPXZNI:	$LCALL	ZNI		;[4002] zero**neg , result=infinity
	JRST	RST234		;[4002] restore registers 2, 3, and 4.


ZERO:	SETZ	T0,		;Result is (0,0). Clear T0.
				;  T1 already has 0.
	JRST	RST234		;Restore registers 2, 3, and 4.
;
;	>>>>  End of X = Y = 0  <<<<
;
;      >>>>  Main Flow continues at NXTST  <<<<
;
NXTST:	JUMPN	T3,POLAR	;If A and B are not both zero,
	JUMPN	T4,POLAR	;  continue with main flow.
				;  in BLOCK 2.
;
;	>>>>  Special Case: A = B = 0 (X,Y not both 0)  <<<<
;
	MOVE	T0,SPONE	;A = B = 0; X, Y not both zero
	SETZ	T1,		;  so result is (1,0)
RST234:	MOVE	T2,SAVT2	;Restore registers 2
	MOVE	T3,SAVT3	;  and 3
	MOVE	T4,SAVT4	;  and 4
	POPJ	P,		;and return.
;
;	>>>>  End of A = B = 0  <<<<
;
;End of BLOCK 1. The cases corresponding to X = Y = 0 and/or
;A = B = 0 are complete.

;
;			*****************
;			*    BLOCK 2    *
;			*****************
;
;In this block of code X and Y are not both zero. A and B
;are not used. (X + iY) is rewritten in the form
;EXP(CDLOG(X + iY)), and the CDLOG routine is invoked. The
;real part of CDLOG is stored in LNRHO, and its imaginary
;part in THETA. Either, but not both of, THETA and LNRHO
;can underflow, and code is inserted to return scaled values
;instead of zero for these cases. On entry to this block X
;and Y are in T1 and T2. T3 and T4, which contain A and B,
;are used here, but A and B will be recovered later, when
;they are needed,from REXP and IEXP, respectively.
;
POLAR:	MOVEM	T5,SAVT5	;Save another register.
	SETZ	T5,		;Clear T5 for initialization
	MOVEM	T5,RSIGN	;  of flags for signs of the
	MOVEM	T5,ISIGN	;  real and imaginary parts
				;  of the final result.
	MOVEM	T5,PHIFLG	;Initialize flag for
				;  underflow of PHI.
	DMOVEM	T4,IEXP		;B is now zero extended.
				;  Store B in double precision.
	SETZ	T4,		;Zero extend A, and store
	DMOVEM	T3,REXP		;  it in double precision.
	MOVE	T4,T1		;Y in T2, copy X to T4
	SETZ	T3,		;Zero extend Y to double.
				;  Already done for X in T4,T5.
;
;Separate out cases leading to underflow of THETA and
;LNRHO, so as to calculate scaled values. All branches follow
;completion of the main flow for this BLOCK.
;
	MOVE	T0,T2		;Copy of Y to T0.
JOIN3:	FDV	T0,T1		;[3247]S.P. chopped DIV shows
	  JFCL	EXCEP		;  possible exceptions.
	MOVM	T0,T0		;Test |Y/X|.
	CAML	T0,TWO63	;If |Y/X| huge: |THETA| = PI/2.
	  JRST	YMGTX		;  Go to special cases.
	CAMG	T0,TWOM63	;If |Y/X| tiny: |THETA| = 0 or
	 JRST	 XMGTY		;   PI. Go to special cases.
;
;
	DMOVEM	T4,ARG1		;Prepare to get CDLOG
	DMOVEM	T2,ARG1+2	;  of (X,Y). There will be
	FUNCT	CDLOG.,<ARG1,LNRHO>	;  no exceptions.
	JRST	PHI0		;Go to BLOCK 3.0 for PHI. Neither
				;  THETA nor LNRHO is scaled.

;
;Main Flow for BLOCK 2 is complete. Now take care of the
;special cases.
;
EXCEP:	JUMPN	T0,YMGTX	;On overflow, |Y| >> |X|.
;
;	>>>>  Treatment of underflow on Y/X.  <<<<
;
	JUMPL	T4,THETPI	;Underflow and X < 0: THETA is PI.
				;Underflow implies |X| not 1.
				;Also |X| >> |Y|; hence LNRHO
	DMOVEM	T4,LNARG	;  is LN(|X|) which is not 0.
				;  X > 0, so |X| = X. Get
	FUNCT	DLOG.,<LNARG>	;  LN(X) = LN(RHO)

	DMOVEM	T0,LNRHO	;  and store it.
;
				;Since Y/X underflows and X > 0
				;  Theta underflows. Get and
				;  store a scaled value.
	FSC	T2,140		;[3247]Scale Y by 2**96 and X by
	FSC	T4,-140		;[3247]  2**(-96). No exceptions on
	DFDV	T2,T4		;Y/X scaled up by 2**192
	DMOVEM	T2,THETA	;Store scaled THETA.
	JRST	PHI1		;Continue with main flow in
				;  BLOCK 3.2; THETA is scaled
				;  and LNRHO is not.
;
THETPI:	DMOVE	T0,PI		;Underflow on Y/X, and X < 0
	JUMPGE	T2,ABSX		;THETA is PI for Y > 0
	  DMOVN	T0,T0		;  or -PI for Y < 0.
ABSX:	MOVN	T2,T4		;[3247]  T2 <-- |X|.
	JRST	GETLN		;Go store THETA and get LNRHO.
;
;	>>>>  End of underflow on Y/X  <<<<
;
;	>>>>  Treatment of very large |Y/X|  <<<<
;
YMGTX:	MOVM	T1,T2		;Get |Y|. If it is not 1, then
	CAME	T1,SPONE	;  neither THETA nor LNRHO will
	  JRST	XTINY		;[3247]  underflow; Merge with XTINY.
	DMOVE	T0,PIOV2	;If |Y| = 1, LNRHO may underflow.
	JUMPGE	T2,THPOS	;THETA is PI/2
	  DMOVN	T0,T0		;  or -PI/2 if Y < 0. |Y| = 1 so
THPOS:	MOVE	T2,T4		;LNRHO = (X**2)/2. Get X in
	JRST	GETLOG		;  T2 and T4. To GETLOG to store
				;  THETA and get scaled LNRHO.
;
;	>>>>  End of treatment of very large |Y/X|  <<<<

;
;	>>>> Treatment of |Y/X| very small; no underflow  <<<<
;
XMGTY:	JUMPGE	T4,TINYTH	;If X > 0, THETA tiny but not 0.
	MOVN	T4,T4		;If X < 0, change its sign.
	DMOVE	T0,PI		;  THETA is PI
	JUMPGE	T2,ISX1		;  if Y > 0, or
	  DMOVN	T0,T0		;  -PI, if Y < 0.
	JRST	ISX1		;Go test whether |X| = 1.
;
TINYTH:	DMOVE	T0,T2		;Get Y to T0,T1, and calculate
	DFDV	T0,T4		;  DP value for small THETA.
				;  No exceptions on DFDV.
;
ISX1:	MOVM	T2,T2		;|Y| to T2.
	CAMN	T4,SPONE	;If |X| is 1, LNRHO may underflow.
	  JRST	XIS1		;  Go find out.
	DMOVEM	T4,LNARG	;X not 1, so no underflow on LNRHO.
	JRST	STHETA		;Go store THETA and get LNRHO.
;
XIS1:	MOVE	T4,T2		;THETA is not scaled.
	JRST	GETLOG		;Get LNRHO = (Y**2)/2.
;
;	>>>>  END of treatment for very small |Y/X|  <<<<
;
;	>>>>  Special case X = 0 from BLOCK 1  <<<<
;	>>>>  Also branch from |Y/X| very large  <<<<
;
XIS0:	MOVEM	T5,SAVT5	;[3247]Save another register
	SETZ	T5,		;[3247]Zero extend T4, and save DP value
	DMOVEM	T4,IEXP		;[3247]  of imaginary part of exponent.
	SETZ	T4,		;[3247]Zero extend T3, and save DP
	DMOVEM	T3,REXP		;[3247]  value of real part of exponent.
	SETZ	T3,		;[3247]Zero extend Y in T2
	SETZM	PHIFLG		;Initialize flag for underflow of PHI
	SETZM	RSIGN		;Clear flags for signs of 
	SETZM	ISIGN		;  final result
XTINY:	DMOVE	T0,PIOV2	;[3247]If X = 0 or |Y| >> X
	JUMPGE	T2,GETLN	;THETA is PI/2
	  DMOVN	T0,T0		;  or -PI/2 if Y < 0.
	  MOVN	T2,T2		;  T2 <-- |Y|. Now LN(RHO)
				;  will be LN(|Y|). Go store
				;  THETA and get LNRHO. No
				;  underflow on LNRHO because
				;  if X not 0, |Y| is not 1 and
				;  if X = 0, LN(|Y|) cannot underflow.
;
;	>>>>  End of special case X = 0  <<<<
;	>>>>  LNRHO from CALL to LOG routine  <<<<
;
GETLN:	DMOVEM	T2,LNARG	;|X| not 1 and >> |Y|
				;  or |Y| not 1 and >> |X|
STHETA:	DMOVEM	T0,THETA	;Store THETA.
	FUNCT	DLOG.,<LNARG>	;Get LN of larger
	DMOVEM	T0,LNRHO	;  and store it.
	JRST	PHI0		;Continue with main flow in
				;  BLOCK 3.0; neither THETA
				;  not LNRHO is scaled.
;
;	>>>>  End of LNRHO by CALL to LOG routine  <<<<
;
;	>>>> LNRHO for RHO = 1 + (X**2 or Y**2)  <<<<
;
GETLOG:	DMOVEM	T0,THETA	;Store unscaled THETA.
	DFMP	T4,T4		;Get X**2 or Y**2. Copy of X or Y
	  JFCL	LNUND		;  in T2, in case of underflow.
	FSC	T4,-1		;DIV by 2.
	  JFCL	LNUND		;  Underflow may occur.
	DMOVEM	T4,LNRHO	;Store unscaled LNRHO if no
				;  underflows occurred.
	JRST	PHI0		;Continue with main flow in
				;  BLOCK 3.0; neither THETA
				;  not LNRHO is scaled.
;
LNUND:	FSC	T2,140		;[3247]Scale up saved X or Y.
	DFMP	T2,T2		;Square. No exceptions possible.
	FSC	T2,-1		;DIV by 2. No exceptions possible.
	DMOVEM	T2,LNRHO	;Store scaled LNRHO
	JRST	PHI2		;Rejoin main flow at BLOCK 3.1 with
				;  THETA unscaled, LNRHO scaled.

;
;	>>>>  End of LNRHO for RHO near 1  <<<<
;
;End of BLOCK 2. The CDLOG(X + iY) has real part LNRHO and
;imaginary part THETA, which are stored in locations with
;these names. If no exceptions occur, the main flow continues
;in BLOCK 3.0 at PHI0. No overflows occur. Either, but
;not both of THETA and LNRHO may underflow. For these
;cases values scaled up by 2**192 decimal, or 2**300
;octal, have been stored, and the main flow continues
;in BLOCK 3.2 at PHI1 for underflow of THETA, and in
;BLOCK 3.1 at PHI2 for underflow of LNRHO.
;			*****************
;			*    BLOCK 3    *
;			*****************
;
;The complex power function is to be evaluated as the
;complex exponential of (A + iB)*(LNRHO + iTHETA), which
;has real part ALPHA = A*LNRHO - B*THETA, and imaginary
;part PHI = A*THETA + B*LNRHO. The calculation of PHI is carried
;out in BLOCK 3. There are three sub-blocks 3.0, 3.2, and
;3.1 for the cases of no underflow, underflow of THETA and
;underflow of LNRHO, respectively.
;
;Since it is ultimately EXP(i*PHI) that is needed, it would
;appear that SIN and COS of PHI are the final parameters
;needed. However, these functions will get multiplied by
;EXP(ALPHA), and the handling of exception boundaries on
;the product will be expedited by using instead LN(SIN(PHI))
;and LN(COS(PHI)), which will be added to ALPHA before
;invoking the EXP function. Actually is is the absolute
;values of SIN and COS which will be used as arguments to
;the LN function; the signs of SIN and COS will be stored
;in flags for use in determining the signs for the real and
;imaginary parts of the complex exponential. Recall that these
;flags were initialized to zero at the beginning of BLOCK 2.
;
;			*****************
;			*   BLOCK 3.0   *
;			*****************
;
;Neither THETA nor LNRHO is scaled.
;
PHI0:	DMOVE	T0,LNRHO	;Get LNRHO and multiply
	DFMP	T0,IEXP		;  by B.
	  JFCL	PHI0A		;  Branch on exception.
	DMOVE	T2,THETA	;Get THETA and multiply
	DFMP	T2,REXP		;  by A.
	  JFCL	PHI0B		;  Branch on exception
	DMOVE	T4,T0		;Get copy of B*LNRHO in case
				;  of exception on next step.
;
;Rejoin the main flow of BLOCK 3.0 here for underflows
;  that do no harm. No exception on DFAD for such cases.
;
PHISUM:	DFAD	T0,T2		;PHI = A*THETA + B*LNRHO.
	  JFCL	PHI0C		;  Branch on exception.
	DMOVE	T4,T0		;Get copy of PHI and
	JUMPGE	T4,ABSARG	;  test its sign.
	  DMOVN	T4,T4		;  Get absolute value
ABSARG:	CAMGE	T4,ARGHI	;Test |PHI| against limit
	  JRST	ARGOK		;  for argument reduction
	CAMG	T4,ARGHI	;  in SIN/COS routines.
	CAML	T5,ARGLO	;  If >/= INT(PI*2**31)
	  JRST	PHIBIG		;  |PHI| is too big.
ARGOK:	DMOVEM	T0,TRGARG	;Else, store PHI.
	FUNCT	DSIN.,<TRGARG>	;Get SIN(PHI) and
	DMOVEM	T0,SPHI		;  store for future use.
;
	FUNCT	DCOS.,<TRGARG>	;Get COS(PHI) and
	DMOVEM	T0,CPHI		;  store for future use.
;
	JRST	ALF0		;PHI calculation is complete.
				;Go to BLOCK 4.0 for ALPHA.
;Main flow of BLOCK 3.0 is complete. Now take care of
;special cases.
;
;	>>>>  B*LNRHO and A*THETA OK; exception on sum  <<<<
;
PHI0C:	JUMPN	T0,PHIBIG	;If overflow, PHI is too big.
	DMOVE	T0,T4		;Get B*LNRHO back to T0.
	FSC	T0,300		;Scale both B*LNRHO and A*THETA
	FSC	T2,300		;  up by 2**192.
	JRST	SCLPHI		;Join with other cases to
				;  get scaled PHI, etc.
;
;	>>>>  End of handling exception on sum  <<<<
;
;	>>>>  B*LNRHO OK, exception on A*THETA  <<<<
;
PHI0B:	JUMPN	T2,PHIBIG	;If overflow, PHI is too big.
	JUMPE	T0,SCTH		;[3247]Don't "FSC" 0!
	MOVM	T4,T0		;Get |B*LNRHO|. If it is
	CAML	T4,TWOM66	;  > 2**(-66), underflow does
	  JRST	PHISUM		;  no harm. Go to main flow.
	FSC	T0,300		;Else scale B*LNRHO up by 2**192
SCTH:	DMOVE	T2,THETA	;Also get THETA and A and
	DMOVE	T4,REXP		;  scale each up by 2**96
	FSC	T2,140		;[3247]  so that resulting product
	FSC	T4,140		;[3247]  has same scaling as B*LNRHO.
	DFMP	T2,T4		;No exceptions on either the
	JRST	SCLPHI		;  MUL or subsequent operations
				;  in getting scaled PHI, etc.
;
;	>>>>  End of B*LNRHO OK; exception on A*THETA  <<<<
;
;	>>>>  Exception on B*LNRHO  <<<<
;
PHI0A:	JUMPN	T0,PHIBIG	;If overflow, |PHI| too big.
	DMOVE	T2,THETA	;B*LNRHO underflows. Get
	DFMP	T2,REXP		;  A*THETA before scaling.
	  JFCL	PHI0D		;Branch on exception.
	MOVM	T4,T2		;If |A*THETA| is big enough
	CAML	T4,TWOM66	;  main flow can be rejoined
	  JRST	PHISUM		;  with 0 for negligible B*LNRHO
	JUMPE	T2,SCLBLN	;[3247]Don't "FSC" 0!
	FSC	T2,300		;Else scale A*THETA up by 2**192
	JRST	SCLBLN		;  and go scale B*LNRHO
;
;	>>>>  End of exception on B*LNRHO  <<<<
;	>>>>  Underflow on B*LNRHO and exception on A*THETA  <<<<
;
PHI0D:	JUMPN	T2,PHIBIG	;If overflow, PHI is too big.
	DMOVE	T2,THETA	;Both products underflow. Get
	DMOVE	T4,REXP		;THETA and A and scale each
	FSC	T2,140		;[3247]  up by 2**96, so that product
	FSC	T4,140		;[3247]  will be scaled up by 2**192
	DFMP	T2,T4		;No exception on product.
SCLBLN:	DMOVE	T0,LNRHO	;Get LNRHO and B and scale
	DMOVE	T4,IEXP		;  each up by 2**96 so
	FSC	T0,140		;[3247]  that product will be
	FSC	T4,140		;[3247]  scaled up by 2**192.
	DFMP	T0,T4		;No exception on product.
SCLPHI:	DFAD	T0,T2		;No exception on sum.
	HRLZI	T5,400000	;Set sign bit of T5.
	MOVEM	T5,PHIFLG	;Set PHIFLG to indicate
				;  underflow of PHI.
	JUMPGE	T0,PHIPL0	;If scaled PHI < 0
	  MOVEM	T5,ISIGN	;  Use T5 to indicate that
				;  imaginary part is negative.
	  DMOVN	T0,T0		;  Get |PHI| in T0.
PHIPL0:	DMOVE	T0,LNARG	;Store scaled |PHI|
	FUNCT	DLOG.,<LNARG>	;  and get its LOG
	DFAD	T0,LN300	;  and add -192*LN(2) to
	DMOVEM	T0,SPHI		;  get LN(SIN(|PHI|))
	SETZ	T0,		;Clear T0 and T1 in order to
	SETZ	T1,		;  store 0 for LN(COS(|PHI|))
	DMOVEM	T0,CPHI		;  Since COS(TINY) > 0
				;  real part will be > 0.
	JRST	ALF0		;Continue with main flow
				;  at BLOCK 4.0.
;	>>>>  End of underflow on both products  <<<<
;	>>>>  PHI is too large  <<<<
;
PHIBIG:	$LCALL	BPI		;Both parts indeterminate
	HRLOI	T0,377777	;Store biggest number in
	MOVE	T1,T0		;  both REAL and IMAG.
	JRST	RESTOR		;Go to end of BLOCK 6
				;  to restore registers.
;
;****If |PHI| is too large something should be stored
;	in T0,T1, and an error message returned saying that
;	neither the real nor the imaginary parts of the
;	result are determinate. If |PHI| > int (2**31 * pi), the
;	argument reduction for SIN and COS breaks down.
;
;End of BLOCK 3.0. Calculation continues at BLOCK 4.0,
;which also applies to the case in which neither LNRHO
;nor THETA underflows. At this point SPHI contains SIN(PHI)
;or, if PHI underflows, LN(SIN(|PHI|)) = LN(|PHI|). CPHI
;contains COS(PHI) or, if PHI underflows, LN(COS(|PHI|)) = 0.
;If PHI underflows, RSIGN and ISIGN contain the signs of the
;real and imaginary parts of the final result, respectively.
;PHIFLG has 0 or 1 in its sign bit, according as PHI has not
;or has underflowed.
;
;			*****************
;			*   BLOCK 3.1   *
;			*****************
;
;LNRHO is scaled up by 2**192. THETA is not scaled. The
;following transformation makes it possible to use the code
;in BLOCK 3.2 for this case also.
;
PHI2:	HRLZI	T0,400000	;Must exchange LNRHO,THETA
	MOVE	T0,XCHFLG	;  and A,B. Set sign bit of
				;  XCHFLG to indicate that
				;  this has been done.
	DMOVE	T0,LNRHO	;Get LNRHO
	DMOVE	T2,THETA	;  and THETA
	DMOVEM	T0,THETA	;  and exchange
	DMOVEM	T2,LNRHO	;  them.
	MOVE	T0,REXP		;Get A
	MOVE	T2,IEXP		;  and B
	MOVEM	T0,IEXP		;  and exchange
	MOVEM	T2,REXP		;  them.
;
;Continue at PHI1 in BLOCK 3.2, which will now carry out
;the PHI calculation correctly for this case.
;
;			*****************
;			*   BLOCK 3.2   *
;			*****************
;
;THETA is scaled up by 2**192. LNRHO is not scaled.
;
;However, the same code is used for LNRHO scaled and THETA
;not scaled on JRST from BLOCK 3.1. For this case read
;THETA for LNRHO, LNRHO for THETA, IEXP for REXP and REXP
;for IEXP throughout.
;
PHI1:	DMOVE	T0,LNRHO	;Get LNRHO and multiply
	DFMP	T0,IEXP		;  by B.
	  JFCL	PHI1A		;  Branch on exception.
	DMOVE	T2,THETA	;Get THETA and multiply
	DFMP	T2,REXP		;  by A.
	  JFCL	PHI1B		;  Branch on exception.
	JUMPE	T2,PHIOK1	;[3247]Don't "FSC" 0!
	DMOVE	T4,T2		;Save copy of A*THETA
	FSC	T2,-300		;Try to unscale A*THETA.
	  JFCL	T2,PHI1C	;  On failure, branch.
	DMOVE	T4,T0		;Save copy of B*LNRHO.
	DFAD	T0,T2		;PHI = A*THETA + B*LNRHO.
	  JFCL	PHI1D		;  Branch on exception.
;
;Rejoin the main flow of BLOCK 3.2 here for underflows
;  that do no harm. PHI is unscaled and ready for final
;  processing.
;
PHIOK1:	DMOVE	T4,T0		;Get copy of PHI and
	JUMPGE	T4,ABARG1	;  test its sign.
	  DMOVN	T4,T4		;  Get absolute value
ABARG1:	CAMGE	T4,ARGHI	;Test |PHI| against limit
	  JRST	ARGOK1		;  for argument reduction
	CAMG	T4,ARGHI	;  in SIN/COS routines.
	CAML	T5,ARGLO	;  If >/= INT(PI*2**31)
	  JRST	PHIBIG		;  |PHI| is too big.
ARGOK1:	DMOVEM	T0,TRGARG	;Store unscaled PHI
	FUNCT	DSIN.,<TRGARG>	;Get its SIN
	DMOVEM	T0,SPHI		;  and store it.
	FUNCT	DCOS.,<TRGARG>	;Get COS(PHI)
	DMOVEM	T0,CPHI		;  and store it.
	JRST	ALF1		;PHI calculation is complete.
				;Go to BLOCK 4.1 for ALPHA.
;Main Flow of BLOCK 3.2 is complete. Now take care of
;special cases.
;
;	>>>>  B*LNRHO OK, unscaled A*THETA OK  <<<<
;
;	>>>>  Exception on their sum  <<<<
;There was no exception on A*(scaled THETA). Hence the unscaled
;product < 2**(-100), hence sum could not overflow, hence sum
;underflowed. T4 has unscaled B*LNRHO and T2 has scaled A*THETA.
;
PHI1D:	DMOVE	T0,T4		;Get B*LNRHO to T0
	FSC	T0,300		;Scale it up by 2**192.
	FSC	T2,300		;Rescale A*THETA
	DFAD	T0,T2		;No exceptions on scaled PHI.
	JRST	SCPHI1		;Go finish up for scaled PHI.
;
;	>>>>  End of exception on sum of unscaled products  <<<<
;

;	>>>>  B*LNRHO OK, A*(scaled THETA) OK  <<<<
;	>>>>  Unscaled A*THETA underflows  <<<<
;T0 has B*LNRHO and T4 has A*(scaled THETA).
;
PHI1C:	MOVM	T2,T0		;Get |B*LNRHO| in T2
	CAML	T2,TWOM66	;If it is big enough
	  JRST	PHIOK1		;  underflow of A*THETA is
				;  harmless: A*THETA is
				;  negligible. Rejoin main
				;  flow at PHIOK1 above.
	DMOVE	T2,T4		;A*(Scaled THETA) back to T2.
	JUMPE	T0,ADD0		;[3247]Avoid an "FSC" of 0
	JRST	SCBLN1		;Go scale B*LNRHO up and
				;  add to A*(scaled THETA)
				;  in T2 to get scaled PHI.
;
;	>>>>  End of small B*LNRHO + underflowed A*THETA  <<<<
;
;	>>>>  B*LNRHO OK, exception on A*(scaled THETA)  <<<<
;
PHI1B:	JUMPE	T2,PHIOK1	;A*(scaled THETA) underflows,
				;  hence A*THETA is negligible.
				;  Rejoin main flow at PHIOK1 above.
	JRST	REMOV1		;Overflow on A*(scaled THETA) can
				;  be removed by unscaling product.
;
;	>>>>  End of B*LNRHO OK, exception on A*(scaled THETA)  <<<<
;	>>>> Exception on B*LNRHO  <<<<
;
PHI1A: JUMPN	T0,PHIBIG	;If overflow, PHI is too big.
	DMOVE	T2,THETA	;B*LNRHO underflows. Get
	DFMP	T2,REXP		;  THETA and multiply by A.
	  JFCL	PHI1E		;  Branch on exception.
	JRST	SCBLN1		;Go scale B*LNRHO up and add
				;  to scaled A*THETA in T2 to
				;  get scaled PHI.
;
PHI1E:	JUMPE	T2,SCBLN1	;B*LNRHO and A*(scaled THETA)
				;  both underflow. A*THETA is
				;  negligible. Go scale B*LNRHO
				;  to get scaled PHI, with 0 in
				;  T2 for negligible A*THETA.
;
;If A*(scaled THETA) overflows, removal of scaling brings
;unscaled |A*THETA| into range between 2**(-100) to 2**100.
;Hence underflowed B*LNRHO is negligible. Code in REMOV1
;will add in 0 (for negligible B*LNRHO) to unscaled A*THETA
;to get unscaled PHI.
;
;Branch from PHI1B to remove overflow from A*(scaled THETA)
;joins here, with B*LNRHO OK and in T0 to be added to
;unscaled A*THETA.
;
REMOV1:	DMOVE	T2,THETA	;Overflow can be removed by
	DMOVE	T4,REXP		;  scaling THETA and A down by
	FSC	T2,-140		;[3247]  2**(-96) each, to get
	FSC	T4,-140		;[3247]  unscaled product.
	DFMP	T2,T4		;No exceptions.
	DFAD	T0,T2		;No exceptions since |A*THETA| is
				;  between 2**(-100) and 2**100
	JRST	PHIOK1		;Rejoin main flow at PHIOK1 to get
				;  SIN and COS of unscaled PHI.
;
;	>>>>  End handling exception on B*LNRHO  <<<<
;
;	>>>>  Completion of calculation for scaled PHI  <<<<
;
SCBLN1:	DMOVE	T0,LNRHO	;Get LNRHO and B and
	DMOVE	T4,IEXP		;  scale each up by 2**96
	FSC	T0,140		;[3247]  so as to get product
	FSC	T4,140		;[3247]  scaled up by 2**192.
	DFMP	T0,T4		;No exception on product or
ADD0:	DFAD	T0,T2		;  on ADD to scaled A*THETA.
;
;T0 now has scaled PHI. Branches from other special cases
;join in here to get signs of the final real and
;imaginary parts, LN(SIN(|PHI|)) and LN(COS(|PHI|)).
;
SCPHI1:	HRLZI	T4,400000	;Set sign bit of T4
	MOVEM	T4,PHIFLG	;Set flag for PHI scaled.
	JUMPGE	T0,PHIPL1	;If scaled PHI is negative
	DMOVN	T0,T0		;  get |scaled PHI| in T0.
	MOVEM	T4,ISIGN	;  Imaginary part negative.
PHIPL1:	DMOVEM	T0,LNARG	;Store scaled PHI
	FUNCT	DLOG.,<LNARG>	;  and get its LOG.
	DFAD	T0,LN300	;  and subtract 192*LN(2).
	DMOVEM	T0,SPHI		;SPHI gets LN(|PHI|) which
				;  = LN(SIN(|PHI|)).
	SETZ	T2,		;Clear T2 and T3 so as to
	SETZ	T3,		;  store 0 for LN(COS|PHI|))
	DMOVEM	T2,CPHI		;  since COS(TINY) = 1. Real
				;  part of result is positive.
	JRST	ALF1		;Now go get ALPHA for LNRHO
				;  unscaled, THETA scaled.
;
;	>>>>  End of calculation for scaled PHI  <<<<
;
;End of BLOCK 3.2. All paths lead either to the error return
;PHIBIG or to ALF1 for the calculation of ALPHA for the case
;LNRHO unscaled and THETA scaled up by 2**192. If PHI is
;unscaled, PHIFLG is clear and SPHI and CPHI contain SIN(PHI)
;and COS(PHI), respectively. IF PHI is scaled, SPHI contains
;LN(SIN(|PHI|)), CPHI contains LN(COS(|PHI|)) = 0, PHIFLG
;has its sign bit set, and RSIGN and ISIGN contain the signs
;of the real and imaginary parts of the final result.
;
;			*****************
;			*    BLOCK 4    *
;			*****************
;
;The real and imaginary parts of the final answer are
;given by EXP(ALPHA) times COS(PHI) and SIN(PHI),
;respectively. ALPHA = A*LNRHO - B*THETA is calculated
;in BLOCK 4.0 for no scaling of THETA or LNRHO. It is
;calculated in BLOCK 4.1 otherwise. It can be shown that
;if ALPHA overflows, neither the SIN nor the COS of PHI
;can underflow sufficiently to bring EXP(ALPHA) back into
;range. If ALPHA overflows and is positive both parts
;of the result overflow: Code providing a message and
;signed largest numbers for the real and imaginary parts
;is given in EXPOV at the end of BLOCK 4.1. The signs of
;the largest numbers are determined by COS(PHI) and SIN(PHI),
;respectively. If ALPHA overflows and is negative, EXP(ALPHA)
;underflows. The code providing a message and returning 0 for
;both parts is deferred to BLOCK 5, in which similar cases
;occur.
;
;			*****************
;			*   BLOCK 4.0   *
;			*****************
;
;Neither THETA nor LNRHO is scaled. ALPHA = A*LNRHO - B*THETA.
;
ALF0:	DMOVE	T0,LNRHO	;Get LNRHO and multiply
	DFMP	T0,REXP		;  by A.
	  JFCL	ALF0A		;  Branch on exception.
	DMOVE	T2,THETA	;Get THETA and multiply
	DFMP	T2,IEXP		;  by B.
	  JFCL	ALF0B		;  Branch on exception
	DFSB	T0,T2		;ALPHA = A*LNRHO - B*THETA.
	  JFCL	ALF0C		;  Branch on exception.
	DMOVEM	T0,ALPHA	;Store ALPHA for future use.
;
	JRST	ANS0		;Calculation of ALPHA complete.
				;Go to BLOCK 5.0 for final answer.
;
;Main flow of BLOCK 4.0 is complete. Now take care of
;special cases.
;
;	>>>>  A*LNRHO and B*THETA OK; exception on  difference  <<<<
;
ALF0C:	JUMPN	T0,DIFFOV	;Exception on final results.
	DMOVEM	T0,ALPHA	;Zero is acceptable if
				;  ALPHA underflows.
	JRST	ANS0		;Calculation of ALPHA complete.
				;Go to BLOCK 5.0 for final answer.
;
DIFFOV:	JUMPL	T2,EXPOV	;Overflow with negative B*THETA
				;  means overflow on final results.
	JRST	EXPUND		;Underflow on final results
				;  for B*THETA positive.
;
;	>>>>  End of handling exception on difference  <<<<
;
;	>>>>  Exception on A*LNRHO  <<<<
;
ALF0A:	JUMPN	T0,ALNOV	;Exception on final reaults.
	DMOVE	T2,THETA	;A*LNRHO underflows. Get
	DFMP	T2,REXP		;  B*THETA.
	  JFCL	ALF0B		;Branch on exception.
	DMOVN	T2,T2		;Underflow on A*LNRHO means it
	DMOVEM	T2,ALPHA	;  is negligible. ALPHA = -B*THETA.
	JRST	ANS0		;Go to BLOCK 5 for final answer.
;
ALNOV:	DMOVE	T2,THETA	;A*LNRHO overflows. What about
	DFMP	T2,IEXP		;  B*THETA?
	  JFCL	BOTHEX		;  On exception to BOTHEX.
	JRST	ALNOV1		;A*LNRHO overflows, B*THETA doesn't.
;
;A*LNRHO overflows and B*THETA has exception.
;
BOTHEX:	JUMPE	T2,ALNOV1	;B*THETA underflows.
	DMOVE	T0,LNRHO	;Both products overflow.
	DMOVE	T2,REXP		;  Get A and LNRHO and scale
	FSC	T0,-140		;[3247]  each down to get scaled
	FSC	T2,-140		;[3247]  A*LNRHO.
	DFMP	T0,T2		;No exception
	DMOVE	T2,THETA	;Get B and THETA
	DMOVE	T4,IEXP		;  and scale each down
	FSC	T2,-140		;[3247]  so get B*THETA
	FSC	T4,-140		;[3247]  scaled like B*LNRHO
	DFMP	T2,T4		;No exception.
	DFSB	T0,T2		;No exception on scaled ALPHA.
	JUMPL	T0,EXPUND	;Results underflow if ALPHA < 0.
	JRST	EXPOV		;Results overflow if ALPHA > 0.
;
;End of overflow on both A*LNRHO and B*THETA.
;
;A*LNRHO overflows, B*THETA does not. Sign of ALPHA is
;determined by sign of A*LNRHO. Neither A nor LNRHO is 0.
;
;NOTE: A similar case from BLOCK 4.1 joins in here.
;
ALNOV1:	SKIPL	T0,LNRHO	;If LNRHO > 0
	  JRST	LNPOS		;  go test B.
	SKIPL	T0,REXP		;LNRHO < 0. If A
	  JRST	ALNNEG		;  > 0, A*LNRHO < 0.
	JRST	ALNPOS		;Else A*LNRHO > 0.
LNPOS:	SKIPL	T0,REXP		;LNRHO > 0. If A
	  JRST	ALNPOS		;  > 0, A*LNRHO > 0.
ALNNEG:	JRST	EXPUND		;A*LNRHO < 0 means
				;  final results underflow.
ALNPOS:	JRST	EXPOV		;A*LNRHO > 0 means
				;  final results overflow.
;
;	>>>>  End of exception on A*LNRHO  <<<<
;
;	>>>>  A*LNRHO OK or underflows, exception on B*THETA  <<<<
;
ALF0B:	JUMPN	T2,BTHOV	;Exception on final results.
	DMOVEM	T0,ALPHA	;Underflow on B*THETA means it
				;  is negligible. ALPHA = A*LNRHO.
				;  Note that if A*LNRHO underflows
				;  it is also negligible and 0 is
				;  an acceptable result for ALPHA.
	JRST	ANS0		;Go to BLOCK 5 for final answer.
;
;On overflow of B*THETA, neither B nor THETA is 0. A*LNRHO OK
;or underflows so sign of ALPHA is negative of sign
;of overflowed B*THETA.
;
BTHOV:	SKIPL	T0,THETA	;If THETA > 0
	  JRST	THPOS1		;  go test B.
	SKIPL	T0,IEXP		;THETA < 0. If B
	  JRST	BTHNEG		;  > 0, B*THETA < 0.
	JRST	BTHPOS		;Else B*THETA > 0.
THPOS1:	SKIPL	T0,IEXP		;THETA > 0. If B
	  JRST	BTHPOS		;  > 0, B*THETA > 0.
BTHNEG:	JRST	EXPOV		;B*THETA < 0 means -B*THETA > 0,
				;  so final results overflow.
BTHPOS:	JRST	EXPUND		;B*THETA > 0 means -B*THETA < 0,
				;  so final results underflow.
;
;	>>>>  End of A*LNRHO OK; exception on B*THETA  <<<<
;
;End of BLOCK 4.0. All paths lead to EXPOV, at the end
;of BLOCK 4.1 or to ANS0 or EXPUND, both in BLOCK 5.
;
;			*****************
;			*   BLOCK 4.1   *
;			*****************
;
;If PHIFLG = 0, THETA is scaled up by 2**192 and LNRHO is not
;scaled. If PHIFLG < 0,LNRHO is scaled up by 2**192, and THETA
;is not scaled.
;
;The code, as written, applies to PHIFLG = 0.
;
;For PHIFLG < 0, read  THETA for LNRHO, LNRHO for THETA, (-IEXP)
;for REXP and (-REXP) for IEXP throughout.
;

ALF1:	SKIPL	PHIFLG		;If PHIFLG set,
	  JRST	ALF2		;  to code for THETA scaled.
	MOVN	T0,REXP		;Else, exchange already done
	MOVEM	T0,REXP		;  for THETA,LNRHO and A,B.
	MOVN	T0,IEXP		;  Negate exchanged A,B to use
	MOVEM	T0,IEXP		;  the code for scaled LNRHO,
				;  instead of scaled THETA.
;
ALF2:	DMOVE	T0,LNRHO	;Get LNRHO and multiply
	DFMP	T0,REXP		;  by A.
	  JFCL	ALF2A		;  Branch on exception.
	DMOVE	T2,THETA	;Get THETA and multiply
	DFMP	T2,IEXP		;  by B.
	  JFCL	ALF2B		;  Branch on exception
	JUMPE	T2,SUB0		;[3247]Don't "FSC" 0
	DMOVE	T4,T2		;Save a copy of B*THETA
	FSC	T2,-300		;  and try to remove scaling.
	  JFCL	STALF		;  Didn't work: B*THETA
				;  underflows and hence is
				;  negligible. ALPHA is
				;  A*LNRHO in To.
SUB0:	DFSB	T0,T2		;[3247]Neither product has an exception.
				;ALPHA = A*LNRHO - B*THETA.
	  JFCL	ALF2C		;  Branch on exception.
STALF:	DMOVEM	T0,ALPHA	;Store ALPHA for future use.
	JRST	ANS0		;Calculation of ALPHA complete.
				;Go to BLOCK 5.0 for final answer.
;
;
;Main flow of BLOCK 4.1 is complete. Now take care of
;special cases.
;
;	>>>>  A*LNRHO and B*THETA OK, exception on difference  <<<<
;
ALF2C:	JUMPE	T0,STALF1	;On underflow 0 OK for ALPHA.
	JUMPL	T2,EXPOV	;On overflow ALPHA has sign
	JRST	EXPUND		;  opposite to B*THETA. Hence
				;  results overflow for T2 < 0
				;  and underflow for T2 > 0.
;
;	>>>>  End of exception on difference  <<<<
;
;	>>>>  Exception on A*LNRHO  <<<<
;
ALF2A:	JUMPN	T0,ALNOV1	;A*LNRHO overflows, B*THETA
				;  will be negligible. Final
				;  results will overflow
				;  or underflow. Go merge with
				;  similar case in BLOCK 4.0.

	DMOVE	T2,THETA	;A*LNRHO underflows and is
	DFMP	T2,IEXP		;  negligible. Get B*(scaled
	  JFCL	ALF2D		;  THETA). Branch on exception.
	JUMPE	T2,STALF1	;[3247]Don't "FSC" 0!
	FSC	T2,-300		;Try to unscale B*THETA.
	  JFCL			;  Failed: B*THETA underflows too.
				;  Hence 0 is acceptable for ALPHA.
STALF1:	DMOVEM	T0,ALPHA	;Store 0 for ALPHA and go to
	JRST	ANS0		;  BLOCK 5 for final answer.
;
;Exception on B*(scaled THETA); A*LNRHO underflows.
;
ALF2D:	JUMPN	T2,REMOV2	;Go unscale overflowed
				;  B*(scaled THETA). Negligible
				;  A*LNRHO stored as 0 in T0.
	DMOVEM	T0,ALPHA	;Both products underflow so
	JRST	ANS0		;  0 is acceptable for ALPHA
				;  To BLOCK 5 for final answer.
;
;	>>>>  End of Exception on A*LNRHO  <<<<
;
;	>>>>  A*LNRHO OK, exception on B*(scaled THETA)  <<<<
;
ALF2B:	JUMPN	T2,REMOV2	;Overflow: unscale A*THETa.
	DMOVEM	T0,ALPHA	;Underflow: A*THETA negligible.
				;  ALPHA = A*LNRHO (or 0 if
				;  here from ALF2C).
	JRST	ANS0		;ALPHA done. Go to BLOCK 5.
;
REMOV2:	DMOVE	T2,THETA	;Get scaled THETA and B
	DMOVE	T4,IEXP		;  and scale each down by
	FSC	T2,-140		;[3247]  2**-96 to remove scaling
	FSC	T4,-140		;[3247]  from product.
	DFMP	T2,T4		;|B*THETA| in 2**(-100) to 2**100
	DFSB	T0,T2		;  so no exception on DFSB.
	DMOVEM	T0,ALPHA	;Store ALPHA and go to
	JRST	ANS0		;  to get final answer.
;
;	>>>>  End of exception on A*(scaled THETA)  <<<<
;
;	>>>>  Overflow of EXP(ALPHA)  <<<<
;
EXPOV:	HRLOI	T0,377777	;Get biggest number in T0.
	SKIPN	T0,PHIFLG	;If PHI underflowed get
	  JRST	EXPOV1		;  signs from ISIGN and RSIGN.
	DMOVE	T2,SPHI		;Get SIN(PHI)
	JUMPN	T2,EXPOV2	;If SIN(PHI) is not 0
	JUMPN	T3,EXPOV2	;  both parts overflow.
	SETZ	T1,		;If SIN(PHI) = 0, so does
				;  imaginary part. Real part
	SKIPGE	T0,CPHI		;  overflows. If COS(PHI) < 0
	  MOVN	T0,T0		;  real part is negative.
	$LCALL	RPO		;Real part overflow
	JRST	RESTOR		;Go restore registers and exit.
;
EXPOV2:	MOVE	T1,T0		;Get Biggest in T1 too.
	SKIPGE	T0,CPHI		;If COS(PHI) < 0, so
	  MOVN	T0,T0		;  is real part.
	SKIPGE	T0,SPHI		;If SIN(PHI) < 0, so
	  MOVN	T1,T1		;  so is imaginary part.
	JRST	ERRMSG		;Go get error message.
;
EXPOV1:	MOVE	T1,T0		;Get Biggest in T1 too.
	SKIPGE	T0,RSIGN	;If real part is negative
	  MOVN	T0,T0		;  negate biggest number.
	SKIPGE	T0,ISIGN	;If imaginary part is negative
	  MOVN	T1,T1		;  negate biggest number.
ERRMSG:	$LCALL	BPO		;Both parts overflow
	JRST	RESTOR		;Go restore registers and exit.
;
;RESTOR for restoring registers is at the end of BLOCK 5.
;
;	>>>>  End of overflow of EXP(ALPHA)  <<<<
;
;End of BLOCK 4.1. All paths lead to ANS0 or EXPUND, both
;in BLOCK 5, or to RESTOR at the end of BLOCK 6.

;			*****************
;			*    BLOCK 5    *
;			*****************
;
;The main flow for this BLOCK starts at ANS0. If PHIFLG
;is 0, SPHI and CPHI contain SIN(PHI) and COS(PHI). PHI
;did not underflow. ALPHA is in ALPHA; no tests on its size
;have yet been made. The special cases, including PHIFLG = 1
;start at ANS1, following completion of the main flow.
;
ANS0:	DMOVE	T0,ALPHA	;Get ALPHA. Will EXP(ALPHA)
	CAMLE	T0,AMINHI	;  underflow? If ALFHI > AMINHI
	  JRST	EXPOK		;  EXP(ALPHA) won't underflow.
	CAML	T0,AMINHI	;If ALFHI < AMINHI, EXP(ALPHA)
				;  underflows, and so do both
				;  parts of answer.
	CAMG	T1,AMINLO	;[3247]ALFHI = AMINHI; test LO's
	  JRST	EXPUND		;  ALFLO </= AMINLO: underflow.
EXPOK:	SKIPGE	T0,PHIFLG	;Was PHI scaled? If so
	  JRST	ANS2		;  go to alternate code.
	CAML	T0,AMAXHI	;If ALFHI >/= AMAXHI
	  JRST	ANS1		;  use LOG approach.
	FUNCT	DEXP.,<ALPHA>	;No exception on EXP(ALPHA).
	DMOVE	T2,T0		;Keep copy of EXP(ALPHA)
	DFMP	T0,SPHI		;Get imaginary part; no overflow
	  JFCL	IMUND		;  get message on underflow.
STIMP:	DMOVEM	T0,IMAG		;Store imaginary part.
	DFMP	T2,CPHI		;Get real part, no overflow
	  JFCL	RLUND		;  get message on underflow
STRLP:	DMOVEM	T2,REAL		;Store real part.
	JRST	DPTOSP		;Go convert to single precision.
				;  in BLOCK 6. Main Flow of
				;  BLOCK 5 complete.
;
;	>>>>  Underflow on MUL by SIN(PHI)  <<<<
;
IMUND:	$LCALL	IPU		;Imaginary part underflow
	JRST	STIMP		;Rejoin main flow above.
;
;	>>>>  End of underflow on MUL by SIN(PHI)  <<<<
;
;	>>>>  Underflow on MUL by COS(PHI)  <<<<
;
RLUND:	$LCALL	RPU		;Real part underflow
	JRST	STRLP		;Rejoin main flow above.
;
;	>>>>  End of underflow on MUL by COS(PHI)  <<<<
;
;	>>>>  Underflow on EXP(ALPHA)  <<<<
;
;Multiplication by SIN(PHI) and COS(PHI) can only aggravate
;the underflow of EXP(ALPHA). Hence both real and imaginary
;parts underflow.
;
;NOTE: Other cases for underflow of EXP(ALPHA) from
;BLOCK 4 join here.
;
EXPUND:	;Write a message saying that EXP(ALPHA) underflows
	;and hence both parts of answer also underflow.
	$LCALL	BPU		;Both parts underflow
	SETZ	T0,		;Store zero for real part.
	SETZ	T1,		;Store zero for imaginary part.
	JRST	RESTOR		;Go to end of BLOCK 6 and
				;  restore registers 2 through 5.
;
;	>>>>  End of underflow on EXP(ALPHA)  <<<<
;
;	>>>>  Overflow on EXP(ALPHA)  <<<<
;
;On entry to the following code, SPHI and CPHI contain
;SIN and COS of PHI. Multiplication by SIN and/or COS
;might compensate the overflow expected on EXP(ALPHA).
;Hence get LN(SIN(|PHI|)) and LN(|COS(PHI)|) and set the
;sign flags for the real and imaginary parts of the
;results. Then merge at ANS2 with the case for underflow
;of PHI (PHIFLG = 1), where this has already been done, to
;complete the calculation.
;
ANS1:	SKIPL	T2,SPHI		;Get SPHI, and if >/= 0
	  JRST	PLUSI		;  no changes needed.
	MOVSI	T5,400000	;Set sign bit of T5
	MOVEM	ISIGN		;  and copy to ISIGN.
	DMOVN	T2,T2		;Get |SIN(PHI)|
PLUSI:	DMOVEM	LNARG		;  and copy to memory.
	FUNCT	DLOG.,<LNARG>	;Get LN(|SIN(PHI)|)
	DMOVEM	T0,SPHI		;  and store in SPHI.
;
	SKIPL	T2,CPHI		;Get CPHI, and if >/= 0
	  JRST	PLUSR		;  no changes needed.
	MOVSI	T5,400000	;Else set sign bit of T5
	MOVEM	RSIGN		;  and copy to RSIGN.
	DMOVN	CPHI		;Get |COS(PHI)|)
PLUSR:	DMOVEM	LNARG		;  and copy to memory.
	FUNCT	DLOG.,<LNARG>	;Get LN(|COS(PHI)|)
	DMOVEM	T0,CPHI		;  and store in CPHI.
;
;	>>>>  End of overflow on EXP(ALPHA)  <<<<
;
;	>>>>  Code for PHIFLG negative  <<<<
;
;The code for ANS1 merges here. SPHI and CPHI contain,
;respectively, LN(SIN(|PHI|)) and LN(|COS(PHI)|). ISIGN
;and RSIGN contain the signs for the imaginary and real
;parts, and ALPHA contains ALPHA, which has already had
;the case for underflow of EXP(ALPHA) separated out. No
;test for overflow has yet been made. However, the
;magnitudes of the real and imaginary parts will be
;calculated as
;
;	EXP(ALPHA + SPHI) and EXP(ALPHA + CPHI)
;
;respectively. Tests for overflow on EXP will be postponed
;until after ALPHA + SPHI and ALPHA + CPHI are calculated.
;
;	>>>> Calculate EXP(ALPHA + SPHI)  <<<<
;
ANS2:	DMOVE	T0,ALPHA	;Get ALPHA and make
	DMOVE	T2,ALPHA	;  a copy.
	DFAD	T0,SPHI		;Add in SPHI. On exception
	  JFCL	IMEX		;  branch to special case.
	CAMGE	T0,AMAXHI	;If T0 is not too big,
	  JRST	NOVIM		;  no overflow on imag. part.
	CAMG	T0,AMAXHI	;If T0 > AMAXHI, overflow
				;  on imaginary part.
	CAML	T1,AMAXLO	;T0 = AMAXHI. Test LO's.
	  JRST	IMOV		;  Overflow for T1 >/= AMAXLO.
NOVIM:	CAMLE	T0,AMINHI	;Next test for EXP underflow:
	  JRST	IMOK		;  If T0 > AMINHI, no underflow.
	CAML	T0,AMINHI	;If T0 </= AMINHI, EXP will
				;  underflow, and also imag. part.
	CAMG	T1,AMINLO	;T0 = AMINHI. If T1 </= AMINLO
	  JRST	IMUND1		;  imaginary part underflows.
IMOK:	DMOVE	T0,EXPARG	;Copy argument for EXP
	FUNCT	DEXP.,<EXPARG>	;  and invoke EXP function.
STIMP1:	SKIPGE	T0,ISIGN	;No exceptions. If ISIGN < 0
	  DMOVN	T0,T0		;  so is imaginary part.
	DMOVEM	T0, IMAG	;Store imaginary part.
	JRST	ANS3		;Now get real part
;
;	>>>>  End of main flow for EXP(ALPHA + SPHI)  <<<<
;
;	>>>>  Exceptions on EXP(ALPHA + SPHI)  <<<<
;
IMEX:	JUMPE	T0,IMAG1	;Underflow on sum: |imag. part| = 1.
	JUMPL	T2,IMUND1	;On overflow ALPHA and SPHI have
				;  same signs. EXP underflows
				;  if they are negative.
IMOV:	$LCALL	IPO		;Imaginary part overflow
	HRLOI	T0,377777	;Store zero extended
	SETZ	T1,		;  biggest number.
	JRST	STIMP1		;Return to main flow in ANS2 to
				;  get correct sign, and store.

IMUND1:	$LCALL	IPU		;Imaginary part underflow
	SETZ	T0,		;Store double precision
	SETZ	T1,		;  zero for imaginary part
	DMOVEM	T0,IMAG		;  and return to main flow
	JRST	ANS3		;  in ANS2 to get real part.

IMAG1:	MOVE	T0,SPONE	;Store double precision
	SETZ	T1,		;  unity for imaginary part
	JRST	STIMP1		;Return to main flow in ANS2 to
				;  get correct sign, and store.
;
;	>>>>  End of exceptions on imaginary part  <<<<
;
;	>>>>  Continue main flow: Get real part  <<<<
;
ANS3:	DMOVE	T0,ALPHA	;Get ALPHA and make
	DMOVE	T2,ALPHA	;  a copy.
	DFAD	T0,CPHI		;Add in CPHI. On exception
	  JFCL	RLEX		;  branch to special case.
	CAMGE	T0,AMAXHI	;If T0 is not too big, no
	  JRST	NOVRL		;  overflow on real part.
	CAMG	T0,AMAXHI	;If T0 > AMAXHI, overflow
				;  on real part.
	CAML	T1,AMAXLO	;T0 = AMAXHI. Test LO's.
	  JRST	RLOV		;  Overflow for T0 >/= AMAXLO
NOVRL:	CAMLE	T0,AMINHI	;Next test for EXP underflow:
	  JRST	RLOK		;  If T0 > AMINHI, no underflow.
	CAML	T0,AMINHI	;If T0 < AMINHI, EXP will
		;  underflow, and also real part.
	CAMG	T1,AMINLO	;T0 = AMINHI. If T1 </= AMINLO
	  JRST	RLUND1		;  real part underflows.
RLOK:	DMOVE	T0,EXPARG	;Copy argument for EXP
	FUNCT	DEXP.,<EXPARG>	;  and invoke EXP function.
STRLP1:	SKIPGE	T0,RSIGN	;No exceptions. If RSIGN < 0
	  DMOVN	T0,T0		;  so is real part.
	DMOVEM	T0, REAL	;Store real part.
	JRST	DPTOSP		;Go to BLOCK 6 to convert
				;  to single precision.
;
;	>>>>  Exceptions on EXP(ALPHA + CPHI)  <<<<
;
RLEX:	JUMPE	T0,REAL1	;Underflow on sum: |real part| = 1.
	JUMPL	T2,RLUND1	;On overflow ALPHA and CPHI have
				;  same signs. EXP underflows
				;  if they are negative.
RLOV:	$LCALL	RPO		;Real part overflow
	HRLOI	T0,377777	;Store zero extended
	SETZ	T1,		;  biggest number.
	JRST	STRLP1		;Return to main flow in ANS3 to
				;  get correct sign, and store.
;
RLUND1:	$LCALL	RPU		;Real part underflow
	SETZ	T0,		;Store double precision
	SETZ	T1,		;  zero for real part
	DMOVEM	T0,REAL		;  and continue with main
	JRST	DPTOSP		;  flow in BLOCK 6 to convert
				;  to single precision.
;
REAL1:	MOVE	T0,SPONE	;Store double precision
	SETZ	T1,		;  unity for real part
	JRST	STRLP1		;Return to main flow in ANS3 to
				;  get correct sign, and store.
;
;	>>>>  End of exceptions on real part  <<<<
;
;End of BLOCK 5. All exceptions in the double precision
;calculations have been properly handled. All paths lead to
;DPTOSP, in BLOCK 6, for conversion to single precision.
;
;			*****************
;			*    BLOCK 6    *
;			*****************
;
DPTOSP:	DMOVE	T0,REAL		;Get D.P. real part to T0.
	JUMPL	T0,RNEG		;Branch if real part < 0.
	TLNE	T1,(1B1)	;Is rounding bit 1?
	TRON	T0,1		;Yes. If possible round by
	JRST	IMCVT		;  setting LSB. Go convert
				;  imaginary part.
	MOVE	T1,T0		;Else, get copy of high word.
	AND	T0,[777000,,1]	;Construct unnormalized LSB
				;  with same exponent.
	FAD	T0,T1		;Round DP to SP
	  JFCL	RCVTX1		;Branch on exception.
	JRST	IMCVT		;Go convert imaginary part.
;
RNEG:	DMOVN	T0,T0		;REAL < 0; get magnitude.
	TLNE	T1,(1B1)	;Is rounding bit 1?
	TRON	T0,1		;Yes. If possible round by
				;  setting LSB. Go restore
	  JRST	MINUS1		;  minus sign.
	MOVE	T1,T0		;Else, get copy of high word.
	AND	T0,[777000,,1]	;Construct unnormalized LSB
	FAD	T0,T1		;Round DP to SP
	  JFCL	RCVTX2		;Branch on exception.
;
MINUS1:	MOVN	T0,T0		;Restore minus sign.
	JRST	IMCVT		;Now convert imaginary part.
;
;	>>>>  Exceptions on conversion of real part.  <<<<
;
RCVTX1:	JUMPE	T0,RCVTUN	;Branch on underflow.
;	Write message that real part overflows on conversion
	$LCALL	RPO		;Real part overflow
	HRLOI	T0,377777	;Store positive biggest for real.
	JRST	IMCVT		;  and go convert imaginary part.
;
RCVTX2:	JUMPE	T0,RCVTUN	;Branch on underflow.
;	Write message that real part overflows on conversion.
	$LCALL	RPO		;Real part overflow
	HRLOI	T0,400000	;Store negative biggest
	TRO	T0,(1B1)	;  for real part.
	JRST	IMCVT		;Now convert imaginary part.
;
RCVTUN:	SETZ	T0,		;Zero for underflow of real part.
;	Write message that real part underflows on conversion.
	$LCALL	RPU		;Real part underflow
	JRST	IMCVT		;Now convert imaginary part.
;
;			*****************
;			*   BLOCK 6.1   *
;			*****************
;
IMCVT:	DMOVE	T1,IMAG		;Get D.P. imaginary part to T1.
	JUMPL	T1,INEG		;Branch if imaginary part < 0.
	TLNE	T2,(1B1)	;Is rounding bit 1?
	TRON	T1,1		;Yes. If possible round by
	JRST	RESTOR		;  setting LSB. Go restore
				;  registers and exit.
	MOVE	T2,T1		;Else, get copy of high word.
	AND	T1,[777000,,1]	;Construct unnormalized LSB
				;  with same exponent.
	FAD	T1,T2		;Round DP to SP
	  JFCL	ICVTX1		;Branch on exception.
	JRST	RESTOR		;Go restore registers and exit.
;
INEG:	DMOVN	T1,T1		;IMAG < 0; get magnitude.
	TLNE	T2,(1B1)	;Is rounding bit 1?
	TRON	T1,1		;Yes. If possible round by
				;  setting LSB. Go restore
	  JRST	MINUS2		;  minus sign.
	MOVE	T2,T1		;Else, get copy of high word.
	AND	T1,[777000,,1]	;Construct unnormalized LSB
	FAD	T1,T2		;Round DP to SP
	  JFCL	ICVTX2		;Branch on exception.
;
MINUS2:	MOVN	T1,T1		;Restore minus sign.
	JRST	RESTOR		;Go restore registers and exit.
;
;	>>>>  Exceptions on conversion of imaginary part.  <<<<
;
ICVTX1:	JUMPE	T1,ICVTUN	;Branch on underflow.
;	Write message that imaginary part overflows on conversion
	$LCALL	IPO		;Imaginary part overflow
	HRLOI	T1,377777	;Store positive biggest for
	JRST	RESTOR		;  imaginary part, go restore
				;  registers, and exit.
;
ICVTX2:	JUMPE	T1,ICVTUN	;Branch on underflow.
;	Write message that imaginary part overflows on conversion.
	$LCALL	IPO		;Imaginary part overflow
	HRLOI	T1,400000	;Store negative biggest
	TRO	T1,(1B1)	;  for imaginary part.
	JRST	RESTOR		;Go restore registers and exit.
;
ICVTUN:	SETZ	T1,		;Zero for underflow of
				;  imaginary part.
;	Write message that imaginary part underflows on conversion.
	$LCALL	IPU		;Imaginary part underflow

RESTOR:	MOVE	T2,SAVT2	;Restore registers T2,
	MOVE	T3,SAVT3	;  T3,
	MOVE	T4,SAVT4	;  T4,
	MOVE	T5,SAVT5	;  and T5.
;
	GOODBY	(1)		;Return.

;Constants

PI:	EXP	202622077325,021026430215 ;pi
PIOV2:	EXP	201622077325,021026430215 ;pi/2
LN300:	EXP	210412126217,316725563320 ;192 * ln(2) for LOG of
					  ;  underflowed PHI
SPONE:	EXP	201400000000		  ;1.0

ARGHI:	EXP	241622077325		  ;Double precision threshold
ARGLO:	EXP	020000000000		  ;  for argument reduction in
					  ;  SIN/COS routines.
TWO63:	EXP	300400000000		  ;2**63  Threshold for |THETA| = PI/2
TWOM63:	EXP	102400000000		  ;2**-63 Threshold for |THETA| = PI
					  ;	    or THETA = Y/X
TWOM66:	EXP	077400000000		  ;2**-66 Threshold for underflow
					  ;	    negligible.
AMAXHI:	EXP	207540074636		  ;88. ... double precision
AMAXLO:	EXP	077042573256		  ;   overflow boundary for EXP.
AMINHI:	EXP	570232254036		  ;-88. ... double precision
AMINLO:	EXP	301750634730		  ;   underflow boundary for EXP.


;Data
	SEGMENT	DATA
LNRHO:	BLOCK	2		;LNRHO and THETA must be
THETA:	BLOCK	2		;  contiguous since they 
				;  provide output for CDLOG.
REXP:	BLOCK	2		;Temporary storage for A
IEXP:	BLOCK	2		;Temporary storage for B
CPHI:	BLOCK	2		;For COS(PHI) or LN(COS(|PHI|)) 
SPHI:	BLOCK	2		;For SIN(PHI) or LN(SIN(|PHI|))
ALPHA:	BLOCK	2		;Storage for ALPHA.
ARG1:	BLOCK	4		;Temporary storage for inputs to CDLOG.
LNARG:	BLOCK	2		;Temporary storage for inputs
TRGARG:	BLOCK	2		;  to DLOG, DSIN, and
EXPARG:	BLOCK	2		;  DCOS and DEXP routines.
REAL:	BLOCK	2		;Real part of result
IMAG:	BLOCK	2		;Imaginary part of result

RSIGN:	BLOCK	1		;Storage for signs of real and
ISIGN:	BLOCK	1		;  imaginary parts of result.
PHIFLG:	BLOCK	1		;Storage for flag indicating
				;  underflow of PHI.
XCHFLG:	BLOCK	1		;Storage for flag indicating
				;  exchange of LNRHO,THETA and A,B.

SAVT2:	BLOCK	1		;SAVT2, SAVT3, SAVT4,
SAVT3:	BLOCK	1		;  must be contiguous. The
SAVT4:	BLOCK	1		;  corresponding registers are
SAVT5:	BLOCK	1		;  restored by DMOVEs.


	PRGEND
TITLE	CLOG	COMPLEX NATURAL LOG FUNCTION
;		(SINGLE PRECISION FLOATING POINT)

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 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	CLOG
EXTERN	CLOG.
CLOG=CLOG.
PRGEND
TITLE	CLOG.	COMPLEX NATURAL LOG FUNCTION
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	Mary Payne	16-Oct-81

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

;CLOG(Z) IS CALCULATED AS FOLLOWS

;  LET Z = X + I*Y

;  IF Y IS NOT ZERO, THEN
;    CLOG(Z) = U + I*V
;      U = (1/2)*ALOG(X**2 + Y**2)
;      V = ATAN2(Y,X)

;  IF X IS NOT ZERO AND Y IS ZERO, THEN
;    IF X IS POSITIVE, CLOG(Z) = ALOG(X) + I*0
;    IF X IS NEGATIVE, CLOG(Z) = ALOG(ABS(X)) + I*PI

;  IF X AND Y ARE ZERO, THEN
;    CLOG(Z) = -MACHINE INFINITY + I*0
;    AND AN ERROR MESSAGE IS RETURNED

;  THE APPROXIMATION USED FOR ALOG IS FORMULA 2662 FROM
;    "COMPUTER APPROXIMATIONS" BY HART ET AL.

;THE REAL PART OF THE RESULT IS RETURNED IN REGISTER T0.
;THE IMAGINARY PART OF THE RESULT IS RETURNED IN REGISTER T1.

;REQUIRED (CALLED) ROUTINES:  ALOG, ATAN, AND ATAN2

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

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  XMOVEI	L,ARG
;  PUSHJ	P,CLOG

;THE RESULT IS RETURNED IN ACCUMULATORS T0 AND T1

	EXTERN	SNG.0

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(CLOG,.)	;ENTRY TO CLOG ROUTINE
	DMOVE	T0,@(L)		;REAL PART OF ARGUMENT TO T0
				;  IMAGINARY PART TO T1

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

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

	PUSH	P,T2		;SAVE REGISTERS T2
	PUSH	P,T3		;  AND T3
	MOVM	T2,T0		;GET ABSOLUTE VALUES OF
	MOVM	T3,T1		;  X AND Y
	AND	T2,MASK		;ISOLATE EXPONENT FIELDS
	AND	T3,MASK		;  OF X AND Y
	SUB	T3,T2		;EXPONENT OF Y - EXPONENT OF X
	CAML	T3,P14		;IF DIFFERENCE > 14
	  JRST	BIG		;  |Y/X| IS LARGE.
	CAMGE	T3,M14		;IF DIFFERENCE < -14
	  JRST	SMALL		;  |Y/X| IS SMALL

	PUSH	P,T4		;SAVE REGISTERS T4
	PUSH	P,T5		;  AND T5
	MOVM	T4,T0		;SAVE MAGNITUDES OF X
	MOVM	T5,T1		;  AND Y IN T4,T5
	MOVEM	T0,TT0		;STORE X AND Y FOR CALL
	MOVEM	T1,TT2		;  TO ATAN2.
	FUNCT	ATAN2.,<TT2,TT0>	;NO EXCEPTIONS
	MOVE	T3,T0		;SAVE IMAG. PART IN T3

	CAML	T4,T5		;COMPARE |X| AND |Y|
	  JRST NOXCH		;  AND GET LARGER = A IN T4
	EXCH	T4,T5		;  AND SMALLER = B IN T5

NOXCH:	MOVE	T2,T4		;LARGER TO T2
	LSH	T2,-33		;ISOLATE ITS BIASED EXPONENT
	SUBI	T2,200		;  AND UNBIAS IT
	JUMPLE	T2,LE		;IF UNBIASED EXPONENT > 0,
	  SUBI	T2,1		;  DECREMENT IT
LE:	MOVN	T2,T2		;NEGATE SCALE INDEX
	FSC	T4,(T2)		;SCALE A TO [1/2,2) AND B BY
	FSC	T5,(T2)		;  SAME FACTOR. NO EXCEPTIONS
	LSH	T2,1		;MULTIPLY NEGATED SCALE INDEX BY 2
	DMOVE	T0,T4		;COPY |X| AND |Y| TO T0,T1
	FMPR	T0,T0		;SQUARE LARGER.  NO EXCEPTIONS
	FMPR	T1,T1		;SQUARE SMALLER.  NO OVERFLOW AND
	  JFCL			;  UNDERFLOW WON'T MATTER
	MOVEM	T1,SMALSQ	;SAVE SMALLER SQUARED
	FADR	T1,T0		;GET SCALED SQUARE OF MODULUS
	CAMLE	T1,RT2		;IF T1 .GT. SQRT(2)
	  JRST	MORE		;  GO TO MORE SCALING
	CAMGE	T1,RT2O2	;IF T1 .LT. SQRT(1/2)
	  JRST	MORE		;  GO TO MORE SCALING
	JUMPN	T2,JOIN		;IF SCALE INDEX NOT 0, GO TO JOIN

;At this point the scale index is zero, implying 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 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:

	MOVE	T0,T4		;COPY LARGER TO T0
	FSBR	T0,ONE		;GET A - 1
	FADR	T4,ONE		;AND A + 1
	FMPR	T4,T0		;AND THEIR PRODUCT
	FADR	T4,SMALSQ	;T4 NOW HAS A**2 + B**2 - 1
	MOVE	T1,T4		;  = NUMERATOR OF Z. COPY TO T1
	FADR	T1,TWO		;GET DENOMINATOR OF Z IN T1
	JRST	MERGE		;MERGE FOR POLYNOMIAL 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	T4,T1		;MAKE COPY OF A**2 + B**2
	LSH	T4,-33		;ISOLATE ITS BIASED EXPONENT
	SUBI	T4,200		;UNBIAS THE EXPONENT
	MOVN	T4,T4		;  AND NEGATE IT
	FSC	T1,(T4)		;SCALED MODULUS NOW IN [1/2,1)
	ADD	T2,T4		;TOTAL NEGATED SCALE INDEX TO T2
	CAML	T1,RT2O2	;IF T1 .GT. SQRT(1/2)
	  JRST	JOIN		;  SCALING IS DONE
	FSC	T1,1		;SCALE T1 TO [1,SQRT(2)]
	ADDI	T2,1		;  AND INCREMENT NEGATED INDEX

;At this point the scaled (A**2 + B**2) is in the interval
;[SQRT(1/2),SQRT(2)]. T2 contains the negative of the index
;necessary to compensate later for the scaling.  In the following
;lines of code, a polynomial approximation is used to evaluate
;the natural log of the scaled (A**2 + B**2).

JOIN:	MOVN	T2,T2		;RESTORE SIGN OF SCALE INDEX
	MOVE	T4,T1		;GET COPY OF SCALED (A**2 + B**2)
	FSBR	T4,ONE		;  -1 = NUMERATOR OF Z
	FADR	T1,ONE		;GET DENOMINATOR OF Z IN T1

;The following code is taken from the AlOG routine.  The constants
;L3, L4, and L5 are those used in ALOG.

MERGE:	FDVR	T4,T1		;Z = ZNUM / ZDEN. NO EXCEPTIONS
	MOVE	T5,T4		;SAVE A COPY OF Z
	FMPR	T4,T4		;W = Z**2
	MOVE	T0,T4		;  AND MAKE COPY
	FMPR	T0,L3		;FORM P(Z) = W*L3
	FADR	T0,L4		;  + L4
	FMPR	T0,T4		;  *W
	FADR	T0,L5		;  + L5
	FMPR	T0,T4		;  *W
	FMPR	T0,T5		;  *Z. NO EXCEPTIONS
	FSC	T5,1		;GET 2*Z

;2*Z still needs to be added into P(Z).
	JUMPN	T2,REST		;IF SCALE INDEX = 0,
	FADR	T0,T5		;  COMPLETE P(Z) BY ADDING IN Z
	JRST	FINISH		;GO FINISH UP

;The following lines of code complete the calculation of the real
;part for the index not zero.  The real part is the double
;precision sum of INDEX*LN(2), Z, and P(Z), rounded to single
;precision.

REST:	MOVE	T4,T5		;GET 2*Z IN DOUBLE
	SETZB	T1,T5		;  PRECISION
	DFAD	T0,T4		;GET DP SUM. FLOAT SCALE
	FLTR	T4,T2		;INDEX -- DP SINCE T5 HAS 0
	DFMP	T4,LN2		;GET INDEX*LN(2)
	DFAD	T0,T4		;AND ADD TO PREVIOUS SUM
	PUSHJ	P,SNG.0		;ROUND TO SINGLE PRECISION	

FINISH:	FSC	T0,-1		;DIVIDE BY 2 TO GET LN(MODULUS)
	MOVE	T1,T3		;IMAG PART TO T1

	POP	P,T5		;POP REGISTERS
	POP	P,T4		;  IN REVERSE
	POP	P,T3		;  ORDER IN WHICH
	POP	P,T2		;  THEY WERE PUSHED
	GOODBY	(1)

BIG:	DMOVE	T2,T0		;SAVE X AND Y
	MOVMM	T3,TT0		;STORE |Y| FOR ALOG
	FUNCT	ALOG.,<TT0>
	MOVE	T1,PIO2		;IMAGINARY PART NEAR +/- PI/2
	JUMPGE	T3,PO2POS	;  + FOR Y POSITIVE
	  MOVN	T1,T1		;  - FOR Y NEGATIVE
PO2POS:	FDVR	T2,T3		;GET X/Y; NO OVERFLOW
	  JFCL	UND1		;  UNDERFLOW IS POSSIBLE
	FSBR	T1,T2		;SMALL CORRECTION TO IMAG
	FMPR	T2,T2		;GET (X/Y)**2. NO OVERFLOW
	  JFCL	UND1		;  UNDERFLOW IS POSSIBLE
	FSC	T2,-1		;(1/2)*(X/Y)**2. NO OVERFLOW
	  JFCL	UND1		;  UNDERFLOW IS POSSIBLE
	FADR	T0,T2		;+ LN(|Y|) = REAL. NO EXCEPTIONS
	JRST	REST23		;GO RESTORE T2 AND T3

SMALL:	DMOVE	T2,T0		;SAVE X AND Y
	MOVMM	T2,TT0		;STORE |X| FOR ALOG
	FUNCT	ALOG.,<TT0>
	SETZ	T1,		;IMAG. PART NEAR 0 OR +/- PI
	JUMPGE	T2,IMAGOK	;  0 FOR X POSITIVE
	  MOVE	T1,CPI		;  +/- PI FOR X NEGATIVE
	JUMPGE	T3,IMAGOK	;    + PI FOR Y POSITIVE
	  MOVN	T1,T1		;    - PI FOR Y NEGATIVE
IMAGOK:	FDVR	T3,T2		;GET Y/X; NO OVERFLOW
	  JFCL	UND2		;  UNDERFLOW IS POSSIBLE
	FADR	T1,T3		;SMALL CORRECTION FOR IMAG.
	FMPR	T3,T3		;GET (X/Y)**2. NO OVERFLOW
	  JFCL	UND1		;  UNDERFLOW IS POSSIBLE
	FSC	T3,-1		;(1/2)*(X/Y)**2. NO OVERFLOW
	  JRST	UND1		;  UNDERFLOW IS POSSIBLE
	FADR	T0,T3		;+ LN(|X|) = REAL. NO EXCEPTIONS
	JRST	REST23		;GO RESTORE REGISTERS T2 AND T3

UND2:	JUMPN	T1,REST23	;DONE IF |IMAG| NEAR PI OR PI/2
	$LCALL	IPU
;LERR	(LIB,%,<CLOG: Imaginary part underflow>)

UND1:	JUMPN	T0,REST23	;DONE IF LN(|X| OR |Y|) NOT 0
	$LCALL	RPU
;LERR	(LIB,%,<CLOG: Real part underflow>)

REST23:	POP	P,T3		;RESTORE REGISTERS
	POP	P,T2		;  T2 AND T3

	GOODBY	(1)

YZRO:	JUMPE	T0,XYZRO	;X = Y = 0 IS A SPECIAL CASE
	JUMPGE	T0,RPOS		;Y IS 0, X ISN'T. IF X < 0
	MOVN	T0,T0		;X=ABS(X)
	MOVE	T1,CPI		;V=PI

RPOS:	MOVEM	T1,TT2		;SAVE IMAGINARY PART
	MOVEM	T0,TT0		;MOVE X TO MEMORY
	FUNCT	ALOG.,<TT0>	;COMPUTE U=DLOG(X)
	MOVE	T1,TT2		;RESTORE IMAGINARY PART
	GOODBY	(1)		;RETURN; RESULT IN T0,T1

XZRO:	MOVE	T0,T1		;X IS 0, Y ISN'T. COPY Y TO T0
	MOVE	T1,PIO2		;IMAG PART IS +/- PI/2
	JUMPGE	T0,YPLUS	;IF Y < 0,
	MOVN	T0,T0		; MAKE IT POSITIVE, AND
	MOVN	T1,T1		; MAKE IMAG PART NEGATIVE
YPLUS:	MOVEM	T0,TT0		;PREPARE TO GET LOG
	MOVEM	T1,TT2		;SAVE IMAGINARY PART
	FUNCT	ALOG.,<TT0>	;GET LOG(|Y|), NO EXCEPTIONS
	MOVE	T1,TT2		;RESTORE IMAGINARY PART.
	GOODBY	(1)		;RETURN: RESULTS IN T0, T1

XYZRO:	$LCALL	ZIZ
;LERR	(LIB,%,<entry CLOG; arg is (0.0,0.0); result=(-infinity,0.0)>)
	MOVE	T0,[400000000001] ;REAL ANSWER IS NEGATIVE INFINITY
	GOODBY	(1)		;RETURN; RESULT IN T0,T1

;CONSTANTS

MASK:	EXP 377000000000		;MASK FOR EXPONENT FIELD
P14:	EXP 016000000000		;HIGH 9 BITS = 16 OCTAL
M14:	EXP 762000000000		;HIGH 9 BITS = -16 OCTAL
ONE:	EXP 201400000000		;1.0
TWO:	EXP 202400000000		;2.0
CPI:	EXP 202622077325		;PI
PIO2:	EXP 201622077325		;PI / 2
RT2O2:	EXP 200552023632		;SQRT(2) / 2
RT2:	EXP 201552023632		;SQRT(2)
L3:	EXP 177464164321		;.301003281
L4:	EXP 177631177674		;.39965794919
L5:	EXP 200525253320		;.666669484507

LN2:	EXP 200542710277,276434757153	;LN(2)

;DATA
	SEGMENT	DATA
TT0:	BLOCK	1			;FOR ATAN2 AND ALOG CALLS
TT2:	BLOCK	1			;FOR ATAN2 CALL
SMALSQ:	BLOCK	1

	PRGEND
TITLE	CSIN	COMPLEX SINE FUNCTION  
;		(SINGLE 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	CSIN
EXTERN	CSIN.
CSIN=CSIN.
PRGEND
TITLE	CCOS	COMPLEX COSINE FUNCTION  
;		(SINGLE 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	CCOS
EXTERN	CCOS.
CCOS=CCOS.
PRGEND
TITLE	CSIN.	COMPLEX SINE AND COSINE ROUTINES
;		(SINGLE 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.

;CSIN(Z) AND CCOS(Z), WHERE Z = X + I*Y, ARE CALCULATED AS FOLLOWS

;  CSIN(Z) = SIN(X)*COSH(Y) + I*COS(X)*SINH(Y)
;  CCOS(Z) = COS(X)*COSH(Y) - I*SIN(X)*SINH(Y)

;  THE FUNCTIONS SINH AND COSH 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 CSIN AND CCOS IS AS FOLLOWS
;FOR
;      Z = X+I*Y.  IF ABS(X) > 823549.66 THE RESULT IS SET TO 0.0 AND
;      AN ERROR MESSAGE IS RETURNED.

;FOR ABS(Y) > 88.029692 CALCULATIONS PROCEED AS FOLLOWS:

;      FOR THE REAL PART OF THE RESULT
;          LET T = ABS(SIN(X)), THEN

;          FOR T = 0.0 THE REAL PART OF THE RESULT IS SET TO 0.0
;          FOR LOG(T) + ABS(Y) > 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 = ABS(COS(X)), THEN

;          FOR T = 0.0 THE IMAGINARY PART OF THE RESULT IS SET TO 0.0
;          FOR LOG(T) + ABS(Y) > 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:  SIN,COS,EXP,ALOG

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

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  MOVEI	L,ARG
;  PUSHJ	P,CSIN
;	OR
;  PUSHJ	P,CCOS

;THE COMPLEX ANSWER IS RETURNED IN ACCUMULATORS T0 AND T1.

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(CCOS,.)		;ENTRY TO CCOS ROUTINE
	PUSH	P,T3
	MOVEI	T3,1			;SET FLAG TO 1 FOR CCOS
	JRST	CSIN1			;GO TO CSIN ROUTINE

	HELLO	(CSIN,.)		;ENTRY TO CMPLX SINE ROUTINE
	PUSH	P,T3
	SETZI	T3,			;SET FLAG TO 0 FOR CSIN

CSIN1:	PUSH	P,T2
	PUSH	P,T4
	PUSH	P,T5
	XMOVEI	T1,@(L)			;GET ADDRESS OF ARG
	MOVE	T0,(T1)			;X = REAL(Z)
	MOVM	T2,T0                   ; GET ABS(X)
	CAMG	T2,CUT                  ;IF ABS(X) IS NOT TOO LARGE
	  JRST	XOK			; GO TO GET Y
	$LCALL	ARZ
;LERR	(LIB,%,<CSIN or CCOS: ABS(REAL(arg)) too large; result = (0.0,0.0)>)
	SETZI	T0,			;SET REAL(RESULT) TO 0.0
	SETZI	T1,			;SET IMAG(RESULT) TO 0.0
	JRST 	RET			;RETURN

XOK:	MOVE	T1,1(T1)		;Y = IMAG(Z)
	MOVEM	T1,YSAVE		;SAVE Y
	MOVEM	T0,ARGSAV
	FUNCT	SIN.,<ARGSAV>		;CALCULATE SIN(X)
	MOVEM	T0,SX			;SX = SIN(X)
	FUNCT	COS.,<ARGSAV>		;CALCULATE COS(X)
	MOVEM	T0,CX			;CX = COS(X)
	SKIPN	T3			;IF THIS IS THE CSIN ROUTINE
	  JRST	NOSWT			;THEN GO TO NOSWT
	MOVN	T2,SX			;OTHERWISE, PUT -SIN(X)
	EXCH	T2,CX			;IN CX, AND COS(X)
	MOVEM	T2,SX			;IN SX.

NOSWT:	SKIPN	T1,YSAVE		;SKIP UNLESS Z IS REAL, IN
	  JRST	CSIN2			;WHICH CASE, GO TO THE SIMPLE CASE
	MOVM	T2,T1			;AY = ABS(YSAVE)
	CAMLE 	T2,XMAX			;IF AY IS TOO LARGE FOR SIMPLE
	  JRST	OVFL			;SINH OR COSH, GO TO OVFL.
	CAMG	T2,TWOM14		;IF |Y| <= 2**(-14),
	  JRST	EASY			;  SINH = Y, COSH = 1
	MOVEM	T2,ARGSAV
	FUNCT	EXP.,<ARGSAV>		;CALCULATE EXP(AY)
	MOVE	T3,T0
	MOVE    T5,ONE
	FDVR	T5,T3
	CAMLE	T2,ONE			;IF AY IS GREATER THAN 1
	  JRST	GTET			;THEN GO TO GTET
	MOVE	T1,T2			;GET COPY OF |Y|
	FMPR	T1,T1			;SS = AY**2
	MOVE	T4,T1
	FMPR	T4,R4			;SS = SS*R4
	FADR	T4,R3			;SS = SS+R3
	FMPR	T4,T1			;SS = SS*(AY**2)
	FADR	T4,R2			;SS = SS+R2
	FMPR	T4,T1			;SS=SS*(AY**2)
	FADR	T4,R1			;SS = SS+R1
	FMPR	T1,T4			;SS = SS*(AY**2)
	FMPR	T1,T2			;SS = SS*AY
	FADR	T1,T2			;SS = SS+AY
	JRST	CCCC			;GO TO CCCC

GTET:	MOVN	T1,T5
	FADR	T1,T3			;SS = T+SS
	FSC	T1,-1			;SS = SS/2.
	
CCCC:	MOVE	T0,T5
     	FADR	T0,T1			;CC = SS+(1/T)
	MOVE	T3,YSAVE		;RESTORE Y
	CAIGE	T3,0			;IF Y IS LESS THAN 0.0
	  MOVN	T1,T1			;THEN NEGATE SS
     	FMPR	T0,SX			;GET REAL PART OF RESULT
	FMPR	T1,CX			;GET IMAG PART OF RESULT
	  JFCL	IMUND
	JRST	RET			;RETURN

EASY:	MOVE	T0,SX			;REAL PART = 1 * SX
	FMPR	T1,CX			;IMAG PART = Y * CX
	  JFCL	IMUND			;UNDERFLOW POSSIBLE
	JRST	RET			;RETURN

IMUND:	$LCALL	IPO,RET
;LERR	(LIB,%,<CSIN or CCOS: imaginary part overflow>,RET)

OVFL:	MOVM	T0,SX			;T=ABS(SX)
	CAIN	T0,			;IF T IS EQUAL TO 0.
	  JRST	GOTR			;THEN GO TO GOTR
	MOVEM	T0,ARGSAV		;OTHERWISE,
	FUNCT	ALOG.,<ARGSAV>		;CALCULATE LOG(T)
	FADR	T0,T2			;T = LOG(T)+AY
	FADR	T0,LN2			;T = T-LOG(2.0)
	CAMG	T0,XMAX			;IF T IS .LE. XMAX
	  JRST	CONTIN			;THEN GO TO CONTIN
	$LCALL	RPO
;LERR	(LIB,%,<CSIN or CCOS: real part overflow>)
	HRLOI	T0,377777	        ;REAL PART OF RESULT SET TO INFINITY
	JRST	SKP2			;SKIP NEXT 2 INSTRUCTIONS

CONTIN:	MOVEM	T0,ARGSAV
	FUNCT	EXP.,<ARGSAV>		;RRES = EXP(T)

SKP2:	SKIPGE	SX			;IF SX IS LESS THAN 0.
	  MOVN	T0,T0			;THEN NEGATE RRES

GOTR:	SKIPN	T1,CX		        ;IF CX = 0,
	  JRST  RET			;THEN RETURN
	MOVEM	T0,SX			;OTHERWISE SAVE RRES
	MOVMM	T1,ARGSAV		;T = ABS(CX)
	FUNCT	ALOG.,<ARGSAV>		;CALCULATE T=LOG(T)
	MOVE	T1,T0
	FADR	T1,T2			;T = T+AY
	FADR	T1,LN2			;T = T-LOG(2)
	CAMG	T1,XMAX			;IF T IS .LE. XMAX
	  JRST	CONTN2			; THEN GO TO CONTN2
	$LCALL	IPO
;LERR	(LIB,%,<CSIN or CCOS: imaginary part overflow>)
	HRLOI	T1,377777		;SET IRES TO INFINITY
	JRST	SKP3			;SKIP NEXT 3 INSTRUCTIONS
	
CONTN2:	MOVEM	T1,ARGSAV
	FUNCT	EXP.,<ARGSAV>		;IRES = EXP(T)
	MOVE	T1,T0

SKP3:	MOVE	T0,SX
	MOVE	T2,YSAVE
	XOR	T2,CX			;T2 HAS THE SIGN OF CX*Y
	JUMPGE	T2,RET			;IF SIGNS ARE THEN SAME, DONE
	  MOVN	T1,T1			;OTHERWISE NEGATE IMAG PART OF RESULT
	JRST    RET			;RETURN

CSIN2:	MOVE	T0,SX			;SIMPLE CASE, Z IS REAL

RET:	POP	P,T5
	POP	P,T4
	POP	P,T2
	POP	P,T3
	GOODBY	(1)

TWOM14:	163400000000			;2**(-14)
LN2:	577235067500			;-(NATURAL LOG OF 2.0)
XMAX:	207540074636			;88.029692
R1:	176525252524			;1.666666643E-1
R2:	172421042352			;8.333352593E-3
R3:	164637771324			;1.983581245E-4
R4:	156572227373			;2.818523951E-6
ONE:	201400000000			;1.0
CUT:    234622077325			;2**26 * PI

	SEGMENT	DATA
CX:	0
SX:	0
YSAVE:	0
ARGSAV: 0

	PRGEND
TITLE	CSQRT	COMPLEX SQUARE ROOT FUNCTION  
;		(SINGLE 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	CSQRT
EXTERN	CSQRT.
CSQRT=CSQRT.
PRGEND
TITLE	CSQRT.	COMPLEX SQUARE ROOT FUNCTION
;		(SINGLE 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.

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

;  IF X AND Y ARE BOTH >= 0.0, THEN
;    U = SQRT((ABS(X)+CABS(Z))/2.0)
;    V = Y/(2.0*U)

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

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

;  IF X AND Y ARE BOTH < 0.0, THEN
;    U = Y/(2.0*V)
;    V = -SQRT((ABS(X)+CABS(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:  SQRT

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

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  MOVEI	L,ARG
;  PUSHJ	P,CSQRT

;THE REAL PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0
;THE IMAG PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(CSQRT,.)	;ENTRY TO COMPLEX SQUARE ROOT ROUTINE.
	XMOVEI	T1,@0(L)	;PICK UP ADDRESS OF ARGUMENT
	MOVE	T0,(T1)		;PICK UP X IN AC T0.
	MOVE	T1,1(T1)	;PICK UP Y IN AC T1.
	PUSH	P,T2		;SAVE AC T2.
	JUMPE	T1,ZREAL	;JUMP IF Y=0
	JUMPE	T0,ZIMAG	;JUMP IF X=0
	MOVM	T2,T0		;ABS(X) TO AC T2.

	PUSH	P,T3		;SAVE REGISTER T3.
	PUSH	P,T4		;SAVE REGISTER T4
	PUSH	P,T5		;SAVE REGISTER T5
	PUSH	P,G1		;SAVE REGISTER G1
	MOVE	T4,T0		;SAVE X IN AC T4.
	MOVE	T5,T1		;SAVE Y IN AC T5.
	MOVM	T3,T1		;ABS(Y) TO AC T3.
	SETZ	G1,		;G1 CONTAINS 0 IF ABS(X)
	CAMLE	T2,T3		;<= ABS(Y), ELSE IT CONTAINS 1.  PUT
	AOJA	G1,DWN2		;THE LARGER OF ABS(X),ABS(Y) IN AC T2
	EXCH	T2,T3		;AND THE SMALLER IN AC T3.
DWN2:	FDVR	T3,T2		;CALC S/L.
	JFCL			;NO UNDERFLOW ERROR MESSAGE ALLOWED.
	FMPR	T3,T3		;CALC (S/L)**2.
	JFCL			;NO UNDERFLOW ERROR MESSAGE ALLOWED.
	FADRI	T3,201400	;HAVE (1+(S/L)**2) IN AC T3.
	MOVEM	T3,TEMP
	FUNCT	SQRT.,<TEMP>	;[3220] CALC. THE SQRT OF
				;(1+(S/L)**2)=1+EPS.
	JUMPE	G1,YGTX		;GO TO YGTX IF ABS(Y) > ABS(X).

XGTY:	FADRI	T0,201400	;CALC. 2 + EPS.
	FSC	T0,-1		;CALC. 1+EPS/2.
	MOVMM	T0,TEMP		;PUT IN TEMP FOR
	FUNCT	SQRT.,<TEMP>	;[3220] CALL TO SQRT
	MOVE	T3,T0		;SAVE SQRT(1+EPS/2) IN AC T3.
	MOVEM	T2,TEMP
	FUNCT	SQRT.,<TEMP>	;[3220] CALC.
	FMPR	T0,T3		;CALC. SQRT(ABS(X)*(1+EPS/2)).
	JRST	OUT1		;GO TO REST OF CALC.

YGTX:	CAMGE	T2,[1.0]	;IF ABS(Y)>1, GO TO POSSIBLE OVERFLOW CASE.
	JRST	YXSMAL		;ELSE, GO TO POSSIBLE UNDERFLOW.
	FSC	T0,-3		;CALC. (1+EPS)/8.
	FMPR	T2,T0		;CALC. ABS(Y)*(1+EPS)/8.
	MOVM	T3,T4		;ABS(X) TO AC T3.
	FSC	T3,-3		;CALC. ABS(X)/8.
	JFCL			;SUPPRESS UNDERFLOW ERROR MESSAGE.
	FADR	T3,T2		;CALC.(ABS(X)/8)+(ABS(Y)*(1+EPS)/8).
	MOVEM	T3,TEMP
	FUNCT	SQRT.,<TEMP>	;[3220] CALC.
	FSC	T0,1	

OUT1:	MOVM	T1,T5	 	;ABS(Y) TO AC T1.
	FDVR	T1,T0		;DIVIDE ABS(Y)/2 BY THE
	  JFCL	UNDFL
	FSC	T1,-1		;SQRT TERM.
	  JFCL	UNDFL
	JRST	SIGNS		;GO TO REST OF CALC.

YXSMAL:	FMPR	T0,T2		;ABS(Y)*(1+EPS) = CDABS(Z).
	MOVM	T3,T4		;ABS(X) TO AC T3.
	FADR	T3,T0		;GET |X| + CDABS(Z)
	FSC	T3,1		;GET 2*(|X| + CABS(Z))
	MOVEM	T3,TEMP
	FUNCT	SQRT.,<TEMP>	;[3220] CALC. SQRT OF 2*(|X| + CABS(Z))
	MOVE	T3,T0		;GET SQRT((|X| + CABS(Z))*2)
	FSC	T0,-1		;GET SQRT((|X| + CABS(Z))/2)
	FDVR	T2,T3		;GET |Y| / SQRT(2*(|X| + CABS(Z))).
	MOVE	T1,T2		;PUT A TEMP ANSWER IN AC T1.

SIGNS:  CAIGE   T4,0		;SET UP THE REAL AND
  	  EXCH	T1,T0		;THE IMAGINARY ANSWERS WITH
	CAIGE	T5,0		;THE APPROPRIATE
	  MOVN	T1,T1		;SIGNS.
	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	T2,T0		;X NOT ZERO; SAVE IT.
	MOVMM	T2,TEMP		;GET ABS(X) FOR SQRT.
	FUNCT	SQRT.,<TEMP>	;[3220] T0 GETS SQRT(ABS(X)).
	SETZ	T1,		;T1 smashed by call; Set to 0 again
	JUMPG	T2,DONE		;T0,T1 OK FOR X>0.
	  EXCH	T0,T1		;  INTERCHANGE FOR X<0.
	POP	P,T2		;RESTORE T2.
	GOODBY	(1)		;RETURN.

ZIMAG:	MOVE	T2,T1		;X = 0, SAVE Y.
	FSC	T1,-1		;DIVIDE Y BY 2.
	MOVMM	T1,TEMP		;GET ABS(Y/2) FOR SQRT.
	FUNCT	SQRT.,<TEMP>	;[3220] T0 GETS SQRT(ABS(Y/2)).
	MOVE	T1,T0		;SO DOES T1.
	JUMPG	T2,DONE		;  DONE IF Y>0.
	  MOVN	T1,T1		;NEGATE IMAG PART IF Y<0.
DONE:	POP	P,T2		;RESTORE T2.
	GOODBY	(1)		;RETURN

UNDFL:	CAIGE	T4,0		;SKIP IF REAL & IMAG SWITCHED
	 $LCALL	RPU,SIGNS
;	 LERR	(LIB,%,<CEXP: real part underflow>,SIGNS)
	$LCALL	IPU,SIGNS
;	LERR	(LIB,%,<CEXP: imaginary part underflow>,SIGNS)

SQ2:	201552023632		;SQRT(2).

	SEGMENT	DATA
TEMP:	0			;TEMPORARY STORAGE FOR SQRT ARGS

	PRGEND
TITLE	CFM	COMPLEX MULTIPLICATION
SUBTTL	KAREN KOLLING /KK/CKS		28-Oct-81

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 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.

;FROM LIB40 VERSION .027
;FROM	V.021	16-OCT-70

;COMPLEX FLOATING MULTIPLY ROUTINE.

;THIS ROUTINE MULTIPLIES TWO COMPLEX NUMBERS ACCORDING
;TO THE RULE:

; (A+IB)*(C+ID)=(A*C-B*D)+I*(B*C+A*D)

;IF EITHER PRODUCT IN EITHER TERM OVER OR UNDERFLOWS,
;A JUMP IS MADE TO A SPECIAL CASE CALCULATION. OVER
;AND UNDERFLOW ERROR MESSAGES ARE RETURNED ONLY IF THE
;ENTIRE TERM OVER OR UNDERFLOWS.

;THERE ARE TWO SETS OF ENTRY POINTS TO THIS ROUTINE:
;THE MULTIPLY, AND THE MULTIPLY-TO-MEMORY.
;THE CALLING SEQUENCE FOR THE MULTIPLY ROUTINES IS AS
;FOLLOWS:
;	XMOVEI	L,ARG2
;	PUSHJ	P,CFM.N
;WHERE N=0,2,4,6,10,12,14. ARG1 IS ASSUMED TO BE IN
;AC'S N AND N+1, AND THE RESULTS ARE LEFT IN AC'S N AND N+1.
;THE CALLING SEQUENCE FOR THE MULTIPLY TO MEMORY ROUTINES IS AS
;FOLLOWS:
;	XMOVEI	L,ARG2
;	PUSHJ	P,CMM.N
;WHERE N AND ARG1 ARE AS BEFORE AND THE RESULTS ARE
;LEFT IN ARG2 AND ARG2+1.

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

DEFINE SYM (N)<
	ENTRY	CFM.'N,CMM.'N

	SIXBIT	/ CFM.'N/
CFM.'N:	DMOVEM	N,ARG11		;STORE FIRST ARG
IFN N,<	DMOVEM	0,SAVE0 >	;SAVE REGS 0-1
	DMOVE	0,(L)		;GET SECOND ARG
	PUSHJ	P,CFM		;MULTIPLY
IFN N,<	DMOVE	N,0		;STORE RESULT
	DMOVE	0,SAVE0 >	;RESTORE REGS 0-1
	POPJ	P,		;RETURN

	SIXBIT	/ CMM.'N/
CMM.'N:	DMOVEM	N,ARG11		;STORE FIRST ARG
	DMOVEM	0,SAVE0		;SAVE REGS 0-1
	DMOVE	0,(L)		;GET SECOND ARG
	PUSHJ	P,CFM		;MULTIPLY
	DMOVEM	0,(L)		;STORE RESULT
	DMOVE	0,SAVE0		;RESTORE REGS 0-1
	POPJ	P,		;RETURN
>

XLIST	;GENERATE CFM.0 ... CFM.14 AND CMM.0 ... CMM.14
.N.=16
REPEAT 7,<SYM(\<.N.==.N.-2>)>
LIST
CFM:	DMOVEM	T0,ARG21	;STORE SECOND ARG
	DMOVEM	T2,SAVE2	;SAVE SCRATCH REGISTERS
	PUSHJ	P,CALC		;CALC (A*C-B*D)
	MOVEM	T0,REAL		;SAVE REAL(ANS) IN REAL.
	MOVE	T0,ARG12	;SET UP
	EXCH	T0,ARG11	;AND CALCULATE
	MOVNM	T0,ARG12	;(B*C+A*D)
	PUSHJ	P,CALC		;AND LEAVE IT IN T0.
	MOVE	T1,T0		;STORE IMAG(ANS) IN T1
	MOVE	T0,REAL		;GET REAL(ANS) IN T0
	DMOVE	T2,SAVE2	;RESTORE SCRATCH REGISTERS
	POPJ	P,		;RETURN

CALC:	MOVE	T0,ARG11	;"A" TO AC0.
	MOVE	T2,ARG12	;"B" TO AC2.
	FMPR	T0,ARG21	;CALC A*C AND IF
	JFCL	OUFLO1		;IT OVER/UNDERFLOWS, GO TO OUFLO1.

SECOND:	FMPR	T2,ARG22	;CALC B*D AND IF
	JFCL	OUFLO2		;IT OVER/UNDERFLOWS, GO TO OUFLO2.
	FSBR	T0,T2		;CALC A*C-B*D AND
	POPJ	P,		;RETURN.

OUFLO1:	MOVE	T3,ARG22	;"D" TO AC3.
	JUMPE	T0,UNDER1	;UNDERFLOW, GO TO UNDER1
	FMPRB	T2,T3		;CALC B*D AND
	JFCL	OFL1OU		;IF OVER/UNDERFLOW GO TO OFL1OU.
	XOR	T3,T0		;OVERFLOW + NORMAL CASE
	JUMPGE	T3,SROVNM	;IF SIGNS ARE THE SAME, GO TO SROVNM
INFIN:	MOVEM	T0,ANSTMP	;STORE ANSWER.
INFIN1:	$LCALL	999
;LERR	(999,%,Complex overflow at $1L,<T1>)
	MOVE	T0,ANSTMP	;PICK UP ANSWER.
	POPJ	P,		;RETURN.
OFL1OU:	JUMPE	T2,INFIN	;OVERFLOW + UNDERFLOW CASE GOES TO INFIN.
	XOR	T3,T0		;OVERFLOW + OVERFLOW CASE.
	JUMPGE	T3,SROVOV	;IF SIGNS ARE THE SAME, GO TO SROVOV
	JRST	INFIN		;O'E, GO TO INFIN.

SROVNM:	JUMPE	T2,INFIN	;OVERFLOW + ZERO, GO TO INFIN.
	MOVE	T0,ARG21	;C TO AC0.
	FSC	T0,-1		;CALC. C/2.
	FMPR	T0,ARG11	;CALC. (C/2)*A AND IF
	JFCL	SROV1		;IT OVERFLOWS, IMMEDIATELY RETURN.
	FSBRM	T0,T2		;CALC ((C/2*A)-(B*D).
	FADR	T0,T2		;CALC (A*C-B*D) AND
	POPJ	P,		;RETURN.
SROV1:	MOVEM	T0,ANSTMP	;STORE ANSWER
	JRST	INFIN1		;GO TO ERROR

SROVOV:	MOVE	T0,ARG21	;C TO AC0.
	FDVR	T0,ARG22	;CALC. (C/D).
	MOVE	T1,ARG12	;B TO AC1.
	FDVR	T1,ARG11	;CALC. (B/A).
	FSBR	T0,T1		;CALC. ((C/D)-(B/A)) AND
	JFCL	FIXUP		;IF IT UNDERFLOWS, GO TO FIXUP.
	FMPR	T0,ARG22	;CALC. ((C)-(B/A)*D)) AND
	JFCL			;SUPPRESS OVERFLOW MESSAGE.
	FMPR	T0,ARG11	;CALC (A*C-B*D) AND
	POPJ	P,		;RETURN.
FIXUP:	FMPR	T1,ARG22	;CALC. ((B/A)*D).
	MOVE	T0,ARG21	;C TO AC0.
	FSBR	T0,T1		;CALC. (C-(B/A)*D).
	FMPR	T0,ARG11	;CALC. (A*C-B*D) AND
	POPJ	P,		;RETURN.

UNDER1:	FMPR	T2,T3		;CALC. B*D AND GO
	JFCL	U1OU		;TO U1OU IF OVER/UNDERFLOW.
	MOVM	T3,T2		;IF B*D IS <2**-105,
	CAMGE	T3,[030400000000] ;GO
	JRST	UNDR0		;TO .+3.
	MOVN	T0,T2		;NO BITS CAN BE SAVED,
	POPJ	P,		;FIX SIGN AND RETURN.
UNDR0:	FSC	T2,32		;CALC B*D*(2**27).
UNDR1:	MOVE	T0,ARG11	;CALC
	FSC	T0,32		;A*C*(2**27)
	FMPR	T0,ARG21	;AND SUPPRESS AN
	JFCL			;ERROR MESSAGE.
UNDR4:	FSBR	T0,T2		;CALC (2**27)(A*C-B*D),
UNDR5:	FSC	T0,-32		;CALC (A*C-B*D) AND
	POPJ	P,		;RETURN.
U1OU:	JUMPE	T2,U1OU1	;IF OVERFLOW + UNDERFLOW, GO
	MOVNM	T2,ANSTMP	;TO ERROR
	JRST	INFIN1		;RETURN.
U1OU1:	MOVE	T2,ARG12	;O'E, TRY TO SAVE BITS.
	FSC	T2,32		;CALC. B*(2**27).
	FMPR	T2,ARG22	;CALC B*D*(2**27)
	JFCL			;AND SUPPRESS UNDERFLOW MESSAGE
	MOVE	T0,ARG11	;CALC. A*(2**27) AND
	FSC	T0,32		;A*(2**27)*C
	FMPR	T0,ARG21	;AND SUPPRESS
	JFCL			;UNDERFLOW MESSAGE.
	JUMPN	T2,UNDR4	;IF BOTH UNDERFLOW, RETURN
	JUMPN	T0,UNDR5	;AN ERROR MESSAGE. O'E,
	$LCALL	888
;LERR	(888,%,Complex underflow at $1L,<T1>)
	SETZ	T0,		;CLEAR AC0
	POPJ	P,		;CALC.

OUFLO2:	JUMPN	T2,OUFLO3	;JUMP AHEAD IF (B*D) OVERFLOWED.
	CAML	T0,[030400000000] ;IF NO BITS CAN BE SAVED,
	POPJ	P,		;RETURN.
	FSC	T0,32		;CALC. (2**27)*(A*C)
	MOVE	T2,ARG12	;AND
	FSC	T2,32		;THEN
	FMPR	T2,ARG22	;(2**27)*(B*D)
	JFCL			;AND GO TO THE REST
	JRST	UNDR4		;OF THE CALC.
OUFLO3:	MOVEM	T2,T1		;OVERFLOW + NORMAL.
	XOR	T1,T0		;IF THE SIGNS ARE THE SAME,
	JUMPGE	T1,SROVN	;GO TO SROVN
	MOVNM	T2,ANSTMP	;THE ANSWER AND
	JRST	INFIN1		;GO TO ERROR EXIT.
SROVN:	MOVE	T3,ARG22	;"D" TO AC3.
	FSC	T3,-2		;CALC (D/2).
	FMPR	T3,ARG12	;CALC (B*(D/2)) AND IF
	JFCL	FIXUP2		;IT OVERFLOWS, GO TO FIXUP2.
	FSBR	T0,T3		;CALC. ((A*C)-(B*(D/2))).
	FSBR	T0,T3		;CALC (A*C-B*D) AND
	POPJ	P,		;RETURN.
FIXUP2:	MOVNM	T3,ANSTMP	;SET UP THE ANSWER
	JRST	INFIN1		;AND GO TO ERROR EXIT.

	SEGMENT	DATA
SAVE0:	BLOCK	2
SAVE2:	BLOCK	2
ARG11:	BLOCK	1
ARG12:	BLOCK	1
ARG21:	BLOCK	1
ARG22:	BLOCK	1
REAL:	BLOCK	1
ANSTMP:	BLOCK	1

	PRGEND
TITLE	CFDV	COMPLEX DIVISION
SUBTTL	KAREN KOLLING /KK/CKS		28-Oct-81

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 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.


;FROM LIB40 VERSION .027
;FROM	V.026	24-MAR-70

;THIS PROGRAM CALCULATES THE QUOTIENT OF TWO FLOATING POINT
;NUMBERS IN THE FOLLOWING MANNER:

;(A+IB)/(C+ID) = [(AC+BD)+I(BC-AD)]/(C**2+D**2)

;IF OVER OR UNDERFLOW OCCURS AT ANY POINT IN THE
;CALCULATION, THE PROGRAM JUMPS TO A MORE
;INVOLVED ROUTINE IN WHICH THE ORDER IN WHICH THE
;TERMS ARE MULTIPLIED OR DIVIDED TOGETHER IS VERY
;IMPORTANT.  OVER OR UNDERFLOW IS NOT PERMITTED
;UNLESS THE FINAL ANSWER OVER OR UNDERFLOWS.
;THERE ARE TWO SETS OF ENTRY POINTS FOR THE PROGRAM - ONE
;FOR THE DIVIDE ROUTINES AND ONE FOR THE DIVIDE-TO-MEMORY
;ROUTINES.  THE CALLING SEQUENCE FOR THE DIVIDE ROUTINES IS
;	XMOVEI	L,ARG2
;	PUSHJ	P,CFD.N
;WHERE N=0,2,4,6,10,12,14. ARG1 IS ASSUMED TO BE IN AC N AND (N+1),
;AND THE RESULTS ARE RETURNED IN AC N AND AC (N+1)

;THE CALLING SEQUENCE FOR THE DIVIDE-TO-MEMORY ROUTINES IS
;AS FOLLOWS:
;	XMOVEI	L,ARG2
;	PUSHJ	P,CDM.N
;WHERE N AND ARG2 ARE AS ABOVE AND THE RESULTS ARE RETURNED IN
;ARG2 AND ARG2+1.

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

DEFINE SYM (N) <
	ENTRY	CFD.'N,CDM.'N

	SIXBIT	/ CFD.'N/
CFD.'N:	DMOVEM	N,SAB		;STORE NUMERATOR
IFN N,<	DMOVEM	0,SAVE0	>	;SAVE REGS 0-1
	DMOVE	0,(L)		;GET SECOND ARG
	PUSHJ	P,CFD		;DIVIDE
IFN N,<	DMOVE	N,0		;STORE RESULT
	DMOVE	0,SAVE0 >	;RESTORE REGS 0-1
	POPJ	P,		;RETURN

	SIXBIT	/ CDM.'N/
CDM.'N:	DMOVEM	N,SAB		;STORE NUMERATOR
	DMOVEM	0,SAVE0		;SAVE 0-1
	DMOVE	0,(L)		;GET SECOND ARG
	PUSHJ	P,CFD		;DIVIDE
	DMOVEM	0,(L)		;STORE RESULT
	DMOVE	0,SAVE0		;RESTORE 0-1
	POPJ	P,
>

XLIST	;GENERATE CFD.0 ... CFD.14 AND CDM.0 ... CDM.14
.N.==16
REPEAT 7,<SYM(\<.N.==.N.-2>)>
LIST
CFD:	DMOVEM	T0,SCD		;STORE DENOMINATOR IN SCD AND LCD
	JUMPE	T0,CZERO	;GO TO CZERO IF C=0, TO
	JUMPE	T1,DZERO	;DZERO IF D=0,
	SKIPE	SAB		;TO NORMAL IF ONLY
	JRST	NORMAL		;ONE OF A AND B
	SKIPE	LAB		;=0, TO ABEQZR
	JRST	NORMAL		;IF
ABEQZR:	SETZB	T0,REAL		;A=B=0
	JRST	OUT1		;RETURN (0,0)

DZERO:	MOVE	T0,SAB		;D IS ZERO, SO REAL (ANS)
	FDVR	T0,SCD		;= A/C AND
	MOVEM	T0,REAL		;IMAG (ANS) =
	MOVE	T0,LAB		;B/C
	FDVR	T0,SCD		;AND THEN
	JRST	OUT1		;GO TO EXIT
CZERO:	MOVE	T0,LAB		;C IS ZERO (POSSIBLY D IS
	FDVR	T0,LCD		;ZERO ALSO). REAL (ANS) =
	MOVEM	T0,REAL		;B/D AND
	MOVN	T0,SAB		;IMAG (ANS) =
	FDVR	T0,LCD		;-A/D AND THEN
	JRST	OUT1		;GO TO EXIT.
NORMAL:	FMPR	T0,SCD		;THIS SIMPLE ROUTINE
	JFCL	NTSMPL		;CALCULATES
	FMPR	T1,LCD		;REAL (ANS) =
	JFCL	NTSMPL		;(A*C+B*D)/(C(2)+D(2))
	FADR	T0,T1		;AND
	JFCL	NTSMPL		;IMAG (ANS) =
	MOVEM	T0,TEMP		;(B*C-A*D)/(C(2)+D(2))
	MOVE	T0,SCD		;BUT
	FMPR	T0,SAB		;IF
	JFCL	NTSMPL		;AT
	MOVE	T1,LAB		;ANY
	FMPR	T1,LCD		;POINT
	JFCL	NTSMPL		;OVER
	FADR	T0,T1		;OR
	JFCL	NTSMPL		;UNDERFW
	FDVR	T0,TEMP		;OCCURS
	JFCL	NTSMPL		;IT
	MOVEM	T0,REAL		;JUMPS
	MOVE	T0,SAB		;TO
	FMPR	T0,LCD		;NTSMPL
	JFCL	NTSMPL		;FOR
	MOVE	T1,LAB		;A
	FMPR	T1,SCD		;DIFFERENT
	JFCL	NTSMPL		;CALCULATION.
	FSBRM	T1,T0		;IF THERE IS
	JFCL	NTSMPL		;OVER OR
	FDVR	T0,TEMP		;UNDERFLOW
	JRST	OUT1		;IT EXITS TO OUT1.

NTSMPL:	DMOVEM	T2,SAVE2	;SAVE THE CONTENTS OF AC2-3.
	MOVEM	T4,SAVE4	;SAVE THE CONTENTS OF AC4.
	MOVE	T2,SAB		;REARRANGE THE REAL
	MOVM	T0,SAB		;AND IMAG PARTS OF
	MOVM	T1,LAB		;THE 1ST ARG
	MOVEI	T4,1		;SO THAT THE SMALLER ONE
	CAMG	T0,T1		;(IN MAGNITUDE) IS IN
	JRST	NTSMP1		;SAB, AND THE OTHER IS IN
	EXCH	T2,LAB		;LAB AND SET UP AC4 AS
	MOVEM	T2,SAB		;A FLAG WORD TO TELL
	MOVN	T4,T4		;WHICH PART IS WHERE.
NTSMP1:	MOVE	T2,SCD		;REARRANGE THE REAL
	MOVM	T0,SCD		;AND IMAG PARTS OF THE OTHER ARG
	MOVM	T1,LCD		;SO THAT THE SMALLER ONE
	CAMG	T0,T1		;(IN MAGNITUDE) IS IN SCD AND
	JRST	NTSMP2		;THE OTHER IS IN LCD AND
	EXCH	T2,LCD		;FIX AC4 TO TELL WHICH
	MOVEM	T2,SCD		;PART
	IMULI	T4,3		;IS WHERE.
NTSMP2:	MOVE	T3,LCD		;CALCULATE
	FDVRB	T2,T3		;THE
	JFCL			;TERM
	FMPR	T2,T2		;1+(SCD/LCD)**2
	JFCL			;AND
	FADRI	T2,(1.0)	;STORE IT IN
	MOVEM	T2,DENOM	;DENOM.
	MOVE	T0,SAB		;CALCULATE
	FDVR	T0,LAB		;(SCD/LCD)*(SAB/LAB)
	JFCL			;SUPPRESSING
	FMPRM	T0,T3		;UNDERFLOW
	JFCL			;IF NECESSARY.
	CAIN	T4,1		;FIX THE AC FLAG WORD
	MOVNI	T4,2		;FOR
	ADDI	T4,1		;EASY COMPARISONS.
	SKIPL	T4		;CALCULATE
	MOVN	T3,T3		;+-1 +-(SCD/LCD)*(SAB/LAB),
	FADRI	T3,(1.0)	;WHERE THE SIGNS
	SKIPN	T4		;DEPEND ON
	MOVN	T3,T3		;THE AC FLAG WORD.
	PUSHJ	P,CALC34	;JUMP TO CALC OF (LAB/LCD)*(AC3/DENOM).
	MOVEM	T0,REAL		;STORE IT IN REAL(ANS) AND
	MOVEM	T0,IMAG		;IMAG(ANS)(THE TEMP. LOCATIONS).
	MOVE	T3,SAB		;CALCULATE
	MOVE	T2,SCD		;+-(SAB/LAB)+-(SCD/LCD)
	CAMN	T4,[-2]		;WHERE THE SIGNS
	MOVN	T2,T2		;AGAIN DEPENDS ON
	CAMN	T4,[-1]		;THE AC FLAG WORD,
	MOVN	T3,T3		;AND IF UNDERFLOW
	MOVE	T0,T2		;OCCURS JUMP TO
	MOVE	T1,T3		;THE
	FDVR	T3,LAB		;SPECIAL
	JFCL	OVER1		;ROUTINES-
	FDVR	T2,LCD		;OVER1,
	JFCL	OVER2		;OVER2,
ADD2:	FADR	T3,T2		;OR
	JFCL	OVER3		;OVER3.
CALCIR:	PUSHJ	P,CALC34	;JUMP TO CALC OF(LAB/LCD)*(AC3/DENOM).
STIR:	JUMPGE	T4,STR		;STORE THIS IN
	MOVEM	T0,IMAG		;THE CORRECT
	JRST	STX		;PART OF THE
STR:	MOVEM	T0,REAL		;ANSWER (TEMP. LOCATION).
STX:	DMOVE	T2,SAVE2	;RESTORE THE CONTENTS OF AC2-3.
	MOVE	T4,SAVE4	;RESTORE THE CONTENTS OF AC4.
OUT:	SKIPA	T1,IMAG		;PICK UP THE IMAG (ANS)
OUT1:	MOVE	T1,T0		;STORE IMAG IN T1
	MOVE	T0,REAL		;GET REAL PART
	POPJ	P,		;DONE
CALC34:	MOVM	T1,LAB		;CALC34 CALCS. (LAB/LCD)*(AC3/DENOM).
	MOVM	T2,LCD		;/LAB/ TO AC T1 AND /LCD/ TO AC 2.
	MOVM	T0,T3		;/AC3/ TO AC 0.
	CAMGE	T2,ONE		;GO TO CASEA IF
	JRST	CASEA		;/LCD/<1.0. O'E, STAY HERE.
	CAMGE	T1,ONE		;GO TO CASEB IF /LCD/>1.0 AND
	JRST	CASEB		;/LAB/<1.0 OR IF
	CAMGE	T0,ONE		;(>1)(<1)/(>1)(>1).
	JRST	CASEB		;STAY HERE IF (>1)(>1)/
	FDVR	T2,T1		;(>1)(>1).
STD:	FDVR	T0,DENOM	;CALCULATE
	FDVR	T0,T2		;IT AND GO
	JRST	SETSGN		;TO SETSGN.
CASEB:	FMPR	T0,T1		;CALCULATE CASE B AND
	JRST	STD		;GO TO SETSGN (EVENTUALLY).
CASEA:	FMPR	T2,DENOM	;CONDENSE THE DENOMINATOR INTO AC 2.
	CAMLE	T1,ONE		;THEN (<1)(<1)/(<1) GOES
	JRST	CHKAGN		;TO SR, AND (>1)(><1)/(<>1)
	CAMLE	T0,ONE		;GOES TO CHKAGN,
	JRST	NOTRVS		;AND EVERYTHING ELSE
	CAMG	T2,ONE		;GOES TO
	JRST	SR		;NOTRVS.

NOTRVS:	FMPRM	T1,T0		;CALCULATE
	JFCL	SETSGN		;NOTRVS AND GO
	FDVR	T0,T2		;TO
	JRST	SETSGN		;SETSGN.
CHKAGN:	CAMG	T0,ONE		;(>1)(<1)/(<>1)
	JRST	NOTRVS		;AND (>1)(>1)/(<1)
	CAMG	T2,ONE		;GO TO
	JRST	NOTRVS		;NOTRVS.
	FDVR	T1,T2		;(>1)(>1)/(>1) IS
	FMPRM	T1,T0		;CALCULATED AND
	JRST	SETSGN		;GOES TO SETSGN.
SR:	MOVEM	T1,TEMP		;SR CALCULATES
	FDVR	T1,T2		;(<1)(<1)/(<1)
	JFCL	OV1		;AND SINCE
	FMPRM	T1,T0		;(<1)/(<1)
	JRST	SETSGN		;CAN
OV1:	MOVE	T1,TEMP		;OVERFLOW, IT
	MOVEM	T0,TEMP		;CALCULATES
	FDVR	T0,T2		;IT
	JFCL	OV2		;WHICHEVER
	FMPRM	T1,T0		;WAY
	JRST	SETSGN		;IS
OV2:	MOVE	T0,TEMP		;NECESSARY
	FMPRM	T1,T0		;AND
	FDVR	T0,T2		;THEN GOES TO SETSGN.
SETSGN:	MOVE	T1,LAB		;GET THE
	XOR	T1,LCD		;SIGN OF THE
	XOR	T1,T3		;RESULT IN AC 1.
	SKIPG	T1		;SET THE SIGN OF
	MOVN	T0,T0		;THE ANSWER
	POPJ	P,		;AND EXIT.
OVER1:	FDVR	T2,LCD		;IF ONE
	JFCL	UU		;TERM
	MOVE	T3,T2		;UNDERFLOWS
	MOVE	T0,T1		;AND THE OTHER
	MOVE	T2,LAB		;TERM IS SO LARGE
	JRST	OVER21		;THAT NO BITS
OVER2:	MOVE	T2,LCD		;CAN BE
OVER21:	MOVEM	T2,SAVE5	;SAVED,
	JUMPE	T3,UU		;THEN
	MOVM	T1,T3		;RETURN
	CAML	T1,[030400000000] ;RIGHT
	JRST	CALCIR		;AWAY.
	FSC	T3,70		;O'E, TRY TO
	FSC	T0,70		;SAVE SOME BITS
	FDVRM	T0,T2		;BY MULTIPLYING
	JFCL			;THE TERMS BY 2**56,
	FADRB	T3,T2		;ADDING THE TERMS, AND THEN / BY 2**56.
	FSC	T3,-70		;IF THE RESULT UNDERFLOWS, GO
	JFCL	SRX		;TO SRX, O'E GO BACK
	JRST	CALCIR		;TO THE MAIN ROUTINE.

OVER3:	MOVE	T3,T1		;SET UP
	FDVR	T3,LAB		;AC 2 FOR
	FSC	T3,70		;SRX. AC 2 WILL
	FSC	T2,70		;CONTAIN THE TERM
	FADR	T2,T3		;*(2**56).
SRX:	MOVEM	T5,SAVE5	;SAVE THE CONTENTS OF AC 5.
	MOVE	T0,LCD		;THIS IS AN
	MOVE	T1,LAB		;ALTERNATE
	MOVM	T5,T1		;CALCULATION TO
	CAMGE	T5,ONE		;CALC34, AND TAKES
	JRST	SRX2		;ADVANTAGE OF
	FMPR	T1,T2		;THE FACT THAT HERE
	MOVM	T5,T0		;DENOM CONTAINS 1.0.
	CAML	T5,ONE		;THE ORDER OF
	JRST	SRY		;CALCULATION
	FSC	T0,70		;DEPENDS
	FDVRM	T1,T0		;ON THE SIZE OF
	JRST	SETSN2		;THE TERMS. AFTER
SRY:	FDVRM	T1,T0		;THE CALCULATION A
	FSC	T0,-70		;TRANSFER
	JRST	SETSN2		;IS
SRX2:	FDVRM	T2,T0		;MADE
	FMPR	T0,T1		;TO
	FSC	T0,-70		;SETSN2
SETSN2:	MOVE	T5,SAVE5	;RESTORE THE CONTENTS OF AC 5.
	JRST	STIR		;GO BACK TO MAIN ROUTINE.
UU:	MOVEM	T0,SAB		;ANOTHER ALTERNATE
	MOVEM	T1,SCD		;CALCULATION TO
	FMPR	T1,LCD		;CALC34
	FMPR	T0,LAB		;USED WHEN
	FADR	T0,T1		;S/L FOR
	JFCL	UND		;BOTH SETS
	FDVR	T0,LCD		;HAS UNDERFLOWED
	FDVR	T0,LCD		;OR FOR
	JRST	STIR		;UNDERFLOW PLUS
UND:	$LCALL	888
;LERR	(888,%,Complex underflow at $1L,<T1>)
	JRST	STIR		;TO ADD2+3.

ONE:	201400000000

	SEGMENT	DATA
SAVE0:	BLOCK	0
SAVE2:	BLOCK	1
SAVE3:	BLOCK	1
SAVE4:	BLOCK	1
SAVE5:	BLOCK	1
SAB:	BLOCK	1
LAB:	BLOCK	1
SCD:	BLOCK	1
LCD:	BLOCK	1
TEMP:	BLOCK	1
REAL:	BLOCK	1
IMAG:	BLOCK	1
DENOM:	BLOCK	1

	PRGEND
TITLE	REAL.C	COMPLEX TO REAL CONVERSION
SUBTTL	H. P. WEISS/HPW/CKS/CDM	17-Nov-81


;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 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.

;FROM V.005
;COMPLEX TO REAL CONVERSION FUNCTION
;THIS ROUTINE RETURNS THE REAL PART OF A COMPLEX NUMBER

;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
;	XMOVEI	L,ARGLST
;	PUSHJ	P,REAL.C	
;THE ANSWER IS RETURNED IN ACCUMULATOR A

	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	SALL

	ENTRY	REAL.C

REAL.C:	MOVE	T0,@(L)		;GET REAL PART OF ARGUMENT
	GOODBY	(1)		;RETURN

	PRGEND
TITLE	AIMAG	COMPLEX TO REAL CONVERSION
SUBTTL	H. P. WEISS/HPW/CKS	28-Oct-81


;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 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	AIMAG
EXTERN	AIMAG.
AIMAG=AIMAG.
PRGEND
TITLE	AIMAG.	COMPLEX TO REAL COVERSION
SUBTTL	H. P. WEISS/HPW/CKS	28-Oct-81


;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 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.


;FROM	LIB40	V.022
;FROM V.005
;COMPLEX TO IMAGINARY CONVERSION FUNCTION
;THIS ROUTINE RETURNS THE IMAGINARY PART OF A COMPLEX NUMBER

;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
;	XMOVEI	L,ARGLST
;	PUSHJ	P,AIMAG.
;THE ANSWER IS RETURNED IN ACCUMULATOR T0

	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	SALL

	HELLO	(AIMAG,.)	;ENTRY TO AIMAG ROUTINE
	XMOVEI	T1,@(L)		;PUT IMAG(ARG)
	MOVE	T0,1(T1)	;IN AC T0
	GOODBY	(1)		;RETURN

	PRGEND
TITLE	CONJG	COMPLEX CONJUGATE FUNCTION
SUBTTL	H. P. WEISS/HPW/CKS	28-Oct-81


;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 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	CONJG
EXTERN	CONJG.
CONJG=CONJG.
PRGEND
TITLE	CONJG.	COMPLEX CONJUGATE FUNCTION
SUBTTL	H. P. WEISS/HPW/CKS	28-Oct-81


;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 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.

;FROM	LIB40	V.022

;COMPLEX CONJUGATE FUNCTION
;THIS ROUTINE CALCULATES THE COMPLEX CONJUGATE OF ITS ARGUMENT

;THE REAL PART OF THE ANSWER IS LEFT IN ACCUMULATOR T0, AND THE
;IMAGINARY PART IN ACCUMULATOR T1

	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	SALL

	HELLO	(CONJG,.)	;ENTRY TO CONJG ROUTINE
	DMOVE	T0,@(L)		;GET COMPLEX ARGUMENT
	MOVN	T1,T1		;NEGATE IMAGINARY PART
	GOODBY			;RETURN

	PRGEND
TITLE	CMPLX	REAL TO COMPLEX CONVERSION
SUBTTL	H. P. WEISS/HPW/CKS	28-Oct-81


;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 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	CMPLX
EXTERN	CMPLX.
CMPLX=CMPLX.
PRGEND
TITLE	CMPLX.	REAL TO COMPLEX CONVERSION
SUBTTL	H. P. WEISS/HPW/CKS	28-Oct-81


;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 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.

;FROM	LIB40	V.022
;REAL TO COMPLEX CONVERSION FUNCTION
;THIS ROUTINE CONVERTS TWO REAL ARGUMENTS TO A COMPLEX NUMBER
;OF THE FORM ARG1+I*ARG2

;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
;	MOVEI	L,ARG
;	PUSHJ	P,CMPLX
;THE REAL PART OF THE ANSWER IS LEFT IN ACCUMULATOR A, AND THE
;IMAGINARY PART IS LEFT IN ACCUMULATOR B

	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	SALL

	HELLO	(CMPLX,.)	;ENTRY TO CMPLX ROUTINE
	MOVE	T0,@(L)		;GET REAL PART OF COMPLEX ANSWER
	MOVE	T1,@1(L)	;GET IMAGINARY PART OF COMPLEX ANS
	GOODBY	(2)		;RETURN


	PRGEND
	TITLE	CMPL.C	Conversion from complexes to complex
	SUBTTL	C. McCutcheon - 28-Oct-81

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 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.

; Fortran Instrinsic Function to return a complex number for
; two complex numbers passed to it.

; From ANSI-77 standard:
;	CMPLX(A1,A2) is the complex value whose real part is REAL(A1)
;	and whose imaginary part is REAL(A2).

; Required (called) routines: SNG.0, SNG.2

	SEARCH	MTHPRM
	SALL
	SEGMENT	CODE
	SALL

	HELLO	(CMPL.C,.)	;Enter routine

	MOVE	T0,@(L)		;Real half.
	MOVE	T1,@1(L)	;Imaginary half

	GOODBYE			;Return

	PRGEND
	TITLE	CMPL.D Conversion from double prec to complex
	SUBTTL	C. McCutcheon - 28-Oct-81

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 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.

; Fortran Instrinsic Function to return a complex number for
; two double precision numbers passed to it.

; Required (called) routines: SNG.0, SNG.2

	SEARCH	MTHPRM
	SALL
	SEGMENT	CODE
	SALL
	EXTERN	SNG.0		;Located in
	EXTERN	SNG.2		; FORDBL.MAC

	HELLO	(CMPL.D,.)	;Enter routine

	PUSH	P,T2		;Save Ac's
	PUSH	P,T3

	DMOVE	T0,@(L)		;Real half.
	PUSHJ	P,SNG.0		;Round DP to Real.

	DMOVE	T2,@1(L)	;Imaginary half.
	PUSHJ	P,SNG.2		;Round DP to Real.
				;(There is no SNG.1 for optimizing ac's)

	MOVE	T1,T2		;Set up complex #.

	POP	P,T3		;Restore AC's.
	POP	P,T2
	GOODBYE			;Return

	PRGEND
	TITLE	CMPL.G	Conversion from gfloating to complex
	SUBTTL	Randy Meyers -- Created by Edit [3201] 17-May-82 

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 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.

; Fortran Instrinsic Function to return a complex number for
; two gfloating numbers passed to it.

	SEARCH	MTHPRM
	SALL
	SEGMENT	CODE
	SALL

	HELLO	(CMPL.G,.)	;Enter routine

	PUSH	P,T2		;Save Ac

	DMOVE	T0,@0(L)	;Real half.
	EXTEND	T0,[GSNGL T0]	;Round gfloating to Real.

	DMOVE	T1,@1(L)	;Imaginary half.
	EXTEND	T1,[GSNGL T1]	;Round gfloating to Real.

	POP	P,T2		;Restore AC.
	GOODBYE			;Return

	PRGEND
	TITLE	CMPL.I	Conversion from integers to complex
	SUBTTL	C. McCutcheon - 28-Oct-81

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 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.

; Fortran Instrinsic Function to return a complex number for
; two integer numbers passed to it.

	SEARCH	MTHPRM
	SALL
	SEGMENT	CODE
	SALL

	HELLO	(CMPL.I,.)	;Enter routine

	MOVE	T0,@(L)		;Real half
	MOVE	T1,@1(L)	;Imaginary half

	FLTR	T0,T0		;Convert to 
	FLTR	T1,T1		; real

	GOODBYE			;Return

	END