HfvNQ LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 5IO.ft4 The person charging this material is re- sponsible for its return to the library from which it was withdrawn on or before the Latest Date stamped below. Theft, mutilation, and underlining of books are reasons for disciplinary action and may result in dismissal from the University. UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN DEC 1 4 1976 *P t o : J82 AUG 2 9 1997 L161 — 0-11W6 ' UIUCDCS-R-73-572 572- rru^ti COO-1469-0223 RATION — A Rational Arithmetic Package by R. L. Brown June 1973 UIUCDCS-R-73-572 RATION — A Rational Arithmetic Package* by R. L. Brown June 1973 DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBAN A- CHAMPAIGN URBANA, ILLINOIS 6l801 Supported in part by the Atomic Energy Commission under grant US AEC AT(ll-l)lU69. ABSTRACT In section A, the theoretical basis for representing a rational number losing the prime factorization of its minimal integer ratio is given, and the resulting algorithms for performing addition, subtraction, multi- plication, and division are developed. The conjecture that these algorithms will normally require fewer machine operations than algorithms based on storing rational numbers as two long integers is supported by a worst case machine operation count estimate for each. These operation counts are applied to parameters derived from test results. Section B is a self-contained user's manual for RATION, a rational arithmetic package written in FORTRAN and based on the above algorithms . Digitized by the Internet Archive in 2013 http://archive.org/details/rationarationala572brow 11 TABLE OF CONTENTS Page A. RATIONAL ARITHMETIC ALGORITHMS 1 1. Theoretical Foundations 1 2. Conversion Algorithms 1+ 3. Arithmetic Algorithms 7 k. Package Organization 11 B. RATION— USER'S MANUAL 17 5 . User Program 17 6. I/O and Utility Instructions 18 7. Arithmetic Instructions 20 LIST OF REFERENCES 22 APPENDICES I. Sample Program 23 II. RATION Listing 25 A. RATIONAL ARITHMETIC ALGORITHMS 1. Theoretical Foundations A package of FORTRAN subroutines has been programmed to perform arithmetic operations on a special representation for rational numbers using algorithms that attempt to minimize computer execution time at the expense of space required for the storage of the operands and the program itself. The algorithms and representations depend on some elementary properties of the sets of integers and rational numbers as developed in number theory. These properties are defined below. See [l] for fuller explanations . Definition 1 : The set I of integers consists of the counting (natural) numbers (1,2,3,...) and their negatives (additive inverses) and . Definition 2 : R is an element of the set R of rational numbers if it is any rati or/ 1 /I _, I e I_, I _ e I_. Note that this defines an equivalence class for each rational number. The rational number is the equivalence class (I *I )/(l *I ) for I e I. If Ip > and there exists no integer I, such that I* divides both I and I , then I /I is the minimal element of the equivalence class and will be referred to as the rational number R = I /I . See [6] for further explanation. Definition 3 : A prime number is any positive integer which is divisible without remainder by exactly two distinct integers — 1 and itself. Definition k : Any integer I can be uniquely represented by the product of some set of prime numbers. This unique set is called the prime factorization of I. Using the above concepts , one can see that any rational number involves two integers — each of which has a unique prime factorization. The prime factorization has a unique representation when placed in the form h = P l P 2 — P n a. > 1 where p is the i-th ordered distinct prime appearing in the factorization 1 and a. is the number of times it appears. i Since a. > 0, a representation for l/i p is ~" 6 1 ^2 T 3 m P-L P 2 ...P m e. > o. i Since the rational number R is I /I , a representation for R is a 1 -g 1 a 2 -§ 2 ft -g P l P --'Pq where {p.} is a reordering to include all elements of {p.}, (p.), and ft. = if p. A {p. } 1 1 T 1 3. = if p. £ {p.}. i i T l This representation can be used in a computer by storing a string of nodes, each node containing the exponent of a prime associated with the node. If we let an indexed array P(l) contain our primes with P(l) equal to the total number of primes defined in the array (this may be less than its dimension) and P(2) = 2 , P( 3) = 3, P(U) = 5 , P( 5) = 7,. . . , then each node will take the form of (#p. ) where y. is the exponent a.-0. of the prime p. with index #p. . l One additional node is needed to specify the sign of the rational number and the number of nodes occurring in this particular rational number's representation since this will seldom be the full length of the array. By using a FORTRAN array dimensioned R(2,K) , 10 <_ K <_ 25, we can set R(l,l) so signum[R(l,l) ] = signum[R] |R(l,l)| = number of nodes including the bookkeeping node. Also set R(2,l) = index of largest prime in the array, R(2,J) = index of prime number taken to the R(l,J) power (J >_ 2). In keeping with our notation, -1, 0, 1 have R(l,l) = -1, 0, 1, respectively, and R(2,l) = 1 by convention. Henceforth, a rational number R(2,K) in this form will be called the indexed prime factorization representation (factored form for short) of the rational number R. It will be written in the form used in the following examples . Ex. 1: Ex. 2 Ex. 3: A representation will also be needed for very long integers if numerators are to be added or subtracted. If a suitable base is chosen such that any number multiplying a digit of the integer will leave the product completely within one computer word , thus retaining the carry, an integer can * r-1 if A < signum[A] = / if A = I 1 if A > 1* be represented as a string of computer words, each word representing an (extremely large) digit. Again, a bookkeeping entry is needed to specify the number of digits defined and the sign of the integer I. Otherwise, all entries of the array would have to be zeroed before a change is made. Therefore, |P(l)| is set to the length of the array, signum[P(l) ] to signum[l], and P(j) to the coefficients of BASE**(J-2) for J >_ 2. The following examples illustrate the representation of long integers. The base is 10000 so the digits 0000 to 9999 are represent able. Ex. k: U7926U3 = (3,261+3,i+T9) Ex. 5: -177133= (-3,7133,17) -1 = (-2,1) Ex. 6: = (2,0) 1 = (2,1) 2. Conversion Algorithms It will be shown in section 3 that during addition, factored forms of two integers must be converted to long integer forms, added, and the resulting long integer converted back to the factored form. Two algorithms are necessary for these conversions. Consider the case of a factored form integer A (all positive exponents) and a list of prime numbers none of which exceeds one computer word in length. This is reasonable since n(x) — the number of prime numbers less than x — is approximated by x/log x [2] and on any IBM 360 computer this allows the program to allocate 10 prime numbers — much more than could ever be stored effectively. Then if long integer II is set to 1 and each element Il(2),... is multiplied by prime P(J) (carrying properly) and each P(J) is multiplied times II a. times as J is stepped through all primes with J positive a. in A, II is converted to a long integer form of A. This is J shown in Algorithm C, which, like all algorithms herein, is expressed in Knuth's notation [3]. ALGORITHM C(A,Il) : Convert factored form integer A to long integer II CI. [Initialize. ] Il(l) +■ 2, Il(2) -*- 1, M *- 2. C2. [Prime factor loop.] J •*■ 2. C3. [Exponent loop. ] K «- 1. Ck. [Long integer multiply. ] L *■ 2, Il(L+l) = 0. C5. [Prime P times II.] IL(L) <- IL(L)*P(A(2 ,J) ) . If I1(L) > BASE, I1(L) «- I1(L) - LI1(L)/BASEJ*BASE, Il(L+l) «- LH(L)/BASEJ. If L 4 M, then L +- L+l , go to C5. If Il(L+l) ^ 0, M^ M+l. C6. [End loop. ] If K 4 A(l,j) , K *■ K+l , go to CU. CT. [End loop. ] It J 4 | A(l,l)| , J «- J+l , go to C3. C8. [Set bookkeeping node.] Il(l) <- M*signum(A(l ,l) ) . END. Going the other way (converting a long integer to factored form) one must check if each prime divides the integer, and if so, add 1 to the exponent of that prime in the factored form result; if not, try the next prime stopping when II = 1. This is shown in Algorithm G. It is assumed that a function IDIV(ll,P) for dividing a long integer II "by a single word integer P is available. This would be similar to the multiply routine in C3-C6 of Algorithm C. ALGORITHM G(I2,A): Convert long integer II to factored form A Gl. [Initialize. ] M^2,I^2. G2. [Check for II = 1. ] If Il(l) = 2 and Il(2) = 1, then A(l,l) #P-,-n or 1A lis IA IJd y ± #p > #p T , or #p T . = ^P-r-r,. In each case, the node in the result, ( #p ) Ixi IA IA 1x5 1L which has its own marker IC , takes on the lower indexed prime and the corresponding exponent and increases the corresponding marker by 1. It also increases its own marker by 1 . If the prime index on both operands is the same, their exponents are added. In us , for the configuration above Y IC = 3 IB' Y IC = a iA' ° r Y IC = a iA + ^IB' res P ectivel y- Tnis leads to Algorithm M. ALGORITHM M(A,B,C) : C *- A*B MO. [Set multiply indicator.] IN * 1. Ml. [Initialize. ] IA * 2, IB*- 2, IC *- 2, II * MAX(A(2,l) , B(2,l)) , I *- 2. M2. [Search for next larger prime.] If A(2,IA) = B(2,IB) = I, go to M5. If A(2,IA) = I, go to M3. If B(2,IB) = I , go to Mh. Otherwise, I * 1+1. M3. [A factors C. ] C(2,IC) * I, C(l,IC) *■ A(l,IA), IA * IA+1, IC * IC+1. Go to M6. Mk.[B factors C. ] C(2,IC) * I, C(l,IC) * B(l,IB)*IN, IB * IB+1, IC * IC+1, go to M6. M5. [A and B factor C. ] C(2,IC) * I, C(l,IC) * A(l,IA) + B(l,IB)*lN. IA * IA+1, IB* IB+1, IC * IC+1. M6. [Loop end.] If I 4 II, I*- 1+1, go to M2; otherwise, if C(1,I) = 0,1* 1-1, go to M6. M7. [Bookkeeping node.] C(l,l) * I*signum(A( l,l)*B(l,l) ) . C(2,l) * C(2,I). END. For division, C = A/B— the same method of advancing through the nodes of A and B and comparing the indices of the primes holds , except 6 IB is re P laced by -3 IB since division is accomplished by subtracting exponents of common factors. This leads to Algorithm D, which sets to -1 a marker to be multiplied by 3 . ALGORITHM D( A, B,C) : C *- A/B Dl. [Set divide indicator.] IN* -1. D2.Go to ML. For addition, a least common denominator (l.c.d. ) of the two rational operands must be found. This is actually the least common multiple (l.c.m. ) of the two denominators. The denominators of the operands can he found easily since they are representable by extracting only the nodes with negative exponents. Their l.c.m. can be found by choosing the more negative of any negative exponents for any primes in either of the two factored form operands. Taking absolute values gives the least common denominator. This result is shown in Algorithm L. ALGORITHM L(A,B,R): R^ least common denominator of A and B LI. [Initialize. ] IA^2,IB^2,IR^2,I^2, II +■ MAX(A(2,1), B(2,l)). 12. [Search for next prime.] If A(l,IA) = B(l,IB) = I, go to L5. If A(1,IA) = I, go to L3. If B(1,IB) = I, go to LU. L3. [A factors R. ] R(l,IR) ■*■ -MIN(A(l,IA) ,0) . IA<- IA+1. If R(1,IR) 4 0, R(2,IR)«- I , IR ■*■ IR+1. Go to L6. LU.[B factors R. ] R(l,IR) ■*■ -MIN (B(l ,IB) ,0 ) . IB *■ IB+1. If R(l ,IR) 4 , R(2 ,IR) «- I , IR «- IR+1. Go to L6. L5. [A and B factor R. ] R(l,IR) = -MIN(A(l,IA) , B(l,IB) ,0). IA «- IA+1, IB * IB+1. If R(1,IR) 4 0, R(2,IR) +• 1 9 IR +■ IR+1. L6. [Bookkeeping node.] If I ^ II , I ^ 1+1 , go to L2 ; otherwise, R(l,l) + IR, R(2,l) *■ R(2,IR). 10 If D represents the l.c.d. , then for A,B rational numbers A*D, B*D are factored form integers. These integers can be converted to long integers "by Algorithm C. If each succeeding element of the corresponding integers IA(K) , K = 2(l)|lA(l)| , IB(L), L = 2(l)|lB(l)| is added with appropriate carries or borrows as in Algorithm I, the numerator can be found. Note that , , v _ /0 if x < step(x) -{ lifx 7 . Converting the numerator to factored form using Algorithm G and then dividing by D, the denominator, using Algorithm D, gives the result. Since the signs of the operands do not matter once they are attached to the long form integers , subtraction can be carried out by changing the sign of the subtrahend before and after addition. Thus, addition and subtraction are performed as in Algorithms A and S. ALGORITHM 1(11,12): II -e 11+12 11. [Initialize. ] IP1 «- | Il(l) | , IP2 «- |l2(l)|, I «- 2, II <- MAX(IP1,IP2) , CNEG +■ 1. 12. [Get first integer to be positive.] If Il(l)*I2(2) > 0, go to 15. 13. [Check node lengths.] If IP1 = IP2 , go to Ik \ otherwise, if IP1 > IP2 and Il(l) > 0, go to IT. Otherwise, go to 16. Ik. [Check integer lengths.] If Il(lPl) > I2(IP2) and 11(1) > 0, go to 17; otherwise, go to 16. 15. [If both < 0, reverse sign.] If Il(l) > 0, go to IT. 16. [Reverse signs.] Il(l)«- -Il(l), 12 (l)-*- -I2(l), CNEG «- -1. IT. [Add integers. ] Il(l+l) •*■ 0. Il(l)+- Il(l)*signum(ll(l))*step(lPl-l) + I2(l) *signum(l2(l))*step(lP2-l) + Il(l). If 11(1) > BASE, 11(1) * 11(1) - BASE, 12(1+1) <- l. If 11(1) < o, 11(1) «- 11(1) + BASE, 11(1+1) «• -1. If I 4 II, I «- 1+1, go to IT. 18. [Restore signs.] If Il(l+l) 4 , I +■ 1+1. 11(1) +■ I*CNEG. 11 ALGORITHM A( A, B,C) : C «- A+B AO. [Initialize. ] BNEG <- 1. Al. [Find least common denominator.] Go to Ll(A,B,R). A2. [A*R, B*R. ] Go to MO(B,R,T) . Go to M0( A,R,C) . A3. [Convert to long integer.] Go to Cl(C,Il). Go to Cl(T,I2). AU.[Add integers.] Go to Il(ll,I2). A5. [Convert back.] Go to Gl(ll,T). A6. [Divide numerator by denominator.] Go to Dl(T,R,C) A7. [Restore B. ] B(l,l) + B(l,l)*BNEG. ALGORITHM S( A, B,C) : C «- A-B SI. [Initialize. ] BNEG <- -1, B(l,l) «■ -B(l,l). S2. [A + (-B). ] Go to Al. k. Package Organization Since the four arithmetic operations to be implemented all contain two operands and one result, setting each operator up as a FORTRAN subroutine with three arguments (essentially a three address instruction) is an acceptable form. The user is also provided with limited I/O capabilities involving other FORTRAN subroutines as described in section 7. The two conversion routines from section 2 are additional routines called by the arithmetic operators. Since prime numbers are needed, there is a routine 12 that has some primes defined "by a DATA statement and it can generate more primes when required. If a large number of primes were available on a sequential access device, this would be primarily an I/O file routine. Since it is possible that one of the operands in an arithmetic operation may also be the result, all code is written so that the result is not finalized until all information is gotten from the operands. An inspection of the subtraction and division algorithms reveals that these need only prepare some indication that the appropriate operation is being done and then use the addition or multiplication algorithm, respectively. In dividing C = A/B , if IN is -1 wherever C(l,IC) = A(l,IA) + B(l,IB)*IN is performed, and IN = +1 when multiplication is being performed, then the division routine need only set IN = -1 and call the multiplication routine. Similarly, when the subtraction routine is called, by setting an indicator BNEG = -1 and changing the sign of B(l,l) for C = A+B , subtraction is performed by calling the addition routine. B(l,l) = B(l,l)*BNEG then changes the operand back to its original form on exiting the addition routine An alternative way to organize a rational arithmetic package would be to store two long integers for each rational number with tags denoting one as a numerator and the other as a denominator. The ordinary rules for arithmetic involving multiple digit fractions would hold. For example, A * 1 _ k B • Y M would be performed by 13 L' = A*X M' = B*Y R = greatest common divisor (g.c.d. ) of L', M' found by Euclid's algorithm [l] L = L'/R M = M'/R all performed using routines for arithmetic in which both operands are long integers . A X L Addition of — + — = — would be accomplished as follows. L* = A*Y + B*X M' = B*Y E = g.c.d. of L' , M' L = L ' /E M = M' /E Again, all arithmetic is performed on long integers. If machine instructions are divided up into sets according to comparable execution times, then it is possible to estimate the worst case 1 execution times of the algorithms needed with the two different representa- tions for rational numbers. Denote the set {multiply , divide} by *, and { add, subtract } by +. A long integer multiply with one integer m digits long and the other m digits long requires iil *m *'s and, at worst, iil *nu , +'s. Using this, consider the rational numbers —with A being n digits D 1 Y long, B being d digits long, and rr with X,Y having n , d digits, respectively. Then rational multiplication and division require, at worst, the number of operations indicated in Table U.l. The value K is the number of divisions lU required to find the greatest common divisor R of L', M 1 . Lame [l] has put the upper bound on K as ( BASE/2) while Knuth [3] estimates K as (12 log 2/iT 2 )log N, for N the minimum of L' , M' . It should he noted that the product of two integers will have, at worst, as many digits as the sum of the number of digits in the operands. Multipli cation Reason "1 "2 d l* d 2 K(n 1 +n 2 )(d 1 +d 2 ) (n 1 +n 2 )*min(n 1 +n 2 ,d 1 +d 2 ) ( d ] + d 2 ) *min (n^+rig A ± + ^ 2 ) 1 2 V d 2 K(n 1 +n 2 )(d 1 +d 2 ) ( n +n 2 ) *min ( n -+n 2 , d + d 2 ) (d L +d 2 )*min(n 1 +n 2 ,d 1 +d 2 ) L' = A*X M' = B*Y R = g.c.d. of L' , M 1 L = L'/R M = M' /R Addition n 1 * d 2 +n 2* d i K*max(n 1 +d 2 ,n +d ) *(d 1 +d 2 ) 3/2*max( ni +d 2 ,n 2 +d^) + n 1 *d 2 +n 2 *d 2 V d 2 K*max(n +d ,n +d ) *(d x +d 2 ) L' = A*X+B*Y M' = BY min(max(n +d ,n +d ) ,d +d ) min(max(n +d ,n +d ) ,d +d ) L= L'/E *max( n + d g ,n 2 + d ) *max( n + d 2 ,n 2 + d 1 ) min(max(n +d ,n +d ) ,d +c^) min(max(n 1 +d 2 ,n +d ) 5 d 1 +d 2 ) M = M' /E *(d 1 +d 2 ) *(^d 2 ) Table U.l Worst case operation count estimate for integer rational representation algorithms. 15 In the indexed prime factorization representation, the only multiplications required are due to conversion between factored form and long integers. Table k.2 shows the worst case operation count estimate for two rational operands A, p nodes longi and B, q nodes long. The value K 1 is the average exponent of the primes, K is the number of digits in the converted integer II, K is that for integer 12, and K "- 3 is the index of the largest prime in the final result. In the sample program in Appendix I (using base 1000 integers), average values for the above parameters were obtained. These appear in Table U.3- Multiplication * + Reason p+q C(1,I)=A(1,I)+B(1,I) Addition 1/2^*^ +q*K ) l/2K 3 *max(K 21 ,K 22 ) Table U.2. find g.c.d. R A*R, B*R Convert A*R-KE1, B*R-*I2 add integers convert to factored form result. Worst case operation count estimate for factored form representation algorithms. p+q p+q+2max(p ,q) l/UK 1 *(p*K 21 +q*K 21 ) 3/2 max(K 21 ,K 22 ) l/UK 3 *max(K 21 ,K 22 ) 16 PARAMETER VALUE n i ■ 2 n 2 2 *1 2 ^2 2 K l 2 K 21 2 K 22 2 P 3 q. 3 K„ 100 Table H.3. Average observed parameter values from test program A considerable savings in execution time for the factored form algorithms over the integer ration algorithms can be predicted on the basis of this test. Table k.k lists the comparative average values computed using these parameters. Knuth's estimate for K, based on the average 6 integer being 10 , is approximately 12. MACHINE OPERATIONS INTEGER RATIO FACTORED FORM * + * + * + 2 32 2 32 2 36 2^2 6 112 77 Table h.k. Computed comparisons using test result parameters B. RATION— USER'S MANUAL 5. User Program The RATION package currently available in the PUBLIC Library and listed in Appendix II was designed to minimize -user responsibility for manipulation of rational variables. The user's primary responsibilities in the FORTRAN calling program are to dimension the variables properly and to direct the sequence of I/O and arithmetic operations through FORTRAN subroutine calls. Several basic I/O facilities and all arithmetic operations are provided and diagnostics are available at two levels. The list of prime numbers is generated by the package automatically although all primes less than 1000 are placed in the array P at compilation time. If subroutine RATION is called, a second level of diagnostics can also be invoked. Other- wise, any routine can be called with no further preliminary calls. All rational variables used by the routine must be specified as A is below: INTEGER*2 A(2,K) where K is at least 10. Rational arrays can be created simply by adding the array dimensions after the two required dimensions. For example, B(2,10,5,3) is a (5,3) array of rational numbers. When these array elements are used as a subroutine calling argument, the (l,l) element is placed in the argument list. For example, CALL ENTER(B(l ,1,2 ,U) ) sets the (2,U) element of B equal to an integer read from cards. The sign of any rational variable is the sign of its (l,l) element; if B(l,l) < 0, B is negative. The program can be initialized by CALL RATION (N) where N is the value of K if all rational variables are dimensioned (2,K). If during execution any rational variable requires K larger than N, a message to that effect is printed. If N is set to 0, no checking is done. Any integers that are not factorable in terms of the available primes are placed at the end of the prime array by the conversion routine PTOG. If these integers (to the base PBASE) are more than one word long, all succeeding words of it are tagged with a minus sign. When converting back to integer representations, PI is initialized to this unfactorable integer instead of to 1. In this way, rational numbers of reasonable complexity can be computed. Default diagnostics are printed when an attempt is made to divide by zero. The quotient is set to zero and the program continues. If a conversion routine GTOP (see section 6) is called to convert a factored non-integer to integer form, a message is printed. 6. I/O and Utility Instructions Two input routines are provided for creating factored form integers. These can then be divided to create rational variable values. One output routine writes a fixed format rational as a signed two integer quotient. User calls for the subroutines are listed below. CALL ENTER(A): reads an integer from FORTRAN FILE 5 (card reader). The format of the required cards follows. The first card contains the least integer greater than or equal to the number of decimal digits di viced by four, right-justified in columns 1-5; -l s +1 ^ n columns 6, 7 indicate the integer's sign. The base 10 digits are then divided in groups of four from the least significant digit up. These groups are placed starting with the most significant group 19 in each four columns starting with column 1. For example, -177133 has input #00002-l# #00177133 #_ column 15 8 80 CALL CREATE (I,A): can be used if the integer I can "be represented in an INTEGER*H computer word. I is a signed integer which becomes the value of A on exiting the subroutine. CALL WRITE(A): outputs A in the form + 1 * _ X _ _ 2 Any additional comments or spacing can be provided by FORTRAN I/O statements in the user program. If the user feels he understands the representations for rational numbers and integers as presented in section 1, he can also use two utility routines that are called by the I/O and arithmetic routines. These provide the conversions described in section 2. If I is a long integer variable and A is a factored form rational variable, then CALL PT0G(l,A) transforms I to ; A and CALL GT0P(A,l) transforms A to I . Thus, the program in Figure 6.1. places the factored form of 177133 into A. The user can save execution time by preparing constants in this and similar ways, but only if he understands what he is doing. 20 INTEGERS A(2,10) INTEGER* U 1(3) DATA 1/3, 7133, 17/ CALL GT0P(I,A) Figure 6.1. Sample initialization of factored form integer 7. Arithmetic Instructions Four arithmetic routines are provided — one for each of the four operations closed on the rational numbers. RADD and RMPY are the principle routines ; RSUB and RDIV are supplementary entries to these using the 360 FORTRAN option of multiple entries [5]. In compilers without this feature, RSUB and RDIV could be very short routines that call RADD and RMPY, respectively, after setting a special variable in a COMMON array. The routines are called as follows for A, B, C factored form rational variables, CALL RMPY(A,B,C): C = A*B CALL RDIV(A,B,C): C = A/B CALL RADD(A,B,C): C = A + B CALL RSUB(A,B,C): C = A - B Any operators can be redundant. CALL RADD(A,A,A) will double A, but much more expensively than CALL RMPY(TW0,A,A) for TWO the factored form of two. Logical comparisons may be made using the fact that the sign of a rational variable is the sign of its (l,l) element. To compare A with B, 21 first compare the sign of A(l,l) with that of B(l,l). If this does not clarify matters , perform CALL RSUB(A,B,TEMP) and compare the sign of TEMP(l,l). Obviously, this is quite time- consuming since RSUB is the most expensive routine to execute; however, comparisons are possible. A sample program using all arithmetic routines and the I/O routines other than ENTER is in Appendix I. 22 LIST OF REFERENCES [lj Ore, 0., NUMBER THEORY AND ITS HISTORY, McGraw-Hill: New York (19U8). [2] Niven, I. and Zuckerman , H. S., AN INTRODUCTION TO THE THEORY OF NUMBERS, John Wiley & Sons, Inc.: New York (1972). [3] Knuth, D. E. , THE ART OF COMPUTER PROGRAMMING, vol. 1, Addison- Wesley: Massachusetts (1968). [k] Peterson, J. A. and Hashisaki , J. , THEORY OF ARITHMETIC, John Wiley & Sons, Inc.: New York (1967) . [5] IBM System/ 360 FORTRAN IV Language, S360-50, New York (1971) . [6] Wolf, F. L. , NUMBER SYSTEMS AND THEIR USES, Xerox College Publishing: Massachusetts (1971). 23 APPENDIX I Sample Program 2k IMPLICIT INTEGEK*2CA-H, J-0,b)-£) IMPLICIT INTEGEK*4C I ,P) COMMON /COM/ PLIM, PHASE, IPuC 100), PC 1 5000) DIMENSION CC 2, 2 5* 78 >* ON EC 2, 25), TMP (2, 25), 1 TW0(2,25> C C************************************************************ c C THIS PKOGKAM COMPUTES THE ELEMENTS IN A LO^EK TKIANGULAK C 12*12 MATKIX. CCK,M) IS KEPKESENTED AS CC K* C K- 1 ) /2 + M ) . C THIS PRESENTS S i'OKAGE OF SPAKSE ELEMENTS. C C* **************************************** ****************** c CALL KATIGNC25) C C C(2,2)= .5 C CALL CKEATEC 1,0NE) CALL CKEATEC 2# TWO) CALL KDW CONE, TtfG.»C( 1 , 1 » 3) ) C C C ) C DO 10 1= J, 12 CALL CKEATEC 1-1 , TMP ) CALL KMPY(TWO*TMP*TMP) CALL KDW(ONE*TMH*TMP) CALL KADDCCC \, \, CI-1 ) *( I -2) /2+2>* TMP*C< 1* 1»I*< 1-1 )/2+2) ) CALL Wril TECCC 1 , 1 , I*C 1- 1 )/2+2) ) 10 CONTINUE DO 20 1 = 3, 12 20 CC 1 , 1 , I*C I- 1 )/2+l )=0 C C C C K , M ) = C C K - 1 , M ) + CM-1)/CM*CK-1))*CCK-1,M-1) C M=3, . . .,K C Ks3#. ..# 12 C DO 30 IK=3,12 DO 30 IM=3,IK CALL CKEATEC IM*( IK- 1 ), TMP ) CALL CKEATEC IM-1 , ONE) CALL KDW CONE, TMP, TMP ) CALL KMPYC CCl,l,CIK-l)*ClK-2)/2+IM-l), TMP, TMP) CALL KADDCCC l,l,CIK-l)*CIK-2)/2+IM),TMP,CCl,l,IK*CIK-l)/2+IM)) CALL aIKITECCC1,1,1K*C1K-1)/2+IM)) IFCPC l ) .GE. 1 7 00) GO TO SI * CCK, 1 ) DO SO IK=2*12 KLM=IK*( IK- 1 )/2+ 1 CALL CKEATEC 1 ,CC 1 , 1 ,KLM) ) DO 40 IJ=2,1K CCl,l,lK*CIK-l)/2+IJ) = CCl,l,IK*CIK-l)/2+IJ)*C-l)**CU+l) CALL KADDCCCl,l,IK*ClK-l)/2+IJ),CCl,l,KLM), 1 CC 1 , 1 ,KLM) ) 40 CONTINUE CALL /JK1TECCC 1 , 1 ,KLM ) ) 50 CONTINUE 5 1 S TOP END c 30 CONTINUE c c c CCK, 1 ) = + S UM C - 1 ) * * 1 1 = 1 , . . . ,K 25 APPENDIX II RATION Listing 26 Q********************************************************************** C C C * * * KATIGiM * * * c c c********************************************************************** c C A PROGRAM BY C R.L. BROWN C MARCH 6, 19 13 C C* ********************************************************************** * c c C THIS PACKAGE READS AND CREATES INTEGERS, FORMS RATIONAL NUMBERS, ADDS C AND SUBTRACTS, MULTIPLIES AND DIVIDES KATIONAL NUMBERS, AND OUTPUTS C RATIONAL NUMBERS. ALL KATIONAL VARIABLES ARE INTEGER*2 AND SHOULD C BE DIMENSIONED CC2,K), K>10 BUT USUALLY LESS THAT 25. INTERNAL C REGISTERS ARE PLACED IN COMMON BLOCKS COM, RT, P12, TO SERVE SUCH C PURPOSES AS STORING INTEGER*4 PRIME NUMBERS AND HOLDING C INTEGEK REPRESENTATIONS IN FINDING THE NUMERATOR DURING RATIONAL C ADDITION. THE SUBROUTINES AVAILABLE ARE: C ENTER(A): PLACES AN INTEGEK INPUT FROM CARDS INTO RATIONAL C VARIABLE A. SEE SUBROUTINE WRITEUP FOR CARD FORMAT. C C CREArE(I,A): PLACES ANY INTEGER OF LESS THAT 1 IF C I J .NE.IM)GO TO 9 7 DO 8 ID=2,IMM 8 PI C ID)=PTC ID) 9 CONTINUE IM= IMM 10 CONTINUE 11 PI C 1 ) = PI C 1 )*IM RETURN C IF A=<-1,0,+1> THIS SECTION FIXES PI TO THE PROPER FORM FOR THESE C SPECIAL VALUES. 12 PI C 1 ) = 2 P 1 C 2 ) = II 1 K ( 1 1 .N E . ) P 1 C 1 ) = P I C 1 ) * A C 1 , 1 ) KE I URN C ERROR MESSAGE IF A IS NOT AN INTEGER CPO/JER OF A PRIME < 0). 1 I WK ITE( 6,999) 999 FORMATCSX, 'iriOP CALLED /JITH NON-INTEGER ARGUMENT') I URN F.ND 28 SUBROUTINE PTOGCPI,A) IMPLICIT INTEGER*2CA-H, J-0,U-Z) IMPLICIT INTEGER*4C I,P) COMMON /COM/ PLIM,PBASE,PTC 1 00 ) > P C 1 0000) DIMENSION PI ( 1 )>A(2j. 1 ) DATA IL/1/ C c C THIS KOJTINE CHANGES INTEGER PI TO ITS INDEXED PRIME REPRESENTATION C FORM BY DIVIDING BY SUCCESSIVE PRIMES AND RECORDING THE FACTORS. C PI IS REDUCED TO 1 BY THE SUCCESSIVE DIVISIONS. C IE THE AVAILAHLE PRIMES DON'T COMPLETELY FACTOR PI, THE UMFACTOKED C PART IS APPENDED TO THE ARRAY OF PRIMES. C C c c C SAVE INTEGER BASE AND INITIALIZE PARAMETERS. C II = IABSCPI CI)) IF(PI(2) .LE. 1 .AND. I 1 . EU . 2 ) GO TO 4 5 3 DO 2 1=1,11 2 PTCI ) = AC \, 1 ) = PI ( 1 )/II AC l,2)=0 AC2,2)=0 IA=2 IM = 2 1 PD=PICII) OUT = PTC 1 ) = PI C 1 ) c C DIVIDE BY PCIM) AND SEE IF IT'S A FACTOR C DO 10 1=2,11 IJ = I I -1 + 2 PE = PD/PCIM) PTC IJ) = PE IFC I J.EU. IABSCPTC 1 ) ) .AND.PE.EU.O) 1 PTC 1 ) = ISIGNC IABSCPTC 1))-1,PTM)) PD = CPD-PE+PC IM) )*PBASE C C SET OUT TO 1 IF PICIM) IS NOT A FACTOR C IFCIJ.EU.2.AND.PD.NE.0) OUT = 1 PD = PD + PI C I J- 1 ) 10 CONTINUE IF COUT.EU. 1 ) GO TO 30 C C PCIM) WAS A FACTOR; ADD IT TO A C IFC IM. EU .AC 2, I A) .OR. AC 2, 1A) .EU.O)GO TO 15 IA=IA+1 IFC IA.GT .PLIM)GO TO 4 1 AC 1, IA) = 15 AC 1, IA)=AC 1, IA)+ 1 AC 2, 1A) = IM I I = IABSCPTC 1 ) ) DO 20 IJ= 1 , I I 20 PI C IJ) = PTC IJ) IF C I I .EU. 2 .AND. PI C2) . EU . 1 ) GO TO 40 GO TO 1 30 IF C IM+ 1 .GT.PC 1 ) .AND. PC 1 ) .LI".«000) CALL GENERP I M = I M + 1 IF C IM.LE.HOOO)GO TO 1 IFCAC2,IA).NE.0)IA=1A+1 AC2, IA)=HOOO+IL AC 1, IA)= 1 DO .)', 1 = 2*11 35 PC8000+I-2+lL)=PICI)*iSlONCl,5-2*I) 29 IL=1L+II-1 PCIL+8000)=-9999999 GO TO 40 41 WRITEC6*999) 999 F0RMATC5X, 'FACTORED 40 AC 1 > 1 )=IA*AC 1, 1 ) AC2, 1 )=AC2, IA) 42 RETURN c c IF PI = <-l,0,+l>, SET c FORM VARIABLE NOT LARGE ENOUGH*) A TO ITS SPECIAL FORM. 4 5 AC 1, 1 ) = AC2,1) = RETURN END PI C2)*PI C 1 )/I I 1 THIS ROUTINE GENERATES THE NEAT PRIME IN THE ARKAY MPCI) JSING KNUTH'S ALGORITHM P, SECTION 1.3-2 MPCM+1) CONTAINS A POSSIBLE PRIME NUMBER. MPU IS THE UUOTIENT OF DIVIDING MP(M+1) BY A KNOWN PRIME; PR IS THE REMAINDER. SUBROUTINE GENERP IMPLICIT INTEGER*2CA-L, N-^) REAL*8 PR COMMON /COM/ MLIM, MBASE, MTC 100), MPC 10000) C C********************************************************* ********* c c c c c c c c c c* * * * * * ************************************************************* c M = MP ( 1 ) MPCM+1 ) = MPCM) + 2 1 DO 10 1*2* M MPU = MPCM+ 1 )/MPC I ) IF CMODCMPCM+ 1 ),MPC I ) ) . EU .0)G0 TO 12 IFCMPQ.LE.MPC I ) ) GO TO 11 10 CONTINUE 11 MPC 1 ) = M+ 1 RETURN 12 MPCM + 1 )=MPCM+ 1 )+2 GO TO 1 END BLOCK DA COMMON / INTEGER* INTEGER* +53*59,61 +137, 139, + 223,227, +307, 31 1, +389, 397, +479, 487, + 577, 587, + 659,661, +757, 761, +857,859, +953,967, END TA COM/ ML I 4 MBASE/ 4 MP/ 169 ,67,71,7 149, 151, 229,233, 31 3, 317, 401, 409, 49 1 , 499, 593*599* 673,677, 769, 773, 863,877, 971 ,977, M, MBASE, 10000/, M ,2, 3, b, / 3,79,83, 157, 16 3, 239,241, 331, 337, 419,421, 503, 509, 601,607, 68 3,6 9 1, 787, 797, 881 ,883, 983,991, MTC 1 LIM/ ,11, 89,9 167, 251, 347, 431, 52 1, 613, 701, 809, 887, 997, 00), MPC 10000 9999999/ 1 3, 1 7, 19,23, 7,101, 103, 10 173, 1 79, 181, 257,263,269, 349, 353, 359, 433, 439, 4$3, 523, 541, 547, 6 17, S 19,6 31, 709, 7 19, 727, 8 1 1,82 1,823, 907,91 1,919, 9831*0/ ) 29,3 7, 10 19 1, 271, 36 7, 449, 557, 641, 733, 827, 929, 1, 37 9, 1 1 193, 277, 37 3, 457, 56 3, 643, 739, 829, 937, ,41, 43, 47, 3^ 127, 131, 197, 199,21 1 j 281,283,293; 379, 383, 46 1, 46 3; 5 69, 57 1, 647,6 53, 743, 751, 839,853, 941,947, 467i 30 SUBROUTINE KMPYCA,B,C) IMPLICIT INTEGEK*2CA-H, J-0,U-Z) IMPLICIT INTEGER*4C I,P) INTEGER* 4 MINCMAXO COMMON /COM/ PLIM,PBASE,PTC 100), PC 10000) DIMENSION C(2* 1)*AC2* 1)»B(2* 1)»D«2»50) C THIS ROUTINE MULTIPLIES OK DIVIDES BY ADDING OR SUBTRACTING C THE EXPONENTS OF THE PRIME NUMBERS IN THE FACTORIZATION OF THE C OPERANDS. IF A PRIME IS DEFINED FOR A BUT NOT FOR B, ITS EXPONENT C IN B IS 0. ENTRIES ARE: C RMPY: A*B C KDIV: A/B C* ************************************ **************************** C INITIALIZE PARAMETERS C IN = 1 1 PA=AC1,1) PB = BC 1, 1 ) PA=IABS(PA) PB=IABSCPB) IF (PA*PB.LE.l ) GO TO 25 ID=2 IA = AC2, 1 ) IB=BC2» 1 ) II=MAX0C IA, IB) IA=MIN0CPA,2) IB=MIN0CPB,2) DC2, 1 )=1 DO 10 1=2,11 IF C AC 2, I A) *1 .EO. I )G0 TO 3 GO TO 4 3 IF CAC2, IA) .EU.BC2, IB) ) GO TO ^ DC 1, ID)=AC 1, IA) DC2, ID)=AC2, IA) IA=MINOC IA+l.PA) GO TO 9 4 IFCBC2, IB)*1 .NE.I ) GO TO 10 DC 1, ID) = BC 1, IB)*IN DC2, ID)=BC2,IB) IB=MIN0CIB+1,PB) GO TO 9 6 DC 1, ID)=AC 1, IA)+IN*BC 1, IB) DC2, ID)=BC2, IB) IA=MINOC IA+ 1,PA) IB=MINOC IB+ 1,PB) 9 IFC I . EU. I I )G0 TO 10 IF CDC 1, ID) .NE.O) I D= I D+ 1 IF CID.GT.PLIM)GO TO 13 10 CONTINUE C D IS TEMPORARY STORAGE SINCE A OR B MAY EUUAL C IN THE CALLING C ROUTINE. IF CDC 1, ID) .EU.O)ID=IU- 1 DC 1, 1> = ISIGNCIDjA< 1, 1 )*BC 1, 1 )* 1 ) DC2, 1 )=DC2, ID) GO TO \A 13 WRITEC6,999) 999 FORMATC 5X, 'FACTORED FORM VARIABLE NOT LARGE ENOUGH') 14 DO 15 1= 1 , I D CC 1, I )=DC 1,1) 1 5 CC2, I )=DC2, I ) RETURN 1,0, + 1>, COMPUTE THE SPECIAL FORM OF C PRINT AN ERROR MESSAGE. C C c c IF BOTH A IF B IS A AND B ARE IN <- ZERO DIVISOR, P 25 CC 1, 1 ) = AC 1, 1 )*BC 1,1) CC2, 1 )=1 IFCPB.EO.U .AND. IN.EU.-I) WK I TE Co , 1 000) 1000 F0KMATC5X, 'DI VISION BY ZEROJ QUOTIENT SET TO ZERO') RETURN ENTRY R[)IVCA,B,C) IN = -1 GO TO 1 END 31 SUBROUTINE KADD(A,B*C) IMPLICIT INTEGER*2(A-H, J-G,u-^) IMPLICIT INTEGER*4( I,P) INTEGER*4 MINO,MAXO COMMON /COM/ PLIM,PBASE,PTC 100), P( lOOOO) COMMON /KT/R(2,50),T(2, 50) DIMENSION AC2, 1),B(2, 1 ),CC2, 1 ) COMMON /P12/ P1C50)*P2(50) C C c c C THIS ROUTINE ADDS OK SUBTRACTS A AND B, PLACING THE RESULTS INTO C C ENTRIES ARE: C KADDCA,B,C) 2 A+B=C C RSU8(A,B,C): A-B=C C C Q+ if if if if ^c ^c if :+: if if if if if ■%. if if if if if if if -if. if if if if if if if if if if if if ^c +%^c*: + ^t + ++5f:^t^;^ + + + ^:+:^;+3f:+^c%^:^c^:+:^a|c + c C PF(I,P) IS USED IN ADDING TWO INTEGEKS OF DIFFERENT LENGTH C BY MERGING THE SHORTER WHEN NECESSARY. C PFC IZ,PC)=MAXO(PC-I^+ 1jO)/MAXO(PC-I^+ l>l) BNEG = l GO TO l ENTRY RSUBCA,B,C) BNEG = -l BC l, l ) - -B( l, l ) 1 PA = AC l, l ) PB = B( l , l ) PA=IABS(PA) PB=IABS(PB) IF( < AC 2, PA- l ) . GT.8000.AN0.PA.NE.2) .OR . ( BC 2, PB- l ) .GT.8000 l .AND.PB.NE.2) )GG TO llU CNEG = l DO 2 I = l , l 00 KC l , I ) = 2 PT(I ) = RC l,2)=0 R(2, l)=l IR = 2 II=MAX0(A(2, l )*1#B<2» 1 )*1) IA=MIN0(2,PA) IB=MIN0(2,PB> IFUI .LE. 1 )G0 TO 13 C C FIND THE LEAST COMMON DENOMINATOR BY STOKING THE MOST NEGATIVE C OF ALL NEGATIVE EXPONENTS OF PRIMES. C DO 10 1=2,11 IF(A(2,IA) .EG. I )G0 TO 3 GO TO A 3 IFCAC2, IA) .EU.BC2, IB))GG TO 6 IF (AC 1, IA) .GE.O)GG TO 5 R( 1, IR)=A( 1* IA) H<2, IR)=A(2, IA) IR = MINO( IR+1, I I ) 5 IA=MINO( IA+1,PA) GO TO 10 A 1F(B(2, IB) .NE. I ) GO TO 10 IF CB( 1, IB) .GE.O) GO TO 7 R( It IK)=B< 1, IB) K<2, IR)=BC2, IB) IR=MINO( IR+ 1*11) 7 IB=MINO< IB+ 1,PB) GO TO 1U 6 R( 1, IK)=MIN0( 1*A< 1, I A), 1*B( 1, IB) ) IF(R( 1, IR) .QE.OGO TO 9 R<2»IH)=A(2#IA) IR = MINO( IR+ 1, I I ) 9 IA=MINO( IA+ 1,PA) IB=MINO( 18+1, PB) 10 CONTINUE 32 C IF(KC 1 , i r ) . GE. 0)IR= IK- ■ 1 KC 1 , 1 ) = IR Hi 2 , 1 ) = K(2, IK) GO ro 1 4 K< 1 * i > = 1 R(2 , i ) = 1 13 C C PUT PARTIAL NUMERATORS INTO T AND C BY A*K, B*R, WHERE R IS THE C LEAST COMMON DENOMINATOR. CONVERT TO INTEGERS P1,P2. C 14 CALL RDIV(B,K,T) CALL RDIV(A,K,C) CALL GT0P(C*P1) CALL GTGP(T*P2) C C ADD P2 rO PI AFTER ADJUSTING SIGNS TO GET A POSITIVE RESULT. C IP1 = IABSCP1 ( 1 ) ) IP2 = IABS(P2( 1 ) ) IF(P1 ( 1 )*P2( 1 ) .GT.O) GO TO 7 IFUP1-IP2) 30,20*40 20 IFCPK IP1 ) .GT.P2( IP2) ) GO TO 40 30 IFCP1 ( 1 ) .LT.O) GO TO 60 35 PI ( 1 ) = -P 1 ( 1 ) P2( 1 ) = -P2( 1 ) CNEG = - 1 GO TO 8 40 IF(P2( 1 ) .LT.O) GO TO 80 GO TO 3 5 70 IF (PI ( 1 ) .LT.O) GO TO 35 80 II = MAXOC IP1, IP2) PT(2)=0 DO 90 1=2, I I PT(I+1 )=0 PTCI) = P 1 ( 1 )/IPl *P1 ( I )*PF( I, IP1 )+PF( 1, IP2)*P2( 1 )/IP2*P2( I ) 1 +PTCI) IF (PTC I ) .GE.O) GO TO 8 5 PT(I) = PT(I) + PBASE P T ( I + 1 ) = PT(I+1) - 1 GO TO 9 85 IF(PT( I ).LT. PBASE) GO TO 90 PT(I) = PT(I) - PBASE PT( 1+ 1 ) = PT( 1+1 ) + 1 90 CONTINUE IF(PT( I 1 + 1 ) .NE.O) 11 = 11 + 1 DO 100 1=2,11 100 PI ( I ) = PT( I ) Pl(l) = II * CNEG C C CONVERT FINAL NUMERATOR TO INDEXED PKIME FACTORISATION FORM, C DIVIDE BY K, AND PLACE ANSWER IN C. CALL PT0G(P1,T) B( 1, 1 ) = B( 1, 1 )*BNEG CALL RMPY(T,K,C) RETURN 1 10 a)RITE(6, 1000) 1000 F0RMAT(5X, 'ADDITION OPERAND HAS MORE THAN ONE UNFACTORED 1, 'COMPONENT; RADD JILL YIELD INCORRECT ANSWER.') RETURN END 33 C* C C C C C C C* SUBROUTINE WKITECA) IMPLICIT INTEGER*2CA-H,J-0,U-/) IMPLICIT INTE6ER*4C I*P) INTEGEK*4 MAXO>MINO COMMON /COM/ PLIM,PBASE.»PTC 100), PC lOOOO) COMMON /RT/RC2,50),TC2, 50) DATA IDASH/' '/ COMMON /P12/P1 C50)>P2C50) DIMENSION AC 2, 1 ), IFMC30), I FORM (4 ) > IF 1 ( 1 2 ) > I F2( 1 2 ) , I F 3(S ) EQUIVALENCE C I FM C 1 ) , IF 1 ( 1 ) ) , ( I FM ( 1 3 ) , IF2C 1 ) ) , C I FM ( 2 5 ) , I F3C 1 ) ) DATA IF1/' C 1H+,23X,27I4) C I H + , 39X , 221 4, /7X, 301 4 ) '/ DATA IF2/' C 1H+,55X, 19I4,/7X, 3014) ( 1 ri+ ,1 IX , 1 51 4, /7X , 301 4 ) V DATA IF 3/' ( 1H + ,87X, 101 4,/7X, 3014) '/ THIS ROUTINE WRITES OUT TWO INTEGERS TO REPRESENT A BASE 10 RATIONAL NUMBER II = AC 1, 1 ) II = IABSCII) IF CII .LE. 1 ) GO TO 26 IA = 2 IR = 2 IT = 2 PLACE THE NUMERATOR (POSITIVE EXPONENTS) (NEGATIVE EXPONENTS) IN K. RC2, 1 )=1 TC2, 1 )=1 CONVERT BOTH TO INTEGERS AND WRITE. 2 IFCAC }j IA) )4,4,6 4 R(2, IR)=A(2, IA) RC 1, IR)=-A( 1 , IA) IK=IR+ 1 GO TO 7 6 TC2, IT)=AC2> IA) TC \, IT)=AC 1, IA) IT=IT+1 7 IA=IA+1 IF ( IA.LE.I I )G0 TO 2 RC 1, 1 ) = IK-1 RC2> 1 )=RC2, IR- 1 ) TC 1, 1 ) = IT-1 TC2, 1 )=TC2, IT-1 ) 25 II = AC 1, 1 )/I I IK=IR-1 I T= IT-1 GO TO 29 26 TC 1, 1 ) = AC 1,1) TC2> 1 ) = 1 R( 1, 1 )=1 RC2, 1 )= 1 29 CALL GTGPCTjPI ) IS = 30 IFCTC2, IT- 1 ) .LE.HOOO.OR.IT. EU . 2 ) GO TO 33 IT=IT-1 IS=IS+ 1 ISS=-1 PTCIS)=99999 31 IS=IS+1 ISS=ISS+1 PTCIS)=IABSCPCTC2, IT) + ISS) ) IFCPCTC2, IT) + ISS+1 ) .L T.O.AND. 1 -9999999)G0 TO 31 GO TO 30 33 IS=IS+1 IN TJ DENOMINATOR PCTC2, ID + ISS+ 1 ) .GT. 3k PTCIS)=99999 I 1 = IABS(PJ C 1 ))-l IFR = MGD< I l,20)/4 DO 9 9 10=1,6 99 IFGKMI = l.Il) IFCIS.GT.l )WKlTEC6> I FORM )C PTC IS+2- I ) > I =2* I S) I'S=0 CALL GTGPCk,P2) 34 IF CRC2, Ir<- 1 ) .LE.8000.GK .IT. Eu .2)GG TO 37 IK=IK-1 IS=IS+1 I S S~- 1 PTOS)=99999 35 IS=I$+1 ISS=ISS+1 PTCIS)=PCRC2, IR)+ISS) IFCPCRC2,IR)+ISS+1) .LT.G.AND.PCKC2, IiO + ISS+1 ) .GT. 1 -9999999) GG TG 35 GG TG 34 37 I S= 1S+ 1 PTCIS)=99999 I2=IABSCP2C 1 ) )-l IJ=MINOCMAXOC I 1 , 12)* 20) WRITE (A* 1002) II, C I DASH, IK= 1, IJ) IFR = M0DCI2,20>/4 DO 98 10=1,6 98 IFOKMC IG)=IFMC IG+6*IFR) WK I TEC6>1001)CP2CI2+2-I), 1=1,12) IFCIS.GT.I) WKITEC6, 1F0KM) C PT C I 5 + 2- I ) , I = 2, I 5 ) RETURN 1001 FOKMATC 5( 6X, SC IX, 201 4/) ) ) 1002 FGKMAT(2X, 12, * * ',20A4) END SUBROUTINE CREATECI,A) IMPLICIT INTEGER*2CA-ri, J-0,U- IN COLUMNS 6,7. C FOLLOWING CARDS: C DIGITS OF THE INTEGER IN GROUPS OF FOUR DIGITS, RIGHT JUSTIFIED C WITHIN GROUPS, GROUPS LEFT JUSTIFIED. C EXAMPLE: -177133 C CARD ONE #00002- 1# C CARD TWO #001771 330000.. .0000# C C C* **************************************************************** READC5, 1000) IN,ISIGN READC5, 1001 ) (PC IN+2-I ), 1= 1, IN) P( 1 ) = CIN+1 )*ISIGN CALL PTOG