PDP-10 Archive: forcpx.mac from AP-D480B-SB

Trailing-Edge - PDP-10 Archives - AP-D480B-SB_1978 - forcpx.mac
There are 3 other files named forcpx.mac in the archive. Click here to see a list.
TITLE	FORCPX	%4.(235) COMPLEX FUNCTION MODULE FOR FOROTS
SUBTTL	D. TODD /DRT/HPW		11-DEC-73



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
;FROM	V.32(425)
;FROM	V.021	26-SEP-69
;FROM	V.32(415)	LIB40

;COMPLEX EXP2 FUNCTION.

;THIS ROUTINE CALCULATES A COMPLEX NUMBER TO
;AN INTEGER POWER.

;THE BASIC ALGORITHM USED IS:

;	Z**B=(R**B)*[COS(B*TH)+I*SIN(B*TH)]

;WHERE Z=X+I*Y, R=/Z/, AND TH=ATAN2(Y,X).

;THE CALLING SEQUENCE IS:

;	MOVEI	16,POWER
;	PUSHJ	17,CEXP.2

;WHERE POWER IS THE ADDRESS OF THE INTEGER POWER,
;AND THE BASE IS IN AC'S 0 AND 1.  THE COMPLEX
;RESULT IS RETURNED IN AC'S 0 AND 1.

;BECAUSE SIN AND COS MUST PERFORM A RANGE REDUCTION
;ON THEIR ARGUMENTS, FOR LARGE B*TH THEY WILL
;RETURN ANSWERS OF ZERO.
	SEARCH	FORPRM

IFN F40LIB,<ENTRY	CEXP.2>
IFN F10LIB,<ENTRY	CEXP2.>


IFN F10LIB,<
	SIXBIT	/CEXP2./
CEXP2.:	DMOVE	T0,@(L)		;GET THE BASE
	MOVEI	L,@1(L)		;POINT TO THE POWER
IFN F40LIB,<
	JRST	CEXP.2		;GO TO THE COMMON ROUTINE
>>
IFN F40LIB,<
	SIXBIT	/CEXP.2/	;F40 ENTRY POINT
CEXP.2:
>
	MOVEM	0,SAVEXL	;SAVE X IN SAVEX.
	MOVEM	1,SAVEYR	;SAVE Y IN SAVEYR.
	JUMPN	0,ZNONZR	;JUMP AHEAD UNLESS
	JUMPN	1,ZNONZR	;Z=(0,0).
	MOVE	1,(16)		;PICK UP B IN AC 1.
	SETZ	0,		;SET REAL (ANS.)=0 UNLESS
	JUMPGE	1,.+4		;B < 0, IN WHICH CASE
	ERROR	(APR,5,1,.+1)	;MESSAGE AND REAL (ANS)=
	HRLOI	0,377777	;+ INFINITY.
	SETZ	1,		;IMAG (ANS)=0.
	POPJ	17,		;RETURN.

ZNONZR:	MOVE	0,(16)		;PICK UP B IN AC 0.
	MOVEM	0,SAVEB		;SAVE B IN SAVEB.
	MOVEM	2,SAVE2		;SAVE THE CONTENTS OF AC 2.
	SETZM	MINFLG		;SET THE -INFINITY FLAG TO 0.
IFE CPU-KI10,<FLTR	0,SAVEB		;USE THE HARDWARE>
IFE CPU-KA10,<CAME	0,[400000000000] ;THESE 4 LINES ARE A FIX FOR A****
	JRST	ZNONZ1		;HARDWARE BUG THAT MAKES FLOAT*****
	HRLZI	0,533400	;RETURN 0 FOR - INTEGER INFINITY.**
	JRST	ZNONZ2		;REMOVE THEM WHEN IT IS FIXED.*****
ZNONZ1:	FUNCT FLOAT.,<SAVEB>	;GET B IN FLOATING POINT
ZNONZ2:>
	MOVEM	0,BFLOAT	;IN BFLOAT
	MOVE	2,0		;AND IN AC 2.
	FUNCT ATAN2.,<SAVEYR,SAVEXL>	;FIND TH=
				;ATAN(X/Y)AND  MULT BY
	FMPR	2,0		;BAND STORE IT IN AC 2.
	FUNCT COS.,<2>	;FIND COS(B*TH)
				;AND STORE
	MOVEM	0,COSBTH	;IT IN COSBTH.
	FUNCT SIN.,<2>		;FIND SIN(B*TH)
	MOVEM	0,SINBTH	;IT IN SINBTH.
	MOVM	2,SAVEXL	;/X/ TO AC 2.
	MOVM	1,SAVEYR	;/Y/ TO AC 1.
	CAMLE	2,1		;PUT THE SMALLER OF /X/ AND
	EXCH	2,1		;/Y/ IN AC 2, THE OTHER IN AC 1.
	FDVR	2,1		;CALC S/L AND
	JFCL			;SUPPRESS THE ERROR MESSAGE.
	FMPR	2,2		;CALC (S/L) **2 AND
	JFCL			;SUPPRESS THE ERROR MESSAGE
	FADRI	2,201400	;CALC. (1.0 + (S/L)**2).


	MOVEM	1,SAVEXL	;SAVE L IN SAVEXL.
	FUNCT	SQRT.,<2>	;CALC. SQRT(1.0 + (S/L)**2)
				;AND LEAVE IT IN AC 0.
	MOVE	1,SAVEB		;PICK UP B IN AC 1.
	MOVE	2,0		;STORE THE SQRT IN AC 2.
	FMPR	0,SAVEXL	;CALC. R.
	JFCL	SEPAR		;IF OVERFLOW, GO TO SEPAR.
	MOVEM	0,SAVEYR	;STORE R IN SAVEYR.
	MOVSI	2,201400	;SET UP 1.0 IN AC 2.
	JUMPGE	1,XP2		;GO TO R**B CALC IF B>=0.
	MOVMS	1,1		;GET /B/. IF B WAS
	JFCL	MININF		;-INFINITY GO TO MININF.
	JRST	XP2		;GO TO R**B CALC.
XP1:	FMPR	0,0		;FORM R**N.
	JFCL	OVER		;IF OVER/UNDERFLOW, GO TO OVER.
	LSH	1,-1		;SHIFT EXP. FOR NEXT BIT.
XP2:	TRZE	1,1		;IS THE BIT ON?
	FMPR	2,0		;YES, MULTIPLY ANSWER BY R**N.
	JFCL	OVER		;IF OVER/UNDERFLOW, GO TO OVER.
	JUMPN	1,XP1		;UPDATE R**N UNLESS ALL FINISHED.
	MOVE	0,2		;PICK UP THE RESULT.
	MOVEM	0,1		;PUT IT IN AC 1 ALSO.
	SKIPL	SAVEB		;IF B WAS >=0,
	JRST	.+4		;GO AHEAD.
	MOVSI	1,201400	;O'E, INVERT THE CALCULATED
	FDVRB	1,0		;QUANTITY AND IF IT
	JFCL	OVER2		;OVER/UNDER FLOWS, GO TO OVER2.
	FMPR	0,COSBTH	;FORM REAL (ANS) AND GO
	JFCL	TEMP1		;TO TEMP1 IF IT UNDERFLOWS.
SECOND:	FMPR	1,SINBTH	;FORM IMAG (ANS) AND GO
	JFCL	TEMP2		;TO TEMP2 IF IT UNDERFLOWS.
OUT:	MOVE	2,SAVE2		;RESTORE THE CONTENTS OF AC 2.
	POPJ	17,		;RETURN.
OVER2:	MOVE	2,.JBTPC		;PICK UP THE FLAGS AND
	TLC	2,(1B11)	;COMPLEMENT THE UNDERFLOW FLAG
	JRST	.+2		;JUMP AHEAD.
OVER:	MOVE	2,.JBTPC		;PICK UP THE FLAGS.
	SKIPE	MINFLG		;IF NECESSARY, RESTORE THE
	SUB	17,[XWD 1,1]	;PUSH DOWN LIST COUNTER.
	SKIPG	BFLOAT		;IF TRUE OVERFLOW
	TLC	2,(1B11)	;OCCURED, JUMP TO
	TLNN	2,(1B11)	;ALOGRT. O'E,
	JRST	ALOGRT		;STAY HERE AND
	SKIPE	COSBTH		;RETURN AN ERROR MESSAGE UNLESS
	PUSHJ	17,ERR		;REAL(ANS) IS IDENTICALLY 0.
	SKIPE	SINBTH		;RETURN AN ERROR MESSAGE UNLESS
	PUSHJ	17,ERR		;IMAG(ANS) IS IDENTICALLY 0.
	SETZB	0,1		;ANS=(0,0)
	JRST	OUT		;GO TO RETURN.


ALOGRT:
	FUNCT	ALOG.,<SAVEYR>	;CALC. THE LOG(R)
ALGRT2:	FMPR	0,BFLOAT	;IT BY B.
	MOVEM	0,TEMP		;STORE IT IN TEMP.
	SKIPN	SINBTH		;IF SIN(B*TH)=0,
	JRST	ZRIMAG		;GO TO ZRIMAG.
	MOVM	2,SINBTH	;/SIN(B*TH)/ TO AC 2.
	FUNCT	ALOG.,<2>	;CALC. THE LOG
	FADR	0,TEMP		;ADD IT TO THE OTHER TERM.
	MOVE	2,0		;PUT IT IN AC NE 0 OR 1.
	FUNCT	EXP.,<2>		;CALC THE IMAGINARY
	SKIPGE	SINBTH		;ANSWER AND GIVE IT
	MOVNS	0,0		;THE CORRECT SIGN.
	MOVEM	0,IMAG		;STORE IMAG(ANS) IN IMAG.


REAL:	SKIPN	COSBTH		;IF COS(B*TH)=0,
	JRST	ZRREAL		;GO TO ZRREAL.
	MOVM	2,COSBTH	;/COS(B*TH)/ TO AC 2
	FUNCT	ALOG.,<2>	;CALC. THE LOG
				;(/COS(B*TH)/) AND
	FADR	0,TEMP		;ADD IT TO THE OTHER TERM.
	MOVE	2,0		;PUT IT IN AC NE 0 OR 1.
	FUNCT	EXP.,<2>		;CALC THE REAL PART OF
	SKIPGE	COSBTH		;IT THE
	MOVNS	0,0		;CORRECT SIGN.
	MOVE	1,IMAG		;SET UP IMAG(ANS).
	JRST	OUT		;GO TO RETURN.

ZRIMAG:	SETZM	IMAG		;IMAG(ANS)=0
	JRST	REAL		;GO BACK TO CALC OF REAL(ANS).

ZRREAL:	SETZ	0,		;REAL(ANS)=0
	MOVE	1,IMAG		;SET UP IMAG(ANS).
	JRST	OUT		;GO TO RETURN


