The Efficient Calculation of Powers of Polynomials by Ellis Horowitz Computer Science Department Cornell University Ithaca, New York 14850 TR 72 --H 129 The Efficient Computation of Powers of Polynomials by Ellis Horowitz Abstract Suppose we are given a polynomial P(x1,. ?Xr) in r > 1 variables, let m bound the degree of P in all variables th Xig ll. In a recent paper which compared the iterative versus the binary method it was shown that their respective computing times 2r r+l 2r were O(m n ) versus O((mn) ) when using single precision arithmetic. In this paper a new algorithm is given whose computing time is shown to beo((inn)r+l? Also if we allow for polynomials with multiprecision integer coefficients, the new algorithm presented here will be faster by a factor of mr?lnr over the binary method and faster by a factor of mr?l over the iterative method. Extensive empirical studies of all three methods show that this new algorithm will be superior for polynomials of even relatively small degree, thus guaranteeing a practical as well as a useful result. 1. Introduction Recently several analyses have appeared dealing with the problem of efficient methods for raising a polynomial to a power. The two major competing methods have been the iterative versus the binary approach. The iterative scheme would compute pfl as P, P2, ..., ?n?l, ?n whereas the binary method would use the binary expansion of n to decide when to square the partial result as well as when to multiply by P. In LKNU69], p. -401, Knuth implies that for the application of raising a polynomial to a power the iterative method is far inferior to the binary method since?the number of multipli- cations required in the latter case can be reduced from n to log2n. However, in a recent paper by Heindel, (HE172), he shows that in fact it is more efficient to use the iterative method rather than the binary method whether one is doing either single precision or multiprecision integer arithmetic on polynomials Moreover, if r is the number of variables of P, and we wish to raise P to the nth power, then Heindel shows that the r-l iterative scheme will be better by at least a factor of n in both cases. In this paper a new method for computing powers bf a polynomial is given which is faster than the binary method by a factor of mr?lnr?l for single precision arithmetic and faster than the iterative method by a factor of mr?l Also, this new algorithm will have at its core the use of the binary method applied to elements of a finite field, thus restoring what is perhaps our intuitive notion of the most efficient way of computing powers. In one paper related to this problem, M. Gentleman in [GEN7l] shows that for his specific model of polynomials the iterative method produces an optimal multiplication chain. In particular his model is for sparse polynomials where the number of nonzero terms remain 5 severly bounded; nor does he concern himself with the work needed for multiprecision coefficients. This paper is concerned with a model which assumes that polynomials are dense so that we can obtain worst case computing time bounds. Besides obtaining better asymptotic bounds on the time needed to compute ?n, the algorithm presented in this paper has been tested within an existing algebraic manipulation system and has been found to give a significant improvement over the iterative method for even reasonably sized problems. In Section 2 I will state some elementary results on computing times for basic operations. Then I will summarize the results of the binary and iterative methods. In Section 3 the new algorithm will be given plus an analysis of its computing time. In Section 4 empirical results comparing these several methods will be given. Finally in Section 5 some observations on even faster but only theoretical methods will be given. 2. Known Results In this section is summarized some of the basic computing times for operations 6n elements of a finite field, integers and multivariate polynomials. All computing times will be given in terms of the commonly used 0-notation. Definition2.1: Let f(x1,... ,Xn) and g(x11. ,Xn) be real-valued --H5--H functions defined on some common domain D. Then f = 0(g) if there exists a positive real number c such that 1f(x1,...,xn)i < cjg(x1,...,x?)I for all (x1,...,xn) ? We shall be concerned both with fixed point, single precision as well as arbitrary multiple precision arithmetic in our later analysis. For the single precision case we will assume that all numbers are elements of the finite field with p elements, designated as GF(p). The number p will generally be chosen to be about the size of a machine word. Then the time for adding, subtracting, multiplying or dividing elements in GF(p) will be 0(1). The computing times for operations on multi-precision integers will be given in terms of the maximum number of digits of the integers. Thus, given two integers, say d and e, the time for addition (subtraction) and multi- plication (division) is O(maxtln d, ln e)) and O((ln d) (ln e)) respectively. A new method for integer multiplication due to s. Schonhage and V. Strassen, [KNU69] p. 269takes time O(ln d(ln ln d) (ln ln ln d)) if d > e. However, their method seems to be of more theoretical interest rather than practically useful and so we will assume the use of a "classical" multiplication algorithm. We now need to have the times for performing operation?s on multivariate polynomials. Definition2.2 Let P(x1,...,xr) be a polynomial in r > 1 --H4--H variables with integer coefficients. Then, norm(P) is the sum of the absolute values of the numerical coefficients of P. Definition2.3 : Let Mr(d,m) = fP(x1,...xr): deg(P) in Xr +1 < m and d = norm(P)). Mr(d,m) is the set of polynomials in r variables whose norm is d and whose total number of terms is r at most m Theorem2.1: Let P(x1,. ,Xr) ? Q(x1,. .. ?xr) f?Nr(d,m). Then an upper bound on the computing time to form P + Q and P.Q using multiprecision arithmetic is o(m?(ln d)) and o(m2?(ln d) 2) respectively. Proof: The maximum number of non--Hzero terms in P and Q is mr. Thus for addition we must add mr terms, each addition taking O(ln d). For the proof of PQ see [11E172i p. 9. We note that if instead P and Q had 6oeffi'cients in GF(p), then all the additions and multiplications would take time 0(1) and the time for addition and multiplication for poly- nomials in r variables with coefficients from GF(p) would be o(m?) and o(m2?). There are some other operations for which we will need the computing times. In particular we need to know the times for applying modular and evaluation homomorphisms to polynomials and for applying interpolation and the Chinese Remainder method to integers and elements of GF(p). These --H5--H methods have already been fully documented in, for example, [COL7lJ, pp. 520-522. Let us reviewthose results here. Theorem2.2: Let A E Mr(d,m) such that A = and p a single precision prime. Then an the time to compute A* = (a.mod p)x1 o+.. xreir is o(m? ln d). l 0 variables with coefficients in some finite field GF(p) and an integer n > 0; Output: R*(x1,.. ?Xr) [P* (x11.. 1Xr)ln; 1) [Initial Case] If r=O, compute R* ??p*?? using the binary method; end. 2) [Initialize] mr ? deg(P) in xr; C(x1,... Ixr?l) 0; D(Xr) 1; [Interpolation] b -1; 3) [Choose next] b = b + 1; if b = p then stop. [point] 4) [Evaluation] P?(x1,...,xr?1) P*(x1,...,xr?1,b); [Homomorphism] If deg(P) in Xr < m then go to (3). 5) [Recursion] R*(x1,. . .,xr?1) ? CPPOWER(p,n); 6) [Interpolate] C(x1,. ,xr) R*(x1,...,xr?1)?C(x1,..,xr?1,b) D(b) 7) [Done?] 8) [end o+1] + D(xr) ? (xrb)D(xr); If deg(D) < mrn then go to (3); R*(x1,...,xr) ? C(xi,...,xr); end. Inorder to greatly simplify the analysis and following the procedure used by w,S. Brown in [BRO7l] let us assume that r (deg(P*) in x. + 1) < m for l log2n. (inn)rT0 (in) O((mn)r+l Now we give the algorithin for polynomials with integer cOefficients. AlgorithmPPOWER Input: P(x1,. .. ,xr) a multivariate polynomial in r>l variables with integer coefficients and n>l. Output: ?n 1) [Initialize) d ? norm(P); m deg(p) in xr; i ? 0; f ? (nd) 2) [Next Prime) i ? i+l; Choose next prime, ??;If none left, stop. 3) [Modular) P*(x1,. ..,xr) ? P(x1,...,xr) mod p1; [Homomorphism) If deg(P*) in xr < m, go to (2). 4) [R*=(P*)n) R*(x1,.. ,xr) ? CPPO?7ER(P*,n) 5) [Chinese) If i=l, R?R*; ????; go to (2); [Remainder) Otherwise find D(x1,.. ,xr) using the iterative Chinese Remainder algorithm such that 6) [Done ?] 7) D(x1,.. . ,xr)=R*(x1,...,xr)mod ??? D(x1,... ,xr)sR(xi,... ,xr)mod q; R ? D; g ? P??; If q < f, go to (2) Pn ? R; end. -? Theorem3.2: Let Pr(d,m) be the maximum computing time for Algorithm PPOWER(P,n) when P E M (d,m). Then ?r(d,m) --H O(mr+lnf+2 (ln nd) 2) Proof: Since P contains at most mr terms the time needed to compute d in step (1) is O(m? ln d). Steps (2) --H(6) are executed at most n(ln nd) times since f = (nd')? bounds the size of the norm of ?n Therefore we get the following table for the computing times: step time (1) O(m?ln d) (2) O(n ln nd) (3) O(m?ln d?n ln nd) < O(m?n(ln nd) 2) (4) (5) (6) O((mn)r+l in nd) = O(m?+l r+2 n n (in nd)) O((mn)r(n in nd) 2) = O(m?nr+2 (in nd) 2) O(n2(in nd) 2 Examining the table, steps (4) and (5) bound the computing time for the entire algorithm. Also 0(inr+?l nr+2 in nd + mr nr+2 (in nd) 2) < O(m?+inr+2 (in nd) 2) Looking at the proof of Theorem 3.2 we see that algorithm CPPOWER plus the time for the Chinese Remainder algorithm is what bounds the total time for this method. This is what we would normally expect. Now let us summarize the asymptotic computing times which have been arrived at for the binary, iterative and homomorphism algorithms. If e = in nd we get the following table: Single Precision Coefficients Multiple Precision Coefficients Binary (mn) 2r m2rn2r+2e2 Iterative m2rnr+i Homomorphism r+l m2rn r+2e2 mr+lnr+2e2 4. Empirical Results In this section we examine several tests which were?made on all three methods to see if we can determine just how practical the new homomorphism algorithm really is. All --H14- testing was done on an IBM 360/65 at Cornell University. The basic routines for the manipula?ion of multivariate polynomials were provided through the use of the SAC-i system for symbolic and algebraic calculation, [C0L69]. This is a collection of Fortran subroutines which allows in part for infinite precision integer and multivariate polynomial arithmetic. In particular the binary, iterative and modular algorithms were tested on univariate, bivariate and trivariate polynort?ials whose coefficients were randomly generated in the range [?2lO,2lO]. Tables I, II and III give the resulting computing times. Table I gives the results obtained for polynomials in one variable. The vertical column represents the degree of the polynomials while the horizontal row represents the successive powers to which the polynomials were raised. For every polynomial there are 3 rows of times which stand for the binary (BIN), iterative (ITR) and modular (NOD) methods. Recalling the computing times obtained in the previous section for one variable (r=l), in this case the time for all 3 methods is the same, O((mn) 2) This is reflected in Talbe I as we see there is no substantial difference in any of these methods for the degree and powers which were tested. Table II gives the results for polynomials in 2 variables. Now the vertical column represents the total degree of the polynomial P(x,y). This means that if m is the total degree of P(x,y) the highest power of both x and y appearing in P 2 is m. In particular P(x,y) c?n have no more than (m+l) terms --H15- and will often have that many since zero is a not--Hso-random random number. Examining Table II we see that in all cases the modular method outstripped the other two rather quickly. In fact, for total degree =4, 5,.and 6 the modular method was initially faster and became more so for successively higher powers. Computation limits were set at 4 minutes = 240,000 milliseconds of computing so we see that often the modular method permitted the calculation of one or two more powers-. This advantage will become even better as greater powers are tried. Table III is set up exactly as Table II, but now we have the results for polynomials in 3 variables. Again the modular method is better than the other 2. However, because of the large growth in the number of terms, practical limits in terms of core size as well- as computing time are soon reached. In any case tables II and III do show that the modular method becomes superior for relatively small powers and the theoretical analysis indicates that this advantage will grow larger as bigger problems are attempted. H ? Ln m 0 ? 0 ?) m ? a' H ? 0 H ? o ? o H H H ? H 0 o 0 0 Ln Ln H H 0 ? H ? n Ln ? m o 0 tr) in 0 Co H ?) H Co H Co 0) 0 (0 0 4J 0) w 0) H ? H ? 0 ? H ? Co a' ? a' cy, 0 Co H a' Co H ? Lo in a' a' Ln a' ? 0 0 ? Lo 0 H H ? cy) ? Co ? H H ? 0 ? H ? H Co `D cy) a' a' H H ? H H H in ? a' `D ? Co H in Co H ? Co m 0 Co in H ? Co H H z ? ? 0 H in Co ? in `D a' in `D ? H H- H Co ? Co H 0 a' H H -(0 H ? ? 0 o o 0) H ? (0 ? ?L? 0 H H ???? H o+H (0 0) a' ? H Co ? Co Co a' 0 ? ? 0 H o a' ? in -? Co H H H fri ? 0 fri --H16--H 0' H ? an ? 0 c? 0' ? n an a' a' r4 0' ? 0 `0 LO a' o ? H a' Cr) H LO ? H 0 ? `0 Lo H H `0 ? 0 ? 0 Lo m H 0 H H H Co H Lo ? H ? Co a' ? a' a' a' o an 0 `0 an. o ? n `0 H ? 0 H cf' --H17--H 0 0 0 0 Co H A 0 o ? o Ln o ? Co an A H 0 o ? o an Co o ? cf' Co `0 0 H t% 0 A `0 ? Co Co Co `0 an H 0 a' H `0 an ? cf' cf' a' a' a' an ? an an cf' ? Co `0 ? an H H 0 `0 `0 ? 0 H H H 0 a' a' Co Co cf' Co ? `0 an `0 ? H a' a' ? an H 0" r?0I o 0 o 0 o o o o o 0 A A o ? o `0 o an H o H `0 o H ? H `0 A H ? a' Co ? Co C? `0 Co Ln cf' a' an H Co `0 an o cf' Co cf' 0 `0 H cf' cf' ? ? an ? cf' cf' ? a' cf' H H H `0. `0 LO o cf' a' Co 0 `0 an cf' o ? cf' ? H a' ? cf' H H H 0 0 0 0 0 A an cf' H cf' a' o o o o o o an o 0 ? H 0 cf' A A cf' Co ? a' an 0 Co 0 an ? `0 H cf' 0 H ? a' H a' Co ? an an cf' `0' `n H 0 c) H 0 ? H Z (`) >` H ?? H ? ? H ? `n 0) H o+H o o o 0 o o o o o 0 A A 0 o ? o ? o ? `0 o r?i `0 Co A H a' o `0 0 H cf' H a' H ? an Co ? cf' a' an o Co H `0 a' Co H cf' cf' ? ? a' ? an a' H H a' cf' 0 H H H H H cf' H H H H H cf' cf' ? an Co `0 `0 ? `0 cf'j a' Co Co `0 a' 0 `0 ? 0 cf' cf' ? an 0) 0 0 0) 0 m H z ? ? 0 ? H a ? H ? 0 H ? 0 H ? z ? 0) H 0) 0 0) a' an cf' Co Co ? cf' an `0 --H18--H `0' Lnj 0 0 0 0 0 A 0 o cr? o ? o o H A H 0 0 o `0 `0 o 0 H ? 0 Co c'? 0 `0 A a' LO Lo ? Co a' Co `0 Ln Co Co a' a' ? a' `0 H 0 H H m Ln `0 0 `0 ? `0 Ln 0 `n 4) 0 0 4) 4) 4) 4) H 0 z ? H ? 0 0 0 0 0 A 0 o H o ? o `0 ? a' A H 0 o ? o c?) ? o Co ? ? a' m ? a' A H Co H tn ? 0 `0 Co Co Cfl Co Ln ? `0 `0 ? o H Co Co a' m H ? H H `0 H ? 0 0 0 0 0 0 A o o o o o o o ? o o Co o o c??? `0 A A H Co 0 ? `0 ? Co Lo LO Co o a' ? H H ? 0 0 o u H ? 4) H ? CO H ? -? o H H H -r' CO 4) H ? ? -19- 5. Even Lower Bounds Recently there has been intensive work and success on developing algorithins for basic arithmetic operations whose asymptotic bounds are as low as possible. In our case we are especially interested in an algorithm which multiplies two e-precision numbers in time M(e) = O(e log e(log log e)) which is due to S. Schonhage and V. Strassen, see (KNU69], p. 269. One can directly generalize their method so that one can multiply two m-degree polynomials with e-precision coefficients in time M(me). This is done simply by considering the polynomials' coefficients to be a long integer and inserting a su?ficient number of zeros to prevent carries. This bound of O(M(me)) is considerably better than the bound we have been using for univariate polynomial multiplication, namely O(m2e2). One can now ask what would be the computing time for forming P(x)n using this fast polynomial multiplication algorithm. o+th If we use the iterative method then at the 1 multipli- cation we?would be multiplying a polynomial of degree im and precision ie by a polynomial of degree m and precision e. Thus the total time for the iterative method is M(imie) n3me(log n me) (log log n me) l