SEARCH	FORPRM
	TV	FORCPX	COMPLEX ROUTINES ,6(2031)

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

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

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

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

	SUBTTL	REVISION HISTORY

COMMENT \

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

1351	EDS	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)

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

\

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

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

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

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

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

SEARCH	FORPRM
NOSYM
ENTRY	CABS
EXTERN	CABS.
CABS=CABS.
PRGEND

TITLE	CABS.	COMPLEX ABSOLUTE VALUE FUNCTION
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

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

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

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

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

;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	FORPRM
	TWOSEG	400000
	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.,<IFIW TEMP>	;SQUARE ROOT IS IN AC T0
     	FMPR	T0,T2		;*V 
	  JFCL	OVFL
RET:	POP	P,T2
	GOODBY	(1)		;RETURN

OVFL:	CAIN	T0,0		;OVERFLOW? - NOTE ERROR.  OTHERWISE,
	  JRST	RET		;UNDERFLOW, RETURN
	LERR	(LIB,%,<CABS: result overflow>,RET)

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

TWOM14:	163400000000		;2**-14

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

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

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

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

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

SEARCH	FORPRM
NOSYM
ENTRY	CEXP
EXTERN	CEXP.
CEXP=CEXP.
PRGEND

TITLE	CEXP.	COMPLEX EXPONENTIAL FUNCTION
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

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

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

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

;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. 36394.429 THE RESULT IS SET TO
;    (0.0,0.0) AND AN ERROR MESSAGE IS RETURNED.

;FOR X.LT.-89.415986 THE RESULT IS SET TO (0.0,0.0) AND AN ERROR
;    MESSAGE IS RETURNED.

;FOR X .GT. 88.029692;
;    IF Y = 0.0, THE RESULT IS SET TO (+INFINITY,0.0) AND AN
;    ERROR MESSAGE IS RETURNED.
;    IF 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	FORPRM
	TWOSEG	400000
	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
	LERR	(LIB,%,<CEXP: AIMAG(arg) too large in absolute value; result=(0,0)>)
	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.,<IFIW YSAVE>	;COSY=COS(Y)
	MOVEM	T0,COSY		
	FUNCT	SIN.,<IFIW YSAVE>	;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
	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
	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
	LERR	(LIB,%,<CEXP: imaginary part underflow>)

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

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

YMAX:	36394.4290E0
TBIG:   207540074636		;88.029692
NEGCON: 570232254037		;-89.415986

	RELOC			;DATA
SINY:   0
COSY:   0
XSAVE:  0
EXPX:   0
YSAVE:  0
	RELOC
	PRGEND

TITLE	CEXP2	COMPLEX ** INTEGER

;COPYRIGHT (C) 1981  BY  DIGITAL EQUIPMENT CORPORATION

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

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

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

SEARCH	FORPRM
NOSYM
ENTRY	CEXP2
EXTERN	CEXP2.
CEXP2=CEXP2.
PRGEND
TITLE	CEXP3	COMPLEX ** COMPLEX

;COPYRIGHT (C) 1981  BY  DIGITAL EQUIPMENT CORPORATION

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

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

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

SEARCH	FORPRM
NOSYM
ENTRY	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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(CEXP2,.)	;CEXP2, complex ** integer
	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.
	MOVE	T3,@1(L)	;N to T3.
	PUSHJ	P,DFL.3##	;T3 and T4 = DFLOAT(T3)
	JRST	CJOIN		;Join main routine


	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.
CJOIN:	JUMPN	T1,NXTST	;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).
	HRLOI	T0,377777	;Else results are indeterminate
	HRLOI	T1,377777	;Provide something for storage
	JUMPE	T3,BTHIND	;If A = 0, results indeterminate.
	LERR	(LIB,%,<CEXP3: (0,0)**(0,0) is indeterminate, result = (0,0)>)
	SETZB	T0,T1		;Return (0,0)
	JRST	RST234		;Restore registers 2, 3, and 4.