SEPAR:
	FUNCT	ALOG.,<2>	;CALC. THE LOG (SQRT(
				;1+(S/L)**2)) AND
	MOVEM	0,TEMP		;STORE IT IN TEMP.
	FUNCT	ALOG.,<SAVEXL>	;CALC THE LOG(L)
	FADR	0,TEMP		;OTHER LOG.
	JRST	ALGRT2		;GO TO EXPANDED CALC.

MININF:
	AOS	MINFLG		;SET B=-INFINITY FLAG.
	HRLOI	1,377777	;SET B=+INFINITY AND
	PUSHJ	17,XP2		;GO TO CALC R**B.
	FDVR	0,SAVEYR	;CALC. REAL(ANS).
	FDVR	1,SAVEYR	;CALC. IMAG(ANS).
	JRST	OUT		;GO TO RETURN.

MINFT:	SKIPE	MINFLG		;IF B=-INFINITY,
	JRST	ALOGRT		;GO TO ALOGRT.

ERR:	ERROR	(APR,7,1,.+1)	;MESSAGE AND GO BACK
	POPJ	17,		;TO CALC.

TEMP1:	MOVEM	1,TEMP		;TYPER. DOES NOT SAVE AC 1.
	PUSHJ	17,MINFT	;JUMP TO MINFT.
	MOVE	1,TEMP		;RESTORE AC 1.
	JRST	SECOND		;RETURN TO CALC.

TEMP2:	PUSHJ	17,MINFT	;JUMP TO MINFT.
	SETZM	1		;TYPER. DOES NOT SAVE AC 1.
	JRST	OUT		;GO TO RETURN.

SAVEXL:	0
SAVEYR:	0
SAVEB:	0
SAVE2:	0
MINFLG:	0
BFLOAT:	0
COSBTH:	0
SINBTH:	0
IMAG:	0
TEMP:	0
	PRGEND
TITLE	CEXP.3	%5.(572) ;FROM LIB40 VERSION 	V.023 14-JULY-72
;FROM	V.023	LIB40


;COMPLEX NUMBER TO A COMPLEX POWER FUNCTION.

;THE ALGORITHM USED IS:

;(X+IY)**(A+IB)=

;EXP[A*LOG(R)-B*TH+LOG(TRIG(A*TH+B*LOG(R)))]

;WHERE TRIG=COS FOR THE REAL PART AND SIN
;FOR THE IMAGINARY PART., AND WHERE
;R=(X**2+Y**2)**0.5 AND
;TH=ATAN2(Y,X).

;THIS ROUTINE USES THE FACT THAT EXP(N) WHERE
;N>ABOUT 89.0 IN MAGNITUDE OVER OR UNDERFLOWS,
;AND THAT COS AND SIN OF LARGE ARGUMENTS RETURN
;ZERO BECAUSE OF THEIR RANGE REDUCTION.

;THE CALLING SEQUENCE FOR THIS ROUTINE IS
;	MOVEI	16,ADDR
;	PUSHJ	17,CEXP.3
;THE BASE IS IN AC'S 0 AND 1 AND THE POWER
;IS IN ADDR AND ADDR+1.
;THE ANSWER IS RETURNED IN AC'S 0 AND 1.



;REVISION HISTORY
;EDIT	SPR	COMMENT
;----	---	-------
;572	-----	ALLOW CEXP TO LOAD IF F40LIB TURNED OFF
;****** END OF REVISION HISTORY ******
	SEARCH	FORPRM
IFN F40LIB,<ENTRY CEXP.3>
IFN F10LIB,<ENTRY CEXP3.>

	EXTERNAL EXP3..	;[572] FLOATING TO FLOATING POWER

	MLON
IFN F10LIB,<
	SIXBIT	/CEXP3./	;FORTRAN 10 ENTRY
CEXP3.:	DMOVE	T0,@(L)	;GET THE BASE
	MOVEI	L,@1(L)		;GET THE POWER POINTER
IFN F40LIB,<
	JRST	CEXP.3		;COMMON ROUTINE
>>
IFN F40LIB,<
	SIXBIT	/CEXP.3/	;F40 ENTRY
CEXP.3:>
	JUMPN	0,.+3		;IF THE BASE = (0,0),
	JUMPN	1,.+2		;THEN THE ANSWER IS SET = (0,0),
	POPJ	17,		;AND THE ROUTINE EXITS.
	MOVEM	2,SAVE2		;SAVE THE CONTENTS OF AC 2.
	MOVEM	3,SAVE3		;SAVE THE CONTENTS OF AC 3.
	MOVE	2,(16)		;REAL(POWER) TO AC 2.
	MOVE	3,1(16)		;IMAG(POWER) TO AC 3.
	JUMPN	3,NORMAL	;IF IMAG(BASE) = 0 AND
	JUMPN	2,RLTORL	;IMAG(POWER) = 0, THEN GO TO RLTORL.
	HRLZI	0,201400	;IF THE BASE NE (0,0) AND THE POWER
	SETZM	1		;= (0,0), THE ANSWER =
	JRST	OUT1		;(1.0,0.0).
OUT:	MOVE	4,SAVE4		;RESTORE THE CONTENTS OF AC 4.
OUT1:	MOVE	2,SAVE2		;RESTORE THE CONTENTS OF AC 2.
	MOVE	3,SAVE3		;RESTORE THE CONTENTS OF AC 3
	POPJ	17,		;EXIT.
RLTORL:	JUMPN	1,NORMAL	;RLTORL IS IMAG OF BASE AND POWER = 0.
	JUMPL	0,NORMAL	;BASE MUST BE >=0 HERE.
	MOVE	1,2		;SET UP ARGS FOR EXP3.
;**; [572] CHANGE @ RLTORL+3	CLRH	2-AUG-76
	PUSHJ	17,EXP3..	;[572] CALC REAL(ANS).
	SETZM	1		;IMAG(ANS) = 0.
	JRST	OUT1		;GO TO EXIT.
NORMAL:	MOVEM	4,SAVE4		;SAVE THE CONTENTS OF AC 4.
	MOVEM	0,BASE1		;SAVE THE BASE IN BASE1
	MOVEM	1,BASE2		;AND BASE2.
	MOVEM	2,POWER1	;SAVE THE POWER IN POWER1
	MOVEM	3,POWER2	;AND POWER2.
	SETZM	COSFLG		;ZERO TWO OF
	SETZM	SINFLG		;THE FLAGS.
	FDVR	1,0		;JUMP TO THAUFL IF
	JFCL	THAUFL		;OV/UN/DIV TROUBLE.
N2:	FUNCT	CLOG.,<BASE1>	;CALC THE LOG(R)
	MOVEM	0,LOGR		;AND STORE
	CAML	1,[202622077325] ;THEM
	FSBR	1,[203622077325];IN LOGR
 	MOVEM	1,THETA		;AND THETA.
N1:	PUSHJ	17,ANGLCL	;CALC THE LOG(TRIG).
	SKIPE	OUTFLG		;IF THE ANS. HAS BEEN SETUP, GO TO OUT,
	JRST	OUT		;O'E, STAY HERE.

MAGCLC:	MOVN	0,POWER2	;CALC -B*THETA,
	FMPR	0,THETA		;UNDERFLOW IS OK,
	JFCL	[JUMPE 0,.+1	;GO TO OVER1 IF
		JRST OVER1]	;OVERFLOW OCCURRED.
	MOVE	1,POWER1	;CALC A*LOGR,
MAG1:	FMPR	1,LOGR		;JUMP TO TEMP1
	JFCL	TEMP1		;IF UNDER/OVERFLOW.
MAG2:	FADRB	0,1		;CALC.A*LOGR - B*THETA AND
	JFCL	TEMP2		;GO TO TEMP2 WHEN OVER/UNDERFLOW.
MAG3:	FADRM	0,LOGCOS	;CALC ARG FOR EXP -- REAL
	JFCL			;AND NEVER MIND OVER/UNDERFLOW.
	FADRM	1,LOGSIN	;CALC ARG FOR EXP -- IMAG
	JFCL			;AND NEVER MIND OVER/UNDERFLOW.
TEST:	SKIPE	COSFLG		;COSFLG AND SINFLG ARE 0, UNLESS
	JRST	SETSG3		;COS AND SIN ARE 0, RESPECTIVELY.
	FUNCT	EXP.,<LOGCOS>	;CALC REAL(ANS)
	MOVEM	0,4		;IT IN 4(REAL).
	SKIPE	SINFLG		;SEE NOTE FOR
	JRST	SETSGN		;TEST ABOVE.
SETSG3:	FUNCT	EXP.,<LOGSIN>	;CALC IMAG(ANS)
	MOVEM	0,IMAG		;IT IN IMAG.
SETSGN:	SKIPGE	COSINE		;SET THE SIGN
	MOVNS	4,4		;OF REAL(ANS).
	SKIPGE	SINE		;SET THE SIGN
	MOVNS	IMAG		;OF IMAG(ANS).
STORE:	MOVE	0,4		;RETURN REAL(ANS).
	MOVE	1,IMAG		;RETURN IMAG(ANS).
	JRST	OUT		;GO TO EXIT.

ANGLCL:	SETZM	OUTFLG		;ZERO THE OUT FLAG.
	FMPR	2,1		;CALC A*THETA AND GO TO
	JFCL	ANGL1		;ANGL1 IF OVER/UNDERFLOW OCCURS.
	FMPR	0,POWER2	;CALC B*LOGR AND GO TO
	JFCL	ANGL2		;ANGL2 IF OVER/UNDERFLOW OCCURS.
	FADR	2,0		;CALC A*THETA + B*LOGR AND GO TO
	JFCL	ANGL3		;ANGL3 IF OVER/UNDERFLOW OCCURS.
TRIG:	FUNCT	SIN.,<2>		;CALC SIN(ARG) AND
	MOVEM	0,SINE		;IN SINE.
	FUNCT	COS.,<2>		;CALC COS(ARG) AND
	MOVEM	0,COSINE	;IN COSINE.
	JUMPN	0,LGCOSC	;IF COS=0, THEN REAL=0.
COSZ:	AOS	COSFLG		;SO SET COSFLG NE 0
	SETZM	4		;AND SET UP REAL AND GO
	JRST	CHKSIN		;AHEAD TO CHKSIN.
LGCOSC:	MOVM	2,COSINE	;COS NE 0, SO
	FUNCT	ALOG.,<2>	;CALC LOG(/COS/)
	MOVEM	0,LOGCOS	;IN LOGCOS.
CHKSIN:	SKIPE	SINE		;IF SIN = 0
	JRST	LGSINC		;THEN IMAG = 0,

SINZ:	AOS	SINFLG		;SO SET SINFLG NE 0
	SETZM	IMAG		;AND SET UP IMAG.
	SKIPE	COSFLG		;IF THE ENTIRE ANSWER IS 0,
	JRST	BTHZRO		;GO TO BTHZRO. O'E,
	POPJ	17,		;RETURN.
