Trailing-Edge
-
PDP-10 Archives
-
BB-D480F-BB_1985_short
-
mthcpx.mac
There are 9 other files named mthcpx.mac in the archive. Click here to see a list.
SEARCH MTHPRM
TV MTHCPX COMPLEX ROUTINES ,7(4002)
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SUBTTL REVISION HISTORY
COMMENT \
***** Begin Revision History *****
1351 EDS 16-Mar-81 Q10-04786
Fix TWOSEG and RELOC problems.
1405 DAW 6-Apr-81
Support extended addressing.
1464 DAW 12-May-81
Error messages.
1673 CDM 9-Sept-81
Added complex routines.
(CMPL.I, CMPL.C, CMPL.D)
***** Begin Mathlib *****
3024 CDM 17-Nov-81
Changed REAL. to REAL.C, defining REAL. to have integer
arguments as in ANSI-77 standard.
3200 JLC 10-May-82
Mathlib integration. Change all LERRs to $LCALLs. Change
TWOSEG/RELOC to SEGMENT macros.
3201 RVM 20-May-82
Add the routine CMPL.G (convert two gfloating numbers to
complex).
3217 JLC 8-Oct-82
Modify SIXBIT names of the complex multiply and divide routines
so that TRACE won't look for an argument block.
3220 PLB 12-Oct-82
Remove IFIWs from invocations of FUNCT macro, so that we
can put IFIW into the defn of FUNCT.
***** Begin Version 1A *****
3234 CKS 23-Mar-83
The maximum value of the arg to SIN was given erroneously
as 36394.429 instead of the correct value of 823549.66.
Fix this in the code and any of the comments which it
appears.
3247 10-33795 MRB 24-May-83
Fixed bug in the interface of CEXP2 with CEXP3, and to
fix bug in block 2 (of CEXP3) at XISO.
3251 10-34015 MRB 30-Aug-83
Added a special case in CEXP2. When the imaginary part of
the complex number is zero and the real part in an integer
call EXP2. instead. (This is edit 3344 to FOROT7)
***** End Revision History *****
***** Begin Version 2 *****
4002 JLC 8-Aug-83
Install the beginnings of user-fixable library routine
error results. Fix an overflow msg.
\
PRGEND
TITLE CABS COMPLEX ABSOLUTE VALUE FUNCTION
; (SINGLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CABS
EXTERN CABS.
CABS=CABS.
PRGEND
TITLE CABS. COMPLEX ABSOLUTE VALUE FUNCTION
; (SINGLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;CABS(Z) IS CALCULATED AS FOLLOWS
; LET Z = X + I*Y
; V = MAX(ABS(X),ABS(Y))
; W = MIN(ABS(X),ABS(Y))
; THEN
; CABS = V*SQRT(1.0 + (W/V)**2)
;THE RANGE OF DEFINITION FOR CABS IS THE REPRESENTABLE REAL NUMBERS.
; AN OVERFLOW CAN OCCUR, IN WHICH CASE CABS = + MACHINE INFINITY.
;REQUIRED (CALLED) ROUTINES: SQRT
;REGISTER T2 WAS SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; XMOVEI L,ARG
; PUSHJ P,CABS
;THE ANSWER IS RETURNED IN ACCUMULATOR T0
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (CABS,.) ;ENTRY TO MODULUS ROUTINE
PUSH P,T2
XMOVEI T2,@(L)
MOVM T0,(T2) ;ABS(X)
MOVM T2,1(T2) ;ABS(Y)
; OBTAIN MAX(ABS(X),ABS(Y))IN T2
; AND MIN IN T0
CAML T0,T2
EXCH T2,T0 ;THEN CALCULATE CABS
JUMPE T2,RET ;Z = 0, HENCE CABS = O
FDVR T0,T2 ;U/V
JFCL ANS ;NO OVERFLOW, RATIO NEGLIGIBLE IF UNDERFLOW
CAMG T0,TWOM14 ;IF RATIO .LE. 2**-14,
JRST ANS ;RATIO IS NEGLIGIBLE.
FMPR T0,T0 ;**2
FADRI T0,201400 ;+1.0
MOVEM T0,TEMP
FUNCT SQRT.,<TEMP> ;[3220] SQUARE ROOT IS IN AC T0
FMPR T0,T2 ;*V
JFCL OVFL ;NO UNDERFLOW, GET MESSAGE ON OVERFLOW
RET: POP P,T2
GOODBY (1) ;RETURN
OVFL:
$SSPR ;SAVE FIXED-UP RESULT
$LCALL ROV ;RESULT OVERFLOW
$RSPR ;GET FIXED-UP RESULT
JRST RET ;RETURN
ANS: MOVE T0,T2 ;ANSWER = ABS(LARGER) TO T0
POP P,T2 ;RESTORE T2
GOODBY (1) ;RETURN
TWOM14: 163400000000 ;2**-14
SEGMENT DATA
TEMP: 0 ;TEMPORARY STORAGE USED FOR SQRT CALL
PRGEND
TITLE CEXP COMPLEX EXPONENTIAL FUNCTION
; (SINGLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CEXP
EXTERN CEXP.
CEXP=CEXP.
PRGEND
TITLE CEXP. COMPLEX EXPONENTIAL FUNCTION
; (SINGLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;CEXP(Z), WHERE Z = X + I*Y, IS CALCULATED AS FOLLOWS
; CEXP(Z) = EXP(X)*(COS(Y)+I*SIN(Y))
;THE RANGE OF DEFINITION FOR CEXP IS AS FOLLOWS
;FOR Z = X+I*Y, IF ABS(Y) .GT. 823549.66 THE RESULT IS SET TO
; (0.0,0.0) AND AN ERROR MESSAGE IS RETURNED.
;FOR X.LT.-89.415986 THE RESULT IS SET TO (0.0,0.0) AND AN ERROR
; MESSAGE IS RETURNED.
;FOR X .GT. 88.029692;
; IF Y = 0.0, THE RESULT IS SET TO (+INFINITY,0.0) AND AN
; ERROR MESSAGE IS RETURNED.
; IF X/2. IS .GT. 88.029692, THE ABS(REAL(RESULT)) IS SET TO
; +INFINITY AND ABS(AIMAG(RESULT)) IS SET TO +INFINITY AND
; AN ERROR MESSAGE IS RETURNED.
;REQUIRED (CALLED) ROUTINES: EXP,SIN,COS
;REGISTERS T2, T3, AND T4 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; MOVEI L,ARG
; PUSHJ P,CEXP
;THE REAL PART OF THE SOLUTION IS RETURNED IN ACCUMULATOR T0
;THE IMAG PART OF THE SOLUTION IS RETURNED IN ACCUMULATOR T1
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (CEXP,.) ;ENTRY TO NAME ROUTINE
PUSH P,T2 ;SAVE REGISTER T2
PUSH P,T3 ;SAVE REGISTER T3
PUSH P,T4 ;SAVE REGISTER T4
XMOVEI T1,@(L) ;GET ADDRESS OF Z
MOVE T0,(T1) ;X=REAL(Z)
MOVE T1,1(T1) ;Y=IMAG(Z)
MOVEM T1,YSAVE ;SAVE Y
MOVM T1,T1 ;T1=ABS(T1)
CAMG T1,YMAX ;IF ABS(Y) IS NOT TOO LARGE
JRST YOK ;GO TO YOK
$SCPX ;SAVE FIXED-UP RESULT
$LCALL AIZ
;LERR (LIB,%,<CEXP: AIMAG(arg) too large in absolute value; result=(0,0)>)
$RCPX ;GET FIXED-UP RESULT
SETZ T0, ;REAL PART OF SOLUTION IS ZERO
SETZ T1, ;IMAG PART OF SOLUTION IS ZERO
JRST RET ;RETURN
YOK: MOVEM T0,XSAVE ;SAVE X
CAML T0,NEGCON ;IF X IS NOT TOO SMALL
JRST XOK ;GO TO XOK
SETZ T0, ;REAL PART OF SOLUTION IS ZERO
SETZ T1, ;IMAG PART OF SOLUTION IS ZERO
JRST YCHK ;CHECK Y
XOK: FUNCT COS.,<YSAVE> ;[3220] COSY=COS(Y)
MOVEM T0,COSY
FUNCT SIN.,<YSAVE> ;[3220] SINY=SIN(Y)
MOVEM T0,SINY
MOVE T2,XSAVE ;T2=X
MOVE T1,YSAVE ;T1=Y
CAMG T2,TBIG ;IF X IS NOT TOO LARGE
JRST ALG ;GO TO ALG
CAIN T1,0 ;ELSE, IF Y=0
JRST ERR0 ;GO TO ERR
FSC T2,-1 ;ELSE, S=X/2
CAMLE T2,TBIG ;IF S IS TOO LARGE
JRST STIM ;GO TO STIM
MOVEM T2,XSAVE
FUNCT EXP.,<XSAVE> ;T=EXP(S)
MOVE T1,T0
FMPR T1,SINY ;V=T*SINY
MOVM T1,T1 ;V=ABS(V)
HRLOI T4,377777 ;T4=XMAX
FDVR T4,T0 ;Q=XMAX/T
CAIG T1,T4 ;IF V IS LESS THAN OR EQUAL TO Q
JRST GETI ;GO TO GETI
STIM: HRLOI T1,377777 ;ELSE, SET IMAG SOLUTION TO XMAX
$LCALL IPO
;LERR (LIB,%,<CEXP: imaginary part overflow>)
JRST D2
GETI: FMPR T1,T0 ;IRES = V*T
D2: MOVE T4,SINY
CAIGE T4, ;IF SINY IS LESS THAN 0.0
MOVN T1,T1 ;THEN NEGATE IRES
CAMLE T2,TBIG ;IF S IS TOO LARGE
JRST ERR0 ;GO TO ERR
MOVE T2,T0
FMPR T0,COSY ;V = T*COSY
MOVM T0,T0 ;V = ABS(V)
HRLOI T3,377777 ;T3=XMAX
FDVR T3,T2 ;Q = XMAX/T
CAIG T0,T3 ;IF V IS LESS THAN OR EQUAL TO Q
JRST LAB ; THEN GO TO LAB
ERR0: HRLOI T0,377777 ;RRES=XMAX
$LCALL RPO
;LERR (LIB,%,<CEXP: real part overflow>)
JRST DD2
LAB: FMPR T0,T2 ;RRES = V*T
DD2: MOVE T2,COSY
CAIGE T2, ;IF COSY.LT. 0.0
MOVN T0,T0 ;THEN NEGATE RRES
JRST RET ;RETURN
ALG: MOVEM T2,XSAVE
FUNCT EXP.,<XSAVE> ;EXPX = EXP(X)
MOVEM T0,EXPX
MOVE T1,EXPX
FMPR T1,SINY ;IRES = EXPX*SINY
JFCL
FMPR T0,COSY ;RRES = EXPX*COSY
JFCL
CAIE T1, ;IF IRES .NE. 0.0
JRST RCHK ;THEN GO CHECK RRES
YCHK: MOVE T2,YSAVE
CAIN T2, ;IF Y .NE. 0
JRST RCHK ;THEN GO CHECK RRES
$LCALL IPU
;LERR (LIB,%,<CEXP: imaginary part underflow>)
RCHK: CAIN T0, ;IF R = 0.0
$LCALL RPU
;LERR (LIB,%,<CEXP: real part underflow>)
RET: POP P,T4
POP P,T3
POP P,T2
GOODBY (1) ;RETURN
YMAX: 823549.66 ;[3234] MAX ARG TO SIN/COS
TBIG: 207540074636 ;88.029692
NEGCON: 570232254037 ;-89.415986
SEGMENT DATA
SINY: 0
COSY: 0
XSAVE: 0
EXPX: 0
YSAVE: 0
PRGEND
TITLE CEXP2 COMPLEX ** INTEGER
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CEXP2
EXTERN CEXP2.
CEXP2=CEXP2.
PRGEND
TITLE CEXP3 COMPLEX ** COMPLEX
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CEXP3
EXTERN CEXP3.
CEXP3=CEXP3.
PRGEND
TITLE CEXP3. Complex ** complex exponenentiation
SUBTTL Mary Payne /MHP 29-May-81
; New Routine for CEXP3; Mary Payne, May 29, 1981.
;
;This routine calculates (X + iY)**(A + iB), where i is SQRT(-1).
SEARCH MTHPRM
SEGMENT CODE
SALL
;[3247]THE CEXP2 ENTRY:
;[3247]
;[3247]CEXP2 CALCULATES (X + IY)**N, WHERE N IS A (REAL) INTEGER.
;[3247]THE INITIALIZATION FOR THE CEXP2 ENTRY IS GIVEN BELOW. THIS
;[3247]CODE COMPLETES THE CALCULATION FOR THE EXPONENT = 0 OR 1,
;[3247]AND ALSO COMPLETES THE CASE FOR WHICH THE COMPLEX BASE IS
;[3247]0 + I*0. IF X = 0, THE ROUTINE MERGES WITH CEXP3 AT XTINY;
;[3247]IN BLOCK 2. ALL OTHER CASES MERGE WITH THE CEXP3 ROUTINE AT
;[3247]JOIN3: IN BLOCK 2.
HELLO (CEXP2,.) ;[3247]CEXP2, COMPLEX**INTEGER
DMOVEM T4,SAVT4 ;[3247]SAVE REGISTERS T4,T5
DMOVE T0,@(L) ;[3247]X + IY TO T0,T1
MOVE T4,@1(L) ;[3247]INTEGER EXPONENT TO T4
JUMPE T4,I0 ;[3247]SEPARATE OUT EXPONENT = 0
MOVE T5,T4 ;[3247]COPY EXPONENT TO T5
SOJE T5,I1 ;[3247]SEPARATE OUT EXPONENT = 1
JFCL ;[3247] SUPPRESS ERROR MESSAGE
JUMPE T1,CX1 ;[3251]IMAGINARY PART = 0
JUMPN T0,CX2 ;[3251]ARGS (<>0,<>0) = NORMAL CASE
JRST NOREAL ;[3251]NO REAL PART = SPECIAL CASE
CX1: JUMPL T0,NOIMAG ;[3251]ARGS (-REAL,0) = SPECIAL CASE
JUMPE T0,ARG0 ;[3251]ARGS (0,0) = SPECIAL CASE
;
;[3247]THE CX2: BLOCK OF CODE DOES THE PROPER INITIALIZATION
;[3247]FOR MERGING WITH CEXP3 AT JOIN3: OR XTINY: IN BLOCK 2.
CX2: DMOVEM T2,SAVT2 ;[3247]SAVE REGISTERS T2,T3
PUSHJ P,DFL.4## ;[3247]D. P. FLOAT EXPONENT
DMOVEM T4,REXP ;[3247]SAVE REAL PART OF EXPONENT
SETZB T4,ISIGN ;[3247]INITIALIZE IMAGINARY SIGN FLAG
SETZB T5,RSIGN ;[3247] REAL SIGN FLAG AND PHI
MOVEM T4,PHIFLG ;[3247] FLAG TO ZERO
DMOVEM T4,IEXP ;[3247]Save 0 IMAG. PART OF EXPONENT
MOVE T4,T0 ;[3247]COPY X TO T4 WITH 0 EXTEND
MOVE T2,T1 ;[3247] IN T5. COPY Y TO T2
SETZ T3, ;[3247] AND ZERO EXTEND
EXCH T0,T1 ;[3247]Y TO T0, X TO T1
JUMPE T1,XTINY ;[3247]IF X = 0, MERGE WITH CEXP3 AT XTINY
JRST JOIN3 ;[3247] ELSE AT JOIN3, BOTH IN BLOCK 2.
;[3247]THE REST OF THE CODE DEALS WITH THE SPECIAL CASES OF
;[3247]EXPONENT = 0 OR 1, AND X = Y = 0.
I1: DMOVE T4,SAVT4 ;[3247]EXPONENT IS 1. RESTORE REGISTERS
;[3247]T4, T5, AND RETURN THE COMPLEX
;[3247]ARGUMENT AS THE RESULT
POPJ P, ;[3247]RETURN
I0: JUMPN T0,ANSRVM ;[3247]EXPONENT IS 0
JUMPE T1,INDET ;[3247]IF X = Y = 0, RESULTS INDETERMINATE
ANSRVM: MOVE T0,SPONE ;[3247]IF X AND Y NOT BOTH 0,
SETZ T1, ;[3247] (X+IY)**0 IS 1.0 + I**0
DMOVE T4,SAVT4 ;[3247]RESTORE T4,T5 AND
POPJ P, ;[3247]RETURN
INDET: SETZB T0,T1 ;[3247]FOR 0**0, RETURN 0 + I*0
$SCPX ;SAVE FIXED-UP RESULT
;[3247] LERR (LIB,%,<CEXP2:(0,0)**0 IS INDETERMINATE, RESULT IS (0,0)>)
$LCALL ZZZ ;[3247]ZERO**ZERO IS INDETERMINATE
$RCPX ;GET FIXED-UP RESULT
DMOVE T4,SAVT4 ;[3247]RESTORE REGISTERS T4, T5
POPJ P, ;[3247]RETURN
;[3251]
;[3251] >>>>> Special case (0,iY)**n <<<<<
;[3251]
NOREAL: MOVEM T1,EXPARG ;[3251]COPY OF THE IMAGINARY PART
MOVEM T4,EXPARG+1 ;[3251]COPY OF THE EXPONENT
FUNCT EXP2.,<EXPARG,EXPARG+1>;[3251]DO THE EXPONENT THE EASY WAY
TRNE T4,2 ;[3251]NUMBER NEED TO BE NEGATED?
MOVN T0,T0 ;[3251]YES; NEGATE THE RESULT.
SETZ T1, ;[3251]ASSUME EXPONENT EVEN;
;[3251]RESULT = (IY^N,0)
TRNE T4,1 ;[3251]IS EXPONENT ODD?
EXCH T0,T1 ;[3251]YES; RESULT = (0,IY^N)
DMOVE T4,SAVT4 ;[3251]RESTORE REGISTERS T4, T5
POPJ P, ;[3251]RETURN
;[3251]
;[3251] >>>>> SPECIAL CASE (-X,0)**N <<<<<
;[3251]
NOIMAG: MOVEM T0,EXPARG ;[3251]COPY OF THE REAL PART
MOVEM T4,EXPARG+1 ;[3251]COPY OF THE EXPONENT
FUNCT EXP2.,<EXPARG,EXPARG+1>;[3251]DO THE EXPONENT THE EASY WAY
SETZ T1, ;[3251]NO IMAGINARY PART
DMOVE T4,SAVT4 ;[3251]RESTORE REGISTERS T4, T5
POPJ P, ;[3251]RETURN
ARG0: JUMPL T4,INFIN ;[3247]ARG = (0,0). JUMP IF EXPONENT < 0
DMOVE T4,SAVT4 ;[3247]ELSE RESULT IS 0 + I*0. RESTORE
POPJ P, ;[3247]REGISTERS T4,T5 AND RETURN
INFIN: HRLOI T0,377777 ;[3247](0 + I*0)**NEGATIVE IS IMPLIED DIV
MOVE T1,T0 ;[3247] BY 0. BOTH PARTS = 377777777777
;[3247] LERR (LIB,%,<CEXP2: (0,0)**NEGATIVE: RESULT = (BIGGEST,BIGGEST)>)
$LCALL ZNI ;[3247]ZERO**NEGATIVE IS INFINITY
DMOVE T4,SAVT4 ;[3247]RESTORE REGISTERS T4,T5
POPJ P, ;[3247]Return
;[3247]
;[3247] CEXP3 ROUTINE STARTS HERE:
;[3247]
HELLO (CEXP3,.) ;CEXP3, complex ** complex
;
; *****************
; * BLOCK 1 *
; *****************
;
;Separate out the special cases: X = Y = 0 and A = B = 0:
MOVEM T2,SAVT2 ;Save Registers T2
MOVEM T3,SAVT3 ; and T3
MOVEM T4,SAVT4 ; and T4
DMOVE T1,@(L) ;X and Y to T1 and T2.
DMOVE T3,@1(L) ;A and B to T3 and T4.
JUMPN T1,NXTST ;[3247]If X not zero, continue
; with main flow below.
JUMPN T2,XIS0 ;Special case: X = 0, Y not.
; Treated in BLOCK 2 after
; its main flow is complete.
; >>>> Special Case: X = Y = 0 <<<<
;
;In the following block of code X = Y = 0. No detailed
;calculations are needed. If A > 0, both the real and
;imaginary parts of the result are zero. If A = 0
;the complex result has modulus one, and indeterminate
;direction, so both its real and imaginary parts should
;be identified as indeterminate. If A < 0, the complex
;result has modulus infinity, and indeterminate direction.
;The real and imaginary parts should probably be identified
;as indeterminate. A is in T3.
;
JUMPG T3,ZERO ;If A > 0, result is (0,0).
JUMPL T3,BTHIND ;If A = 0, results indeterminate.
$LCALL ZZZ ;zero**zero is indeterminate, result = zero
SETZB T0,T1 ;Return (0,0)
JRST RST234 ;Restore registers 2, 3, and 4.
BTHIND: HRLOI T0,377777 ;Else results are indeterminate
HRLOI T1,377777 ;Provide something for storage
JUMPE T4,CPXZNI ;[4002]If no complex part, treat like
;[4002] zero**neg real
$LCALL ZCI ;(0,0)**(negative,non-zero) is indeterminate
JRST RST234 ;Restore registers 2, 3, and 4.
CPXZNI: $LCALL ZNI ;[4002] zero**neg , result=infinity
JRST RST234 ;[4002] restore registers 2, 3, and 4.
ZERO: SETZ T0, ;Result is (0,0). Clear T0.
; T1 already has 0.
JRST RST234 ;Restore registers 2, 3, and 4.
;
; >>>> End of X = Y = 0 <<<<
;
; >>>> Main Flow continues at NXTST <<<<
;
NXTST: JUMPN T3,POLAR ;If A and B are not both zero,
JUMPN T4,POLAR ; continue with main flow.
; in BLOCK 2.
;
; >>>> Special Case: A = B = 0 (X,Y not both 0) <<<<
;
MOVE T0,SPONE ;A = B = 0; X, Y not both zero
SETZ T1, ; so result is (1,0)
RST234: MOVE T2,SAVT2 ;Restore registers 2
MOVE T3,SAVT3 ; and 3
MOVE T4,SAVT4 ; and 4
POPJ P, ;and return.
;
; >>>> End of A = B = 0 <<<<
;
;End of BLOCK 1. The cases corresponding to X = Y = 0 and/or
;A = B = 0 are complete.
;
; *****************
; * BLOCK 2 *
; *****************
;
;In this block of code X and Y are not both zero. A and B
;are not used. (X + iY) is rewritten in the form
;EXP(CDLOG(X + iY)), and the CDLOG routine is invoked. The
;real part of CDLOG is stored in LNRHO, and its imaginary
;part in THETA. Either, but not both of, THETA and LNRHO
;can underflow, and code is inserted to return scaled values
;instead of zero for these cases. On entry to this block X
;and Y are in T1 and T2. T3 and T4, which contain A and B,
;are used here, but A and B will be recovered later, when
;they are needed,from REXP and IEXP, respectively.
;
POLAR: MOVEM T5,SAVT5 ;Save another register.
SETZ T5, ;Clear T5 for initialization
MOVEM T5,RSIGN ; of flags for signs of the
MOVEM T5,ISIGN ; real and imaginary parts
; of the final result.
MOVEM T5,PHIFLG ;Initialize flag for
; underflow of PHI.
DMOVEM T4,IEXP ;B is now zero extended.
; Store B in double precision.
SETZ T4, ;Zero extend A, and store
DMOVEM T3,REXP ; it in double precision.
MOVE T4,T1 ;Y in T2, copy X to T4
SETZ T3, ;Zero extend Y to double.
; Already done for X in T4,T5.
;
;Separate out cases leading to underflow of THETA and
;LNRHO, so as to calculate scaled values. All branches follow
;completion of the main flow for this BLOCK.
;
MOVE T0,T2 ;Copy of Y to T0.
JOIN3: FDV T0,T1 ;[3247]S.P. chopped DIV shows
JFCL EXCEP ; possible exceptions.
MOVM T0,T0 ;Test |Y/X|.
CAML T0,TWO63 ;If |Y/X| huge: |THETA| = PI/2.
JRST YMGTX ; Go to special cases.
CAMG T0,TWOM63 ;If |Y/X| tiny: |THETA| = 0 or
JRST XMGTY ; PI. Go to special cases.
;
;
DMOVEM T4,ARG1 ;Prepare to get CDLOG
DMOVEM T2,ARG1+2 ; of (X,Y). There will be
FUNCT CDLOG.,<ARG1,LNRHO> ; no exceptions.
JRST PHI0 ;Go to BLOCK 3.0 for PHI. Neither
; THETA nor LNRHO is scaled.
;
;Main Flow for BLOCK 2 is complete. Now take care of the
;special cases.
;
EXCEP: JUMPN T0,YMGTX ;On overflow, |Y| >> |X|.
;
; >>>> Treatment of underflow on Y/X. <<<<
;
JUMPL T4,THETPI ;Underflow and X < 0: THETA is PI.
;Underflow implies |X| not 1.
;Also |X| >> |Y|; hence LNRHO
DMOVEM T4,LNARG ; is LN(|X|) which is not 0.
; X > 0, so |X| = X. Get
FUNCT DLOG.,<LNARG> ; LN(X) = LN(RHO)
DMOVEM T0,LNRHO ; and store it.
;
;Since Y/X underflows and X > 0
; Theta underflows. Get and
; store a scaled value.
FSC T2,140 ;[3247]Scale Y by 2**96 and X by
FSC T4,-140 ;[3247] 2**(-96). No exceptions on
DFDV T2,T4 ;Y/X scaled up by 2**192
DMOVEM T2,THETA ;Store scaled THETA.
JRST PHI1 ;Continue with main flow in
; BLOCK 3.2; THETA is scaled
; and LNRHO is not.
;
THETPI: DMOVE T0,PI ;Underflow on Y/X, and X < 0
JUMPGE T2,ABSX ;THETA is PI for Y > 0
DMOVN T0,T0 ; or -PI for Y < 0.
ABSX: MOVN T2,T4 ;[3247] T2 <-- |X|.
JRST GETLN ;Go store THETA and get LNRHO.
;
; >>>> End of underflow on Y/X <<<<
;
; >>>> Treatment of very large |Y/X| <<<<
;
YMGTX: MOVM T1,T2 ;Get |Y|. If it is not 1, then
CAME T1,SPONE ; neither THETA nor LNRHO will
JRST XTINY ;[3247] underflow; Merge with XTINY.
DMOVE T0,PIOV2 ;If |Y| = 1, LNRHO may underflow.
JUMPGE T2,THPOS ;THETA is PI/2
DMOVN T0,T0 ; or -PI/2 if Y < 0. |Y| = 1 so
THPOS: MOVE T2,T4 ;LNRHO = (X**2)/2. Get X in
JRST GETLOG ; T2 and T4. To GETLOG to store
; THETA and get scaled LNRHO.
;
; >>>> End of treatment of very large |Y/X| <<<<
;
; >>>> Treatment of |Y/X| very small; no underflow <<<<
;
XMGTY: JUMPGE T4,TINYTH ;If X > 0, THETA tiny but not 0.
MOVN T4,T4 ;If X < 0, change its sign.
DMOVE T0,PI ; THETA is PI
JUMPGE T2,ISX1 ; if Y > 0, or
DMOVN T0,T0 ; -PI, if Y < 0.
JRST ISX1 ;Go test whether |X| = 1.
;
TINYTH: DMOVE T0,T2 ;Get Y to T0,T1, and calculate
DFDV T0,T4 ; DP value for small THETA.
; No exceptions on DFDV.
;
ISX1: MOVM T2,T2 ;|Y| to T2.
CAMN T4,SPONE ;If |X| is 1, LNRHO may underflow.
JRST XIS1 ; Go find out.
DMOVEM T4,LNARG ;X not 1, so no underflow on LNRHO.
JRST STHETA ;Go store THETA and get LNRHO.
;
XIS1: MOVE T4,T2 ;THETA is not scaled.
JRST GETLOG ;Get LNRHO = (Y**2)/2.
;
; >>>> END of treatment for very small |Y/X| <<<<
;
; >>>> Special case X = 0 from BLOCK 1 <<<<
; >>>> Also branch from |Y/X| very large <<<<
;
XIS0: MOVEM T5,SAVT5 ;[3247]Save another register
SETZ T5, ;[3247]Zero extend T4, and save DP value
DMOVEM T4,IEXP ;[3247] of imaginary part of exponent.
SETZ T4, ;[3247]Zero extend T3, and save DP
DMOVEM T3,REXP ;[3247] value of real part of exponent.
SETZ T3, ;[3247]Zero extend Y in T2
SETZM PHIFLG ;Initialize flag for underflow of PHI
SETZM RSIGN ;Clear flags for signs of
SETZM ISIGN ; final result
XTINY: DMOVE T0,PIOV2 ;[3247]If X = 0 or |Y| >> X
JUMPGE T2,GETLN ;THETA is PI/2
DMOVN T0,T0 ; or -PI/2 if Y < 0.
MOVN T2,T2 ; T2 <-- |Y|. Now LN(RHO)
; will be LN(|Y|). Go store
; THETA and get LNRHO. No
; underflow on LNRHO because
; if X not 0, |Y| is not 1 and
; if X = 0, LN(|Y|) cannot underflow.
;
; >>>> End of special case X = 0 <<<<
; >>>> LNRHO from CALL to LOG routine <<<<
;
GETLN: DMOVEM T2,LNARG ;|X| not 1 and >> |Y|
; or |Y| not 1 and >> |X|
STHETA: DMOVEM T0,THETA ;Store THETA.
FUNCT DLOG.,<LNARG> ;Get LN of larger
DMOVEM T0,LNRHO ; and store it.
JRST PHI0 ;Continue with main flow in
; BLOCK 3.0; neither THETA
; not LNRHO is scaled.
;
; >>>> End of LNRHO by CALL to LOG routine <<<<
;
; >>>> LNRHO for RHO = 1 + (X**2 or Y**2) <<<<
;
GETLOG: DMOVEM T0,THETA ;Store unscaled THETA.
DFMP T4,T4 ;Get X**2 or Y**2. Copy of X or Y
JFCL LNUND ; in T2, in case of underflow.
FSC T4,-1 ;DIV by 2.
JFCL LNUND ; Underflow may occur.
DMOVEM T4,LNRHO ;Store unscaled LNRHO if no
; underflows occurred.
JRST PHI0 ;Continue with main flow in
; BLOCK 3.0; neither THETA
; not LNRHO is scaled.
;
LNUND: FSC T2,140 ;[3247]Scale up saved X or Y.
DFMP T2,T2 ;Square. No exceptions possible.
FSC T2,-1 ;DIV by 2. No exceptions possible.
DMOVEM T2,LNRHO ;Store scaled LNRHO
JRST PHI2 ;Rejoin main flow at BLOCK 3.1 with
; THETA unscaled, LNRHO scaled.
;
; >>>> End of LNRHO for RHO near 1 <<<<
;
;End of BLOCK 2. The CDLOG(X + iY) has real part LNRHO and
;imaginary part THETA, which are stored in locations with
;these names. If no exceptions occur, the main flow continues
;in BLOCK 3.0 at PHI0. No overflows occur. Either, but
;not both of THETA and LNRHO may underflow. For these
;cases values scaled up by 2**192 decimal, or 2**300
;octal, have been stored, and the main flow continues
;in BLOCK 3.2 at PHI1 for underflow of THETA, and in
;BLOCK 3.1 at PHI2 for underflow of LNRHO.
; *****************
; * BLOCK 3 *
; *****************
;
;The complex power function is to be evaluated as the
;complex exponential of (A + iB)*(LNRHO + iTHETA), which
;has real part ALPHA = A*LNRHO - B*THETA, and imaginary
;part PHI = A*THETA + B*LNRHO. The calculation of PHI is carried
;out in BLOCK 3. There are three sub-blocks 3.0, 3.2, and
;3.1 for the cases of no underflow, underflow of THETA and
;underflow of LNRHO, respectively.
;
;Since it is ultimately EXP(i*PHI) that is needed, it would
;appear that SIN and COS of PHI are the final parameters
;needed. However, these functions will get multiplied by
;EXP(ALPHA), and the handling of exception boundaries on
;the product will be expedited by using instead LN(SIN(PHI))
;and LN(COS(PHI)), which will be added to ALPHA before
;invoking the EXP function. Actually is is the absolute
;values of SIN and COS which will be used as arguments to
;the LN function; the signs of SIN and COS will be stored
;in flags for use in determining the signs for the real and
;imaginary parts of the complex exponential. Recall that these
;flags were initialized to zero at the beginning of BLOCK 2.
;
; *****************
; * BLOCK 3.0 *
; *****************
;
;Neither THETA nor LNRHO is scaled.
;
PHI0: DMOVE T0,LNRHO ;Get LNRHO and multiply
DFMP T0,IEXP ; by B.
JFCL PHI0A ; Branch on exception.
DMOVE T2,THETA ;Get THETA and multiply
DFMP T2,REXP ; by A.
JFCL PHI0B ; Branch on exception
DMOVE T4,T0 ;Get copy of B*LNRHO in case
; of exception on next step.
;
;Rejoin the main flow of BLOCK 3.0 here for underflows
; that do no harm. No exception on DFAD for such cases.
;
PHISUM: DFAD T0,T2 ;PHI = A*THETA + B*LNRHO.
JFCL PHI0C ; Branch on exception.
DMOVE T4,T0 ;Get copy of PHI and
JUMPGE T4,ABSARG ; test its sign.
DMOVN T4,T4 ; Get absolute value
ABSARG: CAMGE T4,ARGHI ;Test |PHI| against limit
JRST ARGOK ; for argument reduction
CAMG T4,ARGHI ; in SIN/COS routines.
CAML T5,ARGLO ; If >/= INT(PI*2**31)
JRST PHIBIG ; |PHI| is too big.
ARGOK: DMOVEM T0,TRGARG ;Else, store PHI.
FUNCT DSIN.,<TRGARG> ;Get SIN(PHI) and
DMOVEM T0,SPHI ; store for future use.
;
FUNCT DCOS.,<TRGARG> ;Get COS(PHI) and
DMOVEM T0,CPHI ; store for future use.
;
JRST ALF0 ;PHI calculation is complete.
;Go to BLOCK 4.0 for ALPHA.
;Main flow of BLOCK 3.0 is complete. Now take care of
;special cases.
;
; >>>> B*LNRHO and A*THETA OK; exception on sum <<<<
;
PHI0C: JUMPN T0,PHIBIG ;If overflow, PHI is too big.
DMOVE T0,T4 ;Get B*LNRHO back to T0.
FSC T0,300 ;Scale both B*LNRHO and A*THETA
FSC T2,300 ; up by 2**192.
JRST SCLPHI ;Join with other cases to
; get scaled PHI, etc.
;
; >>>> End of handling exception on sum <<<<
;
; >>>> B*LNRHO OK, exception on A*THETA <<<<
;
PHI0B: JUMPN T2,PHIBIG ;If overflow, PHI is too big.
JUMPE T0,SCTH ;[3247]Don't "FSC" 0!
MOVM T4,T0 ;Get |B*LNRHO|. If it is
CAML T4,TWOM66 ; > 2**(-66), underflow does
JRST PHISUM ; no harm. Go to main flow.
FSC T0,300 ;Else scale B*LNRHO up by 2**192
SCTH: DMOVE T2,THETA ;Also get THETA and A and
DMOVE T4,REXP ; scale each up by 2**96
FSC T2,140 ;[3247] so that resulting product
FSC T4,140 ;[3247] has same scaling as B*LNRHO.
DFMP T2,T4 ;No exceptions on either the
JRST SCLPHI ; MUL or subsequent operations
; in getting scaled PHI, etc.
;
; >>>> End of B*LNRHO OK; exception on A*THETA <<<<
;
; >>>> Exception on B*LNRHO <<<<
;
PHI0A: JUMPN T0,PHIBIG ;If overflow, |PHI| too big.
DMOVE T2,THETA ;B*LNRHO underflows. Get
DFMP T2,REXP ; A*THETA before scaling.
JFCL PHI0D ;Branch on exception.
MOVM T4,T2 ;If |A*THETA| is big enough
CAML T4,TWOM66 ; main flow can be rejoined
JRST PHISUM ; with 0 for negligible B*LNRHO
JUMPE T2,SCLBLN ;[3247]Don't "FSC" 0!
FSC T2,300 ;Else scale A*THETA up by 2**192
JRST SCLBLN ; and go scale B*LNRHO
;
; >>>> End of exception on B*LNRHO <<<<
; >>>> Underflow on B*LNRHO and exception on A*THETA <<<<
;
PHI0D: JUMPN T2,PHIBIG ;If overflow, PHI is too big.
DMOVE T2,THETA ;Both products underflow. Get
DMOVE T4,REXP ;THETA and A and scale each
FSC T2,140 ;[3247] up by 2**96, so that product
FSC T4,140 ;[3247] will be scaled up by 2**192
DFMP T2,T4 ;No exception on product.
SCLBLN: DMOVE T0,LNRHO ;Get LNRHO and B and scale
DMOVE T4,IEXP ; each up by 2**96 so
FSC T0,140 ;[3247] that product will be
FSC T4,140 ;[3247] scaled up by 2**192.
DFMP T0,T4 ;No exception on product.
SCLPHI: DFAD T0,T2 ;No exception on sum.
HRLZI T5,400000 ;Set sign bit of T5.
MOVEM T5,PHIFLG ;Set PHIFLG to indicate
; underflow of PHI.
JUMPGE T0,PHIPL0 ;If scaled PHI < 0
MOVEM T5,ISIGN ; Use T5 to indicate that
; imaginary part is negative.
DMOVN T0,T0 ; Get |PHI| in T0.
PHIPL0: DMOVE T0,LNARG ;Store scaled |PHI|
FUNCT DLOG.,<LNARG> ; and get its LOG
DFAD T0,LN300 ; and add -192*LN(2) to
DMOVEM T0,SPHI ; get LN(SIN(|PHI|))
SETZ T0, ;Clear T0 and T1 in order to
SETZ T1, ; store 0 for LN(COS(|PHI|))
DMOVEM T0,CPHI ; Since COS(TINY) > 0
; real part will be > 0.
JRST ALF0 ;Continue with main flow
; at BLOCK 4.0.
; >>>> End of underflow on both products <<<<
; >>>> PHI is too large <<<<
;
PHIBIG: $LCALL BPI ;Both parts indeterminate
HRLOI T0,377777 ;Store biggest number in
MOVE T1,T0 ; both REAL and IMAG.
JRST RESTOR ;Go to end of BLOCK 6
; to restore registers.
;
;****If |PHI| is too large something should be stored
; in T0,T1, and an error message returned saying that
; neither the real nor the imaginary parts of the
; result are determinate. If |PHI| > int (2**31 * pi), the
; argument reduction for SIN and COS breaks down.
;
;End of BLOCK 3.0. Calculation continues at BLOCK 4.0,
;which also applies to the case in which neither LNRHO
;nor THETA underflows. At this point SPHI contains SIN(PHI)
;or, if PHI underflows, LN(SIN(|PHI|)) = LN(|PHI|). CPHI
;contains COS(PHI) or, if PHI underflows, LN(COS(|PHI|)) = 0.
;If PHI underflows, RSIGN and ISIGN contain the signs of the
;real and imaginary parts of the final result, respectively.
;PHIFLG has 0 or 1 in its sign bit, according as PHI has not
;or has underflowed.
;
; *****************
; * BLOCK 3.1 *
; *****************
;
;LNRHO is scaled up by 2**192. THETA is not scaled. The
;following transformation makes it possible to use the code
;in BLOCK 3.2 for this case also.
;
PHI2: HRLZI T0,400000 ;Must exchange LNRHO,THETA
MOVE T0,XCHFLG ; and A,B. Set sign bit of
; XCHFLG to indicate that
; this has been done.
DMOVE T0,LNRHO ;Get LNRHO
DMOVE T2,THETA ; and THETA
DMOVEM T0,THETA ; and exchange
DMOVEM T2,LNRHO ; them.
MOVE T0,REXP ;Get A
MOVE T2,IEXP ; and B
MOVEM T0,IEXP ; and exchange
MOVEM T2,REXP ; them.
;
;Continue at PHI1 in BLOCK 3.2, which will now carry out
;the PHI calculation correctly for this case.
;
; *****************
; * BLOCK 3.2 *
; *****************
;
;THETA is scaled up by 2**192. LNRHO is not scaled.
;
;However, the same code is used for LNRHO scaled and THETA
;not scaled on JRST from BLOCK 3.1. For this case read
;THETA for LNRHO, LNRHO for THETA, IEXP for REXP and REXP
;for IEXP throughout.
;
PHI1: DMOVE T0,LNRHO ;Get LNRHO and multiply
DFMP T0,IEXP ; by B.
JFCL PHI1A ; Branch on exception.
DMOVE T2,THETA ;Get THETA and multiply
DFMP T2,REXP ; by A.
JFCL PHI1B ; Branch on exception.
JUMPE T2,PHIOK1 ;[3247]Don't "FSC" 0!
DMOVE T4,T2 ;Save copy of A*THETA
FSC T2,-300 ;Try to unscale A*THETA.
JFCL T2,PHI1C ; On failure, branch.
DMOVE T4,T0 ;Save copy of B*LNRHO.
DFAD T0,T2 ;PHI = A*THETA + B*LNRHO.
JFCL PHI1D ; Branch on exception.
;
;Rejoin the main flow of BLOCK 3.2 here for underflows
; that do no harm. PHI is unscaled and ready for final
; processing.
;
PHIOK1: DMOVE T4,T0 ;Get copy of PHI and
JUMPGE T4,ABARG1 ; test its sign.
DMOVN T4,T4 ; Get absolute value
ABARG1: CAMGE T4,ARGHI ;Test |PHI| against limit
JRST ARGOK1 ; for argument reduction
CAMG T4,ARGHI ; in SIN/COS routines.
CAML T5,ARGLO ; If >/= INT(PI*2**31)
JRST PHIBIG ; |PHI| is too big.
ARGOK1: DMOVEM T0,TRGARG ;Store unscaled PHI
FUNCT DSIN.,<TRGARG> ;Get its SIN
DMOVEM T0,SPHI ; and store it.
FUNCT DCOS.,<TRGARG> ;Get COS(PHI)
DMOVEM T0,CPHI ; and store it.
JRST ALF1 ;PHI calculation is complete.
;Go to BLOCK 4.1 for ALPHA.
;Main Flow of BLOCK 3.2 is complete. Now take care of
;special cases.
;
; >>>> B*LNRHO OK, unscaled A*THETA OK <<<<
;
; >>>> Exception on their sum <<<<
;There was no exception on A*(scaled THETA). Hence the unscaled
;product < 2**(-100), hence sum could not overflow, hence sum
;underflowed. T4 has unscaled B*LNRHO and T2 has scaled A*THETA.
;
PHI1D: DMOVE T0,T4 ;Get B*LNRHO to T0
FSC T0,300 ;Scale it up by 2**192.
FSC T2,300 ;Rescale A*THETA
DFAD T0,T2 ;No exceptions on scaled PHI.
JRST SCPHI1 ;Go finish up for scaled PHI.
;
; >>>> End of exception on sum of unscaled products <<<<
;
; >>>> B*LNRHO OK, A*(scaled THETA) OK <<<<
; >>>> Unscaled A*THETA underflows <<<<
;T0 has B*LNRHO and T4 has A*(scaled THETA).
;
PHI1C: MOVM T2,T0 ;Get |B*LNRHO| in T2
CAML T2,TWOM66 ;If it is big enough
JRST PHIOK1 ; underflow of A*THETA is
; harmless: A*THETA is
; negligible. Rejoin main
; flow at PHIOK1 above.
DMOVE T2,T4 ;A*(Scaled THETA) back to T2.
JUMPE T0,ADD0 ;[3247]Avoid an "FSC" of 0
JRST SCBLN1 ;Go scale B*LNRHO up and
; add to A*(scaled THETA)
; in T2 to get scaled PHI.
;
; >>>> End of small B*LNRHO + underflowed A*THETA <<<<
;
; >>>> B*LNRHO OK, exception on A*(scaled THETA) <<<<
;
PHI1B: JUMPE T2,PHIOK1 ;A*(scaled THETA) underflows,
; hence A*THETA is negligible.
; Rejoin main flow at PHIOK1 above.
JRST REMOV1 ;Overflow on A*(scaled THETA) can
; be removed by unscaling product.
;
; >>>> End of B*LNRHO OK, exception on A*(scaled THETA) <<<<
; >>>> Exception on B*LNRHO <<<<
;
PHI1A: JUMPN T0,PHIBIG ;If overflow, PHI is too big.
DMOVE T2,THETA ;B*LNRHO underflows. Get
DFMP T2,REXP ; THETA and multiply by A.
JFCL PHI1E ; Branch on exception.
JRST SCBLN1 ;Go scale B*LNRHO up and add
; to scaled A*THETA in T2 to
; get scaled PHI.
;
PHI1E: JUMPE T2,SCBLN1 ;B*LNRHO and A*(scaled THETA)
; both underflow. A*THETA is
; negligible. Go scale B*LNRHO
; to get scaled PHI, with 0 in
; T2 for negligible A*THETA.
;
;If A*(scaled THETA) overflows, removal of scaling brings
;unscaled |A*THETA| into range between 2**(-100) to 2**100.
;Hence underflowed B*LNRHO is negligible. Code in REMOV1
;will add in 0 (for negligible B*LNRHO) to unscaled A*THETA
;to get unscaled PHI.
;
;Branch from PHI1B to remove overflow from A*(scaled THETA)
;joins here, with B*LNRHO OK and in T0 to be added to
;unscaled A*THETA.
;
REMOV1: DMOVE T2,THETA ;Overflow can be removed by
DMOVE T4,REXP ; scaling THETA and A down by
FSC T2,-140 ;[3247] 2**(-96) each, to get
FSC T4,-140 ;[3247] unscaled product.
DFMP T2,T4 ;No exceptions.
DFAD T0,T2 ;No exceptions since |A*THETA| is
; between 2**(-100) and 2**100
JRST PHIOK1 ;Rejoin main flow at PHIOK1 to get
; SIN and COS of unscaled PHI.
;
; >>>> End handling exception on B*LNRHO <<<<
;
; >>>> Completion of calculation for scaled PHI <<<<
;
SCBLN1: DMOVE T0,LNRHO ;Get LNRHO and B and
DMOVE T4,IEXP ; scale each up by 2**96
FSC T0,140 ;[3247] so as to get product
FSC T4,140 ;[3247] scaled up by 2**192.
DFMP T0,T4 ;No exception on product or
ADD0: DFAD T0,T2 ; on ADD to scaled A*THETA.
;
;T0 now has scaled PHI. Branches from other special cases
;join in here to get signs of the final real and
;imaginary parts, LN(SIN(|PHI|)) and LN(COS(|PHI|)).
;
SCPHI1: HRLZI T4,400000 ;Set sign bit of T4
MOVEM T4,PHIFLG ;Set flag for PHI scaled.
JUMPGE T0,PHIPL1 ;If scaled PHI is negative
DMOVN T0,T0 ; get |scaled PHI| in T0.
MOVEM T4,ISIGN ; Imaginary part negative.
PHIPL1: DMOVEM T0,LNARG ;Store scaled PHI
FUNCT DLOG.,<LNARG> ; and get its LOG.
DFAD T0,LN300 ; and subtract 192*LN(2).
DMOVEM T0,SPHI ;SPHI gets LN(|PHI|) which
; = LN(SIN(|PHI|)).
SETZ T2, ;Clear T2 and T3 so as to
SETZ T3, ; store 0 for LN(COS|PHI|))
DMOVEM T2,CPHI ; since COS(TINY) = 1. Real
; part of result is positive.
JRST ALF1 ;Now go get ALPHA for LNRHO
; unscaled, THETA scaled.
;
; >>>> End of calculation for scaled PHI <<<<
;
;End of BLOCK 3.2. All paths lead either to the error return
;PHIBIG or to ALF1 for the calculation of ALPHA for the case
;LNRHO unscaled and THETA scaled up by 2**192. If PHI is
;unscaled, PHIFLG is clear and SPHI and CPHI contain SIN(PHI)
;and COS(PHI), respectively. IF PHI is scaled, SPHI contains
;LN(SIN(|PHI|)), CPHI contains LN(COS(|PHI|)) = 0, PHIFLG
;has its sign bit set, and RSIGN and ISIGN contain the signs
;of the real and imaginary parts of the final result.
;
; *****************
; * BLOCK 4 *
; *****************
;
;The real and imaginary parts of the final answer are
;given by EXP(ALPHA) times COS(PHI) and SIN(PHI),
;respectively. ALPHA = A*LNRHO - B*THETA is calculated
;in BLOCK 4.0 for no scaling of THETA or LNRHO. It is
;calculated in BLOCK 4.1 otherwise. It can be shown that
;if ALPHA overflows, neither the SIN nor the COS of PHI
;can underflow sufficiently to bring EXP(ALPHA) back into
;range. If ALPHA overflows and is positive both parts
;of the result overflow: Code providing a message and
;signed largest numbers for the real and imaginary parts
;is given in EXPOV at the end of BLOCK 4.1. The signs of
;the largest numbers are determined by COS(PHI) and SIN(PHI),
;respectively. If ALPHA overflows and is negative, EXP(ALPHA)
;underflows. The code providing a message and returning 0 for
;both parts is deferred to BLOCK 5, in which similar cases
;occur.
;
; *****************
; * BLOCK 4.0 *
; *****************
;
;Neither THETA nor LNRHO is scaled. ALPHA = A*LNRHO - B*THETA.
;
ALF0: DMOVE T0,LNRHO ;Get LNRHO and multiply
DFMP T0,REXP ; by A.
JFCL ALF0A ; Branch on exception.
DMOVE T2,THETA ;Get THETA and multiply
DFMP T2,IEXP ; by B.
JFCL ALF0B ; Branch on exception
DFSB T0,T2 ;ALPHA = A*LNRHO - B*THETA.
JFCL ALF0C ; Branch on exception.
DMOVEM T0,ALPHA ;Store ALPHA for future use.
;
JRST ANS0 ;Calculation of ALPHA complete.
;Go to BLOCK 5.0 for final answer.
;
;Main flow of BLOCK 4.0 is complete. Now take care of
;special cases.
;
; >>>> A*LNRHO and B*THETA OK; exception on difference <<<<
;
ALF0C: JUMPN T0,DIFFOV ;Exception on final results.
DMOVEM T0,ALPHA ;Zero is acceptable if
; ALPHA underflows.
JRST ANS0 ;Calculation of ALPHA complete.
;Go to BLOCK 5.0 for final answer.
;
DIFFOV: JUMPL T2,EXPOV ;Overflow with negative B*THETA
; means overflow on final results.
JRST EXPUND ;Underflow on final results
; for B*THETA positive.
;
; >>>> End of handling exception on difference <<<<
;
; >>>> Exception on A*LNRHO <<<<
;
ALF0A: JUMPN T0,ALNOV ;Exception on final reaults.
DMOVE T2,THETA ;A*LNRHO underflows. Get
DFMP T2,REXP ; B*THETA.
JFCL ALF0B ;Branch on exception.
DMOVN T2,T2 ;Underflow on A*LNRHO means it
DMOVEM T2,ALPHA ; is negligible. ALPHA = -B*THETA.
JRST ANS0 ;Go to BLOCK 5 for final answer.
;
ALNOV: DMOVE T2,THETA ;A*LNRHO overflows. What about
DFMP T2,IEXP ; B*THETA?
JFCL BOTHEX ; On exception to BOTHEX.
JRST ALNOV1 ;A*LNRHO overflows, B*THETA doesn't.
;
;A*LNRHO overflows and B*THETA has exception.
;
BOTHEX: JUMPE T2,ALNOV1 ;B*THETA underflows.
DMOVE T0,LNRHO ;Both products overflow.
DMOVE T2,REXP ; Get A and LNRHO and scale
FSC T0,-140 ;[3247] each down to get scaled
FSC T2,-140 ;[3247] A*LNRHO.
DFMP T0,T2 ;No exception
DMOVE T2,THETA ;Get B and THETA
DMOVE T4,IEXP ; and scale each down
FSC T2,-140 ;[3247] so get B*THETA
FSC T4,-140 ;[3247] scaled like B*LNRHO
DFMP T2,T4 ;No exception.
DFSB T0,T2 ;No exception on scaled ALPHA.
JUMPL T0,EXPUND ;Results underflow if ALPHA < 0.
JRST EXPOV ;Results overflow if ALPHA > 0.
;
;End of overflow on both A*LNRHO and B*THETA.
;
;A*LNRHO overflows, B*THETA does not. Sign of ALPHA is
;determined by sign of A*LNRHO. Neither A nor LNRHO is 0.
;
;NOTE: A similar case from BLOCK 4.1 joins in here.
;
ALNOV1: SKIPL T0,LNRHO ;If LNRHO > 0
JRST LNPOS ; go test B.
SKIPL T0,REXP ;LNRHO < 0. If A
JRST ALNNEG ; > 0, A*LNRHO < 0.
JRST ALNPOS ;Else A*LNRHO > 0.
LNPOS: SKIPL T0,REXP ;LNRHO > 0. If A
JRST ALNPOS ; > 0, A*LNRHO > 0.
ALNNEG: JRST EXPUND ;A*LNRHO < 0 means
; final results underflow.
ALNPOS: JRST EXPOV ;A*LNRHO > 0 means
; final results overflow.
;
; >>>> End of exception on A*LNRHO <<<<
;
; >>>> A*LNRHO OK or underflows, exception on B*THETA <<<<
;
ALF0B: JUMPN T2,BTHOV ;Exception on final results.
DMOVEM T0,ALPHA ;Underflow on B*THETA means it
; is negligible. ALPHA = A*LNRHO.
; Note that if A*LNRHO underflows
; it is also negligible and 0 is
; an acceptable result for ALPHA.
JRST ANS0 ;Go to BLOCK 5 for final answer.
;
;On overflow of B*THETA, neither B nor THETA is 0. A*LNRHO OK
;or underflows so sign of ALPHA is negative of sign
;of overflowed B*THETA.
;
BTHOV: SKIPL T0,THETA ;If THETA > 0
JRST THPOS1 ; go test B.
SKIPL T0,IEXP ;THETA < 0. If B
JRST BTHNEG ; > 0, B*THETA < 0.
JRST BTHPOS ;Else B*THETA > 0.
THPOS1: SKIPL T0,IEXP ;THETA > 0. If B
JRST BTHPOS ; > 0, B*THETA > 0.
BTHNEG: JRST EXPOV ;B*THETA < 0 means -B*THETA > 0,
; so final results overflow.
BTHPOS: JRST EXPUND ;B*THETA > 0 means -B*THETA < 0,
; so final results underflow.
;
; >>>> End of A*LNRHO OK; exception on B*THETA <<<<
;
;End of BLOCK 4.0. All paths lead to EXPOV, at the end
;of BLOCK 4.1 or to ANS0 or EXPUND, both in BLOCK 5.
;
; *****************
; * BLOCK 4.1 *
; *****************
;
;If PHIFLG = 0, THETA is scaled up by 2**192 and LNRHO is not
;scaled. If PHIFLG < 0,LNRHO is scaled up by 2**192, and THETA
;is not scaled.
;
;The code, as written, applies to PHIFLG = 0.
;
;For PHIFLG < 0, read THETA for LNRHO, LNRHO for THETA, (-IEXP)
;for REXP and (-REXP) for IEXP throughout.
;
ALF1: SKIPL PHIFLG ;If PHIFLG set,
JRST ALF2 ; to code for THETA scaled.
MOVN T0,REXP ;Else, exchange already done
MOVEM T0,REXP ; for THETA,LNRHO and A,B.
MOVN T0,IEXP ; Negate exchanged A,B to use
MOVEM T0,IEXP ; the code for scaled LNRHO,
; instead of scaled THETA.
;
ALF2: DMOVE T0,LNRHO ;Get LNRHO and multiply
DFMP T0,REXP ; by A.
JFCL ALF2A ; Branch on exception.
DMOVE T2,THETA ;Get THETA and multiply
DFMP T2,IEXP ; by B.
JFCL ALF2B ; Branch on exception
JUMPE T2,SUB0 ;[3247]Don't "FSC" 0
DMOVE T4,T2 ;Save a copy of B*THETA
FSC T2,-300 ; and try to remove scaling.
JFCL STALF ; Didn't work: B*THETA
; underflows and hence is
; negligible. ALPHA is
; A*LNRHO in To.
SUB0: DFSB T0,T2 ;[3247]Neither product has an exception.
;ALPHA = A*LNRHO - B*THETA.
JFCL ALF2C ; Branch on exception.
STALF: DMOVEM T0,ALPHA ;Store ALPHA for future use.
JRST ANS0 ;Calculation of ALPHA complete.
;Go to BLOCK 5.0 for final answer.
;
;
;Main flow of BLOCK 4.1 is complete. Now take care of
;special cases.
;
; >>>> A*LNRHO and B*THETA OK, exception on difference <<<<
;
ALF2C: JUMPE T0,STALF1 ;On underflow 0 OK for ALPHA.
JUMPL T2,EXPOV ;On overflow ALPHA has sign
JRST EXPUND ; opposite to B*THETA. Hence
; results overflow for T2 < 0
; and underflow for T2 > 0.
;
; >>>> End of exception on difference <<<<
;
; >>>> Exception on A*LNRHO <<<<
;
ALF2A: JUMPN T0,ALNOV1 ;A*LNRHO overflows, B*THETA
; will be negligible. Final
; results will overflow
; or underflow. Go merge with
; similar case in BLOCK 4.0.
DMOVE T2,THETA ;A*LNRHO underflows and is
DFMP T2,IEXP ; negligible. Get B*(scaled
JFCL ALF2D ; THETA). Branch on exception.
JUMPE T2,STALF1 ;[3247]Don't "FSC" 0!
FSC T2,-300 ;Try to unscale B*THETA.
JFCL ; Failed: B*THETA underflows too.
; Hence 0 is acceptable for ALPHA.
STALF1: DMOVEM T0,ALPHA ;Store 0 for ALPHA and go to
JRST ANS0 ; BLOCK 5 for final answer.
;
;Exception on B*(scaled THETA); A*LNRHO underflows.
;
ALF2D: JUMPN T2,REMOV2 ;Go unscale overflowed
; B*(scaled THETA). Negligible
; A*LNRHO stored as 0 in T0.
DMOVEM T0,ALPHA ;Both products underflow so
JRST ANS0 ; 0 is acceptable for ALPHA
; To BLOCK 5 for final answer.
;
; >>>> End of Exception on A*LNRHO <<<<
;
; >>>> A*LNRHO OK, exception on B*(scaled THETA) <<<<
;
ALF2B: JUMPN T2,REMOV2 ;Overflow: unscale A*THETa.
DMOVEM T0,ALPHA ;Underflow: A*THETA negligible.
; ALPHA = A*LNRHO (or 0 if
; here from ALF2C).
JRST ANS0 ;ALPHA done. Go to BLOCK 5.
;
REMOV2: DMOVE T2,THETA ;Get scaled THETA and B
DMOVE T4,IEXP ; and scale each down by
FSC T2,-140 ;[3247] 2**-96 to remove scaling
FSC T4,-140 ;[3247] from product.
DFMP T2,T4 ;|B*THETA| in 2**(-100) to 2**100
DFSB T0,T2 ; so no exception on DFSB.
DMOVEM T0,ALPHA ;Store ALPHA and go to
JRST ANS0 ; to get final answer.
;
; >>>> End of exception on A*(scaled THETA) <<<<
;
; >>>> Overflow of EXP(ALPHA) <<<<
;
EXPOV: HRLOI T0,377777 ;Get biggest number in T0.
SKIPN T0,PHIFLG ;If PHI underflowed get
JRST EXPOV1 ; signs from ISIGN and RSIGN.
DMOVE T2,SPHI ;Get SIN(PHI)
JUMPN T2,EXPOV2 ;If SIN(PHI) is not 0
JUMPN T3,EXPOV2 ; both parts overflow.
SETZ T1, ;If SIN(PHI) = 0, so does
; imaginary part. Real part
SKIPGE T0,CPHI ; overflows. If COS(PHI) < 0
MOVN T0,T0 ; real part is negative.
$LCALL RPO ;Real part overflow
JRST RESTOR ;Go restore registers and exit.
;
EXPOV2: MOVE T1,T0 ;Get Biggest in T1 too.
SKIPGE T0,CPHI ;If COS(PHI) < 0, so
MOVN T0,T0 ; is real part.
SKIPGE T0,SPHI ;If SIN(PHI) < 0, so
MOVN T1,T1 ; so is imaginary part.
JRST ERRMSG ;Go get error message.
;
EXPOV1: MOVE T1,T0 ;Get Biggest in T1 too.
SKIPGE T0,RSIGN ;If real part is negative
MOVN T0,T0 ; negate biggest number.
SKIPGE T0,ISIGN ;If imaginary part is negative
MOVN T1,T1 ; negate biggest number.
ERRMSG: $LCALL BPO ;Both parts overflow
JRST RESTOR ;Go restore registers and exit.
;
;RESTOR for restoring registers is at the end of BLOCK 5.
;
; >>>> End of overflow of EXP(ALPHA) <<<<
;
;End of BLOCK 4.1. All paths lead to ANS0 or EXPUND, both
;in BLOCK 5, or to RESTOR at the end of BLOCK 6.
; *****************
; * BLOCK 5 *
; *****************
;
;The main flow for this BLOCK starts at ANS0. If PHIFLG
;is 0, SPHI and CPHI contain SIN(PHI) and COS(PHI). PHI
;did not underflow. ALPHA is in ALPHA; no tests on its size
;have yet been made. The special cases, including PHIFLG = 1
;start at ANS1, following completion of the main flow.
;
ANS0: DMOVE T0,ALPHA ;Get ALPHA. Will EXP(ALPHA)
CAMLE T0,AMINHI ; underflow? If ALFHI > AMINHI
JRST EXPOK ; EXP(ALPHA) won't underflow.
CAML T0,AMINHI ;If ALFHI < AMINHI, EXP(ALPHA)
; underflows, and so do both
; parts of answer.
CAMG T1,AMINLO ;[3247]ALFHI = AMINHI; test LO's
JRST EXPUND ; ALFLO </= AMINLO: underflow.
EXPOK: SKIPGE T0,PHIFLG ;Was PHI scaled? If so
JRST ANS2 ; go to alternate code.
CAML T0,AMAXHI ;If ALFHI >/= AMAXHI
JRST ANS1 ; use LOG approach.
FUNCT DEXP.,<ALPHA> ;No exception on EXP(ALPHA).
DMOVE T2,T0 ;Keep copy of EXP(ALPHA)
DFMP T0,SPHI ;Get imaginary part; no overflow
JFCL IMUND ; get message on underflow.
STIMP: DMOVEM T0,IMAG ;Store imaginary part.
DFMP T2,CPHI ;Get real part, no overflow
JFCL RLUND ; get message on underflow
STRLP: DMOVEM T2,REAL ;Store real part.
JRST DPTOSP ;Go convert to single precision.
; in BLOCK 6. Main Flow of
; BLOCK 5 complete.
;
; >>>> Underflow on MUL by SIN(PHI) <<<<
;
IMUND: $LCALL IPU ;Imaginary part underflow
JRST STIMP ;Rejoin main flow above.
;
; >>>> End of underflow on MUL by SIN(PHI) <<<<
;
; >>>> Underflow on MUL by COS(PHI) <<<<
;
RLUND: $LCALL RPU ;Real part underflow
JRST STRLP ;Rejoin main flow above.
;
; >>>> End of underflow on MUL by COS(PHI) <<<<
;
; >>>> Underflow on EXP(ALPHA) <<<<
;
;Multiplication by SIN(PHI) and COS(PHI) can only aggravate
;the underflow of EXP(ALPHA). Hence both real and imaginary
;parts underflow.
;
;NOTE: Other cases for underflow of EXP(ALPHA) from
;BLOCK 4 join here.
;
EXPUND: ;Write a message saying that EXP(ALPHA) underflows
;and hence both parts of answer also underflow.
$LCALL BPU ;Both parts underflow
SETZ T0, ;Store zero for real part.
SETZ T1, ;Store zero for imaginary part.
JRST RESTOR ;Go to end of BLOCK 6 and
; restore registers 2 through 5.
;
; >>>> End of underflow on EXP(ALPHA) <<<<
;
; >>>> Overflow on EXP(ALPHA) <<<<
;
;On entry to the following code, SPHI and CPHI contain
;SIN and COS of PHI. Multiplication by SIN and/or COS
;might compensate the overflow expected on EXP(ALPHA).
;Hence get LN(SIN(|PHI|)) and LN(|COS(PHI)|) and set the
;sign flags for the real and imaginary parts of the
;results. Then merge at ANS2 with the case for underflow
;of PHI (PHIFLG = 1), where this has already been done, to
;complete the calculation.
;
ANS1: SKIPL T2,SPHI ;Get SPHI, and if >/= 0
JRST PLUSI ; no changes needed.
MOVSI T5,400000 ;Set sign bit of T5
MOVEM ISIGN ; and copy to ISIGN.
DMOVN T2,T2 ;Get |SIN(PHI)|
PLUSI: DMOVEM LNARG ; and copy to memory.
FUNCT DLOG.,<LNARG> ;Get LN(|SIN(PHI)|)
DMOVEM T0,SPHI ; and store in SPHI.
;
SKIPL T2,CPHI ;Get CPHI, and if >/= 0
JRST PLUSR ; no changes needed.
MOVSI T5,400000 ;Else set sign bit of T5
MOVEM RSIGN ; and copy to RSIGN.
DMOVN CPHI ;Get |COS(PHI)|)
PLUSR: DMOVEM LNARG ; and copy to memory.
FUNCT DLOG.,<LNARG> ;Get LN(|COS(PHI)|)
DMOVEM T0,CPHI ; and store in CPHI.
;
; >>>> End of overflow on EXP(ALPHA) <<<<
;
; >>>> Code for PHIFLG negative <<<<
;
;The code for ANS1 merges here. SPHI and CPHI contain,
;respectively, LN(SIN(|PHI|)) and LN(|COS(PHI)|). ISIGN
;and RSIGN contain the signs for the imaginary and real
;parts, and ALPHA contains ALPHA, which has already had
;the case for underflow of EXP(ALPHA) separated out. No
;test for overflow has yet been made. However, the
;magnitudes of the real and imaginary parts will be
;calculated as
;
; EXP(ALPHA + SPHI) and EXP(ALPHA + CPHI)
;
;respectively. Tests for overflow on EXP will be postponed
;until after ALPHA + SPHI and ALPHA + CPHI are calculated.
;
; >>>> Calculate EXP(ALPHA + SPHI) <<<<
;
ANS2: DMOVE T0,ALPHA ;Get ALPHA and make
DMOVE T2,ALPHA ; a copy.
DFAD T0,SPHI ;Add in SPHI. On exception
JFCL IMEX ; branch to special case.
CAMGE T0,AMAXHI ;If T0 is not too big,
JRST NOVIM ; no overflow on imag. part.
CAMG T0,AMAXHI ;If T0 > AMAXHI, overflow
; on imaginary part.
CAML T1,AMAXLO ;T0 = AMAXHI. Test LO's.
JRST IMOV ; Overflow for T1 >/= AMAXLO.
NOVIM: CAMLE T0,AMINHI ;Next test for EXP underflow:
JRST IMOK ; If T0 > AMINHI, no underflow.
CAML T0,AMINHI ;If T0 </= AMINHI, EXP will
; underflow, and also imag. part.
CAMG T1,AMINLO ;T0 = AMINHI. If T1 </= AMINLO
JRST IMUND1 ; imaginary part underflows.
IMOK: DMOVE T0,EXPARG ;Copy argument for EXP
FUNCT DEXP.,<EXPARG> ; and invoke EXP function.
STIMP1: SKIPGE T0,ISIGN ;No exceptions. If ISIGN < 0
DMOVN T0,T0 ; so is imaginary part.
DMOVEM T0, IMAG ;Store imaginary part.
JRST ANS3 ;Now get real part
;
; >>>> End of main flow for EXP(ALPHA + SPHI) <<<<
;
; >>>> Exceptions on EXP(ALPHA + SPHI) <<<<
;
IMEX: JUMPE T0,IMAG1 ;Underflow on sum: |imag. part| = 1.
JUMPL T2,IMUND1 ;On overflow ALPHA and SPHI have
; same signs. EXP underflows
; if they are negative.
IMOV: $LCALL IPO ;Imaginary part overflow
HRLOI T0,377777 ;Store zero extended
SETZ T1, ; biggest number.
JRST STIMP1 ;Return to main flow in ANS2 to
; get correct sign, and store.
IMUND1: $LCALL IPU ;Imaginary part underflow
SETZ T0, ;Store double precision
SETZ T1, ; zero for imaginary part
DMOVEM T0,IMAG ; and return to main flow
JRST ANS3 ; in ANS2 to get real part.
IMAG1: MOVE T0,SPONE ;Store double precision
SETZ T1, ; unity for imaginary part
JRST STIMP1 ;Return to main flow in ANS2 to
; get correct sign, and store.
;
; >>>> End of exceptions on imaginary part <<<<
;
; >>>> Continue main flow: Get real part <<<<
;
ANS3: DMOVE T0,ALPHA ;Get ALPHA and make
DMOVE T2,ALPHA ; a copy.
DFAD T0,CPHI ;Add in CPHI. On exception
JFCL RLEX ; branch to special case.
CAMGE T0,AMAXHI ;If T0 is not too big, no
JRST NOVRL ; overflow on real part.
CAMG T0,AMAXHI ;If T0 > AMAXHI, overflow
; on real part.
CAML T1,AMAXLO ;T0 = AMAXHI. Test LO's.
JRST RLOV ; Overflow for T0 >/= AMAXLO
NOVRL: CAMLE T0,AMINHI ;Next test for EXP underflow:
JRST RLOK ; If T0 > AMINHI, no underflow.
CAML T0,AMINHI ;If T0 < AMINHI, EXP will
; underflow, and also real part.
CAMG T1,AMINLO ;T0 = AMINHI. If T1 </= AMINLO
JRST RLUND1 ; real part underflows.
RLOK: DMOVE T0,EXPARG ;Copy argument for EXP
FUNCT DEXP.,<EXPARG> ; and invoke EXP function.
STRLP1: SKIPGE T0,RSIGN ;No exceptions. If RSIGN < 0
DMOVN T0,T0 ; so is real part.
DMOVEM T0, REAL ;Store real part.
JRST DPTOSP ;Go to BLOCK 6 to convert
; to single precision.
;
; >>>> Exceptions on EXP(ALPHA + CPHI) <<<<
;
RLEX: JUMPE T0,REAL1 ;Underflow on sum: |real part| = 1.
JUMPL T2,RLUND1 ;On overflow ALPHA and CPHI have
; same signs. EXP underflows
; if they are negative.
RLOV: $LCALL RPO ;Real part overflow
HRLOI T0,377777 ;Store zero extended
SETZ T1, ; biggest number.
JRST STRLP1 ;Return to main flow in ANS3 to
; get correct sign, and store.
;
RLUND1: $LCALL RPU ;Real part underflow
SETZ T0, ;Store double precision
SETZ T1, ; zero for real part
DMOVEM T0,REAL ; and continue with main
JRST DPTOSP ; flow in BLOCK 6 to convert
; to single precision.
;
REAL1: MOVE T0,SPONE ;Store double precision
SETZ T1, ; unity for real part
JRST STRLP1 ;Return to main flow in ANS3 to
; get correct sign, and store.
;
; >>>> End of exceptions on real part <<<<
;
;End of BLOCK 5. All exceptions in the double precision
;calculations have been properly handled. All paths lead to
;DPTOSP, in BLOCK 6, for conversion to single precision.
;
; *****************
; * BLOCK 6 *
; *****************
;
DPTOSP: DMOVE T0,REAL ;Get D.P. real part to T0.
JUMPL T0,RNEG ;Branch if real part < 0.
TLNE T1,(1B1) ;Is rounding bit 1?
TRON T0,1 ;Yes. If possible round by
JRST IMCVT ; setting LSB. Go convert
; imaginary part.
MOVE T1,T0 ;Else, get copy of high word.
AND T0,[777000,,1] ;Construct unnormalized LSB
; with same exponent.
FAD T0,T1 ;Round DP to SP
JFCL RCVTX1 ;Branch on exception.
JRST IMCVT ;Go convert imaginary part.
;
RNEG: DMOVN T0,T0 ;REAL < 0; get magnitude.
TLNE T1,(1B1) ;Is rounding bit 1?
TRON T0,1 ;Yes. If possible round by
; setting LSB. Go restore
JRST MINUS1 ; minus sign.
MOVE T1,T0 ;Else, get copy of high word.
AND T0,[777000,,1] ;Construct unnormalized LSB
FAD T0,T1 ;Round DP to SP
JFCL RCVTX2 ;Branch on exception.
;
MINUS1: MOVN T0,T0 ;Restore minus sign.
JRST IMCVT ;Now convert imaginary part.
;
; >>>> Exceptions on conversion of real part. <<<<
;
RCVTX1: JUMPE T0,RCVTUN ;Branch on underflow.
; Write message that real part overflows on conversion
$LCALL RPO ;Real part overflow
HRLOI T0,377777 ;Store positive biggest for real.
JRST IMCVT ; and go convert imaginary part.
;
RCVTX2: JUMPE T0,RCVTUN ;Branch on underflow.
; Write message that real part overflows on conversion.
$LCALL RPO ;Real part overflow
HRLOI T0,400000 ;Store negative biggest
TRO T0,(1B1) ; for real part.
JRST IMCVT ;Now convert imaginary part.
;
RCVTUN: SETZ T0, ;Zero for underflow of real part.
; Write message that real part underflows on conversion.
$LCALL RPU ;Real part underflow
JRST IMCVT ;Now convert imaginary part.
;
; *****************
; * BLOCK 6.1 *
; *****************
;
IMCVT: DMOVE T1,IMAG ;Get D.P. imaginary part to T1.
JUMPL T1,INEG ;Branch if imaginary part < 0.
TLNE T2,(1B1) ;Is rounding bit 1?
TRON T1,1 ;Yes. If possible round by
JRST RESTOR ; setting LSB. Go restore
; registers and exit.
MOVE T2,T1 ;Else, get copy of high word.
AND T1,[777000,,1] ;Construct unnormalized LSB
; with same exponent.
FAD T1,T2 ;Round DP to SP
JFCL ICVTX1 ;Branch on exception.
JRST RESTOR ;Go restore registers and exit.
;
INEG: DMOVN T1,T1 ;IMAG < 0; get magnitude.
TLNE T2,(1B1) ;Is rounding bit 1?
TRON T1,1 ;Yes. If possible round by
; setting LSB. Go restore
JRST MINUS2 ; minus sign.
MOVE T2,T1 ;Else, get copy of high word.
AND T1,[777000,,1] ;Construct unnormalized LSB
FAD T1,T2 ;Round DP to SP
JFCL ICVTX2 ;Branch on exception.
;
MINUS2: MOVN T1,T1 ;Restore minus sign.
JRST RESTOR ;Go restore registers and exit.
;
; >>>> Exceptions on conversion of imaginary part. <<<<
;
ICVTX1: JUMPE T1,ICVTUN ;Branch on underflow.
; Write message that imaginary part overflows on conversion
$LCALL IPO ;Imaginary part overflow
HRLOI T1,377777 ;Store positive biggest for
JRST RESTOR ; imaginary part, go restore
; registers, and exit.
;
ICVTX2: JUMPE T1,ICVTUN ;Branch on underflow.
; Write message that imaginary part overflows on conversion.
$LCALL IPO ;Imaginary part overflow
HRLOI T1,400000 ;Store negative biggest
TRO T1,(1B1) ; for imaginary part.
JRST RESTOR ;Go restore registers and exit.
;
ICVTUN: SETZ T1, ;Zero for underflow of
; imaginary part.
; Write message that imaginary part underflows on conversion.
$LCALL IPU ;Imaginary part underflow
RESTOR: MOVE T2,SAVT2 ;Restore registers T2,
MOVE T3,SAVT3 ; T3,
MOVE T4,SAVT4 ; T4,
MOVE T5,SAVT5 ; and T5.
;
GOODBY (1) ;Return.
;Constants
PI: EXP 202622077325,021026430215 ;pi
PIOV2: EXP 201622077325,021026430215 ;pi/2
LN300: EXP 210412126217,316725563320 ;192 * ln(2) for LOG of
; underflowed PHI
SPONE: EXP 201400000000 ;1.0
ARGHI: EXP 241622077325 ;Double precision threshold
ARGLO: EXP 020000000000 ; for argument reduction in
; SIN/COS routines.
TWO63: EXP 300400000000 ;2**63 Threshold for |THETA| = PI/2
TWOM63: EXP 102400000000 ;2**-63 Threshold for |THETA| = PI
; or THETA = Y/X
TWOM66: EXP 077400000000 ;2**-66 Threshold for underflow
; negligible.
AMAXHI: EXP 207540074636 ;88. ... double precision
AMAXLO: EXP 077042573256 ; overflow boundary for EXP.
AMINHI: EXP 570232254036 ;-88. ... double precision
AMINLO: EXP 301750634730 ; underflow boundary for EXP.
;Data
SEGMENT DATA
LNRHO: BLOCK 2 ;LNRHO and THETA must be
THETA: BLOCK 2 ; contiguous since they
; provide output for CDLOG.
REXP: BLOCK 2 ;Temporary storage for A
IEXP: BLOCK 2 ;Temporary storage for B
CPHI: BLOCK 2 ;For COS(PHI) or LN(COS(|PHI|))
SPHI: BLOCK 2 ;For SIN(PHI) or LN(SIN(|PHI|))
ALPHA: BLOCK 2 ;Storage for ALPHA.
ARG1: BLOCK 4 ;Temporary storage for inputs to CDLOG.
LNARG: BLOCK 2 ;Temporary storage for inputs
TRGARG: BLOCK 2 ; to DLOG, DSIN, and
EXPARG: BLOCK 2 ; DCOS and DEXP routines.
REAL: BLOCK 2 ;Real part of result
IMAG: BLOCK 2 ;Imaginary part of result
RSIGN: BLOCK 1 ;Storage for signs of real and
ISIGN: BLOCK 1 ; imaginary parts of result.
PHIFLG: BLOCK 1 ;Storage for flag indicating
; underflow of PHI.
XCHFLG: BLOCK 1 ;Storage for flag indicating
; exchange of LNRHO,THETA and A,B.
SAVT2: BLOCK 1 ;SAVT2, SAVT3, SAVT4,
SAVT3: BLOCK 1 ; must be contiguous. The
SAVT4: BLOCK 1 ; corresponding registers are
SAVT5: BLOCK 1 ; restored by DMOVEs.
PRGEND
TITLE CLOG COMPLEX NATURAL LOG FUNCTION
; (SINGLE PRECISION FLOATING POINT)
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CLOG
EXTERN CLOG.
CLOG=CLOG.
PRGEND
TITLE CLOG. COMPLEX NATURAL LOG FUNCTION
; (SINGLE PRECISION FLOATING POINT)
SUBTTL Mary Payne 16-Oct-81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;CLOG(Z) IS CALCULATED AS FOLLOWS
; LET Z = X + I*Y
; IF Y IS NOT ZERO, THEN
; CLOG(Z) = U + I*V
; U = (1/2)*ALOG(X**2 + Y**2)
; V = ATAN2(Y,X)
; IF X IS NOT ZERO AND Y IS ZERO, THEN
; IF X IS POSITIVE, CLOG(Z) = ALOG(X) + I*0
; IF X IS NEGATIVE, CLOG(Z) = ALOG(ABS(X)) + I*PI
; IF X AND Y ARE ZERO, THEN
; CLOG(Z) = -MACHINE INFINITY + I*0
; AND AN ERROR MESSAGE IS RETURNED
; THE APPROXIMATION USED FOR ALOG IS FORMULA 2662 FROM
; "COMPUTER APPROXIMATIONS" BY HART ET AL.
;THE REAL PART OF THE RESULT IS RETURNED IN REGISTER T0.
;THE IMAGINARY PART OF THE RESULT IS RETURNED IN REGISTER T1.
;REQUIRED (CALLED) ROUTINES: ALOG, ATAN, AND ATAN2
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; XMOVEI L,ARG
; PUSHJ P,CLOG
;THE RESULT IS RETURNED IN ACCUMULATORS T0 AND T1
EXTERN SNG.0
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (CLOG,.) ;ENTRY TO CLOG ROUTINE
DMOVE T0,@(L) ;REAL PART OF ARGUMENT TO T0
; IMAGINARY PART TO T1
JUMPE T1,YZRO ;Y = 0 IS A SPECIAL CASE
JUMPE T0,XZRO ;X = 0 IS A SPECIAL CASE
PUSH P,T2 ;SAVE REGISTERS T2
PUSH P,T3 ; AND T3
MOVM T2,T0 ;GET ABSOLUTE VALUES OF
MOVM T3,T1 ; X AND Y
AND T2,MASK ;ISOLATE EXPONENT FIELDS
AND T3,MASK ; OF X AND Y
SUB T3,T2 ;EXPONENT OF Y - EXPONENT OF X
CAML T3,P14 ;IF DIFFERENCE > 14
JRST BIG ; |Y/X| IS LARGE.
CAMGE T3,M14 ;IF DIFFERENCE < -14
JRST SMALL ; |Y/X| IS SMALL
PUSH P,T4 ;SAVE REGISTERS T4
PUSH P,T5 ; AND T5
MOVM T4,T0 ;SAVE MAGNITUDES OF X
MOVM T5,T1 ; AND Y IN T4,T5
MOVEM T0,TT0 ;STORE X AND Y FOR CALL
MOVEM T1,TT2 ; TO ATAN2.
FUNCT ATAN2.,<TT2,TT0> ;NO EXCEPTIONS
MOVE T3,T0 ;SAVE IMAG. PART IN T3
CAML T4,T5 ;COMPARE |X| AND |Y|
JRST NOXCH ; AND GET LARGER = A IN T4
EXCH T4,T5 ; AND SMALLER = B IN T5
NOXCH: MOVE T2,T4 ;LARGER TO T2
LSH T2,-33 ;ISOLATE ITS BIASED EXPONENT
SUBI T2,200 ; AND UNBIAS IT
JUMPLE T2,LE ;IF UNBIASED EXPONENT > 0,
SUBI T2,1 ; DECREMENT IT
LE: MOVN T2,T2 ;NEGATE SCALE INDEX
FSC T4,(T2) ;SCALE A TO [1/2,2) AND B BY
FSC T5,(T2) ; SAME FACTOR. NO EXCEPTIONS
LSH T2,1 ;MULTIPLY NEGATED SCALE INDEX BY 2
DMOVE T0,T4 ;COPY |X| AND |Y| TO T0,T1
FMPR T0,T0 ;SQUARE LARGER. NO EXCEPTIONS
FMPR T1,T1 ;SQUARE SMALLER. NO OVERFLOW AND
JFCL ; UNDERFLOW WON'T MATTER
MOVEM T1,SMALSQ ;SAVE SMALLER SQUARED
FADR T1,T0 ;GET SCALED SQUARE OF MODULUS
CAMLE T1,RT2 ;IF T1 .GT. SQRT(2)
JRST MORE ; GO TO MORE SCALING
CAMGE T1,RT2O2 ;IF T1 .LT. SQRT(1/2)
JRST MORE ; GO TO MORE SCALING
JUMPN T2,JOIN ;IF SCALE INDEX NOT 0, GO TO JOIN
;At this point the scale index is zero, implying that (A**2 + B**2)
;is unscaled. Also (A**2 + B**2) is between SQRT(1/2) and SQRT(2),
;which implies that a damaging loss of significance can occur in the
;calculation of (A**2 + B**2 - 1), which is one factor of Z, used
;in the approximation. Growth of error due to loss of significance
;can be avoided by calculating instead (A - 1) * (A + 1) + B**2,
;which is done in the following lines of code:
MOVE T0,T4 ;COPY LARGER TO T0
FSBR T0,ONE ;GET A - 1
FADR T4,ONE ;AND A + 1
FMPR T4,T0 ;AND THEIR PRODUCT
FADR T4,SMALSQ ;T4 NOW HAS A**2 + B**2 - 1
MOVE T1,T4 ; = NUMERATOR OF Z. COPY TO T1
FADR T1,TWO ;GET DENOMINATOR OF Z IN T1
JRST MERGE ;MERGE FOR POLYNOMIAL APPROX
;The following code applies to scaled values of (A**2 + B**2) lying
;outside the interval (SQRT(1/2),SQRT(2)); it carries out more
;scaling to get the scaled (A**2 + B**2) inside this interval.
MORE: MOVE T4,T1 ;MAKE COPY OF A**2 + B**2
LSH T4,-33 ;ISOLATE ITS BIASED EXPONENT
SUBI T4,200 ;UNBIAS THE EXPONENT
MOVN T4,T4 ; AND NEGATE IT
FSC T1,(T4) ;SCALED MODULUS NOW IN [1/2,1)
ADD T2,T4 ;TOTAL NEGATED SCALE INDEX TO T2
CAML T1,RT2O2 ;IF T1 .GT. SQRT(1/2)
JRST JOIN ; SCALING IS DONE
FSC T1,1 ;SCALE T1 TO [1,SQRT(2)]
ADDI T2,1 ; AND INCREMENT NEGATED INDEX
;At this point the scaled (A**2 + B**2) is in the interval
;[SQRT(1/2),SQRT(2)]. T2 contains the negative of the index
;necessary to compensate later for the scaling. In the following
;lines of code, a polynomial approximation is used to evaluate
;the natural log of the scaled (A**2 + B**2).
JOIN: MOVN T2,T2 ;RESTORE SIGN OF SCALE INDEX
MOVE T4,T1 ;GET COPY OF SCALED (A**2 + B**2)
FSBR T4,ONE ; -1 = NUMERATOR OF Z
FADR T1,ONE ;GET DENOMINATOR OF Z IN T1
;The following code is taken from the AlOG routine. The constants
;L3, L4, and L5 are those used in ALOG.
MERGE: FDVR T4,T1 ;Z = ZNUM / ZDEN. NO EXCEPTIONS
MOVE T5,T4 ;SAVE A COPY OF Z
FMPR T4,T4 ;W = Z**2
MOVE T0,T4 ; AND MAKE COPY
FMPR T0,L3 ;FORM P(Z) = W*L3
FADR T0,L4 ; + L4
FMPR T0,T4 ; *W
FADR T0,L5 ; + L5
FMPR T0,T4 ; *W
FMPR T0,T5 ; *Z. NO EXCEPTIONS
FSC T5,1 ;GET 2*Z
;2*Z still needs to be added into P(Z).
JUMPN T2,REST ;IF SCALE INDEX = 0,
FADR T0,T5 ; COMPLETE P(Z) BY ADDING IN Z
JRST FINISH ;GO FINISH UP
;The following lines of code complete the calculation of the real
;part for the index not zero. The real part is the double
;precision sum of INDEX*LN(2), Z, and P(Z), rounded to single
;precision.
REST: MOVE T4,T5 ;GET 2*Z IN DOUBLE
SETZB T1,T5 ; PRECISION
DFAD T0,T4 ;GET DP SUM. FLOAT SCALE
FLTR T4,T2 ;INDEX -- DP SINCE T5 HAS 0
DFMP T4,LN2 ;GET INDEX*LN(2)
DFAD T0,T4 ;AND ADD TO PREVIOUS SUM
PUSHJ P,SNG.0 ;ROUND TO SINGLE PRECISION
FINISH: FSC T0,-1 ;DIVIDE BY 2 TO GET LN(MODULUS)
MOVE T1,T3 ;IMAG PART TO T1
POP P,T5 ;POP REGISTERS
POP P,T4 ; IN REVERSE
POP P,T3 ; ORDER IN WHICH
POP P,T2 ; THEY WERE PUSHED
GOODBY (1)
BIG: DMOVE T2,T0 ;SAVE X AND Y
MOVMM T3,TT0 ;STORE |Y| FOR ALOG
FUNCT ALOG.,<TT0>
MOVE T1,PIO2 ;IMAGINARY PART NEAR +/- PI/2
JUMPGE T3,PO2POS ; + FOR Y POSITIVE
MOVN T1,T1 ; - FOR Y NEGATIVE
PO2POS: FDVR T2,T3 ;GET X/Y; NO OVERFLOW
JFCL UND1 ; UNDERFLOW IS POSSIBLE
FSBR T1,T2 ;SMALL CORRECTION TO IMAG
FMPR T2,T2 ;GET (X/Y)**2. NO OVERFLOW
JFCL UND1 ; UNDERFLOW IS POSSIBLE
FSC T2,-1 ;(1/2)*(X/Y)**2. NO OVERFLOW
JFCL UND1 ; UNDERFLOW IS POSSIBLE
FADR T0,T2 ;+ LN(|Y|) = REAL. NO EXCEPTIONS
JRST REST23 ;GO RESTORE T2 AND T3
SMALL: DMOVE T2,T0 ;SAVE X AND Y
MOVMM T2,TT0 ;STORE |X| FOR ALOG
FUNCT ALOG.,<TT0>
SETZ T1, ;IMAG. PART NEAR 0 OR +/- PI
JUMPGE T2,IMAGOK ; 0 FOR X POSITIVE
MOVE T1,CPI ; +/- PI FOR X NEGATIVE
JUMPGE T3,IMAGOK ; + PI FOR Y POSITIVE
MOVN T1,T1 ; - PI FOR Y NEGATIVE
IMAGOK: FDVR T3,T2 ;GET Y/X; NO OVERFLOW
JFCL UND2 ; UNDERFLOW IS POSSIBLE
FADR T1,T3 ;SMALL CORRECTION FOR IMAG.
FMPR T3,T3 ;GET (X/Y)**2. NO OVERFLOW
JFCL UND1 ; UNDERFLOW IS POSSIBLE
FSC T3,-1 ;(1/2)*(X/Y)**2. NO OVERFLOW
JRST UND1 ; UNDERFLOW IS POSSIBLE
FADR T0,T3 ;+ LN(|X|) = REAL. NO EXCEPTIONS
JRST REST23 ;GO RESTORE REGISTERS T2 AND T3
UND2: JUMPN T1,REST23 ;DONE IF |IMAG| NEAR PI OR PI/2
$LCALL IPU
;LERR (LIB,%,<CLOG: Imaginary part underflow>)
UND1: JUMPN T0,REST23 ;DONE IF LN(|X| OR |Y|) NOT 0
$LCALL RPU
;LERR (LIB,%,<CLOG: Real part underflow>)
REST23: POP P,T3 ;RESTORE REGISTERS
POP P,T2 ; T2 AND T3
GOODBY (1)
YZRO: JUMPE T0,XYZRO ;X = Y = 0 IS A SPECIAL CASE
JUMPGE T0,RPOS ;Y IS 0, X ISN'T. IF X < 0
MOVN T0,T0 ;X=ABS(X)
MOVE T1,CPI ;V=PI
RPOS: MOVEM T1,TT2 ;SAVE IMAGINARY PART
MOVEM T0,TT0 ;MOVE X TO MEMORY
FUNCT ALOG.,<TT0> ;COMPUTE U=DLOG(X)
MOVE T1,TT2 ;RESTORE IMAGINARY PART
GOODBY (1) ;RETURN; RESULT IN T0,T1
XZRO: MOVE T0,T1 ;X IS 0, Y ISN'T. COPY Y TO T0
MOVE T1,PIO2 ;IMAG PART IS +/- PI/2
JUMPGE T0,YPLUS ;IF Y < 0,
MOVN T0,T0 ; MAKE IT POSITIVE, AND
MOVN T1,T1 ; MAKE IMAG PART NEGATIVE
YPLUS: MOVEM T0,TT0 ;PREPARE TO GET LOG
MOVEM T1,TT2 ;SAVE IMAGINARY PART
FUNCT ALOG.,<TT0> ;GET LOG(|Y|), NO EXCEPTIONS
MOVE T1,TT2 ;RESTORE IMAGINARY PART.
GOODBY (1) ;RETURN: RESULTS IN T0, T1
XYZRO: $LCALL ZIZ
;LERR (LIB,%,<entry CLOG; arg is (0.0,0.0); result=(-infinity,0.0)>)
MOVE T0,[400000000001] ;REAL ANSWER IS NEGATIVE INFINITY
GOODBY (1) ;RETURN; RESULT IN T0,T1
;CONSTANTS
MASK: EXP 377000000000 ;MASK FOR EXPONENT FIELD
P14: EXP 016000000000 ;HIGH 9 BITS = 16 OCTAL
M14: EXP 762000000000 ;HIGH 9 BITS = -16 OCTAL
ONE: EXP 201400000000 ;1.0
TWO: EXP 202400000000 ;2.0
CPI: EXP 202622077325 ;PI
PIO2: EXP 201622077325 ;PI / 2
RT2O2: EXP 200552023632 ;SQRT(2) / 2
RT2: EXP 201552023632 ;SQRT(2)
L3: EXP 177464164321 ;.301003281
L4: EXP 177631177674 ;.39965794919
L5: EXP 200525253320 ;.666669484507
LN2: EXP 200542710277,276434757153 ;LN(2)
;DATA
SEGMENT DATA
TT0: BLOCK 1 ;FOR ATAN2 AND ALOG CALLS
TT2: BLOCK 1 ;FOR ATAN2 CALL
SMALSQ: BLOCK 1
PRGEND
TITLE CSIN COMPLEX SINE FUNCTION
; (SINGLE PRECISION)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CSIN
EXTERN CSIN.
CSIN=CSIN.
PRGEND
TITLE CCOS COMPLEX COSINE FUNCTION
; (SINGLE PRECISION)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CCOS
EXTERN CCOS.
CCOS=CCOS.
PRGEND
TITLE CSIN. COMPLEX SINE AND COSINE ROUTINES
; (SINGLE PRECISION)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;CSIN(Z) AND CCOS(Z), WHERE Z = X + I*Y, ARE CALCULATED AS FOLLOWS
; CSIN(Z) = SIN(X)*COSH(Y) + I*COS(X)*SINH(Y)
; CCOS(Z) = COS(X)*COSH(Y) - I*SIN(X)*SINH(Y)
; THE FUNCTIONS SINH AND COSH ARE CODED IN LINE.
; THE REAL PART OF THE RESULT CANNOT UNDERFLOW BECAUSE
; COSH(Y) IS NEVER LESS THAN 1.
; THE IMAGINARY PART OF THE RESULT MAY UNDERFLOW; IF THIS HAPPENS THE
; IMAGINARY PART IS SET TO 0 AND A MESSAGE IS PROVIDED.
;THE RANGE OF DEFINITION FOR CSIN AND CCOS IS AS FOLLOWS
;FOR
; Z = X+I*Y. IF ABS(X) > 823549.66 THE RESULT IS SET TO 0.0 AND
; AN ERROR MESSAGE IS RETURNED.
;FOR ABS(Y) > 88.029692 CALCULATIONS PROCEED AS FOLLOWS:
; FOR THE REAL PART OF THE RESULT
; LET T = ABS(SIN(X)), THEN
; FOR T = 0.0 THE REAL PART OF THE RESULT IS SET TO 0.0
; FOR LOG(T) + ABS(Y) > 88.722839 THE REAL PART OF THE
; RESULT IS SET TO PLUS OR MINUS MACHINE INFINITY AND
; AN ERROR MESSAGE IS RETURNED. (88.722839=88.029692+LN2)
; FOR THE IMAGINARY PART OF THE RESULT
; LET T = ABS(COS(X)), THEN
; FOR T = 0.0 THE IMAGINARY PART OF THE RESULT IS SET TO 0.0
; FOR LOG(T) + ABS(Y) > 88.722839 THE IMAGINARY PART OF THE
; RESULT IS SET TO PLUS OR MINUS MACHINE INFINITY AND
; AN ERROR MESSAGE IS RETURNED. (88.722839=88.029692+LN2)
;REQUIRED (CALLED) ROUTINES: SIN,COS,EXP,ALOG
;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; MOVEI L,ARG
; PUSHJ P,CSIN
; OR
; PUSHJ P,CCOS
;THE COMPLEX ANSWER IS RETURNED IN ACCUMULATORS T0 AND T1.
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (CCOS,.) ;ENTRY TO CCOS ROUTINE
PUSH P,T3
MOVEI T3,1 ;SET FLAG TO 1 FOR CCOS
JRST CSIN1 ;GO TO CSIN ROUTINE
HELLO (CSIN,.) ;ENTRY TO CMPLX SINE ROUTINE
PUSH P,T3
SETZI T3, ;SET FLAG TO 0 FOR CSIN
CSIN1: PUSH P,T2
PUSH P,T4
PUSH P,T5
XMOVEI T1,@(L) ;GET ADDRESS OF ARG
MOVE T0,(T1) ;X = REAL(Z)
MOVM T2,T0 ; GET ABS(X)
CAMG T2,CUT ;IF ABS(X) IS NOT TOO LARGE
JRST XOK ; GO TO GET Y
$LCALL ARZ
;LERR (LIB,%,<CSIN or CCOS: ABS(REAL(arg)) too large; result = (0.0,0.0)>)
SETZI T0, ;SET REAL(RESULT) TO 0.0
SETZI T1, ;SET IMAG(RESULT) TO 0.0
JRST RET ;RETURN
XOK: MOVE T1,1(T1) ;Y = IMAG(Z)
MOVEM T1,YSAVE ;SAVE Y
MOVEM T0,ARGSAV
FUNCT SIN.,<ARGSAV> ;CALCULATE SIN(X)
MOVEM T0,SX ;SX = SIN(X)
FUNCT COS.,<ARGSAV> ;CALCULATE COS(X)
MOVEM T0,CX ;CX = COS(X)
SKIPN T3 ;IF THIS IS THE CSIN ROUTINE
JRST NOSWT ;THEN GO TO NOSWT
MOVN T2,SX ;OTHERWISE, PUT -SIN(X)
EXCH T2,CX ;IN CX, AND COS(X)
MOVEM T2,SX ;IN SX.
NOSWT: SKIPN T1,YSAVE ;SKIP UNLESS Z IS REAL, IN
JRST CSIN2 ;WHICH CASE, GO TO THE SIMPLE CASE
MOVM T2,T1 ;AY = ABS(YSAVE)
CAMLE T2,XMAX ;IF AY IS TOO LARGE FOR SIMPLE
JRST OVFL ;SINH OR COSH, GO TO OVFL.
CAMG T2,TWOM14 ;IF |Y| <= 2**(-14),
JRST EASY ; SINH = Y, COSH = 1
MOVEM T2,ARGSAV
FUNCT EXP.,<ARGSAV> ;CALCULATE EXP(AY)
MOVE T3,T0
MOVE T5,ONE
FDVR T5,T3
CAMLE T2,ONE ;IF AY IS GREATER THAN 1
JRST GTET ;THEN GO TO GTET
MOVE T1,T2 ;GET COPY OF |Y|
FMPR T1,T1 ;SS = AY**2
MOVE T4,T1
FMPR T4,R4 ;SS = SS*R4
FADR T4,R3 ;SS = SS+R3
FMPR T4,T1 ;SS = SS*(AY**2)
FADR T4,R2 ;SS = SS+R2
FMPR T4,T1 ;SS=SS*(AY**2)
FADR T4,R1 ;SS = SS+R1
FMPR T1,T4 ;SS = SS*(AY**2)
FMPR T1,T2 ;SS = SS*AY
FADR T1,T2 ;SS = SS+AY
JRST CCCC ;GO TO CCCC
GTET: MOVN T1,T5
FADR T1,T3 ;SS = T+SS
FSC T1,-1 ;SS = SS/2.
CCCC: MOVE T0,T5
FADR T0,T1 ;CC = SS+(1/T)
MOVE T3,YSAVE ;RESTORE Y
CAIGE T3,0 ;IF Y IS LESS THAN 0.0
MOVN T1,T1 ;THEN NEGATE SS
FMPR T0,SX ;GET REAL PART OF RESULT
FMPR T1,CX ;GET IMAG PART OF RESULT
JFCL IMUND
JRST RET ;RETURN
EASY: MOVE T0,SX ;REAL PART = 1 * SX
FMPR T1,CX ;IMAG PART = Y * CX
JFCL IMUND ;UNDERFLOW POSSIBLE
JRST RET ;RETURN
IMUND: $LCALL IPO,RET
;LERR (LIB,%,<CSIN or CCOS: imaginary part overflow>,RET)
OVFL: MOVM T0,SX ;T=ABS(SX)
CAIN T0, ;IF T IS EQUAL TO 0.
JRST GOTR ;THEN GO TO GOTR
MOVEM T0,ARGSAV ;OTHERWISE,
FUNCT ALOG.,<ARGSAV> ;CALCULATE LOG(T)
FADR T0,T2 ;T = LOG(T)+AY
FADR T0,LN2 ;T = T-LOG(2.0)
CAMG T0,XMAX ;IF T IS .LE. XMAX
JRST CONTIN ;THEN GO TO CONTIN
$LCALL RPO
;LERR (LIB,%,<CSIN or CCOS: real part overflow>)
HRLOI T0,377777 ;REAL PART OF RESULT SET TO INFINITY
JRST SKP2 ;SKIP NEXT 2 INSTRUCTIONS
CONTIN: MOVEM T0,ARGSAV
FUNCT EXP.,<ARGSAV> ;RRES = EXP(T)
SKP2: SKIPGE SX ;IF SX IS LESS THAN 0.
MOVN T0,T0 ;THEN NEGATE RRES
GOTR: SKIPN T1,CX ;IF CX = 0,
JRST RET ;THEN RETURN
MOVEM T0,SX ;OTHERWISE SAVE RRES
MOVMM T1,ARGSAV ;T = ABS(CX)
FUNCT ALOG.,<ARGSAV> ;CALCULATE T=LOG(T)
MOVE T1,T0
FADR T1,T2 ;T = T+AY
FADR T1,LN2 ;T = T-LOG(2)
CAMG T1,XMAX ;IF T IS .LE. XMAX
JRST CONTN2 ; THEN GO TO CONTN2
$LCALL IPO
;LERR (LIB,%,<CSIN or CCOS: imaginary part overflow>)
HRLOI T1,377777 ;SET IRES TO INFINITY
JRST SKP3 ;SKIP NEXT 3 INSTRUCTIONS
CONTN2: MOVEM T1,ARGSAV
FUNCT EXP.,<ARGSAV> ;IRES = EXP(T)
MOVE T1,T0
SKP3: MOVE T0,SX
MOVE T2,YSAVE
XOR T2,CX ;T2 HAS THE SIGN OF CX*Y
JUMPGE T2,RET ;IF SIGNS ARE THEN SAME, DONE
MOVN T1,T1 ;OTHERWISE NEGATE IMAG PART OF RESULT
JRST RET ;RETURN
CSIN2: MOVE T0,SX ;SIMPLE CASE, Z IS REAL
RET: POP P,T5
POP P,T4
POP P,T2
POP P,T3
GOODBY (1)
TWOM14: 163400000000 ;2**(-14)
LN2: 577235067500 ;-(NATURAL LOG OF 2.0)
XMAX: 207540074636 ;88.029692
R1: 176525252524 ;1.666666643E-1
R2: 172421042352 ;8.333352593E-3
R3: 164637771324 ;1.983581245E-4
R4: 156572227373 ;2.818523951E-6
ONE: 201400000000 ;1.0
CUT: 234622077325 ;2**26 * PI
SEGMENT DATA
CX: 0
SX: 0
YSAVE: 0
ARGSAV: 0
PRGEND
TITLE CSQRT COMPLEX SQUARE ROOT FUNCTION
; (SINGLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CSQRT
EXTERN CSQRT.
CSQRT=CSQRT.
PRGEND
TITLE CSQRT. COMPLEX SQUARE ROOT FUNCTION
; (SINGLE PRECISION FLOATING POINT)
SUBTTL IMSL, INC. JUNE 1, 1979
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
; LET Z = X + I*Y
; THEN THE ANSWER CSQRT(Z) = U + I*V IS DEFINED AS FOLLOWS
; IF X AND Y ARE BOTH >= 0.0, THEN
; U = SQRT((ABS(X)+CABS(Z))/2.0)
; V = Y/(2.0*U)
; IF X >= 0.0 AND Y < 0.0, THEN
; U = SQRT((ABS(X)+CABS(Z))/2.0)
; V = Y/(2.0*U)
; IF X < 0.0 AND Y >= 0.0, THEN
; U = Y/(2.0*V)
; V = SQRT((ABS(X)+CABS(Z))/2.0)
; IF X AND Y ARE BOTH < 0.0, THEN
; U = Y/(2.0*V)
; V = -SQRT((ABS(X)+CABS(Z))/2.0)
; THE RESULT IS IN THE RIGHT HALF PLANE, THAT IS, THE POLAR ANGLE OF THE
; RESULT LIES IN THE INTERVAL [-PI/2,+PI/2].
;REQUIRED (CALLED) ROUTINES: SQRT
;REGISTERS T2, T3, T4, T5, AND G1 WERE SAVED, USED, AND RESTORED
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
; MOVEI L,ARG
; PUSHJ P,CSQRT
;THE REAL PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0
;THE IMAG PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1
SEARCH MTHPRM
SEGMENT CODE
SALL
HELLO (CSQRT,.) ;ENTRY TO COMPLEX SQUARE ROOT ROUTINE.
XMOVEI T1,@0(L) ;PICK UP ADDRESS OF ARGUMENT
MOVE T0,(T1) ;PICK UP X IN AC T0.
MOVE T1,1(T1) ;PICK UP Y IN AC T1.
PUSH P,T2 ;SAVE AC T2.
JUMPE T1,ZREAL ;JUMP IF Y=0
JUMPE T0,ZIMAG ;JUMP IF X=0
MOVM T2,T0 ;ABS(X) TO AC T2.
PUSH P,T3 ;SAVE REGISTER T3.
PUSH P,T4 ;SAVE REGISTER T4
PUSH P,T5 ;SAVE REGISTER T5
PUSH P,G1 ;SAVE REGISTER G1
MOVE T4,T0 ;SAVE X IN AC T4.
MOVE T5,T1 ;SAVE Y IN AC T5.
MOVM T3,T1 ;ABS(Y) TO AC T3.
SETZ G1, ;G1 CONTAINS 0 IF ABS(X)
CAMLE T2,T3 ;<= ABS(Y), ELSE IT CONTAINS 1. PUT
AOJA G1,DWN2 ;THE LARGER OF ABS(X),ABS(Y) IN AC T2
EXCH T2,T3 ;AND THE SMALLER IN AC T3.
DWN2: FDVR T3,T2 ;CALC S/L.
JFCL ;NO UNDERFLOW ERROR MESSAGE ALLOWED.
FMPR T3,T3 ;CALC (S/L)**2.
JFCL ;NO UNDERFLOW ERROR MESSAGE ALLOWED.
FADRI T3,201400 ;HAVE (1+(S/L)**2) IN AC T3.
MOVEM T3,TEMP
FUNCT SQRT.,<TEMP> ;[3220] CALC. THE SQRT OF
;(1+(S/L)**2)=1+EPS.
JUMPE G1,YGTX ;GO TO YGTX IF ABS(Y) > ABS(X).
XGTY: FADRI T0,201400 ;CALC. 2 + EPS.
FSC T0,-1 ;CALC. 1+EPS/2.
MOVMM T0,TEMP ;PUT IN TEMP FOR
FUNCT SQRT.,<TEMP> ;[3220] CALL TO SQRT
MOVE T3,T0 ;SAVE SQRT(1+EPS/2) IN AC T3.
MOVEM T2,TEMP
FUNCT SQRT.,<TEMP> ;[3220] CALC.
FMPR T0,T3 ;CALC. SQRT(ABS(X)*(1+EPS/2)).
JRST OUT1 ;GO TO REST OF CALC.
YGTX: CAMGE T2,[1.0] ;IF ABS(Y)>1, GO TO POSSIBLE OVERFLOW CASE.
JRST YXSMAL ;ELSE, GO TO POSSIBLE UNDERFLOW.
FSC T0,-3 ;CALC. (1+EPS)/8.
FMPR T2,T0 ;CALC. ABS(Y)*(1+EPS)/8.
MOVM T3,T4 ;ABS(X) TO AC T3.
FSC T3,-3 ;CALC. ABS(X)/8.
JFCL ;SUPPRESS UNDERFLOW ERROR MESSAGE.
FADR T3,T2 ;CALC.(ABS(X)/8)+(ABS(Y)*(1+EPS)/8).
MOVEM T3,TEMP
FUNCT SQRT.,<TEMP> ;[3220] CALC.
FSC T0,1
OUT1: MOVM T1,T5 ;ABS(Y) TO AC T1.
FDVR T1,T0 ;DIVIDE ABS(Y)/2 BY THE
JFCL UNDFL
FSC T1,-1 ;SQRT TERM.
JFCL UNDFL
JRST SIGNS ;GO TO REST OF CALC.
YXSMAL: FMPR T0,T2 ;ABS(Y)*(1+EPS) = CDABS(Z).
MOVM T3,T4 ;ABS(X) TO AC T3.
FADR T3,T0 ;GET |X| + CDABS(Z)
FSC T3,1 ;GET 2*(|X| + CABS(Z))
MOVEM T3,TEMP
FUNCT SQRT.,<TEMP> ;[3220] CALC. SQRT OF 2*(|X| + CABS(Z))
MOVE T3,T0 ;GET SQRT((|X| + CABS(Z))*2)
FSC T0,-1 ;GET SQRT((|X| + CABS(Z))/2)
FDVR T2,T3 ;GET |Y| / SQRT(2*(|X| + CABS(Z))).
MOVE T1,T2 ;PUT A TEMP ANSWER IN AC T1.
SIGNS: CAIGE T4,0 ;SET UP THE REAL AND
EXCH T1,T0 ;THE IMAGINARY ANSWERS WITH
CAIGE T5,0 ;THE APPROPRIATE
MOVN T1,T1 ;SIGNS.
POP P,G1 ;RESTORE AC G1
POP P,T5 ;RESTORE AC T5
POP P,T4 ;RESTORE AC T4
POP P,T3 ;RESTORE AC T3
POP P,T2 ;RESTORE AC T2
GOODBY (1) ;RETURN
ZREAL: JUMPE T0,DONE ;Y = 0, HOW ABOUT X?
MOVE T2,T0 ;X NOT ZERO; SAVE IT.
MOVMM T2,TEMP ;GET ABS(X) FOR SQRT.
FUNCT SQRT.,<TEMP> ;[3220] T0 GETS SQRT(ABS(X)).
SETZ T1, ;T1 smashed by call; Set to 0 again
JUMPG T2,DONE ;T0,T1 OK FOR X>0.
EXCH T0,T1 ; INTERCHANGE FOR X<0.
POP P,T2 ;RESTORE T2.
GOODBY (1) ;RETURN.
ZIMAG: MOVE T2,T1 ;X = 0, SAVE Y.
FSC T1,-1 ;DIVIDE Y BY 2.
MOVMM T1,TEMP ;GET ABS(Y/2) FOR SQRT.
FUNCT SQRT.,<TEMP> ;[3220] T0 GETS SQRT(ABS(Y/2)).
MOVE T1,T0 ;SO DOES T1.
JUMPG T2,DONE ; DONE IF Y>0.
MOVN T1,T1 ;NEGATE IMAG PART IF Y<0.
DONE: POP P,T2 ;RESTORE T2.
GOODBY (1) ;RETURN
UNDFL: CAIGE T4,0 ;SKIP IF REAL & IMAG SWITCHED
$LCALL RPU,SIGNS
; LERR (LIB,%,<CEXP: real part underflow>,SIGNS)
$LCALL IPU,SIGNS
; LERR (LIB,%,<CEXP: imaginary part underflow>,SIGNS)
SQ2: 201552023632 ;SQRT(2).
SEGMENT DATA
TEMP: 0 ;TEMPORARY STORAGE FOR SQRT ARGS
PRGEND
TITLE CFM COMPLEX MULTIPLICATION
SUBTTL KAREN KOLLING /KK/CKS 28-Oct-81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;FROM LIB40 VERSION .027
;FROM V.021 16-OCT-70
;COMPLEX FLOATING MULTIPLY ROUTINE.
;THIS ROUTINE MULTIPLIES TWO COMPLEX NUMBERS ACCORDING
;TO THE RULE:
; (A+IB)*(C+ID)=(A*C-B*D)+I*(B*C+A*D)
;IF EITHER PRODUCT IN EITHER TERM OVER OR UNDERFLOWS,
;A JUMP IS MADE TO A SPECIAL CASE CALCULATION. OVER
;AND UNDERFLOW ERROR MESSAGES ARE RETURNED ONLY IF THE
;ENTIRE TERM OVER OR UNDERFLOWS.
;THERE ARE TWO SETS OF ENTRY POINTS TO THIS ROUTINE:
;THE MULTIPLY, AND THE MULTIPLY-TO-MEMORY.
;THE CALLING SEQUENCE FOR THE MULTIPLY ROUTINES IS AS
;FOLLOWS:
; XMOVEI L,ARG2
; PUSHJ P,CFM.N
;WHERE N=0,2,4,6,10,12,14. ARG1 IS ASSUMED TO BE IN
;AC'S N AND N+1, AND THE RESULTS ARE LEFT IN AC'S N AND N+1.
;THE CALLING SEQUENCE FOR THE MULTIPLY TO MEMORY ROUTINES IS AS
;FOLLOWS:
; XMOVEI L,ARG2
; PUSHJ P,CMM.N
;WHERE N AND ARG1 ARE AS BEFORE AND THE RESULTS ARE
;LEFT IN ARG2 AND ARG2+1.
SEARCH MTHPRM
SEGMENT CODE
SALL
DEFINE SYM (N)<
ENTRY CFM.'N,CMM.'N
SIXBIT / CFM.'N/
CFM.'N: DMOVEM N,ARG11 ;STORE FIRST ARG
IFN N,< DMOVEM 0,SAVE0 > ;SAVE REGS 0-1
DMOVE 0,(L) ;GET SECOND ARG
PUSHJ P,CFM ;MULTIPLY
IFN N,< DMOVE N,0 ;STORE RESULT
DMOVE 0,SAVE0 > ;RESTORE REGS 0-1
POPJ P, ;RETURN
SIXBIT / CMM.'N/
CMM.'N: DMOVEM N,ARG11 ;STORE FIRST ARG
DMOVEM 0,SAVE0 ;SAVE REGS 0-1
DMOVE 0,(L) ;GET SECOND ARG
PUSHJ P,CFM ;MULTIPLY
DMOVEM 0,(L) ;STORE RESULT
DMOVE 0,SAVE0 ;RESTORE REGS 0-1
POPJ P, ;RETURN
>
XLIST ;GENERATE CFM.0 ... CFM.14 AND CMM.0 ... CMM.14
.N.=16
REPEAT 7,<SYM(\<.N.==.N.-2>)>
LIST
CFM: DMOVEM T0,ARG21 ;STORE SECOND ARG
DMOVEM T2,SAVE2 ;SAVE SCRATCH REGISTERS
PUSHJ P,CALC ;CALC (A*C-B*D)
MOVEM T0,REAL ;SAVE REAL(ANS) IN REAL.
MOVE T0,ARG12 ;SET UP
EXCH T0,ARG11 ;AND CALCULATE
MOVNM T0,ARG12 ;(B*C+A*D)
PUSHJ P,CALC ;AND LEAVE IT IN T0.
MOVE T1,T0 ;STORE IMAG(ANS) IN T1
MOVE T0,REAL ;GET REAL(ANS) IN T0
DMOVE T2,SAVE2 ;RESTORE SCRATCH REGISTERS
POPJ P, ;RETURN
CALC: MOVE T0,ARG11 ;"A" TO AC0.
MOVE T2,ARG12 ;"B" TO AC2.
FMPR T0,ARG21 ;CALC A*C AND IF
JFCL OUFLO1 ;IT OVER/UNDERFLOWS, GO TO OUFLO1.
SECOND: FMPR T2,ARG22 ;CALC B*D AND IF
JFCL OUFLO2 ;IT OVER/UNDERFLOWS, GO TO OUFLO2.
FSBR T0,T2 ;CALC A*C-B*D AND
POPJ P, ;RETURN.
OUFLO1: MOVE T3,ARG22 ;"D" TO AC3.
JUMPE T0,UNDER1 ;UNDERFLOW, GO TO UNDER1
FMPRB T2,T3 ;CALC B*D AND
JFCL OFL1OU ;IF OVER/UNDERFLOW GO TO OFL1OU.
XOR T3,T0 ;OVERFLOW + NORMAL CASE
JUMPGE T3,SROVNM ;IF SIGNS ARE THE SAME, GO TO SROVNM
INFIN: MOVEM T0,ANSTMP ;STORE ANSWER.
INFIN1: $LCALL 999
;LERR (999,%,Complex overflow at $1L,<T1>)
MOVE T0,ANSTMP ;PICK UP ANSWER.
POPJ P, ;RETURN.
OFL1OU: JUMPE T2,INFIN ;OVERFLOW + UNDERFLOW CASE GOES TO INFIN.
XOR T3,T0 ;OVERFLOW + OVERFLOW CASE.
JUMPGE T3,SROVOV ;IF SIGNS ARE THE SAME, GO TO SROVOV
JRST INFIN ;O'E, GO TO INFIN.
SROVNM: JUMPE T2,INFIN ;OVERFLOW + ZERO, GO TO INFIN.
MOVE T0,ARG21 ;C TO AC0.
FSC T0,-1 ;CALC. C/2.
FMPR T0,ARG11 ;CALC. (C/2)*A AND IF
JFCL SROV1 ;IT OVERFLOWS, IMMEDIATELY RETURN.
FSBRM T0,T2 ;CALC ((C/2*A)-(B*D).
FADR T0,T2 ;CALC (A*C-B*D) AND
POPJ P, ;RETURN.
SROV1: MOVEM T0,ANSTMP ;STORE ANSWER
JRST INFIN1 ;GO TO ERROR
SROVOV: MOVE T0,ARG21 ;C TO AC0.
FDVR T0,ARG22 ;CALC. (C/D).
MOVE T1,ARG12 ;B TO AC1.
FDVR T1,ARG11 ;CALC. (B/A).
FSBR T0,T1 ;CALC. ((C/D)-(B/A)) AND
JFCL FIXUP ;IF IT UNDERFLOWS, GO TO FIXUP.
FMPR T0,ARG22 ;CALC. ((C)-(B/A)*D)) AND
JFCL ;SUPPRESS OVERFLOW MESSAGE.
FMPR T0,ARG11 ;CALC (A*C-B*D) AND
POPJ P, ;RETURN.
FIXUP: FMPR T1,ARG22 ;CALC. ((B/A)*D).
MOVE T0,ARG21 ;C TO AC0.
FSBR T0,T1 ;CALC. (C-(B/A)*D).
FMPR T0,ARG11 ;CALC. (A*C-B*D) AND
POPJ P, ;RETURN.
UNDER1: FMPR T2,T3 ;CALC. B*D AND GO
JFCL U1OU ;TO U1OU IF OVER/UNDERFLOW.
MOVM T3,T2 ;IF B*D IS <2**-105,
CAMGE T3,[030400000000] ;GO
JRST UNDR0 ;TO .+3.
MOVN T0,T2 ;NO BITS CAN BE SAVED,
POPJ P, ;FIX SIGN AND RETURN.
UNDR0: FSC T2,32 ;CALC B*D*(2**27).
UNDR1: MOVE T0,ARG11 ;CALC
FSC T0,32 ;A*C*(2**27)
FMPR T0,ARG21 ;AND SUPPRESS AN
JFCL ;ERROR MESSAGE.
UNDR4: FSBR T0,T2 ;CALC (2**27)(A*C-B*D),
UNDR5: FSC T0,-32 ;CALC (A*C-B*D) AND
POPJ P, ;RETURN.
U1OU: JUMPE T2,U1OU1 ;IF OVERFLOW + UNDERFLOW, GO
MOVNM T2,ANSTMP ;TO ERROR
JRST INFIN1 ;RETURN.
U1OU1: MOVE T2,ARG12 ;O'E, TRY TO SAVE BITS.
FSC T2,32 ;CALC. B*(2**27).
FMPR T2,ARG22 ;CALC B*D*(2**27)
JFCL ;AND SUPPRESS UNDERFLOW MESSAGE
MOVE T0,ARG11 ;CALC. A*(2**27) AND
FSC T0,32 ;A*(2**27)*C
FMPR T0,ARG21 ;AND SUPPRESS
JFCL ;UNDERFLOW MESSAGE.
JUMPN T2,UNDR4 ;IF BOTH UNDERFLOW, RETURN
JUMPN T0,UNDR5 ;AN ERROR MESSAGE. O'E,
$LCALL 888
;LERR (888,%,Complex underflow at $1L,<T1>)
SETZ T0, ;CLEAR AC0
POPJ P, ;CALC.
OUFLO2: JUMPN T2,OUFLO3 ;JUMP AHEAD IF (B*D) OVERFLOWED.
CAML T0,[030400000000] ;IF NO BITS CAN BE SAVED,
POPJ P, ;RETURN.
FSC T0,32 ;CALC. (2**27)*(A*C)
MOVE T2,ARG12 ;AND
FSC T2,32 ;THEN
FMPR T2,ARG22 ;(2**27)*(B*D)
JFCL ;AND GO TO THE REST
JRST UNDR4 ;OF THE CALC.
OUFLO3: MOVEM T2,T1 ;OVERFLOW + NORMAL.
XOR T1,T0 ;IF THE SIGNS ARE THE SAME,
JUMPGE T1,SROVN ;GO TO SROVN
MOVNM T2,ANSTMP ;THE ANSWER AND
JRST INFIN1 ;GO TO ERROR EXIT.
SROVN: MOVE T3,ARG22 ;"D" TO AC3.
FSC T3,-2 ;CALC (D/2).
FMPR T3,ARG12 ;CALC (B*(D/2)) AND IF
JFCL FIXUP2 ;IT OVERFLOWS, GO TO FIXUP2.
FSBR T0,T3 ;CALC. ((A*C)-(B*(D/2))).
FSBR T0,T3 ;CALC (A*C-B*D) AND
POPJ P, ;RETURN.
FIXUP2: MOVNM T3,ANSTMP ;SET UP THE ANSWER
JRST INFIN1 ;AND GO TO ERROR EXIT.
SEGMENT DATA
SAVE0: BLOCK 2
SAVE2: BLOCK 2
ARG11: BLOCK 1
ARG12: BLOCK 1
ARG21: BLOCK 1
ARG22: BLOCK 1
REAL: BLOCK 1
ANSTMP: BLOCK 1
PRGEND
TITLE CFDV COMPLEX DIVISION
SUBTTL KAREN KOLLING /KK/CKS 28-Oct-81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;FROM LIB40 VERSION .027
;FROM V.026 24-MAR-70
;THIS PROGRAM CALCULATES THE QUOTIENT OF TWO FLOATING POINT
;NUMBERS IN THE FOLLOWING MANNER:
;(A+IB)/(C+ID) = [(AC+BD)+I(BC-AD)]/(C**2+D**2)
;IF OVER OR UNDERFLOW OCCURS AT ANY POINT IN THE
;CALCULATION, THE PROGRAM JUMPS TO A MORE
;INVOLVED ROUTINE IN WHICH THE ORDER IN WHICH THE
;TERMS ARE MULTIPLIED OR DIVIDED TOGETHER IS VERY
;IMPORTANT. OVER OR UNDERFLOW IS NOT PERMITTED
;UNLESS THE FINAL ANSWER OVER OR UNDERFLOWS.
;THERE ARE TWO SETS OF ENTRY POINTS FOR THE PROGRAM - ONE
;FOR THE DIVIDE ROUTINES AND ONE FOR THE DIVIDE-TO-MEMORY
;ROUTINES. THE CALLING SEQUENCE FOR THE DIVIDE ROUTINES IS
; XMOVEI L,ARG2
; PUSHJ P,CFD.N
;WHERE N=0,2,4,6,10,12,14. ARG1 IS ASSUMED TO BE IN AC N AND (N+1),
;AND THE RESULTS ARE RETURNED IN AC N AND AC (N+1)
;THE CALLING SEQUENCE FOR THE DIVIDE-TO-MEMORY ROUTINES IS
;AS FOLLOWS:
; XMOVEI L,ARG2
; PUSHJ P,CDM.N
;WHERE N AND ARG2 ARE AS ABOVE AND THE RESULTS ARE RETURNED IN
;ARG2 AND ARG2+1.
SEARCH MTHPRM
SEGMENT CODE
SALL
DEFINE SYM (N) <
ENTRY CFD.'N,CDM.'N
SIXBIT / CFD.'N/
CFD.'N: DMOVEM N,SAB ;STORE NUMERATOR
IFN N,< DMOVEM 0,SAVE0 > ;SAVE REGS 0-1
DMOVE 0,(L) ;GET SECOND ARG
PUSHJ P,CFD ;DIVIDE
IFN N,< DMOVE N,0 ;STORE RESULT
DMOVE 0,SAVE0 > ;RESTORE REGS 0-1
POPJ P, ;RETURN
SIXBIT / CDM.'N/
CDM.'N: DMOVEM N,SAB ;STORE NUMERATOR
DMOVEM 0,SAVE0 ;SAVE 0-1
DMOVE 0,(L) ;GET SECOND ARG
PUSHJ P,CFD ;DIVIDE
DMOVEM 0,(L) ;STORE RESULT
DMOVE 0,SAVE0 ;RESTORE 0-1
POPJ P,
>
XLIST ;GENERATE CFD.0 ... CFD.14 AND CDM.0 ... CDM.14
.N.==16
REPEAT 7,<SYM(\<.N.==.N.-2>)>
LIST
CFD: DMOVEM T0,SCD ;STORE DENOMINATOR IN SCD AND LCD
JUMPE T0,CZERO ;GO TO CZERO IF C=0, TO
JUMPE T1,DZERO ;DZERO IF D=0,
SKIPE SAB ;TO NORMAL IF ONLY
JRST NORMAL ;ONE OF A AND B
SKIPE LAB ;=0, TO ABEQZR
JRST NORMAL ;IF
ABEQZR: SETZB T0,REAL ;A=B=0
JRST OUT1 ;RETURN (0,0)
DZERO: MOVE T0,SAB ;D IS ZERO, SO REAL (ANS)
FDVR T0,SCD ;= A/C AND
MOVEM T0,REAL ;IMAG (ANS) =
MOVE T0,LAB ;B/C
FDVR T0,SCD ;AND THEN
JRST OUT1 ;GO TO EXIT
CZERO: MOVE T0,LAB ;C IS ZERO (POSSIBLY D IS
FDVR T0,LCD ;ZERO ALSO). REAL (ANS) =
MOVEM T0,REAL ;B/D AND
MOVN T0,SAB ;IMAG (ANS) =
FDVR T0,LCD ;-A/D AND THEN
JRST OUT1 ;GO TO EXIT.
NORMAL: FMPR T0,SCD ;THIS SIMPLE ROUTINE
JFCL NTSMPL ;CALCULATES
FMPR T1,LCD ;REAL (ANS) =
JFCL NTSMPL ;(A*C+B*D)/(C(2)+D(2))
FADR T0,T1 ;AND
JFCL NTSMPL ;IMAG (ANS) =
MOVEM T0,TEMP ;(B*C-A*D)/(C(2)+D(2))
MOVE T0,SCD ;BUT
FMPR T0,SAB ;IF
JFCL NTSMPL ;AT
MOVE T1,LAB ;ANY
FMPR T1,LCD ;POINT
JFCL NTSMPL ;OVER
FADR T0,T1 ;OR
JFCL NTSMPL ;UNDERFW
FDVR T0,TEMP ;OCCURS
JFCL NTSMPL ;IT
MOVEM T0,REAL ;JUMPS
MOVE T0,SAB ;TO
FMPR T0,LCD ;NTSMPL
JFCL NTSMPL ;FOR
MOVE T1,LAB ;A
FMPR T1,SCD ;DIFFERENT
JFCL NTSMPL ;CALCULATION.
FSBRM T1,T0 ;IF THERE IS
JFCL NTSMPL ;OVER OR
FDVR T0,TEMP ;UNDERFLOW
JRST OUT1 ;IT EXITS TO OUT1.
NTSMPL: DMOVEM T2,SAVE2 ;SAVE THE CONTENTS OF AC2-3.
MOVEM T4,SAVE4 ;SAVE THE CONTENTS OF AC4.
MOVE T2,SAB ;REARRANGE THE REAL
MOVM T0,SAB ;AND IMAG PARTS OF
MOVM T1,LAB ;THE 1ST ARG
MOVEI T4,1 ;SO THAT THE SMALLER ONE
CAMG T0,T1 ;(IN MAGNITUDE) IS IN
JRST NTSMP1 ;SAB, AND THE OTHER IS IN
EXCH T2,LAB ;LAB AND SET UP AC4 AS
MOVEM T2,SAB ;A FLAG WORD TO TELL
MOVN T4,T4 ;WHICH PART IS WHERE.
NTSMP1: MOVE T2,SCD ;REARRANGE THE REAL
MOVM T0,SCD ;AND IMAG PARTS OF THE OTHER ARG
MOVM T1,LCD ;SO THAT THE SMALLER ONE
CAMG T0,T1 ;(IN MAGNITUDE) IS IN SCD AND
JRST NTSMP2 ;THE OTHER IS IN LCD AND
EXCH T2,LCD ;FIX AC4 TO TELL WHICH
MOVEM T2,SCD ;PART
IMULI T4,3 ;IS WHERE.
NTSMP2: MOVE T3,LCD ;CALCULATE
FDVRB T2,T3 ;THE
JFCL ;TERM
FMPR T2,T2 ;1+(SCD/LCD)**2
JFCL ;AND
FADRI T2,(1.0) ;STORE IT IN
MOVEM T2,DENOM ;DENOM.
MOVE T0,SAB ;CALCULATE
FDVR T0,LAB ;(SCD/LCD)*(SAB/LAB)
JFCL ;SUPPRESSING
FMPRM T0,T3 ;UNDERFLOW
JFCL ;IF NECESSARY.
CAIN T4,1 ;FIX THE AC FLAG WORD
MOVNI T4,2 ;FOR
ADDI T4,1 ;EASY COMPARISONS.
SKIPL T4 ;CALCULATE
MOVN T3,T3 ;+-1 +-(SCD/LCD)*(SAB/LAB),
FADRI T3,(1.0) ;WHERE THE SIGNS
SKIPN T4 ;DEPEND ON
MOVN T3,T3 ;THE AC FLAG WORD.
PUSHJ P,CALC34 ;JUMP TO CALC OF (LAB/LCD)*(AC3/DENOM).
MOVEM T0,REAL ;STORE IT IN REAL(ANS) AND
MOVEM T0,IMAG ;IMAG(ANS)(THE TEMP. LOCATIONS).
MOVE T3,SAB ;CALCULATE
MOVE T2,SCD ;+-(SAB/LAB)+-(SCD/LCD)
CAMN T4,[-2] ;WHERE THE SIGNS
MOVN T2,T2 ;AGAIN DEPENDS ON
CAMN T4,[-1] ;THE AC FLAG WORD,
MOVN T3,T3 ;AND IF UNDERFLOW
MOVE T0,T2 ;OCCURS JUMP TO
MOVE T1,T3 ;THE
FDVR T3,LAB ;SPECIAL
JFCL OVER1 ;ROUTINES-
FDVR T2,LCD ;OVER1,
JFCL OVER2 ;OVER2,
ADD2: FADR T3,T2 ;OR
JFCL OVER3 ;OVER3.
CALCIR: PUSHJ P,CALC34 ;JUMP TO CALC OF(LAB/LCD)*(AC3/DENOM).
STIR: JUMPGE T4,STR ;STORE THIS IN
MOVEM T0,IMAG ;THE CORRECT
JRST STX ;PART OF THE
STR: MOVEM T0,REAL ;ANSWER (TEMP. LOCATION).
STX: DMOVE T2,SAVE2 ;RESTORE THE CONTENTS OF AC2-3.
MOVE T4,SAVE4 ;RESTORE THE CONTENTS OF AC4.
OUT: SKIPA T1,IMAG ;PICK UP THE IMAG (ANS)
OUT1: MOVE T1,T0 ;STORE IMAG IN T1
MOVE T0,REAL ;GET REAL PART
POPJ P, ;DONE
CALC34: MOVM T1,LAB ;CALC34 CALCS. (LAB/LCD)*(AC3/DENOM).
MOVM T2,LCD ;/LAB/ TO AC T1 AND /LCD/ TO AC 2.
MOVM T0,T3 ;/AC3/ TO AC 0.
CAMGE T2,ONE ;GO TO CASEA IF
JRST CASEA ;/LCD/<1.0. O'E, STAY HERE.
CAMGE T1,ONE ;GO TO CASEB IF /LCD/>1.0 AND
JRST CASEB ;/LAB/<1.0 OR IF
CAMGE T0,ONE ;(>1)(<1)/(>1)(>1).
JRST CASEB ;STAY HERE IF (>1)(>1)/
FDVR T2,T1 ;(>1)(>1).
STD: FDVR T0,DENOM ;CALCULATE
FDVR T0,T2 ;IT AND GO
JRST SETSGN ;TO SETSGN.
CASEB: FMPR T0,T1 ;CALCULATE CASE B AND
JRST STD ;GO TO SETSGN (EVENTUALLY).
CASEA: FMPR T2,DENOM ;CONDENSE THE DENOMINATOR INTO AC 2.
CAMLE T1,ONE ;THEN (<1)(<1)/(<1) GOES
JRST CHKAGN ;TO SR, AND (>1)(><1)/(<>1)
CAMLE T0,ONE ;GOES TO CHKAGN,
JRST NOTRVS ;AND EVERYTHING ELSE
CAMG T2,ONE ;GOES TO
JRST SR ;NOTRVS.
NOTRVS: FMPRM T1,T0 ;CALCULATE
JFCL SETSGN ;NOTRVS AND GO
FDVR T0,T2 ;TO
JRST SETSGN ;SETSGN.
CHKAGN: CAMG T0,ONE ;(>1)(<1)/(<>1)
JRST NOTRVS ;AND (>1)(>1)/(<1)
CAMG T2,ONE ;GO TO
JRST NOTRVS ;NOTRVS.
FDVR T1,T2 ;(>1)(>1)/(>1) IS
FMPRM T1,T0 ;CALCULATED AND
JRST SETSGN ;GOES TO SETSGN.
SR: MOVEM T1,TEMP ;SR CALCULATES
FDVR T1,T2 ;(<1)(<1)/(<1)
JFCL OV1 ;AND SINCE
FMPRM T1,T0 ;(<1)/(<1)
JRST SETSGN ;CAN
OV1: MOVE T1,TEMP ;OVERFLOW, IT
MOVEM T0,TEMP ;CALCULATES
FDVR T0,T2 ;IT
JFCL OV2 ;WHICHEVER
FMPRM T1,T0 ;WAY
JRST SETSGN ;IS
OV2: MOVE T0,TEMP ;NECESSARY
FMPRM T1,T0 ;AND
FDVR T0,T2 ;THEN GOES TO SETSGN.
SETSGN: MOVE T1,LAB ;GET THE
XOR T1,LCD ;SIGN OF THE
XOR T1,T3 ;RESULT IN AC 1.
SKIPG T1 ;SET THE SIGN OF
MOVN T0,T0 ;THE ANSWER
POPJ P, ;AND EXIT.
OVER1: FDVR T2,LCD ;IF ONE
JFCL UU ;TERM
MOVE T3,T2 ;UNDERFLOWS
MOVE T0,T1 ;AND THE OTHER
MOVE T2,LAB ;TERM IS SO LARGE
JRST OVER21 ;THAT NO BITS
OVER2: MOVE T2,LCD ;CAN BE
OVER21: MOVEM T2,SAVE5 ;SAVED,
JUMPE T3,UU ;THEN
MOVM T1,T3 ;RETURN
CAML T1,[030400000000] ;RIGHT
JRST CALCIR ;AWAY.
FSC T3,70 ;O'E, TRY TO
FSC T0,70 ;SAVE SOME BITS
FDVRM T0,T2 ;BY MULTIPLYING
JFCL ;THE TERMS BY 2**56,
FADRB T3,T2 ;ADDING THE TERMS, AND THEN / BY 2**56.
FSC T3,-70 ;IF THE RESULT UNDERFLOWS, GO
JFCL SRX ;TO SRX, O'E GO BACK
JRST CALCIR ;TO THE MAIN ROUTINE.
OVER3: MOVE T3,T1 ;SET UP
FDVR T3,LAB ;AC 2 FOR
FSC T3,70 ;SRX. AC 2 WILL
FSC T2,70 ;CONTAIN THE TERM
FADR T2,T3 ;*(2**56).
SRX: MOVEM T5,SAVE5 ;SAVE THE CONTENTS OF AC 5.
MOVE T0,LCD ;THIS IS AN
MOVE T1,LAB ;ALTERNATE
MOVM T5,T1 ;CALCULATION TO
CAMGE T5,ONE ;CALC34, AND TAKES
JRST SRX2 ;ADVANTAGE OF
FMPR T1,T2 ;THE FACT THAT HERE
MOVM T5,T0 ;DENOM CONTAINS 1.0.
CAML T5,ONE ;THE ORDER OF
JRST SRY ;CALCULATION
FSC T0,70 ;DEPENDS
FDVRM T1,T0 ;ON THE SIZE OF
JRST SETSN2 ;THE TERMS. AFTER
SRY: FDVRM T1,T0 ;THE CALCULATION A
FSC T0,-70 ;TRANSFER
JRST SETSN2 ;IS
SRX2: FDVRM T2,T0 ;MADE
FMPR T0,T1 ;TO
FSC T0,-70 ;SETSN2
SETSN2: MOVE T5,SAVE5 ;RESTORE THE CONTENTS OF AC 5.
JRST STIR ;GO BACK TO MAIN ROUTINE.
UU: MOVEM T0,SAB ;ANOTHER ALTERNATE
MOVEM T1,SCD ;CALCULATION TO
FMPR T1,LCD ;CALC34
FMPR T0,LAB ;USED WHEN
FADR T0,T1 ;S/L FOR
JFCL UND ;BOTH SETS
FDVR T0,LCD ;HAS UNDERFLOWED
FDVR T0,LCD ;OR FOR
JRST STIR ;UNDERFLOW PLUS
UND: $LCALL 888
;LERR (888,%,Complex underflow at $1L,<T1>)
JRST STIR ;TO ADD2+3.
ONE: 201400000000
SEGMENT DATA
SAVE0: BLOCK 0
SAVE2: BLOCK 1
SAVE3: BLOCK 1
SAVE4: BLOCK 1
SAVE5: BLOCK 1
SAB: BLOCK 1
LAB: BLOCK 1
SCD: BLOCK 1
LCD: BLOCK 1
TEMP: BLOCK 1
REAL: BLOCK 1
IMAG: BLOCK 1
DENOM: BLOCK 1
PRGEND
TITLE REAL.C COMPLEX TO REAL CONVERSION
SUBTTL H. P. WEISS/HPW/CKS/CDM 17-Nov-81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;FROM V.005
;COMPLEX TO REAL CONVERSION FUNCTION
;THIS ROUTINE RETURNS THE REAL PART OF A COMPLEX NUMBER
;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
; XMOVEI L,ARGLST
; PUSHJ P,REAL.C
;THE ANSWER IS RETURNED IN ACCUMULATOR A
SEARCH MTHPRM
NOSYM
SEGMENT CODE
SALL
ENTRY REAL.C
REAL.C: MOVE T0,@(L) ;GET REAL PART OF ARGUMENT
GOODBY (1) ;RETURN
PRGEND
TITLE AIMAG COMPLEX TO REAL CONVERSION
SUBTTL H. P. WEISS/HPW/CKS 28-Oct-81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY AIMAG
EXTERN AIMAG.
AIMAG=AIMAG.
PRGEND
TITLE AIMAG. COMPLEX TO REAL COVERSION
SUBTTL H. P. WEISS/HPW/CKS 28-Oct-81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;FROM LIB40 V.022
;FROM V.005
;COMPLEX TO IMAGINARY CONVERSION FUNCTION
;THIS ROUTINE RETURNS THE IMAGINARY PART OF A COMPLEX NUMBER
;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
; XMOVEI L,ARGLST
; PUSHJ P,AIMAG.
;THE ANSWER IS RETURNED IN ACCUMULATOR T0
SEARCH MTHPRM
NOSYM
SEGMENT CODE
SALL
HELLO (AIMAG,.) ;ENTRY TO AIMAG ROUTINE
XMOVEI T1,@(L) ;PUT IMAG(ARG)
MOVE T0,1(T1) ;IN AC T0
GOODBY (1) ;RETURN
PRGEND
TITLE CONJG COMPLEX CONJUGATE FUNCTION
SUBTTL H. P. WEISS/HPW/CKS 28-Oct-81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CONJG
EXTERN CONJG.
CONJG=CONJG.
PRGEND
TITLE CONJG. COMPLEX CONJUGATE FUNCTION
SUBTTL H. P. WEISS/HPW/CKS 28-Oct-81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;FROM LIB40 V.022
;COMPLEX CONJUGATE FUNCTION
;THIS ROUTINE CALCULATES THE COMPLEX CONJUGATE OF ITS ARGUMENT
;THE REAL PART OF THE ANSWER IS LEFT IN ACCUMULATOR T0, AND THE
;IMAGINARY PART IN ACCUMULATOR T1
SEARCH MTHPRM
NOSYM
SEGMENT CODE
SALL
HELLO (CONJG,.) ;ENTRY TO CONJG ROUTINE
DMOVE T0,@(L) ;GET COMPLEX ARGUMENT
MOVN T1,T1 ;NEGATE IMAGINARY PART
GOODBY ;RETURN
PRGEND
TITLE CMPLX REAL TO COMPLEX CONVERSION
SUBTTL H. P. WEISS/HPW/CKS 28-Oct-81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
SEARCH MTHPRM
NOSYM
ENTRY CMPLX
EXTERN CMPLX.
CMPLX=CMPLX.
PRGEND
TITLE CMPLX. REAL TO COMPLEX CONVERSION
SUBTTL H. P. WEISS/HPW/CKS 28-Oct-81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1980, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
;FROM LIB40 V.022
;REAL TO COMPLEX CONVERSION FUNCTION
;THIS ROUTINE CONVERTS TWO REAL ARGUMENTS TO A COMPLEX NUMBER
;OF THE FORM ARG1+I*ARG2
;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
; MOVEI L,ARG
; PUSHJ P,CMPLX
;THE REAL PART OF THE ANSWER IS LEFT IN ACCUMULATOR A, AND THE
;IMAGINARY PART IS LEFT IN ACCUMULATOR B
SEARCH MTHPRM
NOSYM
SEGMENT CODE
SALL
HELLO (CMPLX,.) ;ENTRY TO CMPLX ROUTINE
MOVE T0,@(L) ;GET REAL PART OF COMPLEX ANSWER
MOVE T1,@1(L) ;GET IMAGINARY PART OF COMPLEX ANS
GOODBY (2) ;RETURN
PRGEND
TITLE CMPL.C Conversion from complexes to complex
SUBTTL C. McCutcheon - 28-Oct-81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
; Fortran Instrinsic Function to return a complex number for
; two complex numbers passed to it.
; From ANSI-77 standard:
; CMPLX(A1,A2) is the complex value whose real part is REAL(A1)
; and whose imaginary part is REAL(A2).
; Required (called) routines: SNG.0, SNG.2
SEARCH MTHPRM
SALL
SEGMENT CODE
SALL
HELLO (CMPL.C,.) ;Enter routine
MOVE T0,@(L) ;Real half.
MOVE T1,@1(L) ;Imaginary half
GOODBYE ;Return
PRGEND
TITLE CMPL.D Conversion from double prec to complex
SUBTTL C. McCutcheon - 28-Oct-81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
; Fortran Instrinsic Function to return a complex number for
; two double precision numbers passed to it.
; Required (called) routines: SNG.0, SNG.2
SEARCH MTHPRM
SALL
SEGMENT CODE
SALL
EXTERN SNG.0 ;Located in
EXTERN SNG.2 ; FORDBL.MAC
HELLO (CMPL.D,.) ;Enter routine
PUSH P,T2 ;Save Ac's
PUSH P,T3
DMOVE T0,@(L) ;Real half.
PUSHJ P,SNG.0 ;Round DP to Real.
DMOVE T2,@1(L) ;Imaginary half.
PUSHJ P,SNG.2 ;Round DP to Real.
;(There is no SNG.1 for optimizing ac's)
MOVE T1,T2 ;Set up complex #.
POP P,T3 ;Restore AC's.
POP P,T2
GOODBYE ;Return
PRGEND
TITLE CMPL.G Conversion from gfloating to complex
SUBTTL Randy Meyers -- Created by Edit [3201] 17-May-82
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
; Fortran Instrinsic Function to return a complex number for
; two gfloating numbers passed to it.
SEARCH MTHPRM
SALL
SEGMENT CODE
SALL
HELLO (CMPL.G,.) ;Enter routine
PUSH P,T2 ;Save Ac
DMOVE T0,@0(L) ;Real half.
EXTEND T0,[GSNGL T0] ;Round gfloating to Real.
DMOVE T1,@1(L) ;Imaginary half.
EXTEND T1,[GSNGL T1] ;Round gfloating to Real.
POP P,T2 ;Restore AC.
GOODBYE ;Return
PRGEND
TITLE CMPL.I Conversion from integers to complex
SUBTTL C. McCutcheon - 28-Oct-81
;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY
;TRANSFERRED.
;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE
;AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.
;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.
; Fortran Instrinsic Function to return a complex number for
; two integer numbers passed to it.
SEARCH MTHPRM
SALL
SEGMENT CODE
SALL
HELLO (CMPL.I,.) ;Enter routine
MOVE T0,@(L) ;Real half
MOVE T1,@1(L) ;Imaginary half
FLTR T0,T0 ;Convert to
FLTR T1,T1 ; real
GOODBYE ;Return
END