TITLE Polynomial Zerofinder (Real Coeffiecients) 1.6 **miscdoc+************************************************** ************************************************************ ** ** File: ProotR.s Version 1.6, 05/25/93 ** Machine: Alcuin ** ************************************************************ **+miscdoc************************************************** INCLUDE ./PolyN.h ( include macros and equates ) ASSEMBLE ***************************************************************** * Status Bit Equates sCmpArry EQU 0 Flags real-valued complex array sCITER EQU 0 Flags CITER call for ITRUTI ***************************************************************** STITLE proot ***************************************************************** ***************************************************************** ** ***Name: proot ** ***Category: ** ***Abstract: Polynomial zerofinder. ** ***Stack on Entry: arry% | arryC% ** ***Stack on Exit: rootsC% ** ***Error Exits: ***Author: Paul J. McClellan ***Date Written: February 18, 1987 ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** 08/04/92 Paul McClellan handle real-val cmplx arry ** ***Detail: ** ** Assume the polynomial ** ** N ** --- k ** P(Z) = > A *Z A real ** --- k k ** k=0 ** ** The operand array A is interpreted as containing the real- ** valued coefficients of the given polynomial. The coefficients ** A(N), A(N-1), ..., A(0) are contained in this order from left ** to right along each row from the first row to the last. Thus, ** an (N+1)-element operand array represents an Nth degree ** polynomial. The operand array A may be real or complex. ** ** PROOT will compute all N complex zeros of the polynomial P(Z). ** If A contains less than two elements, then PROOT takes an ** error exit. If the given polynomial has degree N the result ** will be a complex N-vector containing the zeros roughly in ** order of increasing magnitude. ** ** Although A has N+1 elements, the degree of the source poly- ** nomial is not necessarily N. This may be due to leading zeros ** in the operand array and this case is trapped. ** ** The implementation is a modified version of the HP-71 Math ** ROM PROOT runtime execution code implemented by Laurence Grodd. ** The major modifications are: ** ** a. Removal of argument Nan and Infinity checking. ** b. Support for complex coefficients. ** c. Explicit dimensioning of the complex destination array. ** d. Replacement of dedicated RAM math scratch space and utilities ** by space in PROOT hex string and other utilities. ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** RPL NULLNAME proot ( arry ) :: DUP ARSIZE ( arry N+1 ) DUPTWO #< casedrop SETDIMERR ( bc ) ( arry N+1 ) DUP#1- ONE{}N C%0 MAKEARRY UNROT ( rootsC% arry N+1 ) OVER RealEl? ITE ( bc ) proot_r proot_c ( rootsC% ) DUP RealEl? case MATRE ( bc ) ( rootsC% ) ; ASSEMBLE STITLE proot_r ***************************************************************** ***************************************************************** ** ***Name: proot_r ** ***Category: ** ***Abstract: Polynomial zerofinder for real-valued coefficients. ** ***Stack on Entry: rootsC% arry% N+1 | ** rootsC% arryC% N+1 ** ***Stack on Exit: rootsC%' ** ***Error Exits: ***Author: Paul J. McClellan ***Date Written: February 18, 1987 ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** 08/04/92 Paul McClellan handle real-val cmplx arry ** ***Detail: ** ** Assume the polynomial ** ** N ** --- k ** P(Z) = > A *Z A real ** --- k k ** k=0 ** ** The operand array A is interpreted as containing the real- ** valued coefficients of the given polynomial. The coefficients ** A(N), A(N-1), ..., A(0) are contained in this order from left ** to right along each row from the first row to the last. Thus, ** an (N+1)-element operand array represents an Nth degree ** polynomial. The operand array A may be real or complex, but ** if it is complex the imaginary parts will be assumed zero. ** ** PROOT will compute all N complex zeros of the polynomial P(Z). ** If A contains less than two elements, then PROOT takes an ** error exit. If the given polynomial has degree N the result ** will be a complex N-vector containing the zeros roughly in ** order of increasing magnitude. ** ** Although A has N+1 elements, the degree of the source poly- ** nomial is not necessarily N. This may be due to leading zeros ** in the operand array and this case is trapped. See PROOT ** Block 1 for a description of the trap action. ** ** The algorithm is Laguerre iteration with stepsize control ** and is based upon the ZERPOL program contained in the thesis ** entitled "A Zero Finding Algorithm Using Laguerre's Method" ** by Brian T. Smith (directed by W. Kahan). ** ** The implementation is a modified version of the HP-71 Math ** ROM PROOT runtime execution code implemented by Laurence Grodd. ** The major modifications are: ** ** a. Removal of argument Nan and Infinity checking. ** b. Handling of real-valued complex coefficient array. ** c. Replacement of dedicated RAM math scratch space and ** utilities by space in PROOT hex string and other utilities. ** ***Test Classes: ** ***************************************************************** ***************************************************************** RPL NULLNAME proot_r ( rootsC% arry N+1 -> rootsC%' ) :: # 2A #* # 272 #+ ( rootsC% arry #hxs=42*[N+1]+626 ) CREATEhxs ( rootsC% arry hxs ) CODE NIBHEX 823 SB=XM=0 STITLE PROOT Variable Structure ***************************************************************** * * PROOT VARIABLE STRUCTURE * ------------------------ * * The PROOT function reserves 42N+668 nibbles of temporary memory * for its execution. This space is used to hold the values of * the PROOT variables described below. Many PROOT variables are * stored in extended precision as follows: * * * Extended | | * precision | SIGN | (1) ^ * real | MANTISSA | (15) high | * variable | EXPONENT | (5) addresses | * (21 nibbles) | | * * | | * Extended | SIGN(Im) | (1) * precision | MANTISSA(Im) | (15) ^ * complex | EXPONENT(Im) | (5) high | * variable | SIGN(Re) | (1) addresses | * (42 nibbles) | MANTISSA(Re) | (15) * | EXPONENT(Re) | (5) * | | * * * What follows is a description of all PROOT variables: * * CD: 5 nibble binary variable containing the current degree of * the polynomial. * * OD: 5 nibble binary variable containing the original degree of * the polynomial. * * STARTL: 1 nibble logical variable. Set TRUE whenever the * Laguerre iterations have successfully begun; false * otherwise. * * SPIRAL: 1 nibble logical variable. Set TRUE whenever a spiral * search for an initial iterate has begun; false * otherwise. * * CRTA: 5 nibble binary variable containing the address of the * current root vector element to store. * * E: Extended precision real variable. Contains the computed * error bound at the current iterate. * * PSCR1: * PSCR2: * PSCR3: * PSCR4: * PSCR5: * PSCR6: Six extended precision real variables used for general * scratch space during PROOT. * * Z0: Extended precision complex variable. Contains the last * accepted iterate of the iteration process. * * LGS0: Extended precision complex variable. Contains the * Laguerre step at the last accepted iterate. * * F0: Extended precision real variable. Contains the magnitude * of the polynomial at the last accepted iterate. * * ALGS0: Extended precision real variable. Contains |LGS0|. * * OUTR0: Extended precision real variable. Contains the radius * of a circle about Z0 known to contain the closest zero * of the polynomial to Z0. * * ZN: Extended precision complex variable. Contains the current * iterate being tested for acceptance. * * LGSN: Extended precision complex variable. Contains the * Laguerre step at the current iterate. * * FN: Extended precision real variable. Contains the magnitude * of the polynomial at the current iterate. * * ALGSN: Extended precision real variable. Contains |LGSN|. * * OUTRN: Extended precision real variable. Contains the radius * of a circle about the current iterate known to contain * the zero of the polynomial closest to that iterate. * * N: Extended precision real variable. Contains FLOAT(CD). * * INR0: Extended precision real variable. Contains the inner * radius of an annulus about the origin known to contain * the smallest zero of the current polynomial. * * QSCR0: * QSCR1: * QSCR2: * QSCR3: * QSCR4: * QSCR5: Six extended precision real variables used to replace * the dedicated scratch registers SCR0-5 found on some * saturn machines. These could be removed and the access * routines changed if dedicated scratch registers were * added to the machine. * * A(): Extended precision real vector. Contains the coefficients * of the current polynomial for which a root is sought. * Dimensioned to OD+1 upon initialization so that its memory * usage is 21*OD+21 nibbles. Also initialized as the * operand array. This array is allocated in the following * manner: (Let N=CD) * * | A(N) | (low memory) * | A(N-1) | * | A(N-2) | * . . * . . * . . * | A(2) | * | A(1) | * | A(0) | (high memory) * * where A(N) denotes the current leading coefficient and A(0) * denotes the current trailing coefficient. As roots are * found, this array shrinks in the direction of low memory. * * DFLTMP(): Extended precision real vector. Contains the * (temporary) coefficients of the deflated polynomial * stored exactly as array A() is stored. Dimensioned * to OD+1 upon initialization so that its memory usage * is 21*OD+21 nibbles. When an iterate is accepted as * a root, the elements of this array are transferred * to A(). As roots are found, this array shrinks in * the direction of low memory. * ***************************************************************** ***************************************************************** * * MEMORY STRUCTURE FOR PROOT VARIABLES * ------------------------------------ * * ------------ * | ->DOHXS | (5) * |----------| * | LENGTH | (5) * |----------| * [RSMEM]-----------> | CD | (5) * |----------| * [RSMEM]+5--------> | OD | (5) * |----------| * [RSMEM]+10-------> | STARTL | (1) * |----------| * [RSMEM]+11-------> | SPIRAL | (1) * |----------| * [RSMEM]+12-------> | CRTA | (5) * |----------| * [RSMEM]+17-------> | E | (21) * |----------| * [RSMEM]+38-------> | PSCRATCH | (126) * |----------| * [RSMEM]+164------> | Z0 | (42) * |----------| * [RSMEM]+206------> | LGS0 | (42) * |----------| * [RSMEM]+248------> | F0 | (21) * |----------| * [RSMEM]+269------> | ALGS0 | (21) * |----------| * [RSMEM]+290------> | OUTR0 | (21) * |----------| * [RSMEM]+311------> | ZN | (42) * |----------| * [RSMEM]+353------> | LGSN | (42) * |----------| * [RSMEM]+395------> | FN | (21) * |----------| * [RSMEM]+416------> | ALGSN | (21) * |----------| * [RSMEM]+437------> | OUTRN | (21) * |----------| * [RSMEM]+458------> | N | (21) * |----------| * [RSMEM]+479------> | R | (21) * |----------| * [RSMEM]+500------> | QSCRATCH | (126) * |----------| * [RSMEM]+626------> | A() | (21*OD+21) * |----------| * [RSMEM]+21OD+647-> | DFLTMP() | (21*OD+21) * ------------ * [RSMEM]+42OD+668-> * ***************************************************************** STITLE PROOT Block Design ***************************************************************** * * PROOT BLOCK DESIGN * ------------------ * * BLOCK1: INITIALIZATION * * BLOCK2: SCALING OF THE COEFFICIENTS * * BLOCK3: REENTRY FOR CURRENT (REDUCED) POLYNOMIAL; * FACTORIZATION OF LINEAR AND QUADRATIC POLYNOMIALS * * BLOCK4: ANNULUS FOR THE SMALLEST ZERO; * ITERATE INITIALIZATION AT THE ORIGIN * * BLOCK5: ACCEPTANCE OR REJECTION OF THE LAGUERRE STEP * * BLOCK6: ACCEPTANCE AND/OR UPDATING OF THE ITERATE * * BLOCK7: EVALUATION AT A REAL ITERATE * * BLOCK8: EVALUATION AT A COMPLEX ITERATE * * BLOCK9: REJECTION OF AN ITERATE * * BLOCK10: SPIRAL SEARCH FOR AN INITIAL ITERATE * * BLOCK11: DECLARATION OF A ZERO * ***************************************************************** STITLE PROOT BLOCK 1 ***************************************************************** * * PROOT BLOCK 1 * ------------- * * INITIALIZATION * * The function of this block is to initialize the PROOT function. * This initialization proceeds as follows: * * (1) Pop top two arguments and get pointers. * * (2) Initialize OD (original degree) variable to N and CD * (current degree) variable to N. Initialize CRTA (current * root element address) variable to the result array base * address. * * (3) Copy the operand array into RSMEM (reserved memory). * We then split and push each operand element, keeping track * of the maximum exponent. This is equivalent to the * following pseudocode: * * 10 MAXEXP = -10000 * 20 X = NEXT_OPERAND_ELEMENT * 30 SPLIT AND PUSH X INTO RSMEM * 40 IF EXPONENT(X)>MAXEXP THEN MAXEXP=EXPONENT(X) * 50 N = N-1 * 60 IF N >= 0 THEN 20 * * (4) For every leading zero in the coefficient array, store a * root at (Inf,0) (no exceptions). For each such overflow * root stored, decrement CD and update CRTA. If ever CD=0 * then exit PROOT. * * (5) For every trailing zero in the coefficient array, store a * root at (0,0). For each such root stored, decrement CD and * update CRTA. If ever CD=0 then exit PROOT. * * ENTRY: Entry to this block is only made from the calling * RPL object. * * EXIT: Exit to Block 2 with hex mode, R2[A]= MAXEXP, and * R4[A]= RSMEM. * * Also, CD>0, all coefficients are finite, and no leading * or trailing zeros are present. * ***************************************************************** ***************************************************************** * * SUB-BLOCK (1) EXECUTION * * data stack: hxs * arry% | arryC% * rootsC% * * RPL variables in CPU * hexmode, P=0 * * Pop top two arguments and get pointers. * Note real-valued complex operand array. * ***************************************************************** D1=D1+ 10 pop top two items D=D+1 A D=D+1 A GOSBVL =SAVPTR D1=D1- 10 A=DAT1 A A[A] -> hxs D0=A D0=D0+ 10 D0 = RSMEM -> CD AD0EX R4=A R4[A] = RSMEM D1=D1+ 10 A=DAT1 A A[A]-> rootsC% D0=A D0=D0+ 15 D0=D0+ 5 A=DAT0 A A[A]= N R1=A R1[A]= N D0=D0+ 5 D0-> roots[1] AD0EX R3=A R3[A] -> roots[1] ***************************************************************** * Note real-valued complex operand array. ***************************************************************** D1=D1- 5 A=DAT1 A A[A]-> arry% or arryC% D0=A D0=D0+ 10 A=DAT0 A A[A]= element type ST=0 sCmpArry assume arry% LC(5) =DOREAL ?A=C A arry type real? GOYES realArry Yes ( bc ) ST=1 sCmpArry No: flag arryC% realArry D0=D0+ 5 A=DAT0 A A[A]= dim count D0=D0+ 10 A=A-1 A ?A=0 A vector? GOYES vector ( bc ) D0=D0+ 5 vector AD0EX R2=A R2[A] -> arry[1] = A(N) ***************************************************************** * * SUB-BLOCK (2) EXECUTION * * Here: R1[A]= N * R2[A]-> arry[1] = A(N) * R3[A]-> roots[1] * R4[A]= RSMEM * P=0, hexmode * sCmpArry set iff arryC% * * Initialize OD (original degree) variable and CD (current * degree) variable to N. Initialize CRTA (current root * element address) variable to the result array base address. * ***************************************************************** A=R4 D1=A D1-> CD A=R1 A[A]= N DAT1=A A CD = N D1=D1+ 5 D1-> OD DAT1=A A OD = N D1=D1+ 7 D1-> CRTA A=R3 A[A]-> roots[1] DAT1=A A CRTA -> roots[1] AD1EX LC(5) 614 P=0 here A=A+C A D1=A D1-> A() position in RSMEM ***************************************************************** * * SUB-BLOCK (3) EXECUTION * * Here: D1-> A() position in RSMEM * R1[A]= OD (= N) * R2[A]-> arry[1] = A(N) * R4[A]= RSMEM * P=0, hexmode * sCmpArry set iff arryC% * * Copy the operand array into RSMEM (reserved memory). * We then split and push each operand element, keeping track * of the maximum exponent. This is equivalent to the * following pseudocode: * * 10 MAXEXP = -10000 * 20 X = NEXT_OPERAND_ELEMENT * 30 SPLIT AND PUSH X INTO RSMEM * 40 IF EXPONENT(X)>MAXEXP THEN MAXEXP=EXPONENT(X) * 50 N = N-1 * 60 IF N >= 0 THEN 20 * * Throughout the sub-block below, * * D1-> A() position in RSMEM * D0-> A() position in operand * R1[A]= count * R2[A]= MAXEXP * sCmpArry set iff arryC% * ***************************************************************** LCHEX 90000 CR2EX R2[A]= MAXEXP := -10000 D0=C D0-> arry[1] = A(N) copary A=DAT0 W A[W]= Re(A(i)) D0=D0+ 16 D0-> A(i-1) or Im(A(i)) ?ST=0 sCmpArry real arry ? GOYES coparyr ( bc ) D0=D0+ 16 D0-> Re(A(i-1)) coparyr GOSBVL =SPLITA (A,B)=X, decmode GOSUBL putab1 Put X into RSMEM, increment D1 C=R2 C[A]= MAXEXP P= 4 Set pointer for exponent compare ?A=C P Same signs? GOYES samsgn Jump if so ( bc ) ?C#0 P MAXEXP < 0? GOYES abig EXP(X)>MAXEXP if so ( bc ) GONC asmall (BET) samsgn ?A<=C A EXP(X)<=MAXEXP? GOYES asmall Jump if so ( bc ) abig R2=A Update MAXEXP asmall ***************************************************************** * Check for attention key abort here. ***************************************************************** GOSBVL =ATTNCHK hex A=R1 A[A]=count A=A-1 A A[A]=count-1 R1=A Update count GONC copary loop if count>0 ( bc ) P= 0 ***************************************************************** * * SUB BLOCKS (4) AND (5) EXECUTION * * * Here: R2[A]= MAXEXP * R4[A]= RSMEM * P=0, hexmode * * (4) For every leading zero in the coefficient array, store a * root at (Inf,0) (no exceptions). For each such overflow * root stored, decrement CD and update CRTA. If ever CD=0 * then exit PROOT. * * (5) For every trailing zero in the coefficient array, store a * root at (0,0). For each such root stored, decrement CD and * update CRTA. If ever CD=0 then exit PROOT. * ***************************************************************** C=0 S Want exact (Inf,0) in result GOSUBL LEAD0 Eliminate leading zeros C=0 S Want exact (0,0) in result GOSUBL TRAL0 Eliminate trailing zeros STITLE PROOT BLOCK 2 ***************************************************************** * * PROOT BLOCK 2 * ------------- * * SCALING OF THE COEFFICIENTS * * The function of this block is to scale the coefficients of * the current polynomial. * * Entry: This block has two entry points: * * (1) SCALE ENTRY: This point is entered only by falling through * from PROOT Block 1 after initialization of the coefficient * array. The function of entry SCALE is to perform initial * scaling of the polynomial. The maximum exponent in the * coefficient array, MAXEXP, is compared to 480 and if MAXEXP is * greater than or equal to 480, then no initial scaling of the * polynomial is performed and this block is exited via branch to * PROOT Block 3 at the REENT entry. Otherwise, the polynomial * is scaled by the factor SC = 480-MAXEXP, which brings the * maximum magnitude coefficient up to the order 1E480. * * (2) RESCAL ENTRY: This entry is entered only from PROOT Blocks * 7 or 8 when extended range overflow is detected in the * evaluation loops. In this case, the polynomial is scaled by * the factor SC= -12 in the hopes of removing the overflow * condition. * * After any scaling, this block then calls LEAD0 to eliminate * any leading zero coefficients. Note that the leading * coefficient may only vanish through rescaling since it does * not change throughout the iteration. Thus, this call is * vacuous if this block was entered at SCALE. In addition, the * leading coefficient would have to be re-scaled to less than * 1E-5000 by factors of 1E-12 so that the chance of it vanishing * is small. However, should it vanish, then our scaling * procedure has broken down and a serious overflow condition * exists. This is indicated by reporting a zero at (1E1000,0) * with overflow indicated. The remaining coefficients are then * slid back down. If all leading coefficients have vanished, then * PROOT is exited. * * Normal exit is only made to PROOT Block 3 at REENT entry to * begin the iteration for the current polynomial. * ***************************************************************** ***************************************************************** * * BLOCK 2 - SCALE ENTRY * * Here: R2[A]= MAXEXP * R4[A]= RSMEM * P=0, hexmode * ***************************************************************** SCALE A=R2 A[A]= MAXEXP LCHEX 00480 C[A]= 480 P= 4 Set pointer for exponent test ?A#0 P MAXEXP<0? GOYES getsc Definitely scale if so ( bc ) ?A>=C A MAXEXP >= 480? GOYES REENT No initial scaling if so ( bc ) getsc SETDEC Mode for subtract C=C-A A SC (in C[A]) = 480-MAXEXP GOTO doscal Go scale coefficients ( c ) ***************************************************************** * * BLOCK 2 - RESCAL ENTRY * * Here: R4[A]= RSMEM * mode undetermined * ***************************************************************** RESCAL P= 0 Set pointer for load LCHEX 99988 C[A]= -12 ***************************************************************** * Here: C[A]= SC * R4[A]= RSMEM * P=0, mode undetermined * * The code below scales all the coefficients by SC, that is, * sets A(I) = SCALE10(A(I),SC) for all I. * * It then eliminates any leading zeros which may have resulted * from re-scaling (reporting overflow zeros at (1E1000,0)) * and, if not all leading coefficients were zero, falls through * to PROOT Block 3. * ***************************************************************** doscal R0=C R0[A]= SC GOSUBL ->CD,A D1-> CD, A[A]-> A(N), hexmode,P=0 C=DAT1 A C[A]= CD R1=C R1[A]= CD D1=A D1-> A() scallp GOSUBL getab1 (A,B) = A(I), D1->A(I-1) ?B=0 W Don't scale 0 (no dirty zeros) GOYES nodir0 ( bc ) SETDEC Mode for scale C=R0 C[A]= SC A=A+C A Scale coefficient GOSUBL prtPAK Underflows to 0 (no overflow) GOSUBL ptab1- Replace A(I) with scaled A(I) nodir0 GOSBVL =ATTNCHK Exit if attention key abort C=R1 C[A]= CD C=C-1 A CD=CD-1 R1=C Store new CD GONC scallp Scale CD+1 elements ( bc ) C=0 S Want leading zeros to report C=C+1 S overflow messages GOSUBL LEAD0 Eliminate leading zeros. STITLE PROOT BLOCK 3 ***************************************************************** * * PROOT BLOCK 3 * ------------- * * REENTRY FOR CURRENT (REDUCED) POLYNOMIAL. FACTORIZATION * OF LINEAR AND QUADRATIC POLYNOMIALS * * * This block is the entry point whenever the degree of the * polynomial HAS changed due to root declaration or MAY have * changed due to leading zero elimination in scaling. * * The function of this block is first to eliminate trailing * zeros and then to test the degree to determine whether or not * linear or quadratic factorization is in order. * * This block first eliminates trailing zeros from the polynomial. * If underflow occurred during re-scaling or if underflow in the * trailing deflated coefficient occurred during the evaluation * blocks, then the trailing coefficient has been made zero and * this indicates an underflow zero. * * After elimination of any trailing zero coefficients, the * degree of the polynomial is then tested. If 1 or 2, then * linear or quadratic factorization is in order and PROOT is * exited. Otherwise, the branch is made to Block 4 at the BEGIN0 * entry to begin iteration for a new root at the origin. * * This block is entered from one of the following: * * (1) PROOT Block 2 after leading zero tests subsequent to * (re-) scaling. * (2) PROOT Block 11 after a root has been declared and stored. * * On entry the leading coefficient is NON-ZERO. * * Exit is made only to PROOT Block 4 with CD>2 and the leading * and trailing coefficients are both non-zero. * ***************************************************************** ***************************************************************** * * Here: R4[A]= RSMEM * * The call to TRAL0 below checks for trailing zeros in the * coefficient array and returns underflow roots for each one * found. * ***************************************************************** REENT C=0 S Want trailing zeros to C=C+1 S report underflow GOSUBL TRAL0 Eliminate trailing zeros ***************************************************************** * * Here: R4[A]= RSMEM * CD>0 * P=0, hex mode * * We next check for CD<=2. If so, we solve a linear or quadratic * equation and exit. Otherwise, we exit this block via branch to * BEGIN0 at Block 4. * ***************************************************************** C=R4 D1=C D1-> CD C=DAT1 A C[A]=CD C=C-1 A C[A]=CD-1 C=C-1 A C[A]=CD-2 GOC linear Jump if CD=1 ( bc ) C=C-1 A C[A]=CD-3 GONC BEGIN0 Jump if CD>2 ( bc ) ***************************************************************** * * Here: R4[A]-> RSMEM * CD=2 * P=0, hex mode * * We solve a quadratic equation and exit PROOT. * ***************************************************************** GOSUBL A2A1A0 Put A(2),A(1),A(0) in scratch GOSUBL QDRTC Solve the quadratic GOSUBL STORE Store smaller root GOSUBL FSCR1 (A,B)= Re(larger root) GOSUBL FTHCM+ (R0,R1)= Im(larger root) GOLONG STORE Store larger root and exit ( c ) ***************************************************************** * * Here: R4[A]-> RSMEM * CD=1 * P=0, hex mode * * We solve a linear equation and exit PROOT. * ***************************************************************** linear GOSUBL ->CD,A D1-> CD, A[A]-> A(N), hexmode,P=0 D1=A D1-> A(CD) which is A(1) GOSUBL getcd1 (C,D)= A(1) GOSUBL getab1 (A,B)= A(0) SETDEC Mode for compliment A=-A-1 S (A,B)= -A(0) GOSUBL divf (A,B)= -A(0)/A(1) C=0 W Im R0=C of this R1=C root is 0 GOLONG STORE Store root and exit proot_r ( c ) STITLE PROOT BLOCK 4 ***************************************************************** * * PROOT BLOCK 4 * ------------- * * ANNULUS FOR THE SMALLEST ZERO; ITERATE INITIALIZATION AT THE * ORIGIN * * The function of this block is to initialize the Laguerre * iteration to begin at the origin. * * This block first computes the radii of an annulus about the * origin known to contain the smallest zero of the current * polynomial. The inner radius of the annulus is assigned to * the PROOT variable INR0. It is the Cauchy Lower Bound (CLB) * for the smallest zero, that is, the unique positive zero of * the polynomial (N=CD) * * N * --- k * S(Z) = -|A(0)| + > |A(k)|*Z * --- * k=1 * * The outer radius of the annulus is the minimum of * * (1) the Fejer bound at the origin, * (2) the Laguerre bound at the origin, * (3) the Geometric Mean of the zeros, and * (4) the Cauchy Upper Bound at the origin. * * This block first sets the PROOT variable N to FLOAT(CD). * The Laguerre step (LGSN), its magnitude (ALGSN), and the Fejer * bound (OUTRN) are then assigned values for the origin. OUTRN * is then set to min{OUTRN,SQRT(N)*ALGSN}, the second quantity * being the Laguerre bound. The Geometric Mean of the zeros is * then computed next ((|A(0)/A(N)|)^(1/N)) and OUTRN is set to * 1.00000000001*min{OUTRN, Geometric Mean}. The multiplication * by 1.00000000001 assures a good upper bound despite roundoff. * The CLB is then computed by Newton-Raphson iteration with * initial guess OUTRN. The iteration stops with two iterates * X(I) and X(I+1) bracketing the CLB with X(I+1)>=X(I). We * assign R to .999999999999*X(I) and OUTRN = * min{OUTRN,1.445NX(I+1)}. The Cauchy Upper Bound is given by * CLB/(2^(1/N)-1) so that 1.445NX(I+1) is an overestimate. * * We therefore have an annulus about the origin with radii * INR0 CD * R4[A]= RSMEM * P=0, hex mode * * Set N=FLOAT(CD). * Assign the variables OUTRN, LGSN, and ALGSN values for the * origin. * ***************************************************************** BEGIN0 A=DAT1 A A[A]=CD GOSBVL =BtoEREL (A,B)= FLOAT(CD), decmode GOSUBL SN N= FLOAT(CD), hex, D1->INR0 GOSUBL A2A1A0 Set up coefficients for LAGUER GOSUBL LAGUER Compute LGSN, ALGSN, and OUTRN, D1->N ***************************************************************** * * Here: D1-> N * R4[A]= RSMEM * * The code below continues the OUTRN assignment by setting * OUTRN= 1.00000000001*min{OUTRN, Geometric Mean (GM)}. * ***************************************************************** GOSUBL getab1 (A,B)=N GOSUBL stab2 (R2,R3)=N C=R4 D1=C D1-> CD C=DAT1 A C[A]= CD GOSUBL FTCHAI (A,B)=A(CD), hex GOSUBL stab0 (R0,R1)=A(CD) C=0 A Want A(0) GOSUBL FTCHAI (A,B)=A(0) SETDEC Mode for DIVF * Compute Nth Root (transcription of old cNthRoot subroutine) GOSUBL rccd0 (C,D)=A(CD) GOSUBL divf (A,B)=A(0)/A(CD) A=0 S (A,B)=|A(0)/A(CD)| GOSBVL =aLNF (A,B)=LN(|A(0)/A(CD)|) GOSUBL rccd2 (C,D)=N GOSUBL divf (A,B)=LN(|A(0)/A(CD)|)/N GOSBVL =aEXP15F (A,B)=GM * end of NthRoot GOSUBL GOUTRN OUTRN= min{OUTRN,GM}, dec GOSUBL gtab1- (A,B)= OUTRN C=0 W D=0 W P= 14 D=D+1 P D=D+1 M (C,D)=1.00000000001 GOSUBL multf (A,B)=1.00000000001*OUTRN GOSUBL ptab1- Store new OUTRN ***************************************************************** * * Here: (A,B)= OUTRN * R4[A]= RSMEM * decmode * * The code below calculates the Cauchy Lower Bound, that is, the * unique positive root of the polynomial (N=CD) * * N * --- k * S(Z) = -|A(0)| + > |A(k)|*Z * --- * k=1 * * The method used is Newton-Raphson iteration with initial guess * OUTRN. The following pseudocode is executed: * * B = OUTRN * nriter C = |A(N)| * D = 0 * FOR K = N-1 TO 0 * D = C+B*D * P = |A(K)| * IF K = 0 THEN P = -P * C = P+B*C * NEXT K * Q = B * B = B-C/D * IF B < Q THEN nriter * * Thus, C = S(B) and D = S'(B) where B is the current iterate * Z(I). Note that the code exits when B>=Q, that is, * Z(I+1)>=Z(I). * * Throughout the iteration we will have PSCR1 - C * PSCR2 - D * PSCR3 - B * * D1 will point into the array A(), D0 will point into PSCRATCH, * and R3[A] will be K. * ***************************************************************** GOSUBL SSCR3 B=OUTRN, D1->PSCR4 D0=C D0-> B ***************************************************************** * * Here: D0-> B * R4[A]= RSMEM * * Initialize evaluation loop. * ***************************************************************** nriter GOSUBL ->CD,A D1= RSMEM, A[A]-> A(N), hex C=DAT1 A C[A]= CD= N C=C-1 A C[A]= N-1 (initial K) D1=A D1-> A(N) A=0 W Initialize B=0 W D=0 GOSUBL ptab0- D0->B GOSUBL getab1 (A,B)= A=0 S |A(N)|, D1->A(N-1) D0=D0- 16 D0-> D0=D0- 5 D GOSUBL ptab0- Initialize C=|A(N)| ***************************************************************** * * Here: D0-> D * D1-> A(K) * C[A]= K * R4[A]= RSMEM * PSCR1= C * PSCR2= D * PSCR3= B * * Evaluation loop -- all real arithmetic. * ***************************************************************** evallp R3=C R3[A]=K GOSUBL getab0 (A,B)=D, D0->B SETDEC Mode for MULTF GOSUBL multd0 (A,B)=B*D, D0->pastB GOSUBL D0-63 D0-> C GOSUBL getcd0 (C,D)=C, D0->D GOSUBL raddf (A,B)=C+B*D GOSUBL putab0 D=C+B*D, D0->B GOSUBL getab0 (A,B)=B, D0->pastB GOSUBL D0-63 D0-> C GOSUBL multd0 (A,B)=B*C, D0->D GOSUBL getcd1 (C,D)=A(K), D1->A(K-1) C=0 S (C,D)=P AR3EX A[A]=K ?A#0 A K=0? GOYES k#0 Jump if not ( bc ) C=-C-1 S P=-P k#0 AR3EX Restore K to R3[A] GOSUBL raddf (A,B)=P+B*C GOSUBL ptab0- C=P+B*C, D0->D GOSBVL =ATTNCHK Exit if attn key abort, hex C=R3 C[A]=K C=C-1 A C[A]= NEXT K GONC evallp Loop through K=0 ( bc ) ***************************************************************** * * Here: D0-> D * (A,B)= C * R4[A]= RSMEM * hexmode * * Test for iteration complete. * ***************************************************************** SETDEC Reset mode GOSUBL getcd0 (C,D)=D, D0->B GOSUBL divf (A,B)=C/D A=-A-1 S (A,B)=-C/D GOSUBL getcd0 (C,D)=B, D0->pastB GOSBVL =STCD0 (R0,R1)=B (Q) GOSUBL raddf (A,B)=B-C/D (new B) GOSUBL ptab0- B=B-C/D, D0->pastB D0=D0- 16 D0-> D0=D0- 5 B GOSUBL rccd0 (C,D)=Q GOSUBL prTEST Compare B and Q, dec GOC nrdone Finished if B>=Q ( bc ) GOTO nriter Iterate ***************************************************************** * * Here: D0-> B * (R0,R1)= Q * R4[A]= RSMEM * decmode * * The Newton-Raphson iteration is complete. * Set INR0= .999999999999*Q and OUTRN= min{OUTRN,1.445*N*B}. * Initialize the iteration to begin at the origin: * * OUTR0 := OUTRN * ZN := 0 * FN := |A(0)| * STARTL,SPIRAL := FALSE * ***************************************************************** nrdone GOSUBL rcab0 (A,B)=Q C=0 W C=C-1 A D=0 W D=D-1 M GOSUBL multf (A,B)=.999999999999*Q GOSUBL SINR0 INR0=.999999999999*Q, D1->QSCR0, hex GOSUBL D1-42 D1->N GOSUBL getab1 (A,B)=N, D1->INR0 SETDEC Mode for MULTF GOSUBL multd0 (A,B)=N*B, D0->pastB D=0 W C=0 W P= 11 LCHEX 1445 CDEX W (C,D)=1.445 GOSUBL multf (A,B)=1.445*N*B GOSUBL GOUTRN OUTRN= min{OUTRN,1.445*N*B}, D1->N GOSUBL gtab1- (A,B)=OUTRN GOSUBL SOUTR0 OUTR0:=OUTRN, D1->Re(ZN), hex A=0 W (A,B)= B=0 W 0 GOSUBL putab1 ZN= 0 GOSUBL putab1 D1->LGSN C=0 A Want A(0) GOSUBL FTCHAI (A,B)=A(0) A=0 S (A,B)=|A(0)| GOSUBL SFN FN:= |A(0)| C=R4 D1=C D1= RSMEM D1=D1+ 10 D1-> STARTL C=0 B STARTL,SPIRAL DAT1=C B = FALSE STITLE PROOT BLOCK 5 ***************************************************************** * * PROOT BLOCK 5 * ------------- * * ACCEPTANCE OR REJECTION OF THE LAGUERRE STEP * * The function of this block is to accept (with possible * modification) or reject the Laguerre step. This block is * entered at its only entry point AMRSTP from * * (1) PROOT Block 4 after initialization of the iteration * process to begin at the origin in order to test the * Laguerre step at the origin for acceptability; * * (2) PROOT Block 7 or 8 (evaluation blocks) when the Laguerre * step at the current iterate ZN is to be tested for * acceptability - this test will occur when either the * magnitude of the polynomial has decreased (FN OUTR0, then the current iterate ZN is rejected since * the Laguerre step LGSN is too large, that is, leads outside a * circle known to contain the closest zero. This block is exited * via a branch to PROOT Block 9 at the REJECT entry. Otherwise, * ZN is an acceptable iterate, but the Laguerre step may be * modified. * * If the iterate ZN is accepted and ALGSN <= OUTR0/2, then the * Laguerre step LGSN is not modified and this block is exited via * a branch to PROOT Block 6 at the ACCEPT entry. * * If the iterate ZN is accepted, ALGSN > OUTR0/2, the iterate * is the origin, and INR0 > OUTR0/2, then no modification of the * Laguerre step will occur since "this indicates a narrow annular * region containing the closest zero"; in this case, this block * is exited via branch to PROOT Block 6 at the ACCEPT entry. * * Otherwise, the Laguerre step at ZN is modified to * (OUTR0/2)*(LGSN/ALGSN) and ALGSN is accordingly changed to OUTR0/2; * PROOT Block 6 at ACCEPT entry is then branched to. * * Note that exit from this block is only to PROOT Block 9 * (REJECT entry) if the step is rejected, or PROOT Block 6 * (ACCEPT entry) if the step is accepted. * ***************************************************************** ***************************************************************** * * Here: R4[A]= RSMEM * * Compare OUTR0 and ALGSN. If ALGSN > OUTR0 then we reject the * iterate and branch to PROOT Block 9 at REJECT entry. * ***************************************************************** AMRSTP GOSUBL FOUTR0 (A,B)=OUTR0, hex, D1->ZN GOSUBL stab0 (R0,R1)=OUTR0 GOSUBL D1+105 D1-> ALGSN GOSUBL getcd1 (C,D)=ALGSN, D1->OUTRN GOSBVL =STCD2 (R2,R3)=ALGSN GOSUBL prTEST Compare OUTR0 and ALGSN, dec GOC amr1 Jump if OUTR0 >= ALGSN ( bc ) GOTO REJECT Else, reject the iterate ***************************************************************** * * Here: (R0,R1)= OUTR0 * (R2,R3)= ALGSN * D1-> OUTRN * R4[A]= RSMEM * decmode * * Compare OUTR0/2 and ALGSN. If ALGSN <= OUTR0/2, then we accept * the step as it is and branch to PROOT Block 6 at the ACCEPT entry. ***************************************************************** amr1 GOSUB rcab0 (A,B)=OUTR0 GOSUB div2 (A,B)=OUTR0/2 GOSUBL stab0 (R0,R1)=OUTR0/2 GOSUB rccd2 (C,D)=ALGSN GOSUBL prTEST Compare OUTR0/2 and ALGSN GOC ACCEPT To block 6 if OUTR0/2 >= ALGSN ( bc ) ***************************************************************** * * Here: (R0,R1)= OUTR0/2 * (R2,R3)= ALGSN * R4[A]= RSMEM * decmode * * We have OUTR0 >= ALGSN > OUTR0/2. * If STARTL=FALSE and SPIRAL=FALSE, that is, if the current * iterate is the origin, we compare INR0 and OUTR0/2. * If INR0 > OUTR0/2, then we accept the step as it is and branch to * PROOT Block 6 at the ACCEPT entry. ***************************************************************** C=R4 D1=C D1-> RSMEM D1=D1+ 10 D1-> STARTL C=DAT1 B C[B]= STARTL/SPIRAL ?C#0 B Jump if current GOYES amr2 iterate is not origin ( bc ) GOSUBL FINR0 (A,B)=INR0, hex GOSUB stabcd (C,D)=INR0 GOSUB rcab0 (A,B)=OUTR0/2 GOSUBL prTEST Compare INR0 and OUTR0/2, dec GONC ACCEPT To Block 6 if INR0 > OUTR0/2 ( bc ) ***************************************************************** * * Here: (R0,R1)= OUTR0/2 * (R2,R3)= ALGSN * R4[A]= RSMEM * * We have OUTR0 >= ALGSN > OUTR0/2 >= INR0. * Modify the Laguerre step by setting LGSN:= (OUTR0/2)*(LGSN/ALGSN) * and ALGSN:= OUTR0/2. * Fall through to PROOT Block 6 at ACCEPT entry. ***************************************************************** amr2 GOSUB rcab0 (A,B)=OUTR0/2 GOSUBL SALGSN ALGSN:= OUTR0/2, hex, D1->OUTRN SETDEC Mode for DIVF GOSUB rccd2 (C,D)=ALGSN (old) GOSUB divf (A,B)= OUTR0/2*ALGSN GOSUBL D1-84 D1-> Re(LGSN) GOSUBL MPYUT+ LGSN:= (OUTR0/2)*(LGSN/ALGSN) STITLE PROOT BLOCK 6 ***************************************************************** * * PROOT BLOCK 6 * ------------- * * ACCEPTANCE AND/OR UPDATING OF THE ITERATE * * The function of this block is to accept an iterate and/or * update an iterate and/or determine which polynomial evaluation * block (7 or 8) is to be used. * * This block has three entry points: ACCEPT UPDATE BRANCH * * The function of entry ACCEPT is to accept a Laguerre iterate * and assign the values of certain variables in accordance with * continuing the iteration procedure. These variables contain * information associated with the current (accepted) iterate * with the view that this iterate now becomes the "previous * iterate". These variables are: * * (1) Z0:= ZN (The last accepted iterate is now the current * iterate). * (2) F0:= FN (The magnitude of the polynomial at the last * accepted iterate is now that at the current iterate). * (3) LGS0:= LGSN (The Laguerre step at the last accepted iterate * is now that at the current iterate). * (4) ALGS0:= ALGSN (The magnitude of the Laguerre step at the * last accepted iterate is now that at the current iterate). * (5) OUTR0:= OUTRN (The radius of a circle about the last * accepted iterate known to contain at least one zero of * the polynomial is now that at the current iterate. * (6) STARTL:= TRUE (The Laguerre iterations have now success- * fully started - may, of course, be redundant). * * PROOT Block 6 is only entered at the ACCEPT entry point from * PROOT Block 5 when the Laguerre step at the current iterate is * accepted (after any necessary modification). ACCEPT entry * code always falls through to the UPDATE entry. * * The function of entry UPDATE is to define the current iterate * ZN as ZN:= Z0+LGS0, that is, the current iterate is the previous * iterate plus the Laguerre step at the previous iterate. * * UPDATE is entered from * (1) falling through from entry ACCEPT code, * (2) PROOT Block 9 when the current iterate is rejected and the * Laguerre step at the previous iterate is halved. * UPDATE code only exits by falling through to BRANCH entry. * * The function of entry BRANCH is to determine what polynomial * evaluation code will be used to evaluate and test the current * iterate ZN. According to ZERPOL documentation and common * sense "the time required to evaluate the polynomial and its * derivatives at a real point is less than that at a complex * point so this test forces the iterate to be real whenever * |Im(ZN)| <= 0.2*|LGS0| in order to save time whenever * possible." It also provides cleaner roots when real roots * are expected. If the above test is positive, then PROOT Block * 7 (RITER entry) is branched to. Otherwise, we test for * |Re(ZN)| <= 0.2*|LGS0| and, if true, set Re(ZN) to zero. * This provides cleaner roots when imaginary roots are expected. * We then branch to PROOT Block 8 (CITER entry). * * BRANCH is entered from * (1) falling through from entry UPDATE code, * (2) PROOT Block 10 when the next point in the spiral search * for an initial iterate has been found. * * Exit from BRANCH code is exit from this block. * * Thus, this block is only exited by branching to PROOT Blocks 7 * or 8 (RITER or CITER entry). * ***************************************************************** ***************************************************************** * * BLOCK 6 - ACCEPT ENTRY * * Here: R4[A]= RSMEM * decmode * * Set the logical variable STARTL to TRUE. * Assign the remainder of the desired variables. * * STARTL := TRUE * Z0 := ZN * LGS0 := LGSN * F0 := FN * ALGS0 := ALGSN * OUTR0 := OUTRN * ***************************************************************** ACCEPT C=R4 D1=C D1-> CD D1=D1+ 10 D1-> STARTL C=0 S STARTL:= 1 C=C+1 S (=> STARTL:= DAT1=C S TRUE) P= 0 LC(3) 164 Get offset for Re(Z0) GOSUBL PVARAD D1, C[A]-> Re(Z0), hex, P=0 D0=C D0-> Re(Z0) GOSUBL D1+126 D1=D1+ 16 D1=D1+ 5 D1-> Re(ZN) LCHEX 6 C[0]=6 (loop counter) varlp A=DAT1 A Read Exponent D1=D1+ 5 DAT0=A A Write Exponent D0=D0+ 5 A=DAT1 W Read Mantissa and Sign D1=D1+ 16 DAT0=A W Write Mantissa and Sign D0=D0+ 16 C=C-1 P Decrement counter GONC varlp Loop seven times ( bc ) ***************************************************************** * * BLOCK 6 - UPDATE ENTRY * * Here: R4[A]= RSMEM * * Set ZN := Z0+LGS0 * *********************************************************** UPDATE GOSUBL FREZ0 (A,B)=Re(Z0), D1->Im(Z0), hex GOSUB stabcd (C,D)=Re(Z0) GOSUBL getab1 (A,B)=Im(Z0), D1->Re(LGS0) GOSUB stab2 (R2,R3)=Im(Z0) GOSUBL CFTCD1 (A,B),(R0,R1)=LGS0, D1-> F0 SETDEC Mode for ADDF GOSUB raddf (A,B)=Re(Z0+LGS0) GOSUBL exab0 (A,B)=Im(LGS0) GOSUB rccd2 (C,D)=Im(Z0) GOSUB raddf (A,B)=Im(Z0+LGS0) GOSUBL exab0 (A,B),(R0,R1)=Z0+LGS0 GOSUBL D1+63 D1->Re(ZN) GOSUBL CSTOD1 ZN:=Z0+LGS0 ***************************************************************** * * BLOCK 6 - BRANCH ENTRY * * Here: R4[A]= RSMEM * * Compare |Im(ZN)| and 0.2*ALGS0. * If |Im(ZN)| <= 0.2*ALGS0, then branch to PROOT Block 7 at RITER. * Otherwise, compare |Re(ZN)| and 0.2*ALGS0. * If |Re(ZN)| <= 0.2*ALGS0, then set Re(ZN) = 0. * Branch to PROOT Block 8 at CITER entry. ***************************************************************** BRANCH GOSUBL FALGS0 (A,B)=ALGS0, D1->G, hex SETDEC Mode for decrement C=0 W C=C-1 A D=0 W P= 14 D=D+1 P D=D+1 P (C,D)=0.2 GOSUB multf (A,B)=0.2*ALGS0 GOSUBL stab0 (R0,R1)=0.2*ALGS0 GOSUBL D1+42 D1->Im(ZN) GOSUBL getcd1 (C,D)=Im(ZN), D1->LGSN C=0 S (C,D)=|Im(ZN)| GOSUB prTEST Compare |Im(ZN)| and 0.2*ALGS0 GOC RITER To RITER if |Im(ZN)|<=0.2*ALGS0 ( bc ) GOSUB rcab0 (A,B)=0.2*ALGS0 GOSUBL D1-42 D1->Re(ZN) GOSUBL getcd1 (C,D)=Re(ZN), D1->Im(ZN) C=0 S (C,D)=|Re(ZN)| GOSUB prTEST Compare |Re(ZN)| and 0.2*ALGS0 GONC citer jump if |Re(ZN)|>0.2*ALGS0 ( bc ) A=0 W (A,B) B=0 W =0 GOSUBL ptab1- Re(ZN)=0, D1->Im(ZN) citer D1=D1+ 16 D1=D1+ 5 D1->LGSN GOTO CITER ( c ) STITLE PROOT BLOCK 7 ***************************************************************** * * PROOT BLOCK 7 * ------------- * * EVALUATION AT A REAL ITERATE * * The function of this block is to evaluate the polynomial, its * derivatives, the Laguerre step, and the Fejer bound at the real * iterate ZN. * * This block is entered (at its only entry point RITER) only * from PROOT Block 6 when the current iterate ZN is declared to * be "real". * * This block executes the following sequence: * * (1) Compute the value of the polynomial, its first and second * derivatives, and the deflated (by a linear factor) * coefficients of the quotient polynomial. If at any time * during the evaluation process the above values overflow * the extended range, then this block is exited via a branch * to PROOT Block 2 at the RESCAL entry to re-scale the * coefficients. This would cause the iterations to restart * at the origin with the re-scaled coefficients. * If an attention key abort condition is detected during the * evaluation process, then this block is exited via ABORT. * * (2) Assign the variable FN the value |P(ZN)|. * * (3) Compute and assign the variable E, the bound for the * rounding error in the polynomial value. * * (4) If FN <= E, then branch to PROOT Block 11 at the ZERO entry * to store a real zero at the iterate ZN. * * (5) Else, if FN >= F0 and STARTL=TRUE, then branch to PROOT * Block 9 at the REJECT entry to reject the iterate. * * (6) Else, compute and assign the variables LGSN (Laguerre step), * ALGSN (Laguerre step magnitude), and OUTRN (Fejer bound) * values for the current iterate. * * (7) If |ZN|+ALGSN = |ZN| then branch to PROOT Block 11 at * the ZERO entry to declare a real zero at ZN. * * (8) Else, branch to PROOT Block 5 at the AMRSTP entry to accept, * modify, or reject the current Laguerre step. * * Note that steps (3)-(8) are executed by the ITRUTI utility, * code which is shared by PROOT Block 8. * * The polynomial evaluation at a real iterate executes the * following pseudocode: * * 5 ! LET N=CD AND X=Re(ZN) * 10 Im(ZN) = 0 ! FORCES A REAL ITERATE * 20 D = 0 * 30 C = A(N) * 40 DFLTMP(N) = A(N) * 50 DFLTMP(N-1) = A(N-1) + A(N)X * 60 FOR K = N-2 TO 0 * 70 D = C + DX * 80 C = DFLTMP(K+1) + CX * 90 DFLTMP(K) = A(K) + DFLTMP(K+1)X * 100 NEXT K * 110 P(ZN) = (DFLTMP(0),0) * 120 P'(ZN) = (C,0) * 130 P''(ZN)/2 = (D,0) * 140 FN = |P(ZN)| * * Throughout the evaluation process (lines 60-100) the following * machine elements are dedicated: * * D1 -> DFLTMP() * D0 -> PSCRATCH * (R0,R1) = X * R2[A] = K * R3[A] -> A(K) * PSCR1 = C * PSCR2 = D * * Here: D1-> Re(LGSN) * * Initialize the evaluation loop (lines 60-100) by executing * lines 10-50 of the above pseudocode. * * X = Re(ZN) * Im(ZN) = 0 * D = 0 * C = A(N) * DFLTMP(N) = A(N) * K = N-1 * DFLTMP(N-1) = A(N-1) + A(N)X * ***************************************************************** RITER GOSUBL D1-42 D1-> Re(ZN) GOSUBL getab1 (A,B)=Re(ZN), D1->Im(ZN) GOSUBL stab0 (R0,R1)=Re(ZN) (X) A=0 W (A,B) B=0 W =0 GOSUBL putab1 Im(ZN)=0 GOSUBL SSCR2 PSCR2(D):=0, C[A]-> PSCR2, hex D0=C D0->PSCR2(D) GOSUB ->CD,A D1-> CD=N, A[A]-> A(N) C=DAT1 A C[A]=N C=C-1 A R2[A]=N-1 R2=C (initial K) C=C+1 A C[A]=N D1=A D1-> A(N) GOSUBL getab1 (A,B)=A(N), D1->A(N-1) CD1EX C[A]->A(N-1), D1=N R3=C R3[A]->A(N-1) CD1EX C[A]=N GOSUBL ptab0- PSCR1(C):= A(N), D0->PSCR2 GOSUBL DFLST DFLTMP(N)=A(N) GOTO rlpent Jump into evaluation loop ***************************************************************** * * Here: D1 -> DFLTMP(K) * D0-> PSCR2 * (R0,R1) = X * R2[A] = K * R3[A] -> A(K) * PSCR1 = C * PSCR2 = D * * Execution the evaluation loop at a real iterate. * * D = C + DX * C = DFLTMP(K+1) + CX * ***************************************************************** ritrlp GOSUBL getab0 (A,B)=D, D0-> PSCR3 SETDEC Mode for MULTF GOSUB multr01 (A,B)=D*X GOSUBL D0-42 D0->C GOSUBL getcd0 (C,D)=C, D0->D GOSUB raddf (A,B)= C+D*X GOSUB prtPAK Check for extended overflow GOC ritscl Exit block if so ( Fc ) GOSUBL putab0 D:= C+D*X, D0-> PSCR3 GOSUBL D0-42 D0->C GOSUBL getab0 (A,B)=C, D0->D GOSUB multr01 (A,B)=C*X D1=D1- 16 D1-> D1=D1- 5 DFLTMP(K+1) GOSUBL getcd1 (C,D)=DFLTMP(K+1), D1->DFLTMP(K) GOSUB raddf (A,B)=DFLTMP(K+1)+C*X GOSUB prtPAK Check for extended overflow GONC ritr+ Jump if not ( Fc ) ritscl GOTO RESCAL Exit block to RESCAL ( nc ) ritr+ GOSUBL ptab0- C:= DFLTMP(K+1)+C*X, D0->D ***************************************************************** * * Here: D1 -> DFLTMP(K) * D0 -> PSCR2 * (R0,R1) = X * R2[A] = K * R3[A] -> A(K) * PSCR1 = C * PSCR2 = D * * Real evaluation loop (continued). * * DFLTMP(K) = A(K) + DFLTMP(K+1)X * ***************************************************************** rlpent GOSUBL gtab1- (A,B)=DFLTMP(K+1), D1->DFLTMP(K) SETDEC Mode for evaluation loop GOSUB multr01 (A,B)=DFLTMP(K+1)*X AR3EX A[A]->A(K) AD1EX D1->A(K), A[A]->DFLTMP(K) GOSUBL getcd1 D1->A(K-1), (C,D)=A(K) AD1EX A[A]->A(K-1), D1->DFLTMP(K) AR3EX R3[A]->A(K-1) GOSUB raddf (A,B)=A(K)+DFLTMP(K+1)*X GOSUB prtPAK Check for extended overflow GOC ritscl Exit if so ( Fc ) GOSUBL putab1 DFLTMP(K):= A(K)+DFLTMP(K+1)*X, D1->DFLTMP(K-1) GOSBVL =ATTNCHK Check attn key abort, hex C=R2 C[A]=K C=C-1 A C[A]=NEXT K R2=C Into R2[A] GOC ritrdn Loop ( bc ) GOTO ritrlp through K=0 ***************************************************************** * * Here: (A,B) = DFLTMP(0) * D0-> PSCR2 * PSCR1 = C * PSCR2 = D * hexmode * * Set up the PROOT scratch area for ITRUTI, that is, compute * P(ZN), P'(ZN), and P''(ZN)/2. ***************************************************************** ritrdn C=0 W All imaginary parts R0=C are 0 so set R1=C (R0,R1)=0 GOSUB stab2 (R2,R3)=DFLTMP(0) GOSUBL getab0 (A,B)=D, D0-> PSCR3 GOSUB D0-42 D0-> C GOSUBL getcd0 (C,D)=C, D0-> PSCR2 D0=D0- 16 D0=D0- 5 D0-> PSCR1 GOSUB CSTOD0 P''(ZN)/2 = (D,0) GOSUB stcdab (A,B)=C GOSUB CSTOD0 P'(ZN) = (C,0) GOSUB rcab2 (A,B) = DFLTMP(0) GOSUB CSTOD0 P(ZN)= (DFLTMP(0),0) ***************************************************************** * * Here: (A,B) = P(ZN) * R4[A] -> RSMEM * PSCR1 = Re(P''(ZN)/2) * PSCR2 = Im(P''(ZN)/2) * PSCR3 = Re(P'(ZN) * PSCR4 = Im(P'(ZN) * PSCR5 = Re(P(ZN) * PSCR6 = Im(P(ZN) * * Set FN=|P(ZN)| and branch to ITRUTI to complete processing for * the real iterate evaluation. ***************************************************************** A=0 S (A,B)=|P(ZN)| GOSUBL SFN FN:= |P(ZN)| , hex, D1->ALGSN ST=0 sCITER Flag RITER call for ITRUTI GOLONG ITRUTI Finish iterate processing STITLE PROOT BLOCK 8 ***************************************************************** * * PROOT BLOCK 8 * ------------- * * EVALUATION AT A COMPLEX ITERATE * * The function of this block is to evaluate the polynomial, its * derivatives, the Laguerre step, and the Fejer bound at the * complex iterate ZN. * * This block is entered (at its only entry point CITER) only * from PROOT Block 6 when the current iterate ZN is declared to * be complex. * * This block executes the following sequence: * * (1) Compute the value of the polynomial, its first and second * derivatives, and the deflated (by a quadratic factor) * coefficients of the quotient polynomial. If at any time * during the evaluation process the above values overflow * the extended range, then this block is exited via a branch * to PROOT Block 2 at the RESCAL entry to re-scale the * coefficients. If an attention key abort condition is * detected during the evaluation process, then this block is * exited via ABORT. * * (2) Assign the variable FN the value |P(ZN)|. * * (3) Compute and assign the variable E, the bound for the * rounding error in the polynomial value. * * (4) If FN <= E, then branch to PROOT Block 11 at the ZERO * entry to store a complex zero at the iterate ZN. * * (5) Else, if FN >= F0 and STARTL=TRUE, then branch to PROOT * Block 9 at the REJECT entry to reject the iterate. * * (6) Else, compute and assign the variables LGSN (Laguerre step), * ALGSN (Laguerre step magnitude), and OUTRN (Fejer bound) * values for the current iterate. * * (7) If |ZN|+ALGSN = |ZN| then branch to PROOT Block 11 at * the ZERO entry to declare a complex zero at ZN. * * (8) Else, branch to PROOT Block 5 at the AMRSTP entry to accept, * modify, or reject the current Laguerre step. * * Note that steps (3)-(8) are executed by the ITRUTI utility, * code which is shared by PROOT Block 7. * * The polynomial evaluation at a complex iterate executes the * following pseudocode: * * 5 ! LET N=CD, X=Re(ZN), AND Y=Im(ZN) * 10 R = X^2+Y^2 * 20 B = 2X * 30 D1 = D0 = C1 = 0 * 40 C0 = A(N) * 50 DFLTMP(N) = A(N) * 60 FOR K = N-1 TO 0 * 70 IF K = 0 THEN B = X * 80 IF K > N-3 THEN 150 * 90 V = D1*R @ D1 = D0 * 100 P = B*D0-V * 110 D0 = C1+P * 120 V = C1*R @ C1 = C0 * 130 P = B*C0-V * 140 C0 = DFLTMP(K+2)+P * 150 P = A(K)+B*DFLTMP(K+1) * 160 IF K = N-1 THEN 180 * 170 P = P-R*DFLTMP(K+2) * 180 DFLTMP(K) = P * 190 NEXT K * 200 P(ZN) = (DFLTMP(0),Y*DFLTMP(1)) * 210 P'(ZN) = (DFLTMP(1)-2*Y^2*C1,2*Y*C0) * 220 P''(ZN)/2 = (C0-4*Y^2*D0,Y*(3*C1-4*Y^2*D1)) * 230 FN = |P(ZN)| * * Throughout the evaluation process (lines 60-190) the following * machine elements are dedicated: * * D1 -> DFLTMP() * D0 -> PSCRATCH * R2[A] = K * R3[A] -> A() * RSTK= N-2 * R4[A] = RSMEM * PSCR1 = D0 * PSCR2 = D1 * PSCR3 = C0 * PSCR4 = C1 * PSCR5 = B * PSCR6 = R * * Here: D1-> Re(LGSN) * R4[A]= RSMEM * decmode * * Initialize the evaluation loop (lines 60-190) by executing * lines 10-50 of the above pseudocode. * * 10 R = X^2+Y^2 * 20 B = 2X * 30 D1 = D0 = C1 = 0 * 40 C0 = A(N) * 50 DFLTMP(N) = A(N) ***************************************************************** CITER GOSUB D1-42 D1->Re(ZN) GOSUB CFTCD1 (A,B),(R0,R1)=ZN, D1->LGSN GOSUB stab2 (R2,R3)=X GOSUBL CpMAG^2 (A,B)=R, trashes R0,R1 GOSUBL stab0 (R0,R1)=R GOSUB rcab2 (A,B)=X GOSUB mult2 (A,B)=2X (B) GOSUBL SSCR5 Put B in scratch, hex, D1->PSCR6 D0=C D0->PSCR5 GOSUB STOCM+ Put R in scratch GOSUB ->CD,A D1= RSMEM, A[A]-> A(N) C=DAT1 A C[A]=N D1=A D1-> A(N) A=C A A[A]=N C=C-1 A C[A]=N-1 R2=C R2[A]=N-1 (initial K) C=C-1 A C[A]=N-2 RSTK=C RSTK: N-2 C=A A C[A]=N GOSUB getab1 (A,B)=A(N), D1->A(N-1) AD1EX A[A]->A(N-1) R3=A Initialize R3[A]->A(N-1) AD1EX Restore A[A] GOSUBL DFLST DFLTMP(N)=A(N) GOSUB D0-84 D0-> PSCR1 C=0 W Prepare to zero out scratch DAT0=C W D0=D0+ 16 DAT0=C W D0=D0+ 16 DAT0=C W D0=D0+ 10 Initialize D0,D1=0 GOSUB putab0 Initialize C0=A(N), D0->C1 DAT0=C W D0=D0+ 16 DAT0=C A D0=D0+ 5 Initialize C1=0 GOSUB D0+42 D0-> past R for evaluation loop ***************************************************************** * * Here: D1 -> DFLTMP(K) * D0 -> past R * R2[A] = K * R3[A] -> A(K) * RSTK= N-2 * R4[A] = RSMEM * PSCR1,2 = D0, D1 * PSCR3,4 = C0, C1 * PSCR5 = B * PSCR6 = R * hexmode * * Execute lines 70-80 of the pseudocode. * * 60 FOR K = N-1 TO 0 * 70 IF K = 0 THEN B = X * 80 IF K > N-3 THEN 150 ***************************************************************** citrlp A=R2 A[A]=K ?A#0 A K=0? GOYES citr1 Jump if not ( bc ) CD1EX Save D1 in R0=C R0[A] GOSUBL FREZN (A,B)=X C=R0 Restore D1=C D1 GOSUB D0-42 D0->B GOSUB putab0 B:=X, D0->R GOTO citr3 Skip next test ( c ) citr1 C=RSTK C[A]=N-2 RSTK=C ?AB D1=D1- 16 D1-> DFLTMP(K+1) D1=D1- 5 GOTO citr7 ( c ) citr2 D0=D0- 16 D0->R D0=D0- 5 ***************************************************************** * * Here: D1 -> DFLTMP(K) * D0 -> R * R2[A] = K * R3[A] -> A(K) * RSTK= N-2 * R4[A] = RSMEM * PSCR1,2 = D0, D1 * PSCR3,4 = C0, C1 * PSCR5 = B * PSCR6 = R * * Execute lines 90-110 of the pseudocode. * * 90 V = D1*R @ D1 = D0 * 100 P = B*D0-V * 110 D0 = C1+P ***************************************************************** citr3 GOSUB getcd0 (C,D)=R, D0-> past R GOSUB D0-126 D0-> D0 (PSCR1) GOSUB CFTCD0 (A,B)=D0, (R0,R1)=D1, D0->C0 GOSUB ptab0- D1:= D0, D0-> C0 SETDEC Mode for MULTF GOSUB exab0 (A,B)=D1, (R0,R1)=D0 GOSUB multf (A,B)=V GOSUB exab0 (A,B)=D0, (R0,R1)=V GOSUB D0+42 D0->B GOSUB multd0 (A,B)=B*D0, D0->R A=-A-1 S GOSUB raddr01 (A,B)=-P A=-A-1 S (A,B)=P GOSUB D0-42 D0->C1 GOSUB getcd0 (C,D)=C1, D0->B GOSUB raddf (A,B)=C1+P GOSUB D0-84 D0->D0 GOSUB prtPAK Check for extended range overflow GONC citr4 Jump if not ( Fc ) GOTO RESCAL Else, rescale citr4 GOSUB putab0 D0:= C1+P ***************************************************************** * * Here: D1 -> DFLTMP(K) * D0 -> D1 * R2[A] = K * R3[A] -> A(K) * RSTK= N-2 * R4[A] = RSMEM * PSCR1,2 = D0, D1 * PSCR3,4 = C0, C1 * PSCR5 = B * PSCR6 = R * decmode * * Execute lines 120-140 of the pseudocode. * * 120 V = C1*R @ C1 = C0 * 130 P = B*C0-V * 140 C0 = DFLTMP(K+2)+P ***************************************************************** D0=D0+ 16 D0-> C0 D0=D0+ 5 GOSUB CFTCD0 (A,B)=C0, (R0,R1)=C1, D0->B GOSUB ptab0- C1:= C0, D0->B GOSUB exab0 (A,B)=C1, (R0,R1)=C0 D0=D0+ 16 D0-> R D0=D0+ 5 GOSUB multd0 (A,B)= C1*R (V), D0-> past R GOSUB exab0 (A,B)=C0, (R0,R1)=V GOSUB D0-42 D0-> B GOSUB multd0 (A,B)=B*C0, D0->R A=-A-1 S GOSUB raddr01 (A,B)=-P A=-A-1 S (A,B)=P GOSUB D1-42 D1-> DFLTMP(K+2) GOSUB getcd1 (C,D)=DFLTMP(K+2), D0->DFLTMP(K+1) GOSUB raddf (A,B)=DFLTMP(K+2)+P GOSUB D0-63 D0->C0 GOSUB prtPAK Check for extended range overflow GONC citr5 Jump if not ( Fc ) GOTO RESCAL Else, rescale citr5 GOSUB putab0 C0:= DFLTMP(K+2)+P, D0->C1 D0=D0+ 16 D0->B D0=D0+ 5 ***************************************************************** * * Here: D1 -> DFLTMP(K+1) * D0 -> B * R2[A] = K * R3[A] -> A(K) * RSTK= N-2 * R4[A] -> RSMEM * PSCR1,2 = D0, D1 * PSCR3,4 = C0, C1 * PSCR5 = B * PSCR6 = R * * Execute lines 150-190 of the pseudocode. * * 150 P = A(K)+B*DFLTMP(K+1) * 160 IF K = N-1 THEN 180 * 170 P = P-R*DFLTMP(K+2) * 180 DFLTMP(K) = P * 190 NEXT K ***************************************************************** citr7 GOSUB getab1 (A,B)=DFLTMP(K+1), D1->DFLTMP(K) SETDEC Mode for rest GOSUB multd0 (A,B)=B*DFLTMP(K+1), D0->R AR3EX R3=saved A, A[A]->A(K) AD1EX A[A]= saved D1, D1-> A(K) GOSUB getcd1 (C,D)=A(K), D1->A(K-1) AD1EX D1 restored, A[A]-> A(K-1) AR3EX A restored, R3[A]->A(K-1) GOSUB raddf (A,B)=P GOSUB getcd0 (C,D)=R, D0-> past R GOSBVL =STCD0 (R0,R1)=R C=RSTK C[A]= N-2 RSTK=C D=C A D[A]=N-2 C=R2 C[A]=K ?C>D A K = N-1? GOYES citr8 Jump if so ( bc ) GOSUB exab0 (R0,R1)=P, (A,B)=R GOSUB D1-42 D1->DFLTMP(K+2) GOSUB getcd1 (C,D)=DFLTMP(K+2), D1-> DFLTMP(K+1) D1=D1+ 16 D1-> DFLTMP(K) D1=D1+ 5 GOSUB multf (A,B)=R*DFLTMP(K+2) A=-A-1 S (A,B)=-R*DFLTMP(K+2) GOSUB raddr01 P:= P-R*DFLTMP(K+2) citr8 GOSUB prtPAK Check for extended range overflow GONC citr9 Jump if not ( Fc ) GOTO RESCAL Else, rescale citr9 GOSUB putab1 DFLTMP(K)=P, D1-> DFLTMP(K-1) GOSBVL =ATTNCHK Exit if attention key abort, hex C=R2 C[A]=K C=C-1 A C[A]=NEXT K R2=C R2[A]=K-1 GOC citr10 Done if decremented zero ( bc ) GOTO citrlp Loop through K=0 ***************************************************************** * * Here: (A,B)= DFLTMP(0) * D1 -> past DFLTMP(0) * D0 -> past R * R4[A] = RSMEM * PSCR1,2 = D0, D1 * PSCR3,4 = C0, C1 * PSCR5 = B * PSCR6 = R * hexmode * * Execute lines 200-220 of the pseudocode, computing P(ZN), * P'(ZN), and P''(ZN), thus setting up the scratch area for * ITRUTI. * * 200 P(ZN) = (DFLTMP(0),Y*DFLTMP(1)) * 210 P'(ZN) = (DFLTMP(1)-2*Y^2*C1,2*Y*C0) * 220 P''(ZN)/2 = (C0-4*Y^2*D0,Y*(3*C1-4*Y^2*D1)) * 230 FN = |P(ZN)| * * P(ZN) = (DFLTMP(0),Y*DFLTMP(1)) ***************************************************************** citr10 C=RSTK unnecessary ??? GOSUB D0-42 D0-> PSCR5 GOSUB putab0 Re(P(ZN))=DFLTMP(0), D0->PSCR6 GOSUB D1-42 D1->DFLTMP(1) GOSUB getab1 (A,B)=DFLTMP(1), D1->DFLTMP(0) GOSUBL stab0 (R0,R1)=DFLTMP(1) GOSUB FIMZN (A,B)=Y GOSUB stab2 (R2,R3)=Y GOSUB SQSCR0 (QSCR0)= Y, hex GOSUB STOCM+ (QSCR1)= DFLTMP(1), D1-> QSCR2 SETDEC Mode for MULTF GOSUB multr01 (A,B)=Y*DFLTMP(1) GOSUB putab0 Im(P(ZN))=Y*DFLTMP(1), D0->past PSCR6 ***************************************************************** * * Here: (R2,R3)= Y * D0 -> past PSCR6 * R4[A] = RSMEM * PSCR1,2 = D0, D1 * PSCR3,4 = C0, C1 * PSCR5 = DFLTMP(0) (Re(P(ZN)) * PSCR6 = Y*DFLTMP(1) (Im(P(ZN)) * QSCR0 = Y * QSCR1 = DFLTMP(1) * decmode * * P''(ZN)/2 = (C0-4*Y^2*D0,Y*(3*C1-4*Y^2*D1)) ***************************************************************** GOSUB rcab2 (A,B)=Y GOSUB x^2 (A,B)=Y^2 GOSUB stab2 (R2,R3)=Y^2 GOSUB mult2 (A,B)=2*Y^2 GOSUB mult2 (A,B)=4*Y^2 A=-A-1 S (A,B)=-4*Y^2 GOSUB D0-126 D0->PSCR1 (D0) CD0EX D0=C D1=C D0,D1->PSCR1 GOSUB MPYUT+ Replace (D0,D1) by -4*Y^2*(D0,D1) GOSUB getab1 (A,B)=C0, D1-> C1 GOSUB getcd0 (C,D)=-4*Y^2*D0, D0-> PSCR2 GOSUB raddf (A,B)=C0-4*Y^2*D0 GOSUB ptab0- Put Re(P''(ZN)/2), D0->PSCR2 GOSUB getab1 (A,B)=C1, D1-> PSCR5 C=0 W D=0 W P= 14 LCHEX 3 CDEX P (C,D)=3 GOSUB multf (A,B)=3*C1 GOSUB getcd0 (C,D)=-4*Y^2*D1, D0->PSCR3 GOSUB raddf (A,B)=3*C1-4*Y^2*D1 GOSUB QRCL0 (C,D)=Y, D0->QSCR1, dec GOSUB multf (A,B)=Y*(3*C1-4*Y^2*D1) CD1EX CD0EX D0->PSCR5 CD1EX D1->QSCR1 GOSUB D0-42 D0->PSCR3 GOSUB ptab0- Put Im(P''(ZN)/2), D0->PSCR3 ***************************************************************** * * Here: (R2,R3)= Y^2 * D0 -> PSCR3 * D1 -> QSCR1 * R4[A] = RSMEM * PSCR1 = C0-4*Y^2*D0 (Re(P''(ZN)/2)) * PSCR2 = Y*(3*C1-4*Y^2*D1) (Im(P''(ZN)/2)) * PSCR3,4 = C0, C1 * PSCR5 = DFLTMP(0) (Re(P(ZN)) * PSCR6 = Y*DFLTMP(1) (Im(P(ZN)) * QSCR0 = Y * QSCR1 = DFLTMP(1) * decmode * * * P'(ZN) = (DFLTMP(1)-2*Y^2*C1,2*Y*C0) ***************************************************************** GOSUB rcab2 (A,B)=Y^2 GOSUB mult2 (A,B)=2*Y^2 A=-A-1 S (A,B)=-2*Y^2 D0=D0+ 16 D0-> PSCR4 D0=D0+ 5 GOSUB multd0 (A,B)= -2*Y^2*C1, D0->PSCR5 GOSUB getcd1 (C,D)= DFLTMP(1), D1->QSCR2 GOSUB raddf (A,B)= DFLTMP(1)-2*Y^2*C1 GOSUB D0-42 D0-> PSCR3 GOSUB getcd0 (C,D)= C0, D0-> PSCR4 GOSUB ptab0- Put Re(P'(ZN)), D0-> PSCR4 GOSUB D1-42 D1->QSCR0 (Y) GOSUB getab1 (A,B)=Y, D1->QSCR1 GOSUB multf (A,B)=Y*C0 GOSUB mult2 (A,B)=2*Y*C0 GOSUB putab0 Push Im(P'(ZN)), D0->PSCR5 ***************************************************************** * * Here: D0 -> PSCR5 * R4[A] = RSMEM * PSCR1 = Re(P''(ZN)/2) * PSCR2 = Im(P''(ZN)/2) * PSCR3 = Re(P'(ZN) * PSCR4 = Im(P'(ZN) * PSCR5 = Re(P(ZN) * PSCR6 = Im(P(ZN) * QSCR0 = Y * QSCR1 = DFLTMP(1) * decmode * * Set FN:= |P(ZN)| and branch to ITRUTI to complete processing * for the complex iterate evaluation. ***************************************************************** GOSUB CFTCD0 (A,B),(R0,R1)= P(ZN), D0->past PSCR6 GOSUBL CpABS (A,B)= |P(ZN)|, trashes R0,R1 GOSUB SFN FN:= |P(ZN)| , hex, D1-> ALGSN ST=1 sCITER Flag CITER call for ITRUTI GOLONG ITRUTI Finish iterate processing ( c ) STITLE PROOT BLOCK 9 ***************************************************************** * * PROOT BLOCK 9 * ------------- * * REJECTION OF AN ITERATE * * The function of this block is to deal with a rejected iterate. * This block is entered at its only entry point REJECT from * (1) PROOT Block 5 when the Laguerre step LGSN at the iterate * ZN is deemed to be too large * (2) The evaluation blocks (7 and 8) when the magnitude of the * polynomial has not decreased after the iterations have * successfully started. * * If the Laguerre iterations have not yet started (STARTL=FALSE), * then the branch is made to PROOT Block 10 at the SEARCH entry * where the annulus INR0 <= |Z| <= OUTR0 containing the smallest * zero about the origin is searched for an acceptable starting point. * * Otherwise, a successful Laguerre step from the origin has occurred * previously, and the previous Laguerre step LGS0 is halved. * According to ZERPOL documentation, "In the case that FN >= F0, * the reduction of the step LGS0 is an attempt to reach a point * along the direction of the Laguerre step from Z0 where the * magnitude of the polynomial decreases. In the case that LGSN * is too large, the reduction of the step LGS0 is an attempt to * make the Laguerre step from ZN shorter by placing ZN closer to * an acceptable iterate Z0." * * After the previous step size is cut in half, it is tested to * determine whether or not it is negligible (perhaps halved too * many times) with respect to ZN. If not, then we branch to * PROOT Block 6 at the UPDATE entry and the new iterate becomes * ZN:= Z0+LGS0/2. FN and LGSN are then computed and tested at * this new iterate. * * If the step size has become negligible, then PROOT has * problems. However, by ZERPOL directive, "we are willing to * tolerate the error in a zero caused by accepting ZN if * FN < E*N^2." In this case, the branch is made to PROOT Block * 11 at the ZERO entry to accept ZN as a root. Otherwise, PROOT * has failed, an error is reported, and PROOT aborts. * * Thus, possible exits are to * (1) PROOT Block 10 at SEARCH entry if STARTL=FALSE, * (2) PROOT Block 6 at UPDATE entry if the halved step size is * not negligible, * (3) PROOT Block 11 at ZERO entry if the halved step size is * negligible but FN < E*N^2, * (4) Fatal error exit otherwise. * ***************************************************************** ***************************************************************** * * Here: R4[A]= RSMEM * * Test the logical variable STARTL. * If STARTL=FALSE, then branch to PROOT Block 10 at SEARCH entry. ***************************************************************** REJECT C=R4 D1=C D1= RSMEM D1=D1+ 10 D1-> STARTL C=DAT1 S C[S]=STARTL ?C#0 S STARTL=TRUE? GOYES rject1 Jump if so ( bc ) GOTO SEARCH Else, branch to SEARCH ( c ) ***************************************************************** * Set LGS0:= (1/2)LGS0 and ALGS0:= (1/2)ALGS0. ***************************************************************** rject1 SETDEC Mode for decrement A=0 W A=A-1 A B=0 W P= 14 LCHEX 5 B=C P (A,B)=1/2 P= 0 LC(3) 206 C[X]= offset for Re(LGS0) GOSUB MPYUTI LGS0:= (1/2)LGS0, (R0,R1) = 1/2 D1=D1+ 16 D1-> ALGS0 D1=D1+ 5 GOSUB getab1 (A,B)=ALGS0 GOSUB multr01 (A,B)=(1/2)ALGS0 GOSUB ptab1- ALGS0:= (1/2)ALGS0 ***************************************************************** * * Here: (A,B)= new ALGS0 * D1-> G * R4[A]= RSMEM * decmode * * Compare |ZN|+ALGS0 and |ZN|; if unequal, then branch to * PROOT Block 6 at UPDATE entry. ***************************************************************** GOSUB stab2 (R2,R3)=ALGS0 D1=D1+ 16 D1-> Re(ZN) D1=D1+ 5 GOSUB CFTCD1 (A,B),(R0,R1)= ZN, D1->LGSN GOSUBL CpABS (A,B)= |ZN|, trashes R0,R1 GOSUB stab0 (R0,R1)= |ZN| GOSUB rccd2 (C,D)= ALGS0 GOSUB raddf (A,B)= |ZN|+ALGS0 GOSUB rccd0 (C,D)= |ZN| GOSUB prTEST Compare |ZN| and |ZN|+ALGS0 ?B=0 WP |ZN| = |ZN|+ALGS0? GOYES rject2 Jump if step size negligible ( Fc ) GOTO UPDATE Else, branch to UPDATE ( c ) ***************************************************************** * * Here: D1-> Re(LGSN) * R4[A]= RSMEM * decmode * * Since the step size has become negligible, compare FN and * E*N^2. If FN < E*N^2, then branch to PROOT Block 11 to accept * ZN as a zero. Otherwise, take an error exit. ***************************************************************** rject2 GOSUB D1+42 D1->FN GOSUB getab1 (A,B)=FN, D1->ALGSN GOSUB stab0 (R0,R1)=FN GOSUB D1+42 D1-> N GOSUB getab1 (A,B)=N, D1->INR0 GOSUB x^2 (A,B)=N^2 GOSUB stab2 (R2,R3)=N^2 GOSUB FE (A,B)=E, hex, D1->PSCR1 SETDEC Mode for multf GOSUB rccd2 (C,D)=N^2 GOSUB multf (A,B)=E*N^2 GOSUB rccd0 (C,D)=FN GOSUB xyex (A,B)=FN, (C,D)=E*N^2 GOSUB prTEST Compare FN and E*N^2 GOC error Jump if FN >= E*N^2 ( nc ) GOTO ZERO Else, accept ZN as a zero error P= 0 LCHEX 0C001 "Unable to find roots" GOVLNG =GPErrjmpC STITLE PROOT BLOCK 10 ***************************************************************** * * PROOT BLOCK 10 * -------------- * * SPIRAL SEARCH FOR AN INITIAL ITERATE * * The function of this block is to search the annulus * {Z: INR0 <= |ZN| <= OUTR0} for an acceptable starting point * for the Laguerre iteration. * * This block is entered at its only entry point SEARCH from * PROOT Block 9 when the current iterate has been rejected but * the Laguerre iterations have not yet successfully begun * (STARTL=FALSE). * * This block first tests the logical variable SPIRAL. If * SPIRAL=TRUE then the spiral search has already begun; * otherwise, we set SPIRAL=TRUE and begin the search "on the * inner radius of the annulus |Z|=INR0, where INR0 is the Cauchy * lower bound, in the direction of the Laguerre step from the * origin", that is, we set ZN:= (INR0/ALGSN)*LGSN. In addition, we * set ALGS0:= INR0/N^2. The only purpose in setting ALGS0 is for * the test at entry BRANCH in PROOT Block 6 for a real or * complex iterate. ALGS0 retains this value until an acceptable * starting point is found, in which case it is set to its true * value of |LGS0|. We then branch to PROOT Block 6 at BRANCH * entry to evaluate the polynomial at the new iterate. * * If, on entry, SPIRAL=TRUE then the spiral search has begun and * "subsequent trial points lie on the spirals traced out by * INR0*(1,-1.25/N)^K for K=0,1,2,... where the angle between * consecutive trial points is -ATAN(N/1.25)... If every fourth * trial point is examined, the locus is a spiral progressing in * a counterclockwise direction. Therefore, the configuration of * trial points resembles four outward going equiangular spirals." * That is, we set ZN:= ZN*(-1.25/N,1) and branch to PROOT Block * 6 at BRANCH entry to evaluate the polynomial at the new * iterate. * * Note that this block only exits * (1) By branching to PROOT Block 6 at the BRANCH entry, * (2) to ABORT if an attention key abort condition exists. * * Here: D1-> STARTL * R4[A]= RSMEM * * Test SPIRAL branch if set. * Set SPIRAL. * Set ALGS0:= INR0/N^2 * ZN:= (INR0/ALGSN)*LGSN * Branch to PROOT Block 6 at BRANCH entry. ***************************************************************** SEARCH D1=D1+ 1 D1->SPIRAL C=DAT1 S C[S]= SPIRAL ?C#0 S SPIRAL=TRUE? GOYES sptrue Jump if so ( bc ) C=C+1 S Else, set DAT1=C S SPIRAL=TRUE GOSUB FN (A,B)=N, hex, D1->INR0 SETDEC Mode for MULTF GOSUB x^2 (A,B)=N^2 GOSUB stabcd (C,D)=N^2 GOSUB getab1 (A,B)=INR0, D1->QSCR0 GOSUB stab2 (R2,R3)=INR0 GOSUB divf (A,B)=INR0/N^2 GOSUB SALGS0 ALGS0:= INR0/N^2, hex, D1->G GOSUB D1+63 D1-> Re(LGSN) GOSUB CFTCD1 (A,B),(R0,R1)= LGSN, D1->FN GOSUB D1-84 D1-> Re(ZN) GOSUB CSTOD1 ZN:= LGSN, D1->LGSN GOSUB D1+63 D1-> ALGSN GOSUB getcd1 (C,D)=ALGSN, D1->OUTRN SETDEC Mode for DIVF GOSUB rcab2 (A,B)=INR0 GOSUB divf (A,B)=INR0/ALGSN GOSUB D1-126 D1-> Re(ZN) GOSUB MPYUT+ ZN:= LGSN*(INR0/ALGSN) GOTO BRANCH Exit to Block 6 ( c ) ***************************************************************** * Set ZN:= ZN*(-1.25/N,1) and branch to PROOT Block 6 at BRANCH * entry. ***************************************************************** sptrue GOSUB FREZN GOSUB FTHCM+ (A,B),(R0,R1)= ZN, hex, D1->LGSN GOSUB C->SCR01 (QSCR0,QSCR1)= ZN, dec GOSUB D1+105 D1->N GOSUB getab1 (A,B)=N, D1->INR0 C=0 W P= 12 LCHEX 125 D=C W C=0 M C=-C-1 S (C,D)=-1.25 GOSUB xyex (A,B)=-1.25, (C,D)=N GOSUB divf (A,B)= -1.25/N GOSUB stab0 (R0,R1)= -1.25/N A=0 W B=0 W P= 14 B=B+1 P (A,B)=1 GOSUB exab0 (A,B),(R0,R1)= (-1.25/N,1) GOSUB CpMPY (A,B),(R0,R1)= ZN*(-1.25/N,1) GOSUB SREZN ZN:= ZN*(-1.25/N,1), hex GOSUB STOCM+ D1->LGSN GOSBVL =ATTNCHK GOTO BRANCH Exit to Block 6 ( c ) STITLE PROOT BLOCK 11 ***************************************************************** * * PROOT BLOCK 11 * -------------- * * DECLARATION OF A ZERO * * The function of this block is to declare and store a root at * ZN. In addition, the coefficients of the quotient polynomial * in DFLTMP() are transferred to array A() and the degree of the * polynomial is reduced. * * This block is entered (at its only entry point ZERO) from * (1) PROOT Block 7 when a zero is to be stored at the real * iterate ZN; * (2) PROOT Block 8 when a zero is to be stored at the complex * iterate ZN; * (3) PROOT Block 9 when a zero is to be stored at ZN due to the * fact that, although the iterate is being rejected, we have * FN < E*N^2. * * When entered from the evaluation blocks (7 and 8), ZN is being * declared as a zero either because FN <= E or because LGSN has * become negligible with respect to ZN. * * Exit from this block is * (1) to PROOT Block 3 at the REENT entry point, or * (2) to ABORT if an attention key abort is detected during the * store code or during movement of DFLTMP() to A(). * * This block first tests IM(ZN). IM(ZN)=0 if and only if ZN is * to be a real zero; this is designed into the PROOT algorithm. * * If ZN is real, then we store a root at (Re(ZN),0), decrease * the degree of the polynomial by 1 and move DFLTMP() to A(). * * If ZN is complex, then we store a root at ZN and CONJ(ZN), * decrease the degree of the polynomial by 2, and move DFLTMP() * to A(). * * The logic for moving DFLTMP() to A() is as follows (the * coefficient of highest order does not change): * * (let N = new CD) * FOR I=N-1 TO 0 @ A(I)=DFLTMP(I) @ NEXT I * * Store a real zero at ZN or a complex pair at ZN and CONJ(ZN). ***************************************************************** ZERO GOSUB FREZN (A,B)=Re(ZN), hex GOSUB FTHCM+ (A,B),(R0,R1)= ZN C=R1 C= mantissa(Im(ZN)) ?C=0 W Is ZN real? GOYES zero1 Jump if so ( bc ) C=R0 C[S]= sign(Im(ZN)) SETDEC Mode for compliment C=-C-1 S (A,B),(R0,R1)= R0=C CONJ(ZN) GOSUB STORE Store root at CONJ(ZN) GOSUB FREZN (A,B),(R0,R1) = ZN GOSUB FTHCM+ zero1 GOSUB STORE Store root at ZN, hex *********************************************************** * * Here: C[A]= new CD>0 * R4[A]= RSMEM * CRTA updated * CD>0 updated * hex mode * * Move DFLTMP() to A(). *********************************************************** C=C-1 A C[A]=CD-1 R0=C R0[A]=CD-1 GOSUB AIADDR D1,C[A]-> A(CD-1) D0=C D0-> A(CD-1) C=R0 C[A]=CD-1 GOSUB DFLADR D1-> DFLTMP(CD-1) C=R0 C[A]=CD-1 B=C A B[A]=CD-1 zero2 A=DAT1 A Read Exponent D1=D1+ 5 DAT0=A A Write Exponent D0=D0+ 5 A=DAT1 W Read Mantissa and Sign D1=D1+ 16 DAT0=A W Write Mantissa and Sign D0=D0+ 16 GOSBVL =ATTNCHK hex B=B-1 A Decrement count GONC zero2 Move CD elements ( bc ) GOLONG REENT Exit block via REENT STITLE REGISTER MOVEMENT ROUTINES *********************************************************** stabcd C=B W D=C W C=A W RTN *********************************************************** * ->CD,A: * * entry: R4[A]= RSMEM * * exit: D1= RSMEM * A[A]-> A(N) in RSMEM * P=0, hexmode *********************************************************** ->CD,A SETHEX A=R4 D1=A D1= RSMEM P= 0 LC(5) 626 A=A+C A RTN STITLE JUMP TABLE *********************************************************** xyex GOVLNG =XYEX stab2 GOVLNG =STAB2 rcab0 GOVLNG =RCAB0 rcab2 GOVLNG =RCAB2 rccd0 GOVLNG =RCCD0 rccd2 GOVLNG =RCCD2 *********************************************************** * Add and subtract utilities *********************************************************** raddr01 GOSUB rccd0 raddf CD0EX RSTK=C CD0EX GOSBVL =RADDF C=RSTK D0=C RTN *********************************************************** * Multiply and divide utilities *********************************************************** multr01 GOSUB rccd0 multf GOVLNG =MULTF multd0 GOSUB getcd0 GOTO multf x^2 GOSUB stabcd GOTO multf divf GOVLNG =DIVF mult2 GOVLNG =MULT2 div2 GOVLNG =X/2 nonzero arguments!! STITLE STORE ***************************************************************** ***************************************************************** ** ***Name: STORE ** ***Category: ** ***Abstract: Stores a complex root Z into the result array of ** the PROOT function. ** ***Entry Conditions: (A,B) = Re(Z) (separated) ** (R0,R1) = Im(Z) (separated) ** R4[A]= RSMEM ** CD>0 assumed ** ***Exit Conditions: via ABORT if attention key abort ** via LOOP if CD becomes zero ** C[A]= new CD ** D1=R4[A]= RSMEM ** CRTA updated ** CD>0 updated ** hex mode, carry set ** ***Stack Levels: 1 ***Error Exits: ***Calls: =PACK,=EXAB0, =ATTNCHK ** ***Uses: A,B,C,D,R0,R1,P,D0,D1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** STORE SETDEC Mode for PACK GOSBVL =PACK A = packed Re(Z) GOSUB exab0 (A,B)= Im(Z), R0= Re(Z) GOSBVL =PACK A = packed Im(Z) AR0EX A= Re, R0= Im C=R4 D1=C D1-> RSMEM D1=D1+ 12 D1-> CRTA C=DAT1 A C[A]= CRTA D0=C D0-> destination DAT0=A W D0=D0+ 16 A=R0 DAT0=A W D0=D0+ 16 GOSBVL =ATTNCHK Exit if attention key abort, hex CD0EX C[A]= updated CRTA DAT1=C A Update CRTA D1=D1- 12 D1-> CD C=DAT1 A C[A]= old CD C=C-1 A C[A] = new degree DAT1=C A Update CD ?C#0 A CD=0? RTNYES Return if not ( bc ) P= 0 GOVLNG =GETPTRLOOP STITLE LEAD0 ***************************************************************** ***************************************************************** ** ***Name: LEAD0 ** ***Category: ** ***Abstract: Eliminates leading zeros from the coefficient array ** A(). ** ***Entry Conditions: C[S] set appropriately ** R4[A]= RSMEM ** CD>0 assumed ** ***Exit Conditions: via ABORT if attention key abort ** via STORE if CD becomes zero ** R4[A]= RSMEM ** CRTA updated ** CD>0 updated ** P=0, hex mode ** ***Stack Levels: 2 ***Error Exits: ***Calls: FTCHAI,=MOVEDOWN,=STAB0,STORE ** ***Uses: A,B,C,D,R0,R1,P,D0,D1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** This routine is called from two points: ** ** (1) PROOT BLOCK 1 to eliminate leading zeros in the original ** coefficient array. In this case, every leading zero ** produces a root at (Inf,0) with no exceptions raised. ** C[S]=0 to flag this entry. ** ** (2) PROOT BLOCK 2 to eliminate leading zeros resulting from ** re-scaling the coefficients. In this case, every leading ** zero produces a root at (1E1000,0) which signals ** overflow. C[S]#0 to flag this entry. ** ** Logic: 10 IF A(CD)#0 THEN RETURN ** 20 STORE ZERO AT (Inf,0) OR (1E1000,0) ** 30 SLIDE UP COEFFICIENTS A(CD-1)...A(0) ** 40 CRTA = CRTA+DL ! DL = result data length ** 50 CD = CD-1 ** 60 IF CD#0 THEN 10 ELSE EXIT PROOT ** ***Test Classes: ** ***************************************************************** ***************************************************************** LEAD0 A=R4 A=C S R4=A R4[S]= entry point flag D1=A D1= RSMEM -> CD C=DAT1 A C[A]= CD leadlp R0=C R0[A] = CD GOSUB FTCHAI (A,B)= A(CD), C[A]-> A(CD), hex, P=0 ?B#0 W A(CD)=0? RTNYES Return if not ( bc ) CD1EX D1-> destination, C[A]-> source D0=C D0-> source C=R0 C[A]= CD GOSBVL =21*N D[A]=21CD (# nibbles to slide down) C=D A Into C[A] GOSBVL =MOVEDOWN Slide down them thar nibbles A=0 W Prepare (A,B) to hold result R0=A (R0,R1)= 0 (imaginary part) R1=A C=R4 C[S]= entry point flag ?C=0 S Want (Inf,0)? GOYES infans Jump if so ( Tc ) A=A+1 M Else, P= 14 (A,B)= B=B+1 P 1E1000 GONC gotl (BET) infans SETDEC GOSBVL =DVZXCP if returns, have maxreal, else error SETHEX gotl GOSUB STORE Go store root GOTO leadlp ( c ) STITLE TRAL0 ***************************************************************** ***************************************************************** ** ***Name: TRAL0 ** ***Category: ** ***Abstract: Eliminates trailing zeros from the coefficient ** array A(). ** ***Entry Conditions: C[S] set appropriately ** R4[A]= RSMEM ** CD>0 assumed ** ***Exit Conditions: via ABORT if attention key abort ** via STORE if CD becomes zero ** R4[A]= RSMEM ** CRTA updated ** CD>0 updated ** P=0, hex mode, carry set ** ***Stack Levels: 2 ***Error Exits: ***Calls: FTCHAI,=STAB0,STORE ** ***Uses: A,B,C,D,R0,R1,P,D0,D1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** This routine is called from two points: ** ** (1) PROOT BLOCK 1 to eliminate trailing zeros in the original ** coefficient array. In this case, every trailing zero ** produces a root at (0,0) with no exceptions raised. ** C[S]=0 to flag this entry. ** ** (2) PROOT BLOCK 3 to eliminate trailing zeros at reentry ** for the current (reduced) polynomial. In this case, ** every trailing zero produces a root at (1E-1000,0) ** which signals underflow. C[S]#0 to flag this entry. ** ** Logic: 10 IF A(0)#0 THEN RETURN ** 20 STORE ZERO AT (0,0) OR (1E-1000,0) ** 30 CRTA = CRTA+DL ! DL = result data length ** 40 CD = CD-1 ** 50 IF CD#0 THEN 10 ELSE EXIT PROOT ** ***Test Classes: ** ***************************************************************** ***************************************************************** TRAL0 A=R4 A=C S R4=A R4[S]= entry point flag trallp C=0 A C[A]=0 => fetch A(0) GOSUB FTCHAI (A,B)=A(0), hex ?B#0 W A(0)=0? ( bc ) RTNYES Return if not A=0 W Prepare (A,B) for root R0=A (R0,R1)= 0 (imaginary part) R1=A C=R4 C[S]= entry point flag ?C=0 S Want (0,0)? GOYES gott Jump if SO ( Tc ) P= 14 B= B=B+1 P mantissa of 1E-1000 SETDEC Mode for decrement A=A-1 A A[A]= A=0 X exponent of 1E-1000 gott GOSUB STORE Go store root GOC trallp (BET) ( c ) STITLE prtPAK ***************************************************************** ***************************************************************** ** ***Name: prtPAK ** ***Category: ** ***Abstract: Tests the separated number X in (A,B) for special ** underflow/overflow. ** ***Entry Conditions: (A,B)= value to test ** ***Exit Conditions: (A,B)= possibly underflowed value ** carry set iff overflow ** P=0, decmode ** ***Stack Levels: 0 ***Error Exits: ***Calls: ** ***Uses: C[A] ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** Logic: If |X| >= 1E5000 then exit with carry SET ** If |X| < 1E-5000 then set X=0 and exit with carry CLEAR ** Else, exit with carry CLEAR ** ***Test Classes: ** ***************************************************************** ***************************************************************** prtPAK SETDEC Mode for rest C=A A C[A]= exponent(X) C=C+C A CS iff exponent(X)<0 P= 0 C[A]=5000 LCHEX 05000 GOC negexp Jump if exponent(X)<0 ( bc ) ?A>=C A Exponent(X) >= 5000? RTNYES Return CS if so ( Fc ) RTN Else, return CC ( c ) negexp C=-C A C[A]= -5000 ?A=Y ** P=14, decmode ** ***Stack Levels: 1 ***Error Exits: ***Calls: ** ***Uses: A,B,C,D ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** X Carry CLEAR ** X>Y => Carry SET, B[14...0]<>0 ** X=Y => Carry SET, B[14...0]=0 ** ***Test Classes: ** ***************************************************************** ***************************************************************** prTEST SETDEC Mode for rest C=-C-1 S (C,D) = -Y GOSUB raddf (A,B) = X-Y P= 14 Set pointer for mantissa test ?B=0 WP X=Y? RTNYES Return carry set if so ( Fc ) ?A=0 S X>Y? RTNYES Return carry set if so ( bc ) RTN Else, return carry clear ( c ) STITLE CSTOD0 ***************************************************************** ***************************************************************** ** ***Name: CSTOD0 ** ***Category: ** ***Abstract: Stores a complex variable at the location pointed ** to by D0 ** ***Entry Conditions: D0-> where to push complex value ** (A,B)= Re ** (R0,R1)= Im ** ***Exit Conditions: D0-> past pushed value ** (A,B)= Re ** (R0,R1)= Im ** ***Stack Levels: 1 ***Error Exits: ***Calls: =PUTAB0, =EXAB0 ** ***Uses: D0 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** CSTOD0 GOSUB putab0 Push Re STOCM0+ GOSUB exab0 (A,B)= Im GOSUB putab0 Push Im exab0 GOVLNG =EXAB0 Restore Re and Im ***************************************************************** * C->SCR01: Store complex number in (QSCR0,QSCR1) * * Exit: decmode ***************************************************************** C->SCR01 GOSUB QSTO0 dec GOTO STOCM0+ ***************************************************************** * C<->SCR01: Exchange complex number in (QSCR0,QSCR1) * * Exit: decmode ***************************************************************** C<->SCR01 GOSUB QRCL0 (C,D)= Re(Z), D0-> QSCR1, dec GOSUB ptab0- (QSCR0)= Re(2C), D0-> QSCR1 GOSUB getab0 (A,B)= Im(Z), D0-> QSCR2 GOSUB exab0 (A,B)= Im(2C), (R0,R1)= Im(Z) GOSUB ptab0- (QSCR1)= Im(2C), D0-> QSCR2 stcdab (A,B),(R0,R1)= Z A=C W C=D W B=C W RTN STITLE CSTOD1 ***************************************************************** ***************************************************************** ** ***Name: CSTOD1 ** ***Category: ** ***Abstract: Stores a complex variable at the location pointed ** to by D1 ** ***Entry Conditions: D1-> where to push complex value ** (A,B)= Re ** (R0,R1)= Im ** ***Exit Conditions: D1-> past pushed value ** (A,B)= Re ** (R0,R1)= Im ** ***Stack Levels: 1 ***Error Exits: ***Calls: =PUTAB1, =EXAB0 ** ***Uses: D1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** Alternative entry: ** ** STOCM+: Stores a complex variable given that the real part is ** stored and D1-> Im location. ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** CSTOD1 GOSUB putab1 Push Re STOCM+ GOSUB exab0 (A,B)=Im, (R0,R1)=Re GOSUB putab1 Push Im GOTO exab0 Restore (A,B) and (R0,R1) STITLE CFTCD0 ***************************************************************** ***************************************************************** ** ***Name: CFTCD0 ** ***Category: ** ***Abstract: Fetches a complex variable at the location pointed ** to by D0 ** ***Entry Conditions: D0-> where to fetch complex value ** ***Exit Conditions: D0-> past fetched value ** (A,B)= Re ** (R0,R1)= Im ** ***Stack Levels: 1 ***Error Exits: ***Calls: =GETAB0, =STAB0, =EXAB0 ** ***Uses: D0 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** CFTCD0 GOSUB getab0 Pop Re GOSUB stab0 (R0,R1)= Re GOSUB getab0 Pop Im GOTO exab0 (A,B)= Re, (R0,R1)= Im STITLE CFTCD1 ***************************************************************** ***************************************************************** ** ***Name: CFTCD1 ** ***Category: ** ***Abstract: Fetches a complex variable at the location pointed ** to by D1 ** ***Entry Conditions: D1-> where to fetch complex value ** ***Exit Conditions: D1-> past fetched value ** (A,B)= Re ** (R0,R1)= Im ** ***Stack Levels: 1 ***Error Exits: ***Calls: =GETAB1, =STAB0, =EXAB0 ** ***Uses: D1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** Alternative entry: ** ** FTHCM+: Fetches a complex variable given that the real part is ** fetched and D1-> Im location. ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** CFTCD1 GOSUB getab1 (A,B)=Re FTHCM+ GOSUB stab0 (R0,R1)=Re GOSUB getab1 (A,B)=Im GOTO exab0 (A,B)=Re, (R0,R1)=Im STITLE MPYUTI ***************************************************************** ***************************************************************** ** ***Name: MPYUTI ** ***Category: ** ***Abstract: Multiplies a complex variable Z by a real number X ** and stores X*Z over Z. ** ***Entry Conditions: (A,B)= X ** C[X]= offset for Z in PRVARS area ** ***Exit Conditions: (A,B)= Im(X*Z) ** (R0,R1)=X ** Z replaced by X*Z ** D1-> past variable ** decimal mode ** ***Stack Levels: 3 ***Error Exits: ***Calls: PVARAD, MPYSUB, =GETAB1, =STAB0, =RCCD0, =MULTF ** exits thru =PUTAB1 ** ***Uses: A,B,C,D,R0,R1,P ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** Alternate entries: ** ** MPYUT+: Same as MPYUTI but assumes variable address in ** D1 on entry; no C assumptions ** ** MPYU2+: Same as MPYUT+ but sets X in (A,B) to 2 ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** MPYU2+ GOSUB X=2 (A,B)=2 GONC MPYUT+ Do MPYUT+ (BET) MPYUTI GOSUB PVARAD D1-> variable MPYUT+ SETDEC Mode for rest GOSUB stab0 (R0,R1)=X GOSUB MPYSUB Write X*Re(Z) over Re(Z) MPYSUB GOSUB getab1 (A,B)=Im(Z) GOSUB multr01 (A,B)=X*Im(Z) GOTO ptab1- Write X*Im(Z) over Im(Z) X=2 A=0 W B=0 W P= 14 (A,B)=2 B=B+1 P B=B+1 P RTN ***************************************************************** * * D1 adjustments. These routines move D1 and are used when it is * necessary to locate a PROOT variable but not use time going * through the variable fetch routines. Any adjustment of D1 by * less than 42 must be done in-line. * ***************************************************************** D1-126 D1=D1- 16 D1=D1- 5 D1-105 D1=D1- 16 D1=D1- 5 D1-84 D1=D1- 16 D1=D1- 5 D1-63 D1=D1- 16 D1=D1- 5 D1-42 D1=D1- 16 D1=D1- 16 D1=D1- 10 RTN D1+126 D1=D1+ 16 D1=D1+ 5 D1+105 D1=D1+ 16 D1=D1+ 5 D1+84 D1=D1+ 16 D1=D1+ 5 D1+63 D1=D1+ 16 D1=D1+ 5 D1+42 D1=D1+ 16 D1=D1+ 16 D1=D1+ 10 RTN ***************************************************************** * * D0 adjustments. These routines move D0 and are used when it is * necessary to locate a PROOT variable but not use time going * through the variable fetch routines. Any adjustment of D0 by * less than 42 must be done in-line. * ***************************************************************** D0-126 D0=D0- 16 D0=D0- 5 D0-105 D0=D0- 16 D0=D0- 5 D0-84 D0=D0- 16 D0=D0- 5 D0-63 D0=D0- 16 D0=D0- 5 D0-42 D0=D0- 16 D0=D0- 16 D0=D0- 10 RTN D0+84 D0=D0+ 16 D0=D0+ 5 D0+63 D0=D0+ 16 D0=D0+ 5 D0+42 D0=D0+ 16 D0=D0+ 16 D0=D0+ 10 RTN ***************************************************************** ptab1- D1=D1- 16 D1=D1- 5 putab1 GOVLNG =PUTAB1 ***************************************************************** ptab0- D0=D0- 16 D0=D0- 5 putab0 GOVLNG =PUTAB0 ***************************************************************** gtab1- D1=D1- 16 D1=D1- 5 getab1 GOVLNG =GETAB1 ***************************************************************** getab0 GOVLNG =GETAB0 ***************************************************************** getcd1 GOVLNG =GETCD1 STITLE QSCRATCH Store ***************************************************************** ***************************************************************** ** ***Name: QSTOX ** ***Category: ** ***Abstract: Store into QSCRATCH. ** ***Entry Conditions: (A,B)= value to store ** ***Exit Conditions: D0-> past value stored ** C[A]-> at value stored ** P=0, dec mode ** ***Stack Levels: 0 ***Error Exits: ***Calls: ** ***Uses: C,D[A],D0 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** QSTO5 P= 0 LC(5) 605 Offset to QSCR5 GOTO qstox QSTO0 P= 0 LC(5) 500 Offset to QSCR0 qstox D=C A D[A]= offset C=R4 SETHEX C=C+D A C[A]= address SETDEC D0=C Into D0 GOTO putab0 Push value STITLE QSCRATCH Fetch ***************************************************************** ***************************************************************** ** ***Name: QRCLX ** ***Category: ** ***Abstract: Fetch from QSCRATCH. ** ***Entry Conditions: ** ***Exit Conditions: (C,D)= value fetched ** D0-> past value fetched ** P=0, dec mode ** ***Stack Levels: 0 ***Error Exits: ***Calls: ** ***Uses: C,D[A],D0 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** QRCL1 P= 0 LC(5) 521 Offset to QSCR1 GOTO qrclx QRCL0 P= 0 LC(5) 500 Offset to QSCR0 qrclx D=C A D[A]= offset C=R4 SETHEX C=C+D A C[A]= address SETDEC D0=C Into D0 getcd0 GOVLNG =GETCD0 [C,D] <--- QSCRx STITLE Variable Store ***************************************************************** ***************************************************************** ** ***Name: SXXXXX ** ***Category: ** ***Abstract: Store PROOT variable. The PROOT variables are ** stored in the PVARS area and are accessed by an offset from ** the start of the PRVARS area. ** ***Entry Conditions: (A,B)= value to store ** ***Exit Conditions: D1-> past value stored ** C[A]-> at value stored ** P=0, hex mode ** ***Stack Levels: 1 ***Error Exits: ***Calls: ** ***Uses: C,D[A],D1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** The offsets are as follows: ** ** E 17 ** PSCR1 38 ** PSCR2 59 ** PSCR3 80 ** PSCR4 101 ** PSCR5 122 ** PSCR6 143 ** RE(Z0) 164 ** IM(Z0) 185 ** RE(LGS0) 206 ** IM(LGS0) 227 ** F0 248 ** ALGS0 269 ** OUTR0 290 ** RE(ZN) 311 ** IM(ZN) 332 ** RE(LGSN) 353 ** IM(LGSN) 374 ** FN 395 ** ALGSN 416 ** OUTRN 437 ** N 458 ** INR0 479 ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** SE P= 0 LC(3) 17 Offset to store E GOTO PVARST SSCR2 P= 0 LC(3) 59 Offset to store PSCR2 GOTO PVARST SSCR3 P= 0 LC(3) 80 Offset to store PSCR3 GOTO PVARST SSCR5 P= 0 LC(3) 122 Offset to store PSCR5 GOTO PVARST SALGS0 P= 0 LC(3) 269 Offset to store ALGS0 GOTO PVARST SOUTR0 P= 0 LC(3) 290 Offset to store OUTR0 GOTO PVARST SREZN P= 0 LC(3) 311 Offset to store RE(ZN) GOTO PVARST SRELGN P= 0 LC(3) 353 Offset to store RE(LGSN) GOTO PVARST SFN P= 0 LC(3) 395 Offset to store FN GOTO PVARST SALGSN P= 0 LC(3) 416 Offset to store ALGSN GOTO PVARST SOUTRN P= 0 LC(3) 437 Offset to store OUTRN GOTO PVARST SN P= 0 LC(3) 458 Offset to store N GOTO PVARST SINR0 P= 0 LC(3) 479 Offset to store INR0 GOTO PVARST SQSCR0 P= 0 LC(3) 500 Offset to store QSCR0 GOTO PVARST STITLE Variable Fetch ***************************************************************** ***************************************************************** ** ***Name: FXXXXX ** ***Category: ** ***Abstract: Fetch PROOT variable. The PROOT variables are ** stored in the PVARS area and are accessed by an offset from ** the start of the PRVARS area. ** ***Entry Conditions: ** ***Exit Conditions: (A,B)= fetched value ** D1-> past value fetched ** C[A]-> at value fetched ** P=0, hex mode ** ***Stack Levels: 1 ***Error Exits: ***Calls: ** ***Uses: A,B,C,D[A],D1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** The offsets are as shown in variable store documentation. ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** FE P= 0 LC(3) 17 Offset to fetch E GOTO PVARFT FSCR1 P= 0 LC(3) 38 Offset to fetch PSCR1 GOTO PVARFT FSCR3 P= 0 LC(3) 80 Offset to fetch PSCR3 GOTO PVARFT FREZ0 P= 0 LC(3) 164 Offset to fetch RE(Z0) GOTO PVARFT FF0 P= 0 LC(3) 248 Offset to fetch F0 GOTO PVARFT FALGS0 P= 0 LC(3) 269 Offset to fetch ALGS0 GOTO PVARFT FOUTR0 P= 0 LC(3) 290 Offset to fetch OUTR0 GOTO PVARFT FREZN P= 0 LC(3) 311 Offset to fetch RE(ZN) GOTO PVARFT FIMZN P= 0 LC(3) 332 Offset to fetch IM(ZN) GOTO PVARFT FRELGN P= 0 LC(3) 353 Offset to fetch RE(LGSN) GOTO PVARFT FFN P= 0 LC(3) 395 Offset to fetch FN GOTO PVARFT FOUTRN P= 0 LC(3) 437 Offset to fetch OUTRN GOTO PVARFT FN P= 0 LC(3) 458 Offset to fetch N GOTO PVARFT FINR0 P= 0 LC(3) 479 Offset to fetch INR0 GOTO PVARFT STITLE PVARST ***************************************************************** ***************************************************************** ** ***Name: PVARST ** ***Category: ** ***Abstract: Stores a PROOT variable given its offset relative ** to the base address of the PRVARS area. ** ***Entry Conditions: (A,B)= value to store ** C[X] = offset ** R4[A]= RSMEM ** ***Exit Conditions: (A,B)= value stored ** C[A]-> at variable ** D1-> past variable ** R4[A]= RSMEM ** P=0, hex mode ** ***Stack Levels: 1 ***Error Exits: ***Calls: PVARAD ** exits thru =PUTAB1 ** ***Uses: A,B,C,D[A],D1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** PVARST GOSUB PVARAD D1-> variable location GOTO putab1 Push value STITLE PVARFT ***************************************************************** ***************************************************************** ** ***Name: PVARFT ** ***Category: ** ***Abstract: Fetches a PROOT variable given its offset relative ** to the base address of the PRVARS area. ** ***Entry Conditions: C[X] = offset ** R4[A]= RSMEM ** ***Exit Conditions: (A,B)= value fetched ** C[A]-> at variable ** D1-> past variable ** R4[A]= RSMEM ** P=0, hex mode ** ***Stack Levels: 1 ***Error Exits: ***Calls: PVARAD ** exits thru =GETAB1 ** ***Uses: A,B,C,D[A],D1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** PVARFT GOSUB PVARAD D1-> variable location GOTO getab1 Pop value STITLE PVARAD ***************************************************************** ***************************************************************** ** ***Name: PVARAD ** ***Category: ** ***Abstract: Computes the address of a PROOT variable given ** its offset relative to the base address of the PRVARS area. ** ***Entry Conditions: C[X] = offset ** R4[A]= RSMEM ** ***Exit Conditions: D1=C[A]-> variable ** R4[A]= RSMEM ** P=0, hex mode ** ***Stack Levels: 0 ***Error Exits: ***Calls: ** ***Uses: C,D[A],D1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** Formula: Address = [RSMEM] + C[X] ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** ***************************************************************** * * PSCRAD: Sets C[X]=38 to fetch the address of the PROOT scratch * area. * ***************************************************************** PSCRAD P= 0 Offset for PROOT scratch area LC(3) 38 PVARAD SETHEX Mode for rest P= 3 C[A]= offset in 5 nibbles LC(2) 0 P= 0 D=C A D[A]= offset C=R4 C=C+D A C[A]= address D1=C Into D1 RTN STITLE FTCHAI ***************************************************************** ***************************************************************** ** ***Name: FTCHAI ** ***Category: ** ***Abstract: Fetches the coefficient A(I) of the current array ** A() where 0<= I <= CD. ** ***Entry Conditions: C[A] = I ** R4[A]= RSMEM ** ***Exit Conditions: (A,B)= A(I) ** C[A]-> at fetch location ** D1-> past fetch location ** P=0, hex mode ** ***Stack Levels: 2 ***Error Exits: ***Calls: AIADDR ** exits thru =GETAB1 ** ***Uses: A,B,C[A],D[A],D1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** FTCHAI GOSUB AIADDR D1-> A(I) location GOTO getab1 Pop value STITLE AIADDR ***************************************************************** ***************************************************************** ** ***Name: AIADDR ** ***Category: ** ***Abstract: Computes the address of the coefficient A(I) of ** the current array A() where 0 <= I <= CD. ** ***Entry Conditions: C[A] = I ** R4[A]= RSMEM ** ***Exit Conditions: D1=C[A]-> A(I) ** R4[A]= RSMEM ** P=0, hex mode ** ***Stack Levels: 1 ***Error Exits: ***Calls: =21*N ** ***Uses: C[A],D[A],D1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** Formula: Address = [RSMEM] + 626 + 21(CD-I) ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** AIADDR SETHEX Mode for rest D=C A D[A]= I C=R4 D1=C D1= RSMEM C=DAT1 A C[A]= CD C=C-D A C[A]= CD-I GOSBVL =21*N D[A]= 21*(CD-I) P= 0 C[A]= LC(5) 626 626 AIADR+ D=C+D A D[A]= 21*(CD-I) + 626 CD1EX C[A]= RSMEM C=C+D A C[A]= A(I) address D1=C Into D1 RTN STITLE DFLST ***************************************************************** ***************************************************************** ** ***Name: DFLST ** ***Category: ** ***Abstract: Stores the coefficient DFLTMP(I) of the current ** array DFLTMP() where 0 <= I <= CD. ** ***Entry Conditions: (A,B)= value to store ** C[A] = I ** R4[A]= RSMEM ** ***Exit Conditions: (A,B)= A(I) ** C[A]-> at store location ** D1-> past store location ** P=0, hex mode ** ***Stack Levels: 2 ***Error Exits: ***Calls: DFLADR ** exits thru =PUTAB1 ** ***Uses: C[A],D[A],D1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** DFLST GOSUB DFLADR D1-> DFLTMP(I) location GOTO putab1 Push value STITLE DFLFT ***************************************************************** ***************************************************************** ** ***Name: DFLFT ** ***Category: ** ***Abstract: Fetches the coefficient DFLTMP(I) from the current ** array DFLTMP() where 0<= I <= CD. ** ***Entry Conditions: C[A] = I ** R4[A]-> RSMEM ** ***Exit Conditions: (A,B)= DFLTMP(I) ** C[A]-> at fetch location ** D1-> past fetch location ** P=0, hex mode ** ***Stack Levels: 2 ***Error Exits: ***Calls: DFLADR ** exits thru =GETAB1 ** ***Uses: A,B,C[A],D[A],D1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** DFLFT GOSUB DFLADR D1-> DFLTMP(I) location GOTO getab1 Pop value STITLE DFLADR ***************************************************************** ***************************************************************** ** ***Name: DFLADR ** ***Category: ** ***Abstract: Computes the address of the coefficient DFLTMP(I) ** of the current array DFLTMP() where 0 <= I <= CD. ** ***Entry Conditions: C[A] = I ** R4[A]= RSMEM ** ***Exit Conditions: D1=C[A]-> DFLTMP(I) ** R4[A]= RSMEM ** P=0, hex mode ** ***Stack Levels: 1 ***Error Exits: ***Calls: =21*N ** ***Uses: C[A],D[A],D1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** Formula: Address = [RSMEM] + 21(OD+CD-I) + 647 ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** DFLADR SETHEX Mode for rest D=C A D[A]= I C=R4 D1=C D1= RSMEM C=DAT1 A C[A]= CD D=C-D A D[A]= CD-I D1=D1+ 5 D1-> OD C=DAT1 A C[A]= OD C=C+D A C[A]= OD+CD-I GOSBVL =21*N D[A]= 21(OD+CD-I) P= 0 C[A]= 642 (recall LC(5) 642 D1 is 5 nibbles up) GOTO AIADR+ C[A], D1 = address STITLE CpADD ***************************************************************** ***************************************************************** ** ***Name: CpADD ** ***Category: ** ***Abstract: Adds two complex numbers Z and W. ** ***Entry Conditions: (A,B)= Re(W) ** (R0,R1)= Im(W) ** QSCR0= Re(Z) ** QSCR1= Im(Z) ** ***Exit Conditions: (A,B)= Re(Z+W) ** (R0,R1)= Im(Z+W) ** D0-> QSCR2 ** QSCR0,1 preserved ** decmode ** ***Stack Levels: 2 ***Error Exits: ***Calls: QRCL0, raddf, =EXAB0 ** ***Uses: A,B,C,D,P,R0,R1,D0 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** CpADD GOSUB QRCL0 (C,D)= Re(Z), D0->QSCR1, dec GOSUB raddf (A,B)=Re(Z+W) GOSUB exab0 (A,B)=Im(W),(R0,R1)=Re(Z+W) GOSUB getcd0 (C,D)= Im(Z), D0-> QSCR2 GOSUB raddf (A,B)=Im(Z+W) GOTO exab0 (A,B),(R0,R1)= Z+W STITLE CpMPY ***************************************************************** ***************************************************************** ** ***Name: CpMPY ** ***Category: ** ***Abstract: Multiplies two complex numbers Z and W. ** ***Entry Conditions: (A,B)= Re(W) ** (R0,R1)= Im(W) ** QSCR0= Re(Z) ** QSCR1= Im(Z) ** ***Exit Conditions: (A,B)= Re(Z*W) ** (R0,R1)= Im(Z*W) ** D0-> past QSCR5 ** QSCR0,1 preserved ** decmode ** ***Stack Levels: 1 ***Error Exits: ***Calls: ** ***Uses: A,B,C,D,P,R0,R1,R2,R3,D0, QSCR5 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** CpMPY GOSUB QSTO5 (QSCR5)= Re(W), dec CpMPYf GOSUB QRCL1 (C,D)=Im(Z), D0->QSCR2, dec GOSUB multf (A,B)=Re(W)*Im(Z) GOSUB stab2 (R2,R3)=Re(W)*Im(Z) GOSUB rcab0 (A,B)=Im(W) GOSUB D0-42 D0-> PSCR0 GOSUB multd0 (A,B)=Re(Z)*Im(W),D0->QSCR1, dec GOSUB rccd2 (C,D)=Re(W)*Im(Z) GOSUB raddf (A,B)=Im(Z*W) GOSUB exab0 (A,B)=Im(W),(R0,R1)=Im(Z*W) A=-A-1 S GOSUB multd0 (A,B)=-Im(Z)*Im(W) GOSUB stab2 (R2,R3)= -Im(Z)*Im(W) GOSUB D0-42 D0-> PSCR0 GOSUB getab0 (A,B)=Re(Z), D0->QSCR1, dec GOSUB D0+84 D0-> QSCR5 GOSUB multd0 (A,B)=Re(W)*Re(Z) GOSUB rccd2 (C,D)= -Im(Z)*Im(W) GOTO raddf (A,B)=Re(Z*W) STITLE CpDIV ***************************************************************** ***************************************************************** ** ***Name: CpDIV ** ***Category: ** ***Abstract: Divides a complex number Z by a complex number W. ** ***Entry Conditions: (A,B)= Re(W) ** (R0,R1)= Im(W) ** QSCR0= Re(Z) ** QSCR1= Im(Z) ** ***Exit Conditions: (A,B)= Re(Z/W) ** (R0,R1) =Im(Z/W) ** QSCR0: Re(Z) ** QSCR1: Im(Z) ** QSCR4: |W|^2 if >0 ** decmode ** ***Stack Levels: 3 ***Error Exits: ***Calls: ** ***Uses: A,B,C,D,P,R0,R1,R2,R3,D0, QSCR4,QSCR5 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** Note: If W=0 then (1E5000,1E5000) is returned. It is ** assumed that Z and W are not both zero. ** ***Test Classes: ** ***************************************************************** ***************************************************************** CpDIV GOSUB QSTO5 (QSCR5)= Re(W), D0->past QSCR5, dec GOSUB exab0 A=-A-1 S GOSUB stab2 (R2,R3)= -Im(W) GOSUB CpMAG^2 (A,B)= |W|^2, trashes R0,R1 P= 14 Set pointer for mantissa test ?B#0 WP W=0? GOYES W#0 Jump if not ( Tc ) B=0 S B=B+1 P C=0 W P= 3 Answer is (1E5000,1E5000) LCHEX 5 A=C W stab0 GOVLNG =STAB0 (R0,R1)= A ( nc ) ***************************************************************** * * Complex division; non-zero denominator. * * (A,B)= |W|^2 * (R2,R3)= -Im(W) * (QSCR0,QSCR1)= Z * (QSCR5)= Re(W) * D0->past QSCR5 * decmode ***************************************************************** W#0 GOSUB D0-42 D0-> QSCR4 GOSUB putab0 QSCR4= |W|^2, D0-> Re(W) GOSUB rcab2 GOSUB stab0 GOSUB getab0 (A,B)=,(R0,R1)=CONJ(W) GOSUB CpMPYf Get Z*CONJ(W), D0-> past QSCR5 GOSUB D0-42 D0-> QSCR4 GOSUB getcd0 (C,D)= |W|^2, D0-> QSCR5 GOSBVL =STCD2 (R2,R3)=|W|^2 GOSUB divf (A,B)=Re(Z/W) GOSUB exab0 Into (R0,R1); (A,B)=Im(Z*CONJ(W)) GOSUB rccd2 (C,D)=|W|^2 GOSUB divf (A,B)=Im(Z/W) GOTO exab0 Store answer ( c ) STITLE A2A1A0 ***************************************************************** ***************************************************************** ** ***Name: A2A1A0 ** ***Category: ** ***Abstract: Places the coefficients A(2), A(1), A(0) into the ** PROOT scratch area as complex quantities. ** ***Entry Conditions: R4[A]= RSMEM ** ***Exit Conditions: R4[A]= RSMEM ** P=0, hex mode ** ** SCRATCH AREA: PSCR1 - A(2) ** PSCR2 - 0 ** PSCR3 - A(1) ** PSCR4 - 0 ** PSCR5 - A(0) ** PSCR6 - 0 ** ***Stack Levels: 2 ***Error Exits: ***Calls: PSCRAD, AIADDR, =GETAB1, CSTOD0 ** ***Uses: A,B,C,D[A],R0,R1,D1,D0, PSCR1-6 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** A2A1A0 GOSUB PSCRAD D1-> PROOT SCRATCH, hex, P=0 D0=C D0-> PROOT SCRATCH C=0 A C[A]=2 to fetch LCHEX 2 A(2) address GOSUB AIADDR D1-> A(2) C=0 W (R0,R1)=0 R0=C (all imaginary R1=C parts are zero) LCHEX 2 C[A]=2 (initial counter) A2A1A+ GOSUB getab1 (A,B)=A(I) GOSUB CSTOD0 Push (A(I),0) at D0 C=C-1 A Decrement count GONC A2A1A+ Loop three times ( bc ) RTN STITLE ITRUTI ***************************************************************** * * Here: (A,B)= FN * D1-> ALGSN * R4[A]= RSMEM * ST[sCITER]=0 iff RITER call * PSCR1 - Re(P''(ZN)/2) * PSCR2 - Im(P''(ZN)/2) * PSCR3 - Re(P'(ZN)) * PSCR4 - Im(P'(ZN)) * PSCR5 - Re(P(ZN)) * PSCR6 - Im(P(ZN)) * * This code finishes processing for the polynomial evaluation * blocks (7 and 8). * * (1) Compute E, the bound for the rounding error in the * polynomial evaluations, by executing the following * pseudocode: (N=CD) * * 10 E = |DFLTMP(N)| * 20 FOR K = N-1 TO 0 * 30 E = |DFLTMP(K)| + |ZN|*E * 40 NEXT K * 50 E = 2E-13*E ! CITER CALL * 60 E = 1E-13*E ! RITER CALL * * Note: The factors 1E-14 and 2E-14, which might seem * better choices, do not generally improve the accuracy of * the polynomials reconstructed from the computed roots. * * (2) If FN <= E then branch to PROOT Block 11 at the ZERO entry * to accept ZN as a root. * * (3) If FN >= F0 and STARTL=TRUE, then branch to PROOT Block 9 * at the REJECT entry to reject the iterate. * * (4) Else, compute the Laguerre step (LGSN), its magnitude * (ALGSN), and the Fejer bound (OUTRN) via the LAGUER call. * * (5) If |ZN|+ALGSN = |ZN|, then branch to PROOT Block 11 at * the ZERO entry to accept ZN as a root. Else, branch to * AMRSTP to accept, modify, or reject the current Laguerre * step. * ***************************************************************** ITRUTI GOSUB D1-105 D1-> Re(ZN) GOSUB CFTCD1 (A,B),(R0,R1) = ZN, D1->LGSN SETDEC Mode for CpABS GOSUB CpABS (A,B)=|ZN|, trashes R0,R1 ***************************************************************** * Compute E * * E = |DFLTMP(N)| * FOR K = N-1 TO 0 * E = |DFLTMP(K)| + |ZN|*E * NEXT K * E = 2E-13*E ! CITER CALL * E = 1E-13*E ! RITER CALL ***************************************************************** GOSUB stab0 (R0,R1)= |ZN| C=R4 D1=C D1-> CD (N) C=DAT1 A C[A]=N SETHEX Mode for decrement C=C-1 A R2=C R2[A]=N-1 (initial K) C=C+1 A GOSUB DFLFT (A,B)= |DFLTMP(N)|, D1->DFLTMP(K) A=0 S (initial E) eloop SETDEC Mode for MULTF GOSUB multr01 (A,B)= |ZN|*E GOSUB getcd1 (C,D)= |DFLTMP(K)|, D1->DFLTMP(K-1) C=0 S GOSUB raddf (A,B)= new E GOSBVL =ATTNCHK hex C=R2 C[A]=K C=C-1 A C[A]= NEXT K R2=C Into R2[A] GONC eloop Loop through K=0 ( bc ) SETDEC Mode for ADDF ?ST=0 sCITER RITER call? GOYES ritcal Jump if so ( bc ) GOSUB mult2 (A,B)= 2*E ritcal P= 0 C[A]= 13 LCHEX 00013 A=A-C A (A,B) = E*1E-13 GOSUB SE Store into E variable, hex ***************************************************************** * * Here: (A,B)= E * (R0,R1)= |ZN| * R4[A]= RSMEM * hexmode * * (2) If FN <= E then branch to PROOT Block 11 at the ZERO entry * to accept ZN as a root. ***************************************************************** GOSUB stab2 (R2,R3)=E GOSUB FFN (A,B)=FN, D1->ALGSN GOSUB stabcd (C,D)=FN GOSBVL =EXAB2 (A,B)=E, (R2,R3)=FN GOSUB prTEST Compare FN and E, dec GONC itrut1 Jump if E < FN ( bc ) GOTO ZERO Else, branch to block 11 ( c ) ***************************************************************** * * Here: (A,B)= E * (R0,R1)= |ZN| * (R2,R3)= FN * R4[A]= RSMEM * D1-> ALGSN * decmode * * (3) If FN >= F0 and STARTL=TRUE, then branch to PROOT Block 9 * at the REJECT entry to reject the iterate. ***************************************************************** itrut1 C=R4 D1=C D1-> RSMEM D1=D1+ 10 D1-> STARTL C=DAT1 S C[S]= STARTL ?C=0 S STARTL=FALSE? GOYES itrut2 Jump if so ( bc ) GOSUB FF0 (A,B)=F0, hex GOSUB rccd2 (C,D)=FN GOSUB xyex (A,B)=FN, (C,D)=F0 GOSUB prTEST Compare FN and F0, dec GONC itrut2 Jump if FNN ***************************************************************** * * Here: R4[A]= RSMEM * D1-> N * decmode * * (5) If |ZN|+ALGSN = |ZN|, then branch to PROOT Block 11 at * the ZERO entry to accept ZN as a root. Else, branch to * AMRSTP to accept, modify, or reject the current Laguerre * step. ***************************************************************** GOSUB D1-42 D1-> ALGSN GOSUB getab1 (A,B)=ALGSN GOSUB stab2 (R2,R3)=ALGSN GOSUB D1-126 D1-> Re(ZN) GOSUB CFTCD1 (A,B),(R0,R1)=ZN GOSUB CpABS (A,B)=|ZN|, trashes R0,R1 GOSUB stabcd (C,D)=|ZN| GOSBVL =EXAB2 (A,B)=ALGSN, (R2,R3)=|ZN| GOSUB raddf (A,B)=|ZN|+ALGSN GOSUB rccd2 (C,D)=|ZN| GOSUB prTEST Compare |ZN|+ALGSN and |ZN| ?B#0 WP Are they equal? GOYES itrut3 Jump if not ( Tc ) GOTO ZERO Branch to block 11 ( nc ) itrut3 GOLONG AMRSTP Branch to block 5 ( c ) STITLE LAGUER ***************************************************************** ***************************************************************** ** ***Name: LAGUER ** ***Category: ** ***Abstract: Computes and assigns the Laguerre step (LGSN), ** its magnitude (ALGSN) and the Fejer bound (OUTRN) at the ** iterate ZN. ** ***Entry Conditions: R4[A]= RSMEM ** ** PSCRATCH AREA: PSCR1 - Re(P''(ZN)/2) ** PSCR2 - Im(P''(ZN)/2) ** PSCR3 - Re(P'(ZN)) ** PSCR4 - Im(P'(ZN)) ** PSCR5 - Re(P(ZN)) ** PSCR6 - Im(P(ZN)) ** ***Exit Conditions: All values assigned ** R4[A]= RSMEM ** D1->N ** D0 trashed ** dec mode ** ***Stack Levels: ***Error Exits: ***Calls: ** ***Uses: ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** Executes the following pseudocode: ** ** LGSN = ((N-2)/N)*(P'(ZN)/P(ZN)) ** B = (N-1)*P'(ZN) ** C = (N*(N-1)/2)*P(ZN) ** B = ** (P''(ZN)/2)*Z^2 + B*Z +C ** OUTRN = |B| ** LGSN = B/(LGSN*B+(N-1)) ** ALGSN = |LGSN| ** OUTRN = MIN{OUTRN,SQRT(N)*ALGSN} ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** ***************************************************************** * * Here: R4[A]= RSMEM * PSCR1 - Re(P''(ZN)/2) * PSCR2 - Im(P''(ZN)/2) * PSCR3 - Re(P'(ZN)) * PSCR4 - Im(P'(ZN)) * PSCR5 - Re(P(ZN)) * PSCR6 - Im(P(ZN)) * * The code below sets LGSN:= ((N-2)/N)*(P'(ZN)/P(ZN)) * ***************************************************************** LAGUER GOSUB FSCR3 (A,B)= Re(P'(ZN)), D1-> Im GOSUB FTHCM+ (A,B),(R0,R1)= P'(ZN), D1-> P(ZN) GOSUB C->SCR01 (QSCR0,QSCR1)= P'(ZN), dec GOSUB CFTCD1 (A,B),(R0,R1)= P(ZN) GOSUB CpDIV (A,B),(R0,R1)= P'(ZN)/P(ZN) GOSUB SRELGN LGSN:= P'(ZN)/P(ZN), hex GOSUB STOCM+ D1-> FN GOSUB D1+63 D1-> N GOSUB getcd1 D1-> INR0, (C,D)= N GOSBVL =STCD2 (R2,R3)= N SETDEC Mode for DIVF GOSUB X=2 (A,B)= 2 GOSUB divf (A,B)= 2/N A=-A-1 S (A,B)= -2/N GOSBVL =RADD1 (A,B)= 1-2/N (=(N-2)/N), trashes D0 GOSUB D1-126 D1-> Re(LGSN) GOSUB MPYUT+ LGSN:= (N-2)*P'(ZN)/N*P(ZN), D1-> FN ***************************************************************** * * Here: D1-> FN * (R2,R3)=N * R4[A]= RSMEM * PSCR1 - Re(P''(ZN)/2) * PSCR2 - Im(P''(ZN)/2) * PSCR3 - Re(P'(ZN)) * PSCR4 - Im(P'(ZN)) * PSCR5 - Re(P(ZN)) * PSCR6 - Im(P(ZN)) * decmode * * Peplace P(ZN) by (N(N-1)/2)*P(ZN) and P'(ZN) by (N-1)*P'(ZN) * in scratch area. * ***************************************************************** GOSUB rcab2 (A,B)= N GOSBVL =RSUB1 (A,B)= N-1, trashes D0 P= 0 C[X]= 80 LC(3) 80 (offset for PSCR3) (OK) GOSUB MPYUTI Replace P'(ZN) by (N-1)P'(ZN) D1-> P(ZN), (R0,R1)= N-1 GOSUB rcab2 (A,B)= N GOSUB multr01 (A,B)= N(N-1) GOSUB div2 (A,B)= N(N-1)/2 GOSUB MPYUT+ Replace P(ZN) by N(N-1)P(ZN)/2 ***************************************************************** * * Here: R4[A]= RSMEM * PSCR1 - Re(P''(ZN)/2) * PSCR2 - Im(P''(ZN)/2) * PSCR3 - Re((N-1)P'(ZN)) * PSCR4 - Im((N-1)P'(ZN)) * PSCR5 - Re(N(N-1)P(ZN)/2) * PSCR6 - Im(N(N-1)P(ZN)/2) * decmode * * Solve the complex quadratic (P''(ZN)/2)*Z^2 + B*Z + C * where B=(N-1)P'(ZN) and C=(N(N-1)/2)P(ZN). * * Set LGSN= B/(LGSN*B+(N-1)), where B= smaller magnitude zero of * above quadratic. * ***************************************************************** GOSUB QDRTC Solve the quadratic, (A,B),(R0,R1)= (PSCR3,4) = B GOSUB C->SCR01 (QSCR0,QSCR1)= B GOSUB FRELGN (A,B)= Re(LGSN), hex GOSUB FTHCM+ (A,B),(R0,R1)= LGSN GOSUB CpMPY (A,B),(R0,R1)= LGSN*B, dec GOSUB stab2 (R2,R3)= Re(LGSN*B) GOSUB FN (A,B)=N, hex SETDEC Mode for subtract GOSBVL =RSUB1 (A,B)= N-1, trashes D0 GOSUB rccd2 (C,D)= Re(LGSN*B) GOSUB raddf (A,B)= Re(LGSN*B+(N-1)) GOSUB CpDIV (A,B),(R0,R1)= B/(LGSN*B+(N-1)) GOSUB D1-126 D1-> Re(LGSN) GOSUB CSTOD1 LGSN:= B/(LGSN*B+(N-1)), D1-> FN ***************************************************************** * * Here: (A,B),(R0,R1)= LGSN * D1-> FN * R4[A]= RSMEM * PSCR3 - Re(B) * PSCR4 - Im(B) * decmode * * Set ALGSN= |LGSN| and OUTRN= min{|B|,ALGSN} * ***************************************************************** GOSUB CpABS (A,B)= |LGSN|, trashes R0,R1 D1=D1+ 16 D1-> D1=D1+ 5 ALGSN GOSUB putab1 ALGSN= |LGSN| GOSUB stab2 (R2,R3)= ALGSN GOSUB FSCR3 (A,B)= Re(B), hex GOSUB FTHCM+ (A,B),(R0,R1)= B SETDEC Mode for CpABS GOSUB CpABS (A,B)= |B|, trashes R0,R1 GOSUB SOUTRN OUTRN= |B|, hex, D1->N GOSUB getab1 (A,B)= N, D1->INR0 SETDEC Mode for SQRF GOSUB sqrf (A,B)= SQRT(N) GOSUB rccd2 (C,D)= ALGSN GOSUB multf (A,B)= SQRT(N)*ALGSN STITLE GOUTRN ***************************************************************** ***************************************************************** ** ***Name: GOUTRN ** ***Category: ** ***Abstract: Sets OUTRN = min{OUTRN,X} ** ***Entry Conditions: (A,B)=X ** ***Exit Conditions: (R0,R1)= X ** R4[A]= RSMEM ** D1->N ** dec mode ** ***Stack Levels: ***Error Exits: ***Calls: ** ***Uses: FOUTRN, prTEST, =STAB0, =RCCD0, =RCAB0 ** may exit thru =PUTAB1 ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** GOUTRN GOSUB stab0 (R0,R1)= X GOSUB FOUTRN (A,B)= OUTRN, hex, D1->N GOSUB rccd0 (C,D)= X GOSUB prTEST Compare X and OUTRN, dec RTNNC Return if OUTRN < X ( bc ) GOSUB rcab0 Else, set GOTO ptab1- OUTRN:= X, D1->N ( c ) STITLE QDRTC ***************************************************************** ***************************************************************** ** ***Name: QDRTC ** ***Category: ** ***Abstract: Solves the quadratic equation A*Z^2 + B*Z +C ** where A, B, and C are complex. ** ***Entry Conditions: R4[A]= RSMEM ** ** SCRATCH AREA: PSCR1 - Re(A) (separated) ** PSCR2 - Im(A) (separated) ** PSCR3 - Re(B) (separated) ** PSCR4 - Im(B) (separated) ** PSCR5 - Re(C) (separated) ** PSCR6 - Im(C) (separated) ** ***Exit Conditions: (A,B),(R0,R1)= root of smaller magnitude ** R4[A]= RSMEM ** D1-> PSCR5 ** dec mode ** ** SCRATCH AREA: ** PSCR1 - Re(root of larger magnitude) (separated) ** PSCR2 - Im(root of larger magnitude) (separated) ** PSCR3 - Re(root of smaller magnitude) (separated) ** PSCR4 - Im(root of smaller magnitude) (separated) ** PSCR5 - XXX ** PSCR6 - XXX ** ***Stack Levels: 5 ??? ***Error Exits: ***Calls: ** ***Uses: ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** The algorithm uses the following fact: ** |Z+W| >= |Z-W| if and only if >= 0, ** where <,> denotes inner product. ** ** Let Z = -B + SGN(<-B,SQRT(B^2-4AC)>)*SQRT(B^2-4AC) where ** SGN(X) is defined as 1 if X is nonnegative and -1 otherwise. ** Then Z/2A is the root of larger magnitude and 2C/Z is the ** root of smaller magnitude. ** ** It is assumed that C#0. ** If A=B=0, then (1E5000,1E5000) is returned for both roots. ** If A=0 and B#0, then (1E5000,1E5000) and -C/B are returned as ** roots. ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** ***************************************************************** * * The following block of code replaces A by 2A, C by 2C, and * computes SQRT(B^2-4AC). * ***************************************************************** QDRTC GOSUB PSCRAD D1-> Re(A), hex GOSUB MPYU2+ Replace A by 2A, dec GOSUB D1+42 D1-> Re(C) GOSUB MPYU2+ Replace C by 2C GOSUB D1-126 D1-> Re(2A) GOSUB CFTCD1 (A,B),(R0,R1)= 2A GOSUB C->SCR01 (QSCR0,QSCR1)= 2A, D0->QSCR2 GOSUB D1+42 D1-> Re(2C) GOSUB CFTCD1 (A,B),(R0,R1)=2C,D1-> past PSCR6 GOSUB CpMPY Compute 4AC, D0-> past QSCR5 GOSUB CpNEG (A,B),(R0,R1)= -4AC GOSUB D0-84 D0-> QSCR2 GOSUB CSTOD0 (QSCR2,QSCR3)= -4AC, D0-> QSCR4 GOSUB D1-84 D1-> Re(B) GOSUB CFTCD1 GOSUB C->SCR01 (QSCR0,QSCR1)= B, D0->QSCR2 GOSUB CpMPY (A,B),(R0,R1) = B^2 GOSUB C->SCR01 (QSCR0,QSCR1)= B^2, D0->QSCR2 GOSUB CFTCD0 (A,B),(R0,R1)= -4AC GOSUB CpADD (A,B),(R0,R1)=B^2-4AC, D0->QSCR2 GOSUB CpSQRT (A,B),(R0,R1)=SQRT(B^2-4AC) trashes QSCR0 ***************************************************************** * * Here: (A,B),(R0,R1)= SQRT(B^2-4AC) * D1->PSCR5 * D0->QSCR1 * decimal mode * * The code below computes Z. * ***************************************************************** D0=D0+ 16 D0=D0+ 5 D0-> QSCR2 GOSUB CSTOD0 (QSCR2,QSCR3)= SQRT(B^2-4AC) GOSUB D1-42 D1-> Re(B) GOSUB CFTCD1 (A,B),(R0,R1) = B GOSUB CpNEG (A,B),(R0,R1) = -B GOSUB C->SCR01 (QSCR0,QSCR1)= -B, D0->QSCR2 GOSUBL multd0 (A,B)= Re(-B)*Re(SQRT(B^2-4AC)) GOSUB exab0 (A,B)= Im(-B) GOSUBL multd0 (A,B)= Im(-B)*Im(SQRT(B^2-4AC)) GOSUBL raddr01 (A,B)=<-B,SQRT(B^2-4AC)> C=A S C[S]=SGN(<-B,SQRT(B^2-4AC)>) GOSUB D0-42 D0-> QSCR2 GOSUB CFTCD0 (A,B),(R0,R1)= SQRT(B^2-4AC) ?C=0 S SGN(<-B,SQRT(B^2-4AC)> = 1? GOYES sgn=1 Jump if so ( bc ) GOSUB CpNEG (A,B),(R0,R1) = - SQRT(B^2-4AC) sgn=1 GOSUB CpADD (A,B),(R0,R1)= Z, D0->QSCR2 ***************************************************************** * * Here: (A,B),(R0,R1)= Z * D1->PSCR5 * decimal mode * * The code below computes the roots and exits. * *********************************************************** GOSUB C->SCR01 (QSCR0,QSCR1)= Z, D0->QSCR2 GOSUB D1-84 D1-> Re(2A) GOSUB CFTCD1 (A,B),(R0,R1) = 2A GOSUB CpDIV (A,B),(R0,R1)= Z/2A GOSUB D1-42 D1-> PSCR1 GOSUB CSTOD1 Store first root GOSUB D1+42 D1-> Re(2C) GOSUB CFTCD1 (A,B),(R0,R1) = 2C GOSUB C<->SCR01 (A,B),(R0,R1)= Z, (QSCR0,QSCR1)= 2C GOSUB CpDIV (A,B),(R0,R1)= 2C/Z GOSUB D1-84 D1-> PSCR3 GOTO CSTOD1 Store second root ( c ) STITLE CpNEG ***************************************************************** ***************************************************************** ** ***Name: CpNEG ** ***Category: ** ***Abstract: Negates a complex number. ** ***Entry Conditions: (A,B)= Re(Z) ** (R0,R1)= Im(Z) ** decmode ** ***Exit Conditions: (A,B)= Re(-Z) ** (R0,R1)= Im(-Z) ** decmode ** ***Stack Levels: 0 ***Error Exits: none ***Calls: ** ***Uses: A,C ** ***Author: ***Date Written: ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** CpNEG A=-A-1 S C=R0 C=-C-1 S R0=C RTN STITLE CpSQRT ***************************************************************** ***************************************************************** ** ***Name: CpSQRT ** ***Category: ** ***Abstract: Computes the square root of a complex number. ** ***Entry Conditions: (A,B)= Re(Z) ** (R0,R1)= Im(Z) ** ***Exit Conditions: (A,B)= Re(SQRT(Z)) ** (R0,R1)= Im(SQRT(Z)) ** D0-> QSCR1 ** decmode ** ***Stack Levels: 3 ***Error Exits: none ***Calls: CpABS,=RADDF,=DIVF,=SQRF,mult2,div2,xyex ** =STAB0,=STAB2,=EXAB0,=RCCD2,QSTO0,QRCL0 ** may exit thru =STAB0 or EXAB0 ** ***Uses: A,B,C,D,R0,R1,P,D0, QSCR0 ** ***Author: Paul J. McClellan ***Date Written: 2/26/87 ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ** Let r = SQRT(X^2+Y^2), a = SQRT((R+ABS(X))/2), b = Y/(2*a). ** Then ** ** { (0,0) (X,Y) = (0,0) ** SQRT((X,Y)) = { (a,b) X >= 0 ** { (ABS(b), SIGN(b)*a) X < 0 ** ***Test Classes: ** ***************************************************************** ***************************************************************** CpSQRT GOSUBL stab2 (A,B), (R2,R3)= X, (R0,R1)= Y GOSUB exab0 (A,B)= Y, (R0,R1)= X GOSUB QSTO0 [QSCR0]= Y, D0-> QSCR1, dec GOSUB CpABS (A,B)=r =|z|, (R2,R3)= X ?B#0 W r#0 ? GOYES cpsqrt2 ( Tc ) GOTO stab0 (A,B),(R0,R1)= 0 cpsqrt2 GOSUBL rccd2 C=0 S (C,D)= |x| GOSUBL raddf (A,B)= r+|X| GOSUBL div2 (A,B)= a^2 = (R+|X|)/2 GOSUB sqrf (A,B)= a GOSUB stab0 (R0,R1)= a GOSUBL mult2 (A,B)= 2a D0=D0- 16 D0=D0- 5 D0-> QSCR0 GOSUB getcd0 (A,B)= 2a, (C,D)= Y GOSUBL xyex (A,B)= Y, (C,D)= 2a GOSUBL divf (A,B)= b = Y/2a C=R2 C[S]=Sg(X) ?C=0 S Sg(x)="+" ? GOYES exab0+ then sqrt(Z)= (a,b) ( bc ) C=R0 C=Sg(a)&ex(a) C=A S C=Sg(b)&ex(a) R0=C (R0)=Sg(b)*a A=0 S sqrt(Z)= (|b|,Sg(b)*a) RTN exab0+ GOTO exab0 ( c ) STITLE CpABS ***************************************************************** ***************************************************************** ** ***Name: CpABS ** ***Category: ** ***Abstract: Computes the modulus of a complex number. ** ***Entry Conditions: (A,B)= Re(Z) ** (R0,R1)= Im(Z) ** decmode ** ***Exit Conditions: (A,B)= |Z| ** decmode ** ***Stack Levels: 2 ***Error Exits: none ***Calls: CpMAG^2,=RCAB0 ** exits thru =SQRF ** ***Uses: A,B,C,D,R0,R1,P ** ***Author: Paul J. McClellan ***Date Written: 2/26/87 ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** CpABS P= 14 C=R1 C= mantissa(Y) ?C#0 WP Y#0? GOYES ynot0 ( bc ) as=0 A=0 S RTN ynot0 ?B#0 WP X#0? GOYES cpabs2 ( bc ) GOSUBL rcab0 GONC as=0 (BET) ( c ) cpabs2 GOSUB CpMAG^2 (A,B)= |Z|^2, trashes R0,R1 sqrf GOVLNG =SQRF ( c ) STITLE CpMAG^2 ***************************************************************** ***************************************************************** ** ***Name: CpMAG^2 ** ***Category: ** ***Abstract: Computes the squared modulus of a complex number. ** ***Entry Conditions: (A,B)= Re(Z) ** (R0,R1)= Im(Z) ** decmode ** ***Exit Conditions: (A,B)= |Z|^2 ** decmode ** ***Stack Levels: 1 ***Error Exits: none ***Calls: x^2,=EXAB0,=RCCD0 ** exits thru =RADDF ** ***Uses: A,B,C,D,R0,R1,P ** ***Author: Paul J. McClellan ***Date Written: 2/26/87 ** ***Date Changed Programmer Reason ***------------ ---------- ------------------- ** ***Detail: ** ***Test Classes: all branches covered !!! ** ***************************************************************** ***************************************************************** CpMAG^2 GOSUBL x^2 (A,B)= X^2 GOSUB exab0 GOSUBL x^2 (A,B)= Y^2 GOLONG raddr01 (A,B)= X^2+Y^2 ENDCODE ( rootsC% ) ;