LGSINC:	MOVM	2,SINE		;SIN NE 0, SO
	FUNCT	ALOG.,<2>	;CALC LOG(/SIN/)
	MOVEM	0,LOGSIN	;THEN JUMP OUT OF THIS SECTION.
	POPJ	17,		;******END OF CODE FOR MAIN CALC******

OVER1:	MOVE	1,POWER1	;CALC A*LOGR
	FMPR	1,LOGR		;AND IF ANYTHING
	JFCL	[JUMPE 1,.+1	;BUT B*THETA AND A*LOGR
		XOR 1,0		;BOTH OVERFLOW AND THEY HAVE
		JUMPGE 1,.+1	;DIFFERENT SIGNS, GO AHEAD.
		JRST	OVER2]	;FOR THAT SPECIAL CASE, GO TO OVER2.
TEMP4:	MOVEM	0,LOGCOS	;ARG IS ALREADY SO LARGE THAT
	MOVEM	0,LOGSIN	;LOG(TRIG) DOESN'T MATTER, SO STORE
	JRST	TEST		;RESULTS AND GO BACK TO CALC.
OVER2:	MOVN	0,POWER2	;SPECIAL CASE, OVERFLOW
	MOVE	3,LOGR		;PLUS OVERFLOW, WITH
	MOVE	2,THETA		;DIFFERENT SIGNS.
	PUSHJ	17,FOUR		;CALC IT
	JUMPN	0,TEMP4		;AND THEN
	JRST	TEST		;RETURN.

FOUR:	FDVR	0,POWER1	;THIS SECTION CALCULATES
	MOVE	1,3		;X1*X2 - X3*X4
	FDVR	1,2		;IN
	MOVEM	0,TEMP		;DIFFERENT
	FADR	0,1		;WAYS
	JFCL	OVER3		;DEPENDING
	FMPR	0,POWER1	;ON
	JFCL			;WHERE
	FMPR	0,2		;OVERFLOW
	JFCL			;AND
	POPJ	17,		;UNDERFLOW
OVER3:	MOVE	0,2		;OCCUR.  IT IS
	FMPR	0,TEMP		;CALLED
	JFCL			;WITH
	FADR	0,3		;PUSHJ AND
	FMPR	0,POWER1	;EXITS
	JFCL			;WITH
	POPJ	17,		;POPJ.

TEMP1:	JUMPE	1,MAG2		;IGNORE UNDERFLOW,
	JRST	TEMP3		;FOR OVERFLOW PLUS
TEMP2:	JUMPE	1,MAG3		;NORMAL, THIS ANSWER
TEMP3:	MOVEM	1,LOGCOS	;IS OVERFLOW, SO
	MOVEM	1,LOGSIN	;SET IT UP, AND
	JRST	TEST		;RETURN.


ANGL1:	FMPR	0,POWER2	;IF O/U + O/U,
	JFCL	ANGL11		;GO TO ANGL11. O'E, STAY HERE.
	JUMPN	2,BTHZRO	;IF O+N, ANS IS O, GO TO BTHZRO.
	MOVE	2,0		;IF U+N WHERE N>EPS,
	MOVE	0,POWER1	;THEN
BTHM10:	MOVM	3,2		;ANS = N, GO TO TRIG.
	CAML	3,EPS		;IF U+N WHERE N< EPS,
	JRST	TRIG		;THEN
	FMPRI	1,271400	;STAY HERE AND TRY TO
	FMPR	1,0		;SAVE
	JFCL			;SOME BITS
	FMPRI	2,271400	;THEN
	FADR	2,1		;GO TO
	JRST	SUMUFL		;SUMUFL.
BTHZRO:	SETZB	0,1		;ANSWER = (0,0), SO
	AOS	OUTFLG		;SET THE OUT FLAG
	POPJ	17,		;AND EXIT.

ANGL11:	JUMPN	2,.+3		;IF U+U, GO TO UU.
	JUMPE	0,UU		;IF U+O, GO TO
	JRST	BTHZRO		;BTHZRO.
	JUMPE	0,BTHZRO	;IF O+O WITH THE
	XOR	2,0		;SAME SIGN, GO TO BTHZRO.
	JUMPGE	2,BTHZRO	;IF O+O WITH DIFFERENT SIGNS,
	MOVE	0,POWER2	;STAY HERE.
	MOVE	3,THETA		;FOUR CALCULATES
	MOVE	2,LOGR		;THE SUM OF THE TWO
	PUSHJ	17,FOUR		;PRODUCTS -- THEN
	MOVE	2,0		;GO BACK TO
	JRST	TRIG		;TRIG.

UU:	MOVE	0,LOGR		;THIS
	MOVE	1,THETA		;IS
	MOVE	4,LG2126	;THE UNDERFLOW
	FMPRI	0,377400	;PLUS
	FMPRI	1,377400	;UNDERFLOW
	MOVEM	0,TEMP		;CASE --
	MOVEM	1,IMAG		;MULTIPLY
	FMPR	0,POWER2	;IT
	JFCL	U1		;UP
	FMPR	1,POWER1	;AND
	JFCL	U2		;TRY
	MOVEM	0,2		;TO
UU20:	FADR	0,1		;SAVE
	JFCL	DOUBL		;IT.
UU21:	MOVMM	0,2		;4
UU2:	MOVEM	0,SINE		;CONTAINS
	FUNCT	ALOG.,<2>	;A
	SETZM	LOGCOS		;FACTOR
	FSBR	0,4		;TO
	MOVEM	0,LOGSIN	;BE
	POPJ	17,		;SUBTRACTED.

U1:	FMPR	1,POWER1	;CONTINUATION
	JFCL	DOUBL		;OF UNDERFLOW PLUS
	MOVE	0,1		;UNDERFLOW CASE --

U2:	MOVM	2,0		;MULTIPLYING
	CAML	2,EPS		;UP AS
	JRST	UU2		;MUCH

DOUBL:	MOVE	0,POWER2	;AS
	MOVE	1,POWER1	;IS
	FMPRI	0,377400	;NECESSARY.
	FMPRI	1,377400	;DOUBL
	FMPR	0,TEMP		;IS
	FMPR	1,IMAG		;ALSO USED
	FADR	4,LG2126	;BY
	JRST	UU20		;SEVERAL SECTIONS.

ANGL2:	JUMPN	0,BTHZRO	;O+N=O, ANS = (0,0).
	MOVE	1,LOGR		;SET UP FOR
	MOVE	0,POWER2	;U+N TESTS
	JRST	BTHM10		;AND JUMP TO THEM.

ANGL3:	JUMPN	2,BTHZRO	;IF THE SUM OVERFLOWS, ANS=(0,0).
	MOVE	2,POWER1	;THIS
	FMPR	2,THETA		;IS
	FMPRI	2,271400	;UNDERFLOW --
	FMPRI	0,271400	;TRY TO
	FADR	2,0		;SAVE IT.

SUMUFL:	MOVEM	2,SINE		;THE FINAL ARG
	MOVM	2,SINE		;HAD TO BE
	FUNCT	ALOG.,<2>	;SAVED.
	FSBR	0,LG256		;SUBTRACT
	MOVEM	0,LOGSIN	;56.0*LOG(2)BASE2
	SETZM	LOGCOS		;AT THE
	POPJ	17,		;END.


THAUFL:	JUMPE	0,DIVCHK	;GO TO DIVCHK IF X = 0.
	JUMPN	1,N2		;RETURN IMMED. IF REALLY OVFLW.
	MOVM	3,BASE1		;GET
	FUNCT	ALOG.,<3>	;LOGR =
				;LOG (/X/)
	MOVEM	0,LOGR		;AND STORE IT.
	MOVM	1,POWER1	;/A/ < 1.0 ?
	CAML	1,[1.0]		;NO, GO
	JRST	AGTHN1		;TO AGTHN1.
	SETZB	1,THETA		;YES, THETA = 0
	JRST	N1		;AND GO BACK TO MAIN. CALC.

AGTHN1:	MOVE	1,BASE2		;SWAP ARGS., SO
	MOVEM	1,THETA		;THAT A = A/X
	MOVEM	2,TEEMP1	;AND
	FDVR	2,BASE1		;THETA = Y
	MOVE	2,POWER1	;AND CALL THE
	PUSHJ	17,ANGLCL	;ANGLE CALCULATION.
	SKIPE	OUTFLG		;IF (0,0), THEN
	JRST	OUT		;EXIT.
	SETZM	THETA		;SET THETA = 0, AND
	MOVE	1,TEEMP1	;RESTORE A = A
	MOVEM	1,POWER1	;AND GO BACK TO
	JRST	MAGCLC		;MAIN CALC.

DIVCHK:	MOVM	3,BASE2		;X = 0, SO
	FUNCT	ALOG.,<3>	;LOGR =
				;LOG(/Y/)
	MOVEM	0,LOGR		;AND THETA =
	MOVE	1,[201622077325] ;+-
	SKIPGE	BASE2		;PI/2.
	MOVNS	1		;SET UP
	MOVEM	1,THETA		;LOCATIONS AND RETURN
	JRST	N1		;TO MAIN CALC.



EPS:	033400000000	
SAVE2:	0		
SAVE3:	0		
SAVE4:	0
BASE1:	0		
BASE2:	0		
POWER1:	0		
POWER2:	0		
COSFLG:	0		
SINFLG:	0		
LOGR:	0		
THETA:	0		
COSINE:	0		
SINE:	0		
IMAG:	0		
LOGCOS:	0		
LOGSIN:	0		
TEMP:	0		
OUTFLG:	0		
TEEMP1:	0		
LG2126:	207535261175
LG256:	206466417250
	PRGEND
TITLE	CSQRT	%4.(235) COMPLEX SQURE ROOT
SUBTTL	H. P. WEISS/HPW		11-DEC-73



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
SEARCH	FORPRM
ENTRY	CSQRT
EXTERN	CSQRT.
CSQRT=CSQRT.
PRGEND
TITLE	CSQRT.	%4.(235) COMPLEX SQURE ROOT
SUBTTL	D. TODD/DRT/HPW	30-NOV-1973



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
;FORM	V.21	LIB40

;COMPLEX SQUARE ROOT FUNCTION

;THE ALGORITHM IS AS FOLLOWS:

;IF X >= 0 AND Y >= 0.,
;	A=SQRT((/X/+/X+IY/)/2)
;	B=Y/(2*A)

;IF X >= 0 AND Y < 0,
;	A=-SQRT((/X/+/X+IY/)/2)
;	B=Y/(2*A)

;IF X < 0 AND Y >= 0,
;	A=Y/(2*B)
;	B=SQRT((/X/+/X+IY/)/2)

;IF X < 0 AND Y < 0,
;	A=Y/(2*B)
;	B=SQRT((/X/+/X+IY/)/2)

;WHERE /Z/ MEANS THE ABSOLUTE VALUE OF Z.

;BECAUSE OF OVER/UNDERFLOW PROBLEMS, THE SQUARE
;ROOT TERM IS CALCULATED IN SEVERAL DIFFERENT WAYS.