BTHIND:	LERR	(LIB,%,<CEXP3: (0,0)**(negative,any) is indeterminate, result = (infinity,infinity)>)
	JRST	RST234		;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.
	FDV	T0,T1		;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,150		;Scale Y by 2**96 and X by
	FSC	T4,-150		;  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,GETLN	;THETA is PI for Y > 0
	  DMOVN	T0,T0		;  or -PI for Y < 0.
	  MOVN	T2,T4		;  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	XIS0		;  underflow; Merge with XIS0.
	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:	DMOVE	T0,PIOV2	;If X = 0, or |Y| >> |X|
	MOVEM	T3,REXP		;Store real part of exponent
	MOVEM	T4,IEXP		;Store imaginary part of exponent
	SETZM	PHIFLG		;Initialize flag for underflow of PHI
	SETZM	RSIGN		;Clear flags for signs of 
	SETZM	ISIGN		;  final result
	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,150		;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.
	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
	DMOVE	T2,THETA	;Also get THETA and A and
	DMOVE	T4,REXP		;  scale each up by 2**96
	FSC	T2,150		;  so that resulting product
	FSC	T4,150		;  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
	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,150		;  up by 2**96, so that product
	FSC	T4,150		;  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,150		;  that product will be
	FSC	T4,150		;  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:	LERR	(LIB,%,<CEXP3: 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.
	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.
	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,-150		;  2**(-96) each, to get
	FSC	T4,-150		;  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	T4,150		;  so as to get product
	FSC	T4,150		;  scaled up by 2**192.
	DFMP	T0,T4		;No exception on product or
	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,-150		;  each down to get scaled
	FSC	T2,-150		;  A*LNRHO.
	DFMP	T0,T2		;No exception
	DMOVE	T2,THETA	;Get B and THETA
	DMOVE	T4,IEXP		;  and scale each down
	FSC	T2,-150		;  so get B*THETA
	FSC	T4,-150		;  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
	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.
	DFSB	T0,T2		;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.
	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,-150		;  2**-96 to remove scaling
	FSC	T4,-150		;  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.
	LERR	(LIB,%,<CEXP3: 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:	LERR	(LIB,%,<CEXP3: 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.
	CAMLE	T1,AMINLO	;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:	LERR	(LIB,%,<CEXP3: imaginary part underflow>)
	JRST	STIMP		;Rejoin main flow above.
;
;	>>>>  End of underflow on MUL by SIN(PHI)  <<<<
;
;	>>>>  Underflow on MUL by COS(PHI)  <<<<
;
RLUND:	LERR	(LIB,%,<CEXP3: 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.
	LERR	(LIB,%,<CEXP3:  real and imaginary 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.
	HLRZI	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.
	HLRZI	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:	LERR	(LIB,%,<CEXP3: imaginary part overflow>)
	HLROI	T0,377777	;Store zero extended
	SETZ	T1,		;  biggest number.
	JRST	STIMP1		;Return to main flow in ANS2 to
				;  get correct sign, and store.
;
IMUND1:	LERR	(LIB,%,<CEXP3: 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:	LERR	(LIB,%,<CEXP3: real part overflow>)
	HLROI	T0,377777	;Store zero extended
	SETZ	T1,		;  biggest number.
	JRST	STRLP1		;Return to main flow in ANS3 to
				;  get correct sign, and store.
;
RLUND1:	LERR	(LIB,%,<CEXP3: 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
	LERR	(LIB,%,<CEXP3: 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.
	LERR	(LIB,%,<CEXP3: 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.
	LERR	(LIB,%,<CEXP3: 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
	LERR	(LIB,%,<CEXP3: 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.
	LERR	(LIB,%,<CEXP3: 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.
	LERR	(LIB,%,<CEXP3: 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
	RELOC
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.

	RELOC			;Back to high segment

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

;COPYRIGHT (C) 1981  BY  DIGITAL EQUIPMENT CORPORATION

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

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

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

SEARCH	FORPRM
NOSYM
ENTRY	CLOG
EXTERN	CLOG.
CLOG=CLOG.
PRGEND

TITLE	CLOG.	COMPLEX NATURAL LOG FUNCTION
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	Mary Payne	16-Oct-81

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

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

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

;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	FORPRM
	TWOSEG	400000
	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
	LERR	(LIB,%,<CLOG: Imaginary part underflow>)

UND1:	JUMPN	T0,REST23	;DONE IF LN(|X| OR |Y|) NOT 0
	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:	LERR	(LIB,%,<entry CLOG; arg is (0.0,0.0); result=(+infinity,0.0)>)
	HRLOI	T0,377777	;REAL ANSWER IS POSITIVE 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
	RELOC		
TT0:	BLOCK	1			;FOR ATAN2 AND ALOG CALLS
TT2:	BLOCK	1			;FOR ATAN2 CALL
SMALSQ:	BLOCK	1
	RELOC
	PRGEND

TITLE	CSIN	COMPLEX SINE FUNCTION  
;		(SINGLE PRECISION)
SUBTTL	IMSL, INC.	JUNE 1, 1979

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

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

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

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

SEARCH	FORPRM
NOSYM
ENTRY	CSIN
EXTERN	CSIN.
CSIN=CSIN.
PRGEND

TITLE	CCOS	COMPLEX COSINE FUNCTION  
;		(SINGLE PRECISION)
SUBTTL	IMSL, INC.	JUNE 1, 1979

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

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

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

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

SEARCH	FORPRM
NOSYM
ENTRY	CCOS
EXTERN	CCOS.
CCOS=CCOS.
PRGEND

TITLE	CSIN.	COMPLEX SINE AND COSINE ROUTINES
;		(SINGLE PRECISION)
SUBTTL	IMSL, INC.	JUNE 1, 1979

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

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

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

;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) > 36394.429 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	FORPRM
	TWOSEG	400000
	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
	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:	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
	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
	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**28 * PI

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

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

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

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

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

SEARCH	FORPRM
NOSYM
ENTRY	CSQRT
EXTERN	CSQRT.
CSQRT=CSQRT.
PRGEND

TITLE	CSQRT.	COMPLEX SQUARE ROOT FUNCTION
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

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

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

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

;  LET Z = X + I*Y
;  THEN THE ANSWER 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	FORPRM
	TWOSEG	400000
	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.,<IFIW TEMP>	;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.,<IFIW TEMP>	;CALL TO SQRT
	MOVE	T3,T0		;SAVE SQRT(1+EPS/2) IN AC T3.
	MOVEM	T2,TEMP
	FUNCT	SQRT.,<IFIW TEMP>	;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.,<IFIW TEMP>	;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.,<IFIW TEMP>	;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.,<IFIW TEMP>	;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.,<IFIW TEMP>	;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
	  LERR	(LIB,%,<CEXP: real part underflow>,SIGNS)
	LERR	(LIB,%,<CEXP: imaginary part underflow>,SIGNS)

SQ2:	201552023632		;SQRT(2).

	RELOC			;DATA
TEMP:	0			;TEMPORARY STORAGE FOR SQRT ARGS
	RELOC
	PRGEND

TITLE	CFM	COMPLEX MULTIPLICATION
SUBTTL	KAREN KOLLING /KK/CKS		28-Oct-81

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

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

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

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

;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	FORPRM
	TWOSEG	400000
	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:	MOVE	T1,-2(P)	;GET RETURN PC
	SUBI	T1,1		;DECREMENT TO POINT TO CALL INST
	LERR	(APR,%,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,
	MOVE	T1,-2(P)	;GET RETURN PC
	SUBI	T1,1		;DECREMENT TO POINT TO CALL INST
	LERR	(APR,%,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.

	RELOC			;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
	RELOC
	PRGEND

TITLE	CFDV	COMPLEX DIVISION
SUBTTL	KAREN KOLLING /KK/CKS		28-Oct-81

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

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

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

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


;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	FORPRM
	TWOSEG	400000
	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:	MOVE	T1,-2(P)	;GET RETURN PC
	SUBI	T1,1		;DECREMENT TO GET ADDRESS OF CALL INST
	LERR	(APR,%,Complex underflow at $1L,<T1>)
	JRST	STIR		;TO ADD2+3.

ONE:	201400000000

	RELOC			;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
	RELOC
	PRGEND

TITLE	REAL	COMPLEX TO REAL CONVERSION
SUBTTL	H. P. WEISS/HPW/CKS	28-Oct-81


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

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

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

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

SEARCH	FORPRM
NOSYM
ENTRY	REAL
EXTERN	REAL.
REAL=REAL.
PRGEND

TITLE	REAL.	COMPLEX TO REAL CONVERSION
SUBTTL	H. P. WEISS/HPW/CKS	28-Oct-81


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

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

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

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

;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.	
;THE ANSWER IS RETURNED IN ACCUMULATOR A

	SEARCH	FORPRM
	NOSYM
	TWOSEG	400000
	SALL

	HELLO	(REAL,.)	;ENTRY TO REAL ROUTINE
	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) 1980,1981  BY  DIGITAL EQUIPMENT CORPORATION

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

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

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

SEARCH	FORPRM
NOSYM
ENTRY	AIMAG
EXTERN	AIMAG.
AIMAG=AIMAG.
PRGEND

TITLE	AIMAG.	COMPLEX TO REAL COVERSION
SUBTTL	H. P. WEISS/HPW/CKS	28-Oct-81


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

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

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

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


;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	FORPRM
	NOSYM
	TWOSEG	400000
	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) 1980,1981  BY  DIGITAL EQUIPMENT CORPORATION

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

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

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

SEARCH	FORPRM
NOSYM
ENTRY	CONJG
EXTERN	CONJG.
CONJG=CONJG.
PRGEND

TITLE	CONJG.	COMPLEX CONJUGATE FUNCTION
SUBTTL	H. P. WEISS/HPW/CKS	28-Oct-81


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

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

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

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

;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	FORPRM
	NOSYM
	TWOSEG	400000
	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) 1980,1981  BY  DIGITAL EQUIPMENT CORPORATION

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

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

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

SEARCH	FORPRM
NOSYM
ENTRY	CMPLX
EXTERN	CMPLX.
CMPLX=CMPLX.
PRGEND

TITLE	CMPLX.	REAL TO COMPLEX CONVERSION
SUBTTL	H. P. WEISS/HPW/CKS	28-Oct-81


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

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

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

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

;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	FORPRM
	NOSYM
	TWOSEG	400000
	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) 1981  BY  DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.

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

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

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

; 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	FORPRM
	SALL
	TWOSEG	400000
	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) 1981  BY  DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.

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

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

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

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

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

	SEARCH	FORPRM
	SALL
	TWOSEG	400000
	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.I	Conversion from integers to complex
	SUBTTL	C. McCutcheon - 28-Oct-81

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

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

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

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

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

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