From:	IN%"asadikou@email.gc.cuny.edu" 25-OCT-1996 12:59:36.62To:	IN%"caviness@louie.udel.edu", IN%"vpan@lcvax.lehman.cuny.edu"CC:	Subj:	JSC Manuscript #1471 (author: V. Pan).Return-path: <asadikou@email.gc.cuny.edu>Received: from vice1.lehman.cuny.edu by lcvax.lehman.cuny.edu (PMDF #12809) id <01IB276UFX9C0003PL@lcvax.lehman.cuny.edu>; Fri, 25 Oct 1996 12:59 ESTReceived: from gaudi.gc.cuny.edu by vice1.lehman.cuny.edu (5.51/SMI-3.0DEV3) id AA02878; Fri, 25 Oct 96 12:37:04 EDTReceived: from broadway.gc.cuny.edu (broadway.gc.cuny.edu) by GAUDI.GC.CUNY.EDU (PMDF V5.0-6 #13310) id <01IB26WC8RHC008MLN@GAUDI.GC.CUNY.EDU>; Fri, 25 Oct 1996 12:51:02 -0400 (EDT)Received: from [146.96.180.65] ([146.96.180.65]) by broadway.gc.cuny.edu (8.7.6/ank-osf96) with SMTP id MAA25452; Fri, 25 Oct 1996 12:50:55 -0400 (EDT)Date: Fri, 25 Oct 1996 12:54:53 -0500From: asadikou@email.gc.cuny.edu (Akimou  Sadikou)Subject: JSC Manuscript #1471 (author: V. Pan).To: caviness@louie.udel.edu, vpan@lcvax.lehman.cuny.eduMessage-id: <v01530502ae96ae17f0cb@[146.96.180.65]>Content-type: multipart/mixed; boundary="============_-1365856403==_============"Content-transfer-encoding: 7BITX-Sender: asadikou@broadway.gc.cuny.edu (Unverified)Mime-Version: 1.0--============_-1365856403==_============Content-Type: text/plain; charset="us-ascii"Prof V. Pan asked me to send you this file. It is the LATEXversion of his accepted paper for  JSC Manuscript # 1471.You should let me know if you have any problem about my file.Thanks.A. Sadikou..--============_-1365856403==_============Content-Type: text/plain; name="fxmsplod.tex"; charset="us-ascii"Content-Disposition: attachment; filename="fxmsplod.tex"\documentstyle[12pt]{article}\oddsidemargin 0in  %0.5in\topmargin     0in\textheight    9in\textwidth     6.5in %6in\headheight    0in\headsep       0in\footskip      0.5in\footheight    0in\renewcommand{\baselinestretch}{1.5}\newcommand{\qed}{\ \rule[-0.2ex]{0.3em}{1.5ex}}\newcommand{\bfm}[1]{\mbox{\boldmath ${#1}$}}\newcommand{\RR}{{\rm R}\!\!\!\!\!{\rm I}\,\,}\newcommand{\bgamma}{\mbox{\boldmath$\gamma$} }%\newcommand{\bGamma}{\mbox{\boldmath$\Gamma$} }\newcommand{\boeta}{\mbox{\boldmath$\eta$} }\newcommand{\CC}{C\!\!\!\!I}\newcommand{\bzeta}{\mbox{\boldmath$\zeta$} }\begin{document}\begin{center}{\Large \bf Faster Computing $x^m \bmod p(x)$ and an Application \\to Splitting  a Polynomial Into Factors Over a  Fixed Disc}\vskip2mm%\centerline{Victor Y. Pan\footnote{Supported by NSF Grant CCR 9020690and PSC CUNY Awards Nos. 664334 and 665301.}\end{center}\centerline{Mathematics and Computer Science Department}\centerline{Lehman College, CUNY}\centerline{Bronx, NY 10468}\vskip2mm\centerline{Internet: VPAN@LCVAX.LEHMAN.CUNY.EDU}\centerline{Supported by NSF Grant CCR 9020690}\centerline{and PSC-CUNY Award 664334 and 665301}\abstract{Koenig's theorem is a well-known basis for fastsplitting a polynomial into factors over a fixed discin the complex plane.  We simplify the computation of such factorsby means of its reduction to solving a banded triangular Toeplitz linearsystem of equations.  The technique used may be of some interestin its own right, in particular, due to its possible extension tocomputing a power modulo a polynomial. }%\footnote{Supported by NSF Grant CCR 9020690 and PSC CUNY Awards Nos. 664334%and 665301.}\vskip1mm\noindent{\bf Key Words:}  complex polynomial factorization, algorithms,computation of a power modulo a polynomial.\vskip1mm\noindent{\bf 1991 Mathematics Subject Classification:} 68Q40, 68Q25, 30C15\vskip1.5mm\noindentConsider the approach to splitting a given polynomial $p(x)$ into 2factors over a complex disc  based on Koenig's theorem. Practicalimplementation of this approach generallyrequires to overcome some problems of numericalstability (see our remark 1). On the other hand, technically, this approachboils down to computing $x^m \bmod p(x)$ for a large integer $m$ (seetheorem 1 and remark 2). We propose a new solution to the latter problem,which is a problem of substantial independent interest(see remark 2).Let$$p (x) = \Pi_{i=1}^n ( x-z_i ) = \sum_{i=0} ^n p_i x^i ,\eqno (1)$$where$$| z_i | \leq r , i=1, \ldots , k ,\eqno (2)$$$$| z_i | \geq R , i = k+1 , \ldots , n ,\eqno (3)$$$k, n , r , R , p_0 , \ldots , p_n $ are given, $ R > r , $ and $ z_1 , \ldots ,z_n $ are unknown. The 2 polynomials,$$F (x) = \Pi_{i=1}^k (x-z_i ) , G (x) = p (x) / F (x),\eqno (4)$$define the  {\it splitting of $ p (x) $ over the discs} $D(0,r)$ or $ D(0,R)$.Computation of such a splitting is a major part of someeffective algorithms for polynomialfactorization in the complex domain [BP96], [H70], [Kir],[NR], [P94], [P95], [P95,a], [P96], [Sc].A well-known approach to its solution is based onthe following theorem, due to Koenig. We also refer the reader to thealternative recursive splitting algorithms of [Sc], [Kir],and of Appendices A and B of [P95,a], which use self-correctionafter each recursive stage.\vskip2mm\noindent{\bf Theorem 1} [H70].  {\it Under the assumptions (1)-(4), let$ p_0 = p (0) \neq 0 , $ let$$q (x) = \sum_{i=0}^{\infty} q_i x^i = 1 / (x^n p (1/x) )\eqno (5)$$denote the reciprocal of the reverse polynomial $ x^n p (1/x)  =\sum_{i=0}^n p_{n-i} x^i , $ and let $H_{m,k} $ denote the$ (k+1) \times (k+1) $ Hankel matrix with its $(i,j)$th entry equal to$q_{i+j + m - k +1 } , i,j=0,1, \ldots , k .$  Then the scaled lastcolumn of the inverse matrix $ H_{m,k}^{-1} $ approximates allthe components of the coefficient vector of the polynomial$ F (x) $ of (4) (of degree $ k < n ) $within the error norm $ O ( ( r/R ) ^m ) .$}Due to theorem 1, one may approximate the factor $ F(x) $ of $ p (x) $ in 2 stages as follows:\begin{description}\item{a)}choose a sufficiently large $m$ and compute the coefficients$ q_{m-k+1} , \ldots , q_{m+k+1} $ of the power series$ q(x)$ of (5);\item{b)}compute the last column of the matrix $ H_{m,k} ^{-1} .$\end{description}Stage (b) amounts to solving a Hankel linear systemof $ k+1 $ equations and can be performed in $ O ( k \log ^2 k )) $arithmetic operations by meansof various well-known algorithms [BP94].  Once $F(x)$ isoutput in stage (b), the computation of $ G (x) $ canbe performed by means of numerical division of $ p(x) $by a polynomial that approximates $ F (x);$  the cost ofthis computation is $ O ( n \log (n-k) $arithmetic operations (see [BP86], [BP96], or [PLS]).Practical implementation of stage (b) may become a hard problem since thematrix $H_{m,k}$, converging to a rank 1 matrix as $m \longrightarrow \infty$,becomes increasingly ill-conditioned as $m$ growslarge (cf. remark 1).Stage (a) can be performed by means of$ O ( \log m ) $ repeated squarings of $ x $ reducedmodulo $ p (x) $ and followedeither by thesubsequent solution of a general Toeplitz linear system [B] or byapplication of Euclidean algorithm [BG].  Our new alternativealgorithm is simpler, for it removes the latter stage, of solving thelinear system or performing Euclidean algorithm.Specifically, we reduce the entire stage (a) to solvingthe following banded and lowertriangular Toeplitz linear system of $m+k+2$ equations:$ T \vec{v}_m = \vec{e}_m^{(1)} , $ where $ \vec{e}_m ^{(1)} =(1,0, \ldots , 0 ) ^T , T = ( t_{i,j} ) , t_{i,j} = 0 $ if$ i < j $ or $ i-j>n ; t _{i,j} = p_{n+ j -i} $ otherwise.It is easy to verify that the desired coefficients $q_{m-k+1} , \ldots ,q_{m+k+1} $ of the power series $ q (x) $ of (5) are given by thelast $2k+1$ components of the vector $ \vec{v}_m .$  With no loss ofgenerality, we may assume that $(m+k+2) / (n+1) = 2 ^s - 1 $ for an integer$ s . $  (Otherwise, we may embed $T$ into an $ M \times M $ bandedtriangular Toeplitz matrix, for $ M = ( 2 ^s - 1 ) (n+1), M < 2 m .) $Then we may represent $T$ as a block bidiagonal block Toeplitzmatrix with $ k \times k $ nonzero blocks $ L $ (diagonal) and$ U $ (subdiagonal), $ U^T $ and $ L $ being lower triangular Toeplitzmatrices (cf. figure 1).  It is easily verified (for instance, by means ofapplyingblock Gaussian elimination) that the first block column of $ T ^{-1} $ isgiven by the block vector $ ( R _0 , R_1 , \ldots , R_{2^s - 1}) , T $ where$$R _0 = L ^{-1} , R _{i+1} = R_i U L ^{-1} , i = 0,1, \ldots ,2 ^s -2 .\eqno (6)$$Since $ T $ is a triangular Toeplitz matrix, so is $ T ^{-1} .$Therefore, all its blocks $ R_i $ are Toeplitz matrices, and we only need tocompute $ R _{2^s -1 } $ in order to obtain the last $ 2 k + 1\leq 2n+1 $ components of$ \vec{v}_m .$  From (6), we deduce that$$R _{2^{h+1}-1} = R _{2^h - 1} U R _{2^h - 1 } , h = 0 , 1 , \ldots , s-1 .\eqno (7)$$This enables us to reduce the computation of $ R _{2^s-1} $ to $ s = \log(m/ ( n+1) )$ multiplications of triplesof $ (n+1) \times ( n+1 ) $ triangular Toeplitz matrices.Since each resulting product of (7) is a Toeplitz matrix, we may compute itin $ O ( n \log n ) $ arithmetic operations, so that $ O ( ( n \log n )\log (m/n) ) $ operations suffice in order to perform the entire  stage (a).\vskip2mm\noindent{\bf Remark 1.}  The recursive process (7), like ones of [B] and of thealgorithm of [BG], generally magnifies the numericalerrors, so that special techniques of computer algebra [such as shiftingto exact rational arithmetic and application of modular (residue)reduction] should be used at stage (a).The same comments, and even more emphatically, apply to stage (b), whichtends to be hard for numerical implementation, due to the numericalcondition and numerical stability problems already cited.\vskip2mm\noindent{\bf Remark 2.}  Stage (a) is closed related to the problemof computing the coefficients of the polynomial$r(x)=  x ^m \bmod p(x)= x^m -q(x)p(x)$ (forsome polynomial $q(x)$ and for $m$ much larger thanthe degree of $ p(x)$). To see the correlation ofthese 2 problems to each other, substitute $x=1/y$ intothe latter polynomial equations, multiply it by$y^m$, and obtain that$$y^{m-n+1}R(y)=1-Q(y)P(y), \eqno (8)$$where $R(y)=y^{n-1} r(1/y)$, $Q(y)=y^{m-n}q(1/y)$, and$P(y)=y^np(1/y)$ are the reverse polynomials of $r(x)$, $q(x)$,and $p(x)$, respectively. Let $S(y) = \sum_{i=0}^{\infty} s_iy^i$denote the reciprocal power series for $P(y)$, satisfying $S(y)P(y)=1$.Substitute the latter equation into (8) and obtain that$$y^{m-n+1} R(y) = T(y) P(y),  \eqno (9)$$where $T(y)=S(y) - Q(y) = \sum_{i=m-n+1}^{\infty} s_iy^i$.Due to (9), the computation of $R(y)$, and, therefore, of $r(x)$,is essentially reduced  to the computation of the several leadingcoefficients of $T(y)$, which can be reduced to our computation atstage (a). Evaluation of $x^m \bmod p(x)$  is important in thestudy of digital transmission, queuing theory, and Markov chaincomputations [T], so that our techniques may be of some interest in theseareas too.\newpage\begin{figure}\begin{center}\vspace{.25in}$$T= \left( \begin{array}{ccccc}\huge  L    &    O    &      & \cdots     &   O    \\  U    &    L    & \ddots    &     &       \\  O    &    U    &  L   & \ddots    & \vdots      \\\vdots &  \ddots &  U   &  L  &     O  \\  O    &  \cdots       &  O   &  U  &  L\end{array} \right)$$\vskip2mm$U=\left( \begin{array}{lccr} u_0    &  u_1      & \cdots    &  u_k   \\  0     & \ddots    & \ddots    & \vdots \\\vdots  & \ddots    & \ddots    &  u_1   \\   0    & \cdots    &      0    &  u_0          \end{array} \right)$, \hspace{.5in}$L=\left( \begin{array}{lccr} l_0     &  0        & \cdots    &  0      \\ l_1     & \ddots    & \ddots    & \vdots  \\\vdots   & \ddots    & \ddots    &  0      \\ l_k     & \cdots    &    l_1    &  l_0          \end{array} \right)$\vskip3mm\caption{ Representation of abanded triangular Toeplitz matrix $T$}as a $5 \times 5$ block bidiagonal matrix.\end{center}\vskip5mm\end{figure}\newpage\vskip2mm\centerline{REFERENCES}\begin{description}\item{[B]}D. Bini, Complexity of Parallel Polynomial Computations, in {\it Proc.Parallel Computing: Methods, Algorithms, Applications} (J. Evans, C.Nodari editors), 115-126, Adam Hilger, Bristol, 1989.\item{[BG]}D. Bini, L. Gemignani, On the Complexity of Polynomial Zeros, {\it SIAM J.on Sci. Stat. Computing,} 8, 21, 781-799, 1992.\item{[BP86]}D. Bini, V.Y. Pan, Polynomial Division and Its Computational Complexity,{\it J. of Complexity,} 2, 179-203, 1986.\item{[BP94]}D. Bini, V.Y. Pan, {\it Polynomial and Matrix Computations, vol. 1: FundamentalAlgorithms}, Birkhauser, Boston, 1994.\item{[BP96]}D. Bini, V.Y. Pan, {\it  Polynomial and Matrix Computations, vol. 2: SelectedTopics},  Birkhauser, Boston, 1996, to appear.\item{[H70]}A.S. Householder, {\it The Numerical Treatment of a SingleNonlinear Equation,}  McGraw-Hill, New York, 1970.\item{[Kir]}P. Kirrinnis, Polynomial Factorization and Partial Fraction Decompositionby Simultaneous Newton's Iteration, Extended Abstract, 1994.\item{[NR]}C.A. Neff, J.H. Reif, An $ O ( n^{1+\epsilon} \log b ) $ Algorithm for theComplex Root Problem, {\it Proc. 35th Ann. IEEE Symp. on Foundations ofComputer Science,} IEEE Computer Society Press, 540-547, 1994.\item{[P94]}V.Y. Pan, New Resultant Inequalities and Complex Polynomial Factorization,{\it SIAM Journal on Computing}, 23, 5, 934-950, 1994.\item{[P95]}V.Y.Pan, Optimum (up to Polylog Factors) Sequential and ParallelAlgorithms for Approximating Polynomial Zeros, {\it 27th Ann. ACM Symp. onTheory of Computing}, 741--750, ACM Press, New York, 1995.\item{[P95,a]}V.Y. Pan, Deterministic Improvement of Complex Polynomial FactorizationBased on the Properties of the Associated Resultant, {\it Computersand Math. (with Applications}), 30, 2, 71--94, 1995.\item{[P96]}V.Y.Pan, Optimal and Nearly Optimal Algorithms for ApproximatingPolynomial Zeros, {\it Computers \& Math. (with Applications)}, 1996, to appear.\item{[PLS]}V.Y. Pan, E. Landowne, A Sadikou, Polynomial Division with a Remainder byMeans of Evaluation and Interpolation, {\it Information Process. Letters,}44, 149-153, 1992, and {\it Proc. 3rd IEEE Conf. on Parallel and DistributedProc.,} 212-218, IEEE Computer Society Press, 1991.\item{[Sc]}A. Sch\"onhage, The Fundamental Theorem of Algebra in Terms of ComputationalComplexity, manuscript, {\it Math. Dept., University ofT\"ubingen, T\"ubingen,} T\"ubingen, Germany, 1982.\item{[T]}W. Turin, {\it Performance Analysis of Digital Transmission Systems,} ComputerScience Press, 1990.\end{description}\end{document} --============_-1365856403==_============--