;THE CALLING SEQUENCE IS:
;	JSA	16,CSQRT
;	EXP Z
;WHERE Z HAS BEEN DECLARED COMPLEX.

;THE REAL PART OF THE ANSWER IS RETURNED IN AC 0,
;AND THE IMAGINARY PART OF THE ANSWER IS RETURNED IN AC 1.

	SEARCH	FORPRM




	HELLO	(CSQRT,.)	;[235] ENTRY TO COMPLEX SQUARE ROOT ROUTINE.
	MOVEI	1,@0(16)	;[227] PICK UP ADDRESS OF Z
	MOVE	0,(1)		;PICK UP X IN AC 0.
	MOVE	1,1(1)		;PICK UP Y IN AC 1.
	MOVEM	2,SAVE2		;SAVE AC 2.
	MOVM	2,0		;/X/ TO AC 2.
				;IF Y IS NOT ZERO,
	JUMPN	1,NORMAL	;GO TO THE STANDARD ROUTINE.
				;IF Y=0 AND X>=0, Z IS REAL,
	JUMPGE	0,REAL		;SO TRANSFER TO SIMPLE CALC.
NORMAL:	MOVEM	3,SAVE3		;SAVE AC 3.
	MOVEM	0,XSAVE		;SAVE X.
	MOVEM	1,YSAVE		;SAVE Y.
	MOVM	3,1		;/Y/ TO AC 3.
	SETZM	FLAGWD		;FLAGWD CONTAINS 0 IF /X/
	CAMLE	2,3		;<= /Y/, O'E IT CONTAINS 1.  PUT
	AOSA	FLAGWD		;THE LARGER OF /X/,/Y/ IN AC 2
	EXCH	2,3		;AND THE SMALLER IN AC 3.
	FDVR	3,2		;CALC S/L.
	JFCL			;NO UNDERFLOW ERROR MESSAGE ALLOWED.
	FMPR	3,3		;CALC (S/L)**2.
	JFCL			;NO UNDERFLOW ERROR MESSAGE ALLOWED.
	FADRI	3,201400	;HAVE (1+(S/L)**2) IN AC 3.
	FUNCT	SQRT.,<3>	;CALC. THE SQRT OF
				;(1+(S/L)**2)=1+EPS.
	SKIPN	FLAGWD		;GO TO YGTX IF
	JRST	YGTX		;/Y/>/X/.
