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