\documentstyle [11pt] {article}\topmargin 	0in\textheight 9in\textwidth 	6in\headheight 0in\headsep				0in\footskip			0.5in\footheight	0in\newtheorem {algorithm}{\sc \bf Algorithm}[section]\newtheorem {problem}{\sc \bf Problem}[section]\newtheorem {theorem}{\sc \bf Theorem}[section]\newtheorem {proposition}{\sc \bf Proposition}[section]\newtheorem {lemma}{\sc \bf Lemma}[section]\newtheorem {corollary}{\sc \bf Corollary}[section]\newtheorem {fact}{\sc \bf Fact}[section]\newtheorem {remark}{\sc \bf Remark}[section]\newtheorem {definition}{{\sc \bf Definition}}[section]\newtheorem {claim}{\sc \bf Claim}[section]\def\theequation{\arabic{section}.\arabic{equation}}\title{\bf New fast algorithms for polynomial interpolation and evaluation on the Chebyshev node set.}\author{ Victor Y. Pan\\Mathematics and Computer\\ Science Department\\Lehman College, CUNY\\ Bronx, NY 10468\\internet: vpan@lcvax.lehman.cuny.edu\\(Supported by NSF Grant CCR 9625344\\and PSC-CUNY Award No. 667340)}\date{}\begin{document}\maketitle\begin{abstract}For a polynomial $p(x)$ of a degree $n$, we study itsinterpolation and evaluation on a set of Chebyshev nodes, $x_k = \cos (\frac{(2k+1)\pi}{2n +2})$, $k=0,1, \ldots, n$. This is easily reduced to applying discrete Fouriertransforms (DFTs) to the auxiliary polynomial $q(w) =w^np(x)$, where $2 x = \alpha w +(\alpha w)^{-1}$,$\alpha = \exp (\pi \sqrt{-1} /(2n))$. We show the back and forth transition between $p(x)$ and $q(w)$ based on the respective back and forth transformations of the variable:$\alpha w =\frac{1-z}{1+z}$, $y=\frac{x-1}{x+1}$, $y=z^2$. All these transformations (like the DFTs) are performed byusing $O(n \log n)$ arithmetic operations, which thussuffice in order to support both interpolation andevaluation of $p(x)$ on the Chebychev set, as well as on some related sets of nodes. This improves, by factor $\log n$, the known arithmetic time bound for Chebyshev interpolation and gives an alternative supporting algorihtm for the record estimate of $O(n \log n)$ for Chebyshev evaluation, obtained by Gerasoulis in 1987 and based on a distinct algorithm. \end{abstract}\vspace{.2in}\noindent{\bf Key words:} Chebyshev nodes, polynomial interpolation,polynomial evaluation, algorithms, computational complexity. \vspace{.2in}\noindent{\bf 1991 Math. Subject Classification:} 41A05, 41A10, 68Q25,68Q40. \section{Introduction.} Interpolation to a polynomial $p(x)$ of a degree at most $n$ on theset of Chebyshev nodes,\begin{equation}\label{eq1.1} x_k = \cos (\frac{2k+1}{2n +2}\pi ), \  k=0,1, \ldots, n, \end{equation}plays a major role in the study of the approximation of realfunctions (cf. e. g. \cite{cdb}, pages 242--244, or\cite{db}, pages 107--108). In this paper, we propose an algorithm that supports thefollowing result: \begin{theorem}\label{th1.1} Interpolation to a polynomial of a degree at most $n$ on the Chebyshev node set (\ref{eq1.1}) can be performed by using $O(n \log n)$ arithmetic operations. \end{theorem} Theorem~\ref{th1.1} improves, by factor $\log n$, theprevious record complexity estimate for this classicalproblem (cf. \cite{p89}). Our method is immediately extendedto yield an algorithm that supports thefollowing result: \begin{theorem}\label{th1.2}  A polynomial of a degree at most $n$ can be evaluated on theChebyshev set of nodes (\ref{eq1.1}) at the cost of $O(n \log n)$ arithmetic operations. \end{theorem}Theorem~\ref{th1.2} has been proved in \cite{g}, based on adistinct algorithm. Furthermore, theorems \ref{th1.1} and \ref{th1.2} can be supported by the algorithms of \cite{p89}, which rely on the following equations (and their recursive extension):$$ p(x)=p_1(y) + xp_2(y), $$where $$ x=x_k, \  y=y_k=1-2x_k^2=\cos\left( \frac{2k +1} {n+1}\pi \right), \ k=0,1, \ldots,n, $$ $$ 2\deg p_i(y) \leq \deg p(x), \  i=1,2.$$ It remains to observe that $x_k=-x_{n-k}$ for all $k$ and then to repeat literally the derivation of the forward and inverse FFTs as presented onpages 9--12 of \cite{bp} but with $x_k$ playing the role of the roots of 1. Our present derivation of theorems~\ref{th1.1} and~\ref{th1.2} will reveal their correlation to FFTs more directly, based on back and forth transformations between a fixed circle $C$ on the complex plane and the real coordinate line $R$. Such transformations  can be also applied to the extension of an effective algorithm for the evaluation or approximation of polynomials and rational functions on the circle $C$ to one for $R$ and vice versa.Both of our theorems are immediately extended to the evaluation andinterpolation on the expanded Chebyshev set, \begin{equation}\label{eq1.2} y_k = ax_k + b, \   \  k=0,1, \ldots, n, \end{equation}for any pair of real (or even complex) $a \neq 0$ and $b$(cf. \cite{cdb}, page 244), as well as on the{\em Chebyshev-like} set, \begin{equation}\label{eq1.3} x_k^*=\cos\frac{\pi k}{n+1}, \  \    k=0, 1, \ldots, n, \end{equation} and on the {\em expanded Chebyshev-like} set, \begin{equation}\label{eq1.4} y_k^* = ax_k^* + b, \   \  k=0,1, \ldots, n. \end{equation}Hereafter, {\em " ops"} stands for "arithmetic operations". We organize our paper as follows. In the next section, we recall the known correlationsbetween the Chebyshev set and the roots of 1. In section~3,we reduce the evaluation and interpolation problems on theChebyshev set for a given polynomial $p(x)$ to discreteFourier transforms (DFTs) applied to an auxiliarypolynomial $q(w)$. In section~5, we showtransformations of $p(x)$ into $q(w)$ and viceversa, based on some elementary transformations of the variable, for which, in section~4, we give algorithms using $O(n\log n)$ ops.  Since thesame asymptotic cost bound is known for DFT, we arrive at theorems \ref{th1.1} and \ref{th1.2}.  \section{Correlations between the Chebyshev nodesand the roots of 1.} Recall Euler's formula: $$ 2 \cos \phi = \exp ( \phi \sqrt{-1}) + \exp (-\phi \sqrt{-1} ). $$ Substitute \begin{equation}\label{eq2.1} \phi =\frac{2 \pi k}{2n +2}, \  \  w_k = \gamma^k, \  \  \gamma =\exp(\frac{2\pi\sqrt{-1}}{2n+2}), \end{equation} $\gamma =\gamma_{2n +2}$ being a primitive $(2n +2)$-nd root of $1$, and arrive at the following expressions for the Chebyshev-like node set of (\ref{eq1.3}): \begin{equation}\label{eq2.2} 2x^*_k = 2 \cos \frac{2\pi k}{2n +2} = w_k + w_{-k} = \gamma^k + \gamma^{-k}, \  \   k=0, 1, \ldots, n.  \end{equation}Similarly, substitute $\phi = \frac{(2k+1)\pi}{2n +2}$ and arrive at the following expressions for the Chebyshev node set of (\ref{eq1.1}): \begin{equation}\label{eq2.3}2x_k =2\cos \frac{(2k +1) \pi}{2n +2}=\alpha w_k +\frac{w_{-k}}{\alpha} = \alpha \gamma^k_{2n +2} +\frac{1}{\alpha \gamma^k_{2n +2}}, \  \    k=0,1, \ldots, n, \end{equation}where \begin{equation}\label{eq2.4} \alpha =  \exp (\pi \sqrt{-1}/(2n+2)) = \sqrt{\gamma}.\end{equation} \section{Reduction of the polynomial evaluation andinterpolation on Chebyshev node set to discrete Fourriertransforms.} We write \begin{equation}\label{eq3.1}p(x) = \sum_{i=0}^n p_i x^i, \end{equation} $$ p(x) = p(\frac{\alpha w + (\alpha w)^{-1}} {2}), $$ \begin{equation}\label{eq3.2} q(w)=w^n p(x)= w^n p(\frac{\alpha w + (\alpha w)^{-1}} {2}), \end{equation} \begin{equation}\label{eq3.3} x=\frac{\alpha w + (\alpha w)^{-1}} {2}\end{equation} for $\alpha$ of (\ref{eq2.4}), $\alpha = \sqrt{\gamma}$. Note that $q(w)$ is a polynomial in $w$ of a degree at most$2n$ and that $$ \gamma^{(2k+1)n} q(w_{-k-1})=q(w_k)= w_k^np(x_k) = \gamma^{kn}p(x_k), \  \ k=0, 1, \ldots, n. $$Therefore, the values of $p(x)$ on the Chebyshev node set of(\ref{eq1.1}), (\ref{eq2.3}) immediately define the values of $q(w)$on the set of all the $(2n +2)$-nd roots of 1 and vice versa. The transition from the values of $q(w)$ on the latter setto the coefficients of $q(w)$ or vice versa amounts toforward or inverse discrete Fourier transforms, respectively, and can beperformed by using $O(n\log n)$ ops (cf. e. g. \cite{bp},pages 9--15). We arrive at the next results: \begin{fact}\label{fa3.1} $$    $$     \begin{description}\item[a)] Given the values $p(x_k)$ of the polynomial $p(x)$of (\ref{eq3.1}) on the Chebyshev set $\{ x_{k+1}, \   k=0, 1,\ldots, n\}$ of (\ref{eq1.1}) and (\ref{eq2.3}), the coefficientsof the polynomial $q(w)$ of (\ref{eq3.2}) can be computed in $O(n \log n)$ ops; \item[b)] vice versa, given the coefficients of $q(w)$, thevalues $p(x_k)$, $k=0, 1, \ldots, n$, can be computed in $O(n\log n)$ ops. \end{description} \end{fact} Due to fact \ref{fa3.1}, we will arrive at theorems\ref{th1.1} and \ref{th1.2} as soon as we show how to maketransition from the coefficients of $p(x)$ to ones of $q(w)$and vice versa, by using $O(n \log n)$ ops. We will show thisin section 5, based on the auxiliary transformations to bestudied in the next section. \section{Transformations of variables and polynomials.} In this section, we will show that $O(n \log n)$ ops sufficein order to compute the transformation of a polynomialimplied by a linear rational transformation of its variable. Clearly, the scaling of the variable, that is, the transitionfrom the coefficients of the polynomial \begin{equation}\label{eq4.1} t(x) = \sum_{i=0}^n t_i x^i \end{equation} to ones of $t(ax) =  \sum_{i=0}^n (t_i a^i) x^i, \  \   a \neq 0$, only requires $2n$ multiplications. Thefollowing extension (including the shift ofthe variable) is well-known (cf. e. g.\cite{bp}, pages 15--16): \begin{fact}\label{fa4.1} Given the coefficients $t_0, t_1, \ldots, t_n$of the polynomial $t(x)$ of (\ref{eq4.1}) andcomplex $a \neq 0$ and $b$, the coefficients $u_0, u_1, \ldots, u_n$ of the polynomial $t(ax+b) =u(x) = \sum_{i=0}^n u_ix^i$ can becomputed by using $O(n \log n)$ ops. \end{fact} Let us extend this result to a linear rational transformation of $x$. \begin{corollary}\label{co4.2} Given the coefficients $t_0, t_1, \ldots, t_n,\ a, \ b,\  c,\  d$ of a polynomial $t(x)$ of(\ref{eq4.1}) and of a linear fraction$\frac{ax +b}{cx+d}$, where $ad \neq bc$, thecoefficients $v_0, v_1, \ldots, v_n; \  v^*_0, v^*_1, \ldots, v^*_n$ of the rational function $ t(\frac{ax+b}{cx+d})= \frac{v(x)}{v^*(x)}$, where $v(x) = \sum_{i=0}^n v_ix^i$, $v^*(x)= (cx+d)^n=\sum_{i=0}^n v^*_ix^i$, can becomputed by using $O(n \log n)$ ops. \end{corollary} {\bf Proof.} Let $c \neq 0$ (otherwise, $d\neq 0$, $v(x)= d^n t((a/d)x + b/d)$, and it suffices to applyfact~\ref{fa4.1}). Then $\frac{ax+b}{cx+d} =\frac{e}{cx+d} + f$, where $f=\frac{a}{c}$, $e =b-\frac{ad}{c} \neq 0$. Let $y = \frac{1}{cx+d}$. Due tofact~\ref{fa4.1}, $O(n \log n )$ ops suffice in order to compute the coefficients of the polynomial$t(f+ey)=u(y)=\sum_{i=0}^nu_iy^i$. Substitute $y^*=cx+d$into the reverse polynomial$u_r(y^*)=\sum_{i=0}^nu_{n-i} y^{*i}= y^{*n} u(1/y^*)$ andcompute the coefficients $v_0, v_1, \ldots, v_n$ of thepolynomial $v(x)=u_r(cx+d)=\sum_{i=0}^nv_ix^i$. Observethat $$ v(x) = u_r(cx+d) = (cx+d)^n u(\frac{1}{cx+d}) =(cx+d)^n t(f+\frac{e}{cx+d}) = (cx+d)^n t(\frac{ax+b}{cx+d}).$$Therefore, $t(\frac{ax+b}{cx+d}) = v(x)/(cx+d)^n$. We already have the coefficients of $v(x)$, and it remains to compute the coefficients of $(cx+d)^n = c^n \sum_{i=0}^n \left( \begin{array}{c} n \\ i \end{array} \right) \left( \frac{d}{c} \right)^{n-1}x^i$, which takes $O(n)$ops. \hfill $\Box$ \section{Back and forth transformations between $q(w)$ and$p(x)$.} Let us show the desired transformations of the polynomials$q(w)$ and $p(x)$ into each other performed by means oflinear rational transformations and squaring the variable. We write $y=z^2$, slightly abusing the notation, \begin {equation}\label{eq5.1} \alpha w = \frac{1-z}{1+z}, \  \  y=\frac{x-1}{x+1} \end{equation} and obtain that \begin{equation}\label{eq5.2} \alpha w + (\alpha w)^{-1}= \frac{1}{2}(\frac{1-z}{1+z} +\frac{1+z}{1-z})= \frac{1+z^2}{1-z^2} = \frac{1+y}{1-y}=x. \end{equation} Now, we recall from (\ref{eq3.2}) that $w^{-n}q(w)=p(x)$ and write \begin{equation}\label{eq5.3} q(w)/w^n =\alpha^{-n} (\frac{1+z}{1-z})^n q(\frac{1-z}{1+z} \cdot\frac{1}{\alpha})=  r(z)=s(z^2)=s(y)=s(\frac{x-1}{x+1}). \end{equation} Then (\ref{eq3.2}) implies that \begin{equation}\label{eq5.4} p(x)=s(\frac{x-1}{x+1}). \end{equation}  The successive transformations from the coefficients of rational functions $q(w)/w^n$ and $s(y)$ to the coefficients of $r(z)$ and $p(x)$,respectively, can be performed in $O(n \log n)$ ops, due tocorollary~\ref{co4.2}.  To complete the proof of theorem~\ref{th1.1}, it remains to show that $r(z)$ is actually a rational function in $z^2$. This, however,immediately follows from (\ref{eq5.2}) and from the equation $p(x) = r(z)$, implied by (\ref{eq3.2}) and (\ref{eq5.3}).  \hfill $\Box$  To prove theorem \ref{th1.2}, we reverse the transformations of (\ref{eq5.3}) and (\ref{eq5.4}), from $q(w)/w^n$ to $p(x)$, so as to proceed from $p(x)$ to $q(w)/w^n$. Theorem~\ref{th1.2} follows, bythe same argument, since the reverse of a linear rational transformation$\frac{ax+b}{cx+d}$ is again linear rational, provided that$ad \neq cd$; in particular, the transformations $z= \frac{1-\alpha w}{1+\alpha w}$ and $x= \frac{1+y}{1-y}$ reverse ones of (\ref{eq5.1}). \hfill  $\Box$  The extension to the case of using Chebyshev-like set (\ref{eq1.3}),(\ref{eq2.2}) is immediate if we replace $\alpha$ by 1 in the proofs oftheorems \ref{th1.1} and \ref{th1.2}. Further extensions to the expandedChebyshev and Chebyshev-like sets of (\ref{eq1.2}),(\ref{eq1.4}) immediately follow from fact~\ref{fa4.1}.\begin{thebibliography}{92} \bibitem[BP]{bp} D.~Bini, V.~Y.~Pan, {\em Polynomial and Matrices Computations, vol. 1: Fundamental Algorithms}, Birkh\"{a}user, Boston, 1994.\bibitem[CdB]{cdb} C. D. Conte, C. de~Boor, {\em Elementary NumericalAnalysis: an Algorithmic Approach}, McGraw-Hill, New York, 1980.\bibitem[DB]{db} G. Dahlquist, \AA. Bj\"{o}rk, {\em Numerical Methods}, Prentice-Hall, New Jersey, 1974.  \bibitem[G]{g} A. Gerasoulis, A Fast Algorithm for the Multiplication of Generalized Hilbert Matrices with Vectors, {\em Math. Comp.}, 50, 181, 179--188, 1987.\bibitem[P89]{p89} V.~Y.~Pan, Fast and Efficient Parallel Evaluation of the Zeros of a Polynomial Having Only Real Zeros, {\em Computers and Math (withApplications)}, 17 (1989) 11, pp.~1475--1480. \end{thebibliography} \end{document}  