XGTY:	FADRI	0,201400	;CALC. 2 + EPS.
	FDVRI	0,202400	;CALC. (2+EPS)/2.
	MOVM	3,0		;PUT IT IN AC NE 0 OR 1 FOR
	FUNCT	SQRT.,<3>		;CALL TO SQRT
	MOVEM	0,3		;SAVE SQRT(1+EPS/2) IN AC 3.
	FUNCT	SQRT.,<2>	;CALC.
	FMPR	0,3		;CALC. SQRT(/X/*(1+EPS/2)).
	JRST	OUT1		;GO TO REST OF CALC.
YGTX:	CAMGE	2,[1.0]		;IF /Y/>1, GO TO POSSIBLE OVERFLOW CASE.
	JRST	YXSMAL		;O'E, GO TO POSSIBLE UNDERFLOW.
	FDVRI	0,203400	;CALC. (1+EPS)/4.
	FMPR	2,0		;CALC. /Y/*(1+EPS)/4.
	MOVM	3,XSAVE		;/X/ TO AC 3.
	FDVRI	3,203400	;CALC. /X//4.
	JFCL			;SUPPRESS UNDERFLOW ERROR MESSAGE.
	FADR	3,2		;CALC.(/X//4)+(/Y/*(1+EPS)/4).
	FUNCT	SQRT.,<3>	;CALC.
	FMPR	0,SQ2		;MULTIPLY IT BY SQRT(2).




OUT1:	MOVM	1,YSAVE		;/Y/ TO AC 1.
	FDVR	1,0		;CALC. /Y//2 TIMES THE
	FDVRI	1,202400	;SQRT TERM.
	JRST	SIGNS		;GO TO REST OF CALC.
YXSMAL:	FMPR	0,2		;/Y/*(1+EPS).
	MOVM	3,XSAVE		;/X/ TO AC 3.
	FADR	3,0		;CALC. /X/+/Y/*(1+EPS).
	FUNCT	SQRT.,<3>	;CALC. SQRT OF 
	MOVEM	0,3		;SAVE IT IN AC 3.
	FDVR	0,SQ2		;DIVIDE IT BY THE SQRT(2).
	FMPR	3,SQ2		;OR MULTIPLY IT BY SQRT(2).
	FDVR	2,3		;THEN CALC /Y//THIS.
	MOVE	1,2		;PUT A TEMP ANSWER IN AC 1.
SIGNS:	SKIPGE	XSAVE		;SET UP THE REAL AND
	EXCH	1,0		;THE IMAGINARY ANSWERS WITH
	SKIPGE	YSAVE		;THE APPROPRIATE
	MOVNS	0,0		;SIGNS.
	MOVE	2,SAVE2		;RESTORE AC 2.
	MOVE	3,SAVE3		;RESTORE AC 3.
	GOODBY	(1)	;RETURN

REAL:	FUNCT	SQRT.,<2>	;CALC. THE SQRT(X)
	SETZ	1,		;IMAG. ANSWER = 0.
	MOVE	2,SAVE2		;RESTORE AC 2.
	GOODBY	(1)	;RETURN

SQ2:	201552023632		;SQRT(2).
XSAVE:	0
YSAVE:	0
SAVE2:	0
SAVE3:	0
FLAGWD:	0
	PRGEND
TITLE	CLOG	%4.(235) COMPLEX LOG FUNCTION
SUBTTL	H. P. WEISS/HPW		11-DEC-73



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
SEARCH	FORPRM
ENTRY	CLOG
EXTERN	CLOG.
CLOG=CLOG.
PRGEND
TITLE	CLOG.	%4.(235) COMPLEX LOG FUNCTION
SUBTTL	D. TODD /DRT 15-FEB-1973



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
;FROM	LIB40 ..(050)
;FROM	V.021	22-SEP-69
;FROM	V.027	LIB40
;COMPLEX LOGARITHM FUNCTION.
;THIS ROUTINE CALCULATES THE LOGARITHM OF A COMPLEX
;ARGUMENT., Z = X + I*Y, WITH THE FOLLOWING ALGORITHM:

;	CLOG(Z) = LOG(ABS(Z)) + I*ATAN2(Y,X)

;WHERE ABS(Z) IS CALCULATED AS L*SQRT(1.0 + (S/L)**2)
;BECAUSE OF OVERFLOW/UNDERFLOW (L IS THE LARGER 
;OF X AND Y, AND S IS THE SMALLER.)


;A SPECIAL CASE CHECK IS MADE FOR X=0 AND Y=0
;AND AN ANSWER WITH REAL = - INFINITY AND IMAG = 0
;IS RETURNED.

;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
;	JSA	Q,CLOG
;	EXP	ARG
;THE REAL PART OF THE ANSWER IS RETURNED IN ACCUMULATOR A
;AND THE IMAGINARY PART IS RETURNED IN ACCUMULATOR B.

	SEARCH	FORPRM

A=0
B=1
C=10
D=11
Q=16
P=17

	HELLO	(CLOG,.)	;[235] ENTRY TO COMPLEX LOG ROUTINE.
	MOVEI	1,@(Q)		;GET ADDRESS OF COMPLEX ARG.
	MOVE	0,(1)		;PICK UP REAL PART OF ARG.
	MOVE	1,1(1)		;PICK UP IMAG. PART OF ARG.
	JUMPN	0,NORMAL	;SPECIAL CASE CHECK
	JUMPN	1,NORMAL	;FOR X=0 AND Y=0.
	ERROR	(APR,5,1,.+1)	;TYPE AN ERROR MESSGAE
	MOVE	A,[400000000001]	;THE ANSWER RETURNED IS
	SETZ	B,		;REAL=-INFINITY, IMAG. = 0.
	GOODBY	(1)	;RETURN

NORMAL:	MOVEM	C,SAVEC		;SAVE AC C.
	MOVEM	D,SAVED		;SAVE AC D.
	MOVE	C,0		;REAL OF ARG. TO AC C.
	MOVE	D,1		;IMAG OF ARG. TO AC D.
	FUNCT	ATAN2.,<D,C>	;CALCULATE IMAG.
	JUMPGE	A,.+2		;ANGLE IS OK IF >= 0,
	FADR	A,TWOPI		;ADD 2*PI IF ANGLE WAS < 0.
	MOVEM	A,IMAG		;STORE ANSWER IN IMAG.
	MOVMS	C,C		;GET ABS(X).
	MOVMS	D,D		;GET ABS(Y).
	CAMGE	D,C		;PUT LARGER OF X AND Y
	EXCH	D,C		;IN AC D AND SMALLER IN AC C.
	FDVR	C,D		;CALC. S/L.
	JFCL			;NO UNDERFLOW ERROR MESSAGE ALLOWED.
	FMPR	C,C		;CALC. (S/L)**2.
	JFCL			;NO UNDERFLOW ERROR MESSAGE ALLOWED.
	FADRI	C,201400	;CALC. 1.0 + (S/L)**2.
	FUNCT	SQRT.,<C>	;CALC. THE SQRT OF THIS
	MOVE	C,A		;THE ANSWER IS IN AC A AND AC C.
	FMPR	C,D		;CALC. L*SQRT.
	JFCL	OVER		;JUMP IF IT OVERFLOWS.
	FUNCT	ALOG.,<C>	;CALC. LOG(ABS(Z)).
OUT:	MOVE	B,IMAG		;SET UP IMAG. PART OF ANSWER.
	MOVE	C,SAVEC		;RESTORE AC C.
	MOVE	D,SAVED		;RESTORE AC D.
	GOODBY	(1)	;RETURN

OVER:	MOVEM	A,SAVEHF	;ARG. TO SAVEHF.
	FUNCT	ALOG.,<SAVEHF>	;CALC. ALOG(SQRT).
	MOVEM	A,SAVEHF	;STORE ANSWER IN SAVEHF.
	FUNCT	ALOG.,<D>	;CALC. ALOG(L).
	FADR	A,SAVEHF	;SUM THE TWO LOGS.
	JRST	OUT		;GO TO RETURN.



TWOPI:	203622077325		;2*PI.
SAVEHF:	0
IMAG:	0
SAVEC:	0
SAVED:	0
	PRGEND
TITLE	CABS	%4.(235) COMPLEX ABSOLUTE VALUE FUNCTION
SUBTTL	H. P. WEISS/HPW		11-DEC-73



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
SEARCH	FORPRM
ENTRY	CABS
EXTERN	CABS.
CABS=CABS.
PRGEND
TITLE	CABS.	%4.(235) COMPLEX ABSOLUTE VALUE FUNCTION
SUBTTL	D. TODD /DRT 15-FEB-1973



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
;FROM	V.21	LIB40

;COMPLEX ABSOLUTE VALUE FUNCTION
;THIS ROUTINE CALCULATES THE ABSOLUTE VALUE OF A COMPLEX NUMBER
;Z = X +I*Y
;THE CALLING SEQUENCE IS AS FOLLOWS:
;	JSA Q.,CABS
;	EXP ARG
;THE ANSWER IS LEFT IN ACCUMULATOR "A"

;CALCULATE SQRT(Y*Y+X*X) BX USING:
;	Y*SQRT[(X/Y)**2+1] WITH X<Y
;IF (X/Y)**2 UNDERFLOWS, THEN RETURN Y


	SEARCH	FORPRM

Q=16
A=0
B=1

MLON	;ALLOW MULTI LINE LITERALS

	HELLO	(CABS,.)
	MOVEI B,@(Q)		;GET ADR OF COMPLEX ARG
	MOVM A,(B)		;GET MAGNITUDE OF REAL PART
	MOVM B,1(B)		;GET MAGNITUDE OF IMAGINARY PART
	CAML A,B		;GET SMALLER ARG IN A AND
	EXCH A,B		; LARGER IN B
	FDVR A,B		;X/Y
	JFCL OUT		;IF UNDERFLOW, GO TO OUT VIA OVTRAP
	FMPR A,A		;(X/Y)**2
	JFCL OUT		;IF UNDERFLOW, GO TO OUT VIA OVTRAP
	FADRI A,201400		;(X/Y)**2+1.0
	MOVEM	A,TT		;STORE ARG. IN NE 0 OR 1.
	MOVEM	B,TTT		;SAVE LARGER.
	FUNCT	SQRT.,<TT>	;TAKE SQUARE ROOT
	FMPR A,TTT		;Y*SQRT[(X/Y)**2+1.0]
RET:	GOODBY	(1)		;RETURN
OUT:	MOVE	A,B		;Y IS THE ANSWER
	GOODBY	(1)	;RETURN

TT:	0
TTT:	0
	PRGEND
TITLE	CEXP	%4.(235) COMPLEX EXPONENTIATION ROUTINE
SUBTTL	H. P. WEISS/HPW		11-DEC-73



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
SEARCH	FORPRM
ENTRY	CEXP
EXTERN	CEXP.
CEXP=CEXP.
PRGEND
TITLE	CEXP.	%4.(235) COMPLEX EXPONENTIATION ROUTINE
SUBTTL	D. TODD/DRT/HPW		30-NOV-1973



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
;FROM	V.021	7-SEP-69
;FROM	V.27	LIB40

;COMPLEX EXPONENTIAL FUNCTION.

;THIS ROUTINE CALCULATES THE EXPONENTIAL OF A
;COMPLEX ARGUMENT Z=X+I*Y IN THE FOLLOWING
;MANNER:

;IF -89.415... <= X <= +88.029....,

;	CEXP(Z)=EXP(X)*[COS(Y)+I*SIN(Y)]

;IF X < -89.415...,

;	CEXP(Z)=0

;AND AN ERROR MESSAGE IS RETURNED FOR BOTH THE
;REAL AND THE IMAGINARY PARTS UNLESS EITHER
;COS(Y) OR SIN(Y) IS IDENTICALLY ZERO.

;IF X>+88.029...,

;  REAL(CEXP(Z))=[EXP((ALOG(/COS(Y)/))+X)]*[SIGN(COS(Y))]

;  IMAG(CEXP(Z))=[EXP((ALOG(/SIN(Y)/))+X)]*[SIGN(SIN(Y))]


;THE CALLING SEQUENCE FOR THE ROUTINE IS AS FOLLOWS:

;	JSA	16,CEXP
;	EXP	ARG

;THE REAL PART OF THE ANSWER IS LEFT IN AC 0.
;THE IMAGINARY PART OF THE ANSWER IS LEFT IN AC 1.


	SEARCH	FORPRM




	HELLO	(CEXP,.)	;[235] ENTRY TO COMPLEX EXP. ROUTINE.
	MOVEI	1,@0(16)	;PUT ADDR OF Z IN AC 1
	MOVE	0,(1)		;X TO AC 0.
	MOVE	1,1(1)		;Y TO AC 1.
	MOVEM	2,SAVE2		;SAVE THE CONTENTS OF AC 2.
	MOVEM	3,SAVE3		;SAVE THE CONTENTS OF AC 3.
	MOVEM	0,2		;X TO AC 2.
	MOVEM	1,3		;X TO AC 3.
	FUNCT	COS.,<3>	;FIND THE
				;COS(Y) AND STORE
	MOVEM	0,COSY		;IT IN COSY.
	FUNCT	SIN.,<3>		;FIND THE
	MOVEM	0,SINY		;IT IN SINY.
	CAMGE	2,MIN89		;IF X < -89.415...,
	JRST	XVRYNG		;THEN GO TO XVRYNG.
	CAMLE	2,PLS88		;IF X >+88.029...,
	JRST	XVRYPL		;THEN GO TO XVRYPL.
	FUNCT	EXP.,<2>		;O'E, FIND E**X
	MOVE	1,SINY		;FORM (E**X)*(SIN(Y)) AND
	FMPR	1,0		;STORE IT IN AC 1.
	FMPR	0,COSY		;STORE (E**X)*(COS(Y)) IN AC 0.
OUT:	MOVE	2,SAVE2		;RESTORE THE CONTENTS OF AC 2.
	MOVE	3,SAVE3		;RESTORE THE CONTENTS OF AC 3.
	GOODBY	(1)	;RETURN


XVRYNG:	SKIPN	COSY		;IF COS(Y)=0, THE
	JRST	.+3		;REAL(ANS.) IS EXACTLY 0.
	ERROR	(APR,7,1,.+1)	;O'E, REAL(ANS.)=0 PLUS
	SKIPN	SINY		;IF SIN(Y)=0, THE
	JRST	.+3		;IMAG(ANS.) IS EXACTLY 0.
	ERROR	(APR,7,1,.+1)	;O'E, IMAG(ANS.)=0 PLUS
	SETZB	0,1		;ANSWER=(0,0)
	JRST	OUT		;GO TO RETURN.

XVRYPL:	SKIPN	SINY		;IF SIN(Y)=0, THEN IMAG(ANS)
	JRST	IZERO		;=0, GO TO IZERO.
	MOVM	3,SINY		;O'E, /SIN(Y)/TO AC NE 0 OR 1.
	FUNCT	ALOG.,<3>	;FIND THE LOG (/SIN(Y)/)
	FADR	0,2		;CALC. (LOG(/SIN(Y)/)) + X
	MOVE	3,0		;AND STORE IT IN AC NE 0 OR 1.
	FUNCT	EXP.,<3>		;CALC. (E**X)*(/SIN(Y)/)
	MOVEM	0,SAVIMG	;IT IN SAVIMG.
	SKIPGE	SINY		;MAKE IT NEGATIVE IF
	MOVNM	0,SAVIMG	;SIN(Y)<0.
REAL:	SKIPN	COSY		;IF COS(Y)=0, THEN REAL(ANS.)
	JRST	RZERO		;=0, SO GO TO REZERO.
	MOVM	3,COSY		;O'E, /COS(Y)/ TO AC NE 0 OR 1.
	FUNCT	ALOG.,<3>		;FIND THE LOG (/COS(Y)/)
	FADRM	0,2		;CALC. (LOG(/COS(Y)/)) + X
	FUNCT	EXP.,<2>		;CALC. (E**X)*(/COS(Y)/)
	SKIPGE	COSY		;MAKE IT NEGATIVE IF
	MOVNS	0		;COS(Y) WAS NEGATIVE.
	MOVE	1,SAVIMG	;RESTORE IMAG(ANS.)
	JRST	OUT		;GO TO RETURN.

IZERO:	SETZM	SAVIMG		;SET LOC OF IMAG(ANS.)=0.
	JRST	REAL		;RETURN TO CALC.

RZERO:	SETZ	0,		;SET REAL(ANS.)=0.
	MOVE	1,SAVIMG	;SET UP IMAG(ANS.).
	JRST	OUT		;GO TO RETURN.

SINY:	0
COSY:	0
SAVIMG:	0
SAVE2:	0
SAVE3:	0
MIN89:	570232254037
PLS88:	207540074636

	PRGEND
TITLE	CSIN	%4.(235) COMPLEX SIN FUNCTION
SUBTTL	H. P. WEISS/HPW		11-DEC-73



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
SEARCH	FORPRM
ENTRY	CSIN
EXTERN	CSIN.
CSIN=CSIN.
PRGEND
TITLE	CCOS	%4.(235) COMPEX COSINE FUNCTION
SUBTTL	H. P. WEISS/HPW		11-DEC-73



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
SEARCH	FORPRM
ENTRY	CCOS
EXTERN	CCOS.
CCOS=CCOS.
PRGEND
TITLE	CSIN.	%4.(235) COMPLEX SIN AND COSINE ROUTINES
SUBTTL	D. TODD/DRT/HPW	30-NOV-1973



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION

;FROM V.21	LIB40

;COMPLEX SINE AND COSINE FUNCTIONS.

;THESE ROUTINES CALCULATE THE SINE AND COSINE
;OF Z., WHERE Z=X+I*Y.  THE ALGORITHM USED IS:
;	SIN(Z)=SIN(X)COSH(Y)+I*COS(X)SINH(Y)
;	COS(Z)=COS(X)COSH(Y)-I*SIN(X)SINH(Y)

;THE ROUTINES ARE CALLED IN THE FOLLOWING MANNER:
;	JSA	16,CSIN
;	EXP	ARG
;WHERE ARG IS THE ADDRESS OF THE REAL PART OF THE
;COMPLEX ARGUMENT, THE IMAGINARY PART BEING IN
;ARG+1.  THE COMPLEX ANSWER IS RETURNED IN AC'S
;0 AND 1.

	SEARCH	FORPRM



	HELLO	(CCOS,.)	;[235] ENTRY TO COMPLEX COSINE ROUTINE.
	MOVEI	1,1		;SET THE FLAGWORD AC = 1.
	JRST	CSIN1		;GO TO SINE ROUTINE.

	HELLO	(CSIN,.)	;[235] ENTRY TO COMPLEX SINE ROUTINE.
	SETZM	1		;SET THE FLAGWORD AC TO 0.

CSIN1:	MOVEM	1,FLGWRD	;SET UP THE FLAGWORD ITSELF.
	MOVEI	1,@0(16)	;[227] PUT ADDR OF Z IN AC 1
	MOVE	0,(1)		;PUT X IN AC 0.
	MOVE	1,1(1)		;PUT Y IN AC 1.
	MOVEM	1,YSAVE		;SAVE Y.
	MOVEM	2,SAVE2		;SAVE THE CONTENTS OF AC 2.
	MOVE	2,0		;PUT X IN AN AC NE 0 OR 1.
	FUNCT	SIN.,<2>		;CALCULATE
				;THE SIN(X).
	MOVEM	0,SINX		;STORE RESULT OF SIN(X)
	FUNCT	COS.,<2>		;CALC. COS(X) AND
	MOVEM	0,COSX		;IN COSX.
	SKIPN	FLGWRD		;IF THIS IS CSIN ROUTINE,
	JRST	.+4		;THEN JUMP AHEAD.
	MOVN	2,SINX		;OTHERWISE, PUT -SIN(X)
	EXCH	2,COSX		;IN COSX, AND COS(X)
	MOVEM	2,SINX		;IN SINX.
	SKIPN	YSAVE		;SKIP UNLESS Z IS REAL, IN WHICH
	JRST	CSIN2		;CASE, GO TO THE SIMPLE CASE.
	MOVM	2,YSAVE		;/Y/ TO AC 2.
	CAMLE	2,EIGHT8	;IF /Y/ TOO LARGE FOR SIMPLE
	JRST	OVFL		;SINH OR COSH, GO TO OVFL.
	MOVE	2,YSAVE		;Y TO AC 2.
	FUNCT	SINH.,<2>	;CALC. SINH(Y) AND
	MOVEM	0,FLGWRD	;IN FLGWRD.
	FUNCT	COSH.,<2>	;CALC. COSH(Y) AND
	FMPR	0,SINX		;CALC. THE REAL PART OF THE ANS.
	MOVE	1,FLGWRD	;SINH(Y) TO AC 1.
	FMPR	1,COSX		;CALC. THE IMAG. PART OF THE ANS.





OUT:	MOVE	2,SAVE2		;RESTORE AC 2.
	GOODBY	(1)	;RETURN

CSIN2:	MOVEI	1,0		;SET THE IMAG. PART OF THE ANSWER = 0.
	MOVE	0,SINX		;SET REAL(ANSWER)=TRIG. FN.
	JRST	OUT		;GO TO RETURN.

OVFL:	MOVM	2,SINX		;/SIN(X)/ TO AC 2.
	JUMPN	2,.+3		;IF ARG FOR ALOG IS
	SETZM	FLGWRD		;ZERO, SET REAL(ANS) = 0,
	JRST	IMAG		;AND JUMP AHEAD.
	FUNCT	ALOG.,<2>	;CALC. THE LOG(/SIN(X)/) AND
	FSBR	0,LOG2BE	;CALC. (LOG(/SIN(X)/))-(LOG(2)).
	MOVM	2,YSAVE		;/Y/ TO AC 2.
	FADR	2,0		;CALC. (LOG(/SIN(X)/))-(LOG(2))+/Y/.
	FUNCT	EXP.,<2>		;CALC. (E**/Y/)*/SIN(X)/)/2.
	MOVEM	0,FLGWRD	;SAVE THE ANSWER IN FLGWRD.
	SKIPGE	SINX		;IF ANS. IS TO BE < 0,SET
	MOVNM	0,FLGWRD	;ANS. < 0.
IMAG:	MOVM	2,COSX		;/COS(X)/ TO AC 2.
	JUMPN	2,.+3		;IF ARG FOR ALOG = 0,
	SETZM	1		;SET IMAG(ANSWER)=0, AND
	JRST	IMAG2		;JUMP AHEAD.
	FUNCT	ALOG.,<2>	;CALC. LOG(/COS(X)/) AND
	FSBR	0,LOG2BE	;CALC. LOG(/COS(X)/))-(LOG(2)).
	MOVM	2,YSAVE		;/Y/ TO AC 2.
	FADR	2,0		;CALC. LOG(/COS(X)/)-(LOG(2))+/Y/.
	FUNCT	EXP.,<2>		;CALC. (/COS(X)/)*
				;(E**/Y/)/2.
	SKIPGE	YSAVE		;SET THE SIGN NEGATIVE IF
	MOVN	0,0		;SINH(Y) WAS NEGATIVE.
	SKIPGE	COSX		;NEGATE THE SIGN IF
	MOVN	0,0		;COS(X) WAS NEGATIVE. PUT
	MOVE	1,0		;THE IMAG. PART OF THE ANS. IN AC 1.
IMAG2:	MOVE	0,FLGWRD	;PUT THE REAL PART OF THE ANS. IN AC 0.
	JRST	OUT		;GO TO RETURN.

EIGHT8:	207540074636
LOG2BE:	200542710300
PIOT:	201622077325
YSAVE:	0
SAVE2:	0
SINX:	0
COSX:	0
FLGWRD:	0
	PRGEND
TITLE	CFM	%4.(207) ;FROM LIB40 VERSION .027
	SUBTTL	D. TODD /DRT/HPW/	29-AUGUST-1973



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
;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:
;	MOVEI	16,ARG2
;	PUSHJ	17,CFM.N
;WHERE N=0,2,4,OR 6. 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:
;	MOVEI	16,ARG2
;	PUSHJ	17,CMM.N
;WHERE N AND ARG1 ARE AS BEFORE AND THE RESULTS ARE
;LEFT IN ARG2 AND ARG2+1.

	SEARCH	FORPRM

	DEFINE SYM(A,B)<
	SIXBIT	/A'B/
A'B:	HRLI	16,B
	JRST	A'0
>
ZZ.=2
REPEAT 6,<
	SYM(CFM.,\ZZ.)
	SYM(CMM.,\ZZ.)
	ZZ.=ZZ.+2
>
	SIXBIT	/CMM.0/
CMM.0:	MOVSS	16		;SET UP THE ADDRESS WORD.
	JRST	CFM.0
	SIXBIT	/CFM.0/
CFM.0:
	PUSH	P,(L)		;SAVE THE LINK
	POP	P,ARG21		;COPY TO ARG21
	PUSH	P,1(L)		;COPY THE ARGUMENTS
	POP	P,ARG22
	MOVSS	L		;GET THE POINTER TO THE FIRST ARG
	PUSH	P,(L)		;SAVE FIRST ARG
	POP	P,ARG11
	PUSH	P,1(L)
	POP	P,ARG12
	PUSH	P,T0		;SAVE FOUR REGISTERS
	PUSH	P,T1
	PUSH	P,T2
	PUSH	P,T3
	PUSHJ	17,CALC		;CALC (A*C-B*D)
	MOVEM	0,REAL		;SAVE REAL (ANS) IN REAL.
	MOVE	0,ARG12		;SET UP
	EXCH	0,ARG11		;AND CALCULATE
	MOVNM	0,ARG12		;(B*C+A*D)
	PUSHJ	17,CALC		;AND LEAVE IT IN AC0.
	POP	P,T3		;RESTORE THE SCRATCH REGISTERS
	POP	P,T2
	POP	P,T1
	MOVEM	T0,1(L)		;STORE THE IMAG ANS
	POP	P,T0
	PUSH	P,REAL		;PUT REAL ANS ON THE STACK
	POP	P,(L)		;STORE THE REAL PART
	POPJ	P,		;RETURN


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

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

OUFLO1:	MOVE	3,ARG22		;"D" TO AC3.
	SKIPN	0		;JUMP IF OVERFLOW.
	JRST	UNDER1		;UNDERFLOW, GO TO UNDER1
	FMPRB	2,3		;CALC B*D AND
	JFCL	OFL1OU		;IF OVER/UNDERFLOW GO TO OFL1OU.
	XOR	3,0		;OVERFLOW + NORMAL CASE
	TLNN	3,400000	;IF THE SIGNS ARE THE SAME, GO
	JRST	SROVNM		;TO SROVNM.O'E, GO TO INFIN.
INFIN:	MOVEM	0,ANSTMP	;STORE ANSWER.
	ERROR	(APR,5,1,.+1)	;RETURN WITH ERROR MESSAGE
	MOVE	0,ANSTMP	;PICK UP ANSWER.
	POPJ	17,		;RETURN.
OFL1OU:	JUMPE	2,INFIN		;OVERFLOW + UNDERFLOW CASE GOES TO INFIN.
	XOR	3,0		;OVERFLOW + OVERFLOW CASE.
	TLNN	3,400000	;IF THE SIGNS ARE THE
	JRST	SROVOV		;SAME, GO TO SROVOV,
	JRST	INFIN		;O'E, GO TO INFIN.
SROVNM:	JUMPE	2,INFIN		;OVERFLOW + ZERO, GO TO INFIN.
	MOVE	0,ARG21		;C TO AC0.
	FDVRI	0,202400	;CALC. C/2.
	FMPR	0,ARG11		;CALC. (C/2)*A AND IF
	JFCL	1,[POPJ 17,]	;IT OVERFLOWS, IMMEDIATELY RETURN.
	FSBRM	0,2		;CALC ((C/2*A)-(B*D).
	FADR	0,2		;CALC (A*C-B*D) AND
	POPJ	17,		;RETURN.
SROVOV:	MOVE	0,ARG21		;C TO AC0.
	FDVR	0,ARG22		;CALC. (C/D).
	MOVE	1,ARG12		;B TO AC1.
	FDVR	1,ARG11		;CALC. (B/A).
	FSBR	0,1		;CALC. ((C/D)-(B/A)) AND
	JFCL	FIXUP		;IF IT UNDERFLOWS, GO TO FIXUP.
	FMPR	0,ARG22		;CALC. ((C)-(B/A)*D)) AND
	JFCL			;SUPPRESS OVERFLOW MESSAGE.
	FMPR	0,ARG11		;CALC (A*C-B*D) AND
	POPJ	17,		;RETURN.
FIXUP:	FMPR	1,ARG22		;CALC. ((B/A)*D).
	MOVE	0,ARG21		;C TO AC0.
	FSBR	0,1		;CALC. (C-(B/A)*D).
	FMPR	0,ARG11		;CALC. (A*C-B*D) AND
	POPJ	17,		;RETURN.

UNDER1:	FMPR	2,3		;CALC. B*D AND GO
	JFCL	U1OU		;TO U1OU IF OVER/UNDERFLOW.
	MOVM	3,2		;IF B*D IS <2**-105,
	CAMGE	3,[030400000000] ;GO
	JRST	.+3		;TO .+3.
	MOVN	0,2		;NO BITS CAN BE SAVED,
	POPJ	17,		;FIX SIGN AND RETURN.
	FMPRI	2,232400	;CALC B*D*(2**27).
UNDR1:	MOVE	0,ARG11		;CALC
	FMPRI	0,232400	;A*C*(2**27)
	FMPR	0,ARG21		;AND SUPPRESS AN
	JFCL			;ERROR MESSAGE.
	FSBR	0,2		;CALC (2**27)(A*C-B*D),
	FDVRI	0,232400	;CALC (A*C-B*D) AND
	POPJ	17,		;RETURN.

U1OU:	JUMPE	2,.+3		;IF OVERFLOW + UNDERFLOW, GO
	MOVNM	2,ANSTMP	;TO ERROR
	JRST	INFIN+1		;RETURN.
	MOVE	2,ARG12		;O'E, TRY TO SAVE BITS.
	FMPRI	2,232400	;CALC. B*(2**27).
	FMPR	2,ARG22		;CALC B*D*(2**27)
	JFCL			;AND SUPPRESS UNDERFLOW MESSAGE
	MOVE	0,ARG11		;CALC. A*(2**27) AND
	FMPRI	0,232400	;A*(2**27)*C
	FMPR	0,ARG21		;AND SUPPRESS
	JFCL			;UNDERFLOW MESSAGE.
	JUMPN	2,UNDR1+4	;IF BOTH UNDERFLOW, RETURN
	JUMPN	0,UNDR1+5	;AN ERROR MESSAGE. O'E,
	ERROR	(APR,7,1,.+1)	;ERROR MESSGE
	SETZ	0,		;CLEAR AC0
	POPJ	17,		;CALC.

OUFLO2:	JUMPN	2,.+11		;JUMP AHEAD IF (B*D) OVERFLOWED.
	CAML	0,[030400000000] ;IF NO BITS CAN BE SAVED,
	POPJ	17,		;RETURN.
	FMPRI	0,232400	;CALC. (2**27)*(A*C)
	MOVE	2,ARG12		;AND
	FMPRI	2,232400	;THEN
	FMPR	2,ARG22		;(2**27)*(B*D)
	JFCL			;AND GO TO THE REST
	JRST	UNDR1+4		;OF THE CALC.
	MOVEM	2,1		;OVERFLOW + NORMAL.
	XOR	1,0		;IF THE SIGNS ARE THE SAME,
	TLNN	1,400000	;GO TO SROVN; O'E
	JRST	SROVN		;SET UP
	MOVNM	2,ANSTMP	;THE ANSWER AND
	JRST	INFIN+1		;GO TO ERROR EXIT.
SROVN:	MOVE	3,ARG22		;"D" TO AC3.
	FDVRI	3,202400	;CALC (D/2).
	FMPR	3,ARG12		;CALC (B*(D/2)) AND IF
	JFCL	FIXUP2		;IT OVERFLOWS, GO TO FIXUP2.
	FSBR	0,3		;CALC. ((A*C)-(B*(D/2))).
	FSBR	0,3		;CALC (A*C-B*D) AND
	POPJ	17,		;RETURN.
FIXUP2:	MOVNM	3,ANSTMP	;SET UP THE ANSWER
	JRST	INFIN+1		;AND GO TO ERROR EXIT.

ARG11:	0
ARG12:	0
ARG21:	0
ARG22:	0
REAL:	0
ANSTMP:	0
	ENTRY	CFM.0,CFM.2,CFM.4,CFM.6,CFM.10,CFM.12,CFM.14
	ENTRY	CMM.0,CMM.2,CMM.4,CMM.6,CMM.10,CMM.12,CMM.14
IFN F40LIB,<
	;DEFINE THE OLD F40 NAMES
	ENTRY	CFMM.0,CFMM.2,CFMM.4,CFMM.6
CFMM.0=:CMM.0
CFMM.2=:CMM.2
CFMM.4=:CMM.4
CFMM.6=:CMM.6
>

	PRGEND
TITLE	CFDV	%4.(207) ;FROM LIB40 VERSION .027
	SUBTTL	D. TODD /DRT/HPW	29-AUGUST-1973



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
;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
;	MOVEI	16, ARG2
;	PUSHJ	17, CFD.N
;WHERE N=0,2,4, OR 6. 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:
;	MOVEI	16,ARG2
;	PUSHJ	17,CFDM.N
;WHERE N=0,2,4, OR 6. ARG1 IS ASSUMED TO BE IN AC N AND (N+1),


;AND THE RESULTS ARE RETURNED IN ARG2 AND ARG2+1.

	SEARCH	FORPRM

	DEFINE SYM(A,B)<
	SIXBIT	/A'B/
A'B:	HRLI	16,B
	JRST	A'0
>
ZZ.=2
REPEAT 6,<
	SYM(CFD.,\ZZ.)
	SYM(CDM.,\ZZ.)
	ZZ.=ZZ.+2
>
	SIXBIT	/CDM.0/
CDM.0:	SETOM	MEMFLG		;MEMORY FLAG.
	JRST	CFD.0+1		;GO TO THE MAIN ROUTINE.

	SIXBIT	/CFD.0/
CFD.0:	SETZM	MEMFLG		;SET UP THE MEMORY FLAG
	MOVSS	16		;SET UP 16.
	PUSH	P,T0		;SAVE TWO SCRATCH REGISTERS
	PUSH	P,T1
	MOVEM	0,LCD		;SAVE AC 0 IN CASE IT IS AN ARG.
	MOVE	0,(16)		;PICK UP REAL OF ONE ARG AND
	MOVEM	0,SAB		;STORE IT IN SAB.
	MOVE	0,1(16)		;PICK UP IMAG OF ONE ARG AND
	MOVEM	0,LAB		;STORE IT IN LAB.
	MOVE	0,LCD		;RESTORE AC0.
	MOVSS	16		;SET UP ADDRESS WORD.
	MOVE	0,(16)		;PICK UP REAL OF THE OTHER ARG
	MOVEM	0,SCD		;AND STORE IT IN SCD
	MOVE	1,1(16)		;PICK UP IMAG. OF THIS ARG AND
	MOVEM	1,LCD		;STORE IT IN LCD.
	JUMPE	0,CZERO		;GO TO CZERO IF C=0, TO
	JUMPE	1,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:	SETZM	REAL		;A=
	SETZM	0		;B=
	JRST	OUT+1		;0

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


NTSMPL:	MOVEM	2,SAVE2		;SAVE THE CONTENTS OF AC2.
	MOVEM	3,SAVE3		;SAVE THE CONTENTS OF AC3.
	MOVEM	4,SAVE4		;SAVE THE CONTENTS OF AC4.
	MOVE	2,SAB		;REARRANGE THE REAL
	MOVM	0,SAB		;AND IMAG PARTS OF
	MOVM	1,LAB		;THE 1ST ARG
	MOVEI	4,1		;SO THAT THE SMALLER ONE
	CAMG	0,1		;(IN MAGNITUDE) IS IN
	JRST	.+4		;SAB, AND THE OTHER IS IN
	EXCH	2,LAB		;LAB AND SET UP AC4 AS
	MOVEM	2,SAB		;A FLAG WORD TO TELL
	MOVNS	4,4		;WHICH PART IS WHERE.
	MOVE	2,SCD		;REARRANGE THE REAL
	MOVM	0,SCD		;AND IMAG PARTS OF THE OTHER ARG
	MOVM	1,LCD		;SO THAT THE SMALLER ONE
	CAMG	0,1		;(IN MAGNITUDE) IS IN SCD AND
	JRST	.+4		;THE OTHER IS IN LCD AND
	EXCH	2,LCD		;FIX AC4 TO TELL WHICH
	MOVEM	2,SCD		;PART
	IMULI	4,3		;IS WHERE.
	MOVE	3,LCD		;CALCULATE
	FDVRB	2,3		;THE
	JFCL			;TERM
	FMPR	2,2		;1+(SCD/LCD)**2
	JFCL			;AND
	FADRI	2,201400	;STORE IT IN
	MOVEM	2,DENOM		;DENOM.
	MOVE	0,SAB		;CALCULATE
	FDVR	0,LAB		;(SCD/LCD)*(SAB/LAB)
	JFCL			;SUPPRESSING
	FMPRM	0,3		;UNDERFLOW
	JFCL			;IF NECESSARY.
	CAIN	4,1		;FIX THE AC FLAG WORD
	MOVNI	4,2		;FOR
	ADDI	4,1		;EASY COMPARISONS.
	SKIPL	4		;CALCULATE
	MOVNS	3,3		;+-1 +-(SCD/LCD)*(SAB/LAB),
	FADRI	3,201400	;WHERE THE SIGNS
	SKIPN	4		;DEPEND ON
	MOVNS	3,3		;THE AC FLAG WORD.
	PUSHJ	17,CALC34	;JUMP TO CALC OF (LAB/LCD)*(AC3/DENOM).
	MOVEM	0,REAL		;STORE IT IN REAL(ANS) AND
	MOVEM	0,IMAG		;IMAG(ANS)(THE TEMP. LOCATIONS).
	MOVE	3,SAB		;CALCULATE
	MOVE	2,SCD		;+-(SAB/LAB)+-(SCD/LCD)
	CAMN	4,[-2]		;WHERE THE SIGNS
	MOVNS	2,2		;AGAIN DEPENDS ON
	CAMN	4,[-1]		;THE AC FLAG WORD,
	MOVNS	3,3		;AND IF UNDERFLOW
	MOVE	0,2		;OCCURS JUMP TO
	MOVE	1,3		;THE
	FDVR	3,LAB		;SPECIAL
	JFCL	OVER1		;ROUTINES-
	FDVR	2,LCD		;OVER1,
	JFCL	OVER2		;OVER2,
ADD2:	FADR	3,2		;OR
	JFCL	OVER3		;OVER3.
	PUSHJ	17,CALC34	;JUMP TO CALC OF(LAB/LCD)*(AC3/DENOM).
	JUMPGE	4,.+3		;STORE THIS IN
	MOVEM	0,IMAG		;THE CORRECT
	JRST	.+2		;PART OF THE
	MOVEM	0,REAL		;ANSWER (TEMP. LOCATION).
	MOVE	2,SAVE2		;RESTORE THE CONTENTS OF AC2.
	MOVE	3,SAVE3		;RESTORE THE CONTENTS OF AC3.
	MOVE	4,SAVE4		;RESTORE THE CONTENTS OF AC4.
OUT:	MOVE	0,IMAG		;PICK UP THE IMAG (ANS)
	SKIPN	MEMFLG		;SET UP AC L TO
	MOVSS	L		;CONTAIN THE ANSWER ADDRESS.
	POP	P,T1		;RESTORE T1
	MOVEM	0,1(L)		;STORE THE IMAG (ANS) IN ITS FINAL LOCATION, AND
	POP	P,T0		;RESTORE T0
	PUSH	P,REAL		;GET THE REAL PART
	POP	P,(L)		;STORE THE REAL PAR
	POPJ	17,		;EXIT.

CALC34:	MOVM	1,LAB		;CALC34 CALCS. (LAB/LCD)*(AC3/DENOM).
	MOVM	2,LCD		;/LAB/ TO AC 1 AND /LCD/ TO AC 2.
	MOVM	0,3		;/AC3/ TO AC 0.
	CAMGE	2,ONE		;GO TO CASEA IF
	JRST	CASEA		;/LCD/<1.0. O'E, STAY HERE.
	CAMGE	1,ONE		;GO TO CASEB IF /LCD/>1.0 AND
	JRST	CASEB		;/LAB/<1.0 OR IF
	CAMGE	0,ONE		;(>1)(<1)/(>1)(>1).
	JRST	CASEB		;STAY HERE IF (>1)(>1)/
	FDVR	2,1		;(>1)(>1).
	FDVR	0,DENOM		;CALCULATE
	FDVR	0,2		;IT AND GO
	JRST	SETSGN		;TO SETSGN.
CASEB:	FMPR	0,1		;CALCULATE CASE B AND
	JRST	.-4		;GO TO SETSGN (EVENTUALLY).
CASEA:	FMPR	2,DENOM		;CONDENSE THE DENOMINATOR INTO AC 2.
	CAMLE	1,ONE		;THEN (<1)(<1)/(<1) GOES
	JRST	CHKAGN		;TO SR, AND (>1)(><1)/(<>1)
	CAMLE	0,ONE		;GOES TO CHKAGN,
	JRST	NOTRVS		;AND EVERYTHING ELSE
	CAMG	2,ONE		;GOES TO
	JRST	SR		;NOTRVS.
NOTRVS:	FMPRM	1,0		;CALCULATE
	JFCL	1,SETSGN	;NOTRVS AND GO
	FDVR	0,2		;TO
	JRST	SETSGN		;SETSGN.
CHKAGN:	CAMG	0,ONE		;(>1)(<1)/(<>1)
	JRST	NOTRVS		;AND (>1)(>1)/(<1)
	CAMG	2,ONE		;GO TO
	JRST	NOTRVS		;NOTRVS.
	FDVR	1,2		;(>1)(>1)/(>1) IS
	FMPRM	1,0		;CALCULATED AND
	JRST	SETSGN		;GOES TO SETSGN.
SR:	MOVEM	1,TEMP		;SR CALCULATES
	FDVR	1,2		;(<1)(<1)/(<1)
	JFCL	OV1		;AND SINCE
	FMPRM	1,0		;(<1)/(<1)
	JRST	SETSGN		;CAN
OV1:	MOVE	1,TEMP		;OVERFLOW, IT
	MOVEM	0,TEMP		;CALCULATES
	FDVR	0,2		;IT
	JFCL	OV2		;WHICHEVER
	FMPRM	1,0		;WAY
	JRST	SETSGN		;IS
OV2:	MOVE	0,TEMP		;NECESSARY
	FMPRM	1,0		;AND
	FDVR	0,2		;THEN GOES TO SETSGN.
SETSGN:	MOVE	1,LAB		;GET THE
	XOR	1,LCD		;SIGN OF THE
	XOR	1,3		;RESULT IN AC 1.
	SKIPG	1		;SET THE SIGN OF
	MOVNS	0,0		;THE ANSWER
	POPJ	17,		;AND EXIT.

OVER1:	FDVR	2,LCD		;IF ONE
	JFCL	UU		;TERM
	MOVE	3,2		;UNDERFLOWS
	MOVE	0,1		;AND THE OTHER
	MOVE	2,LAB		;TERM IS SO LARGE
	JRST	OVER2+1		;THAT NO BITS
OVER2:	MOVE	2,LCD		;CAN BE
	MOVEM	2,SAVE5		;SAVED,
	JUMPE	3,UU		;THEN
	MOVM	1,3		;RETURN
	CAML	1,[030400000000] ;RIGHT
	JRST	ADD2+2		;AWAY.
	FMPRI	3,270400	;O'E, TRY TO
	FMPRI	0,270400	;SAVE SOME BITS
	FDVRM	0,2		;BY MULTIPLYING
	JFCL			;THE TERMS BY 2**56,
	FADRB	3,2		;ADDING THE TERMS, AND THEN / BY 2**56.
	FDVRI	3,270400	;IF THE RESULT UNDERFLOWS, GO
	JFCL	SRX		;TO SRX, O'E GO BACK
	JRST	ADD2+2		;TO THE MAIN ROUTINE.
OVER3:	MOVE	3,1		;SET UP
	FDVR	3,LAB		;AC 2 FOR
	FMPRI	3,270400	;SRX. AC 2 WILL
	FMPRI	2,270400	;CONTAIN THE TERM
	FADR	2,3		;*(2**56).
SRX:	MOVEM	5,SAVE5		;SAVE THE CONTENTS OF AC 5.
	MOVE	0,LCD		;THIS IS AN
	MOVE	1,LAB		;ALTERNATE
	MOVM	5,1		;CALCULATION TO
	CAMGE	5,ONE		;CALC34, AND TAKES
	JRST	SRX2		;ADVANTAGE OF
	FMPR	1,2		;THE FACT THAT HERE
	MOVM	5,0		;DENOM CONTAINS 1.0.
	CAML	5,ONE		;THE ORDER OF
	JRST	.+4		;CALCULATION
	FMPRI	0,270400	;DEPENDS
	FDVRM	1,0		;ON THE SIZE OF
	JRST	SETSN2		;THE TERMS. AFTER
	FDVRM	1,0		;THE CALCULATION A
	FDVRI	0,270400	;TRANSFER
	JRST	SETSN2		;IS
SRX2:	FDVRM	2,0		;MADE
	FMPR	0,1		;TO
	FDVRI	0,270400	;SETSN2
SETSN2:	MOVE	5,SAVE5		;RESTORE THE CONTENTS OF AC 5.
	JRST	ADD2+3		;GO BACK TO MAIN ROUTINE.


UU:	MOVEM	0,SAB		;ANOTHER ALTERNATE
	MOVEM	1,SCD		;CALCULATION TO
	FMPR	1,LCD		;CALC34
	FMPR	0,LAB		;USED WHEN
	FADR	0,1		;S/L FOR
	JFCL	UND		;BOTH SETS
	FDVR	0,LCD		;HAS UNDERFLOWED
	FDVR	0,LCD		;OR FOR
	JRST	ADD2+3		;UNDERFLOW PLUS
UND:	ERROR	(APR,7,1,.+1)
	JRST	ADD2+3		;TO ADD2+3.


SAVE2:	0
SAVE3:	0
SAVE4:	0
SAVE5:	0
SAB:	0
LAB:	0
SCD:	0
LCD:	0
TEMP:	0
REAL:	0
IMAG:	0
DENOM:	0
MEMFLG:	0
ONE:	201400000000

	ENTRY	CFD.0,CFD.2,CFD.4,CFD.6,CFD.10,CFD.12,CFD.14
	ENTRY	CDM.0,CDM.2,CDM.4,CDM.6,CDM.10,CDM.12,CDM.14

IFN F40LIB,<
	;EQUATE THE OLD F40 NAMES
	ENTRY CFDM.0,CFDM.2,CFDM.4,CFDM.6
CFDM.0=:CDM.0
CFDM.2=:CDM.2
CFDM.4=:CDM.4
CFDM.6=:CDM.6

>
	PRGEND
TITLE	REAL	%4.(235) COMPLEX TO REAL CONVERSION
SUBTTL	H. P. WEISS/HPW		11-DEC-73



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
SEARCH	FORPRM
ENTRY	REAL
EXTERN	REAL.
REAL=REAL.
PRGEND
TITLE	REAL.	%4.(120) ;FROM LIB40 VERSION .022
SUBTTL	D. TODD /DRT 15-FEB-1973



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION

;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:
;	JSA	Q., REAL
;	EXP	ARG
;THE ANSWER IS RETURNED IN ACCUMULATOR A

	SEARCH	FORPRM

	A=	0
	Q=	16

	HELLO	(REAL,.)	;[235] ENTRY TO REAL ROUTINE
	MOVE	A, @(Q)		;GET REAL PART OF ARGUMENT
	GOODBY	(1)	;RETURN

	PRGEND
TITLE	AIMAG	%4.(235) COMPLEX TO REAL CONVERSION
SUBTTL	H. P. WEISS/HPW		11-DEC-73



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
SEARCH	FORPRM
ENTRY	AIMAG
EXTERN	AIMAG.
AIMAG=AIMAG.
PRGEND
TITLE	AIMAG.	%4.(235) COMPLEX TO REAL COVERSION
SUBTTL	D. TODD /DRT 25-MAY-1973



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION

;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:
;	JSA	Q., AIMAG
;	EXP	ARG
;THE ANSWER IS RETURNED IN ACCUMULATOR A

	SEARCH	FORPRM
	A=	0
	B=	1
	Q=	16

	HELLO	(AIMAG,.)	;[235]ENTRY TO AIMAG ROUTINE
	MOVEI	B,@(Q)		;PUT IMAG(ARG)
	MOVE	A,1(B)		;IN AC A.
	GOODBY	(1)	;RETURN

	PRGEND
TITLE	CONJG	%4.(235) COMPLEX CONJUGATE FUNCTION
SUBTTL	H. P. WEISS/HPW		11-DEC-73



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
SEARCH	FORPRM
ENTRY	CONJG
EXTERN	CONJG.
CONJG=CONJG.
PRGEND
TITLE	CONJG.	%4.(235) COMPLEX CONJUGATE FUNCTION
SUBTTL	D. TODD /DRT 15-FEB-1973



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION

;FROM	LIB40	V.022

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

;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
;	JSA	Q., CONJG
;	EXP	ARG
;THE REAL PART OF THE ANSWER IS LEFT IN ACCUMULATOR A, AND THE
;IMAGINARY PART IN ACCUMULATOR B.

	SEARCH	FORPRM

	A=	0
	B=	1
	Q=	16

	HELLO	(CONJG,.)	;[235] ENTRY TO CONJG ROUTINE
	MOVEI	B, @(Q)		;GET ADDRESS OF COMPLEX ARGUMENT
	MOVE	A, (B)		;GET REAL PART OF ARGUMENT
	MOVN	B, 1(B)		;GET NEGATIVE OF IMAGINARY PART
	GOODBY	(1)	;RETURN

	PRGEND
TITLE	CMPLX	%4.(235) REAL TO COMPLEX CONVERSION
SUBTTL	H. P. WEISS/HPW		11-DEC-73



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION
SEARCH	FORPRM
ENTRY	CMPLX
EXTERN	CMPLX.
CMPLX=CMPLX.
PRGEND
TITLE	CMPLX.	%4.(235) REAL TO COMPLEX CONVERSION
SUBTTL	D. TODD /DRT 15-FEB-1973



;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.

;COPYRIGHT (C) 1972,1977 BY DIGITAL EQUIPMENT CORPORATION

;FROM	LIB40	V.022
;FROM V.005
;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:
;	JSA	Q.,CMPLX
;	EXP	ARG1
;	EXP	ARG2
;THE REAL PART OF THE ANSWER IS LEFT IN ACCUMULATOR A, AND THE
;IMAGINARY PART IS LEFT IN ACCUMULATOR B

	SEARCH	FORPRM

	A=0
	B=1
	Q=16

	HELLO	(CMPLX,.)	;[235] ENTRY TO CMPLX ROUTINE
	MOVE	A,@(Q)		;GET REAL PART OF COMPLEX ANSWER
	MOVE	B,@1(Q)		;GET IMAGINARY PART OF COMPLEX ANS
	GOODBY	(2)	;RETURN

	END