\documentclass{article}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\begin{document}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Zn}[1]{(\Z/{#1}\Z)^{*}}
\newcommand{\tf}{\tilde{f}}
\newcommand{\tg}{\tilde{g}}
\newcommand{\rev}{\textrm{rev}}

\title{Multiplication of reciprocal Laurent polynomials with a discrete 
       weighted transform}
\author{A. Kruppa}
\maketitle

\begin{abstract}
This note attempts to develop a method for multiplying two reciprocal
Laurent polynomials of degrees $\leq n - 2$ using two regular polynomial
multiplications degree $\leq n/2 - 1$ each, or with a weighted DFT of 
length $n$.
\end{abstract}

\section{Multiplying reciprocal Laurent polynomials}

Let $f(x)$ be a reciprocal Laurent polynomial
\begin{displaymath}
f(x) = f_0 + \sum_{i=1}^{d_f} f_i (x^i + x^{-i})
\end{displaymath}
of degree $2d_f$ and $\tf(x) = \tf_0 + \sum_{i=1}^{d_f} \tf_i x^i$ with
$\tf_0 = f_0 / 2, \tf_i = f_i \textrm{ for } i>0$ a polynomial so that
$f(x) = \tf(x) + \tf(1/x)$. Likewise for $g(x)$ and $\tg(x)$.
Their product $h(x) = f(x)g(x)$ is a reciprocal Laurent polynomial of degree
$2(d_f + d_g)$, $h(x) = h_0 + \sum_{i=1}^{d_f + d_g} h_i (x^i + x^{-i})$.

\subsection{Without DWT}

Let $\rev(\tf(x)) = x^{d_f} \tf(1/x)$ denote the polynomial with 
reversed sequence of coefficients. We have $\rev(\rev(\tf(x))) = \tf(x)$ and
$\rev(\tf(x) \tg(x)) = \rev(\tf(x)) \rev(\tg(x))$.
Let $\lfloor f(x) \rfloor$ denote a 
polynomial whose coefficients at non-negative exponents of $x$ are equal to 
those in $f(x)$, and whose coefficients at negative exponents of $x$ are 
$0$. We have 
$\lfloor f(x) + g(x) \rfloor = \lfloor f(x) \rfloor + \lfloor g(x) \rfloor$.

Now we can compute the product
\begin{eqnarray*}
   &   & f(x)g(x) \\
   & = & (\tf(x) + \tf(1/x)) (\tg(x) + \tg(1/x)) \\
   & = & \tf(x) \tg(x) + \tf(x) \tg(1/x) + \tf(1/x)\tg(x) + \tf(1/x)\tg(1/x) \\
%   & = & \tf(x) \tg(x) + x^{-d_g} \tf(x) \rev(\tg(x)) + x^{-d_f} \rev(\tf(x))\tg(x) + \tf(1/x)\tg(1/x) \\
   & = & \tf(x) \tg(x) + x^{-d_g} \tf(x) \rev(\tg(x)) + x^{-d_f} \rev(\tf(x) \rev(\tg(x))) + \tf(1/x)\tg(1/x)
\end{eqnarray*}
but we only want to store the coefficients at non-negative exponents in the 
product, so
\begin{eqnarray*}
   &   & \lfloor f(x)g(x) \rfloor \\
   & = & \tf(x) \tg(x) + \lfloor x^{-d_g} \tf(x) \rev(\tg(x)) \rfloor + \lfloor x^{-d_f} \rev(\tf(x) \rev(\tg(x))) \rfloor + \tf_0 \tg_0.
\end{eqnarray*}

We can compute this with two multiplications of a degree $d_f$ and a 
degree $d_g$ polynomial, as well as $d_f + d_g + 2$ additions and one doubling
of the $\tf_0 \tg_0$ term in $\tf(x) \tg(x)$. If an FFT based 
multiplication routine is used, the forward transform of $\tf$ can be re-used
directly, and the forward transform of $\tg$ can be reused by observing that
the $i$-th coefficient in the length $l$ DFT of $\tg(x)$ is 
$\tg(\omega^i)$, $\omega^l=1$, and in the DFT of $rev(\tg(x))$ is
$\omega^{id_g} \tg(\omega^{l-i})$, so the coefficients at indices $1 \ldots l-1$ 
only need to be reversed in order and suitably weighted. (I haven't acrtually 
tried this idea of re-using the forward transform)

\subsection{With DWT}

The product RLP $h(x)$ of degree $2d_h = 2(d_f+d_g)$ has $d_h + 1$ 
possibly distinct coefficients in standard basis which we can obtain from
a discrete weighted FFT convolution of length $l \geq d_h + 1$.

Suppose we compute $\hat{h}(x) = f(x) \cdot g(x) \bmod (x^l - 1/a)$, then
\begin{eqnarray*}
  \hat{h}(x) & = & h_0 + \sum_{i=1}^{d_h} h_i (x^i + ax^{l-i}) \\
             & = & \sum_{i=0}^{d_h} h_i x^i + \sum_{i=1}^{d_h} ah_i x^{l-i}\\
             & = & \sum_{i=0}^{d_h} h_i x^i + \sum_{i=l-d_h}^{l-1} ah_{l-i} x^i\\
             & = & \sum_{i=0}^{l-d_h-1} h_i x^i + 
                   \sum_{i=l-d_h}^{d_h} (h_i + ah_{l-i}) x^i +
                   \sum_{i=d_h+1}^{l-1} ah_{l-i} x^{i} 
\end{eqnarray*}

For $l-d_h \leq i < l/2$, $\hat{h}_i = h_i + ah_{l-i}$ and 
$\hat{h}_{l-i} = ah_i + h_{l-i}$. With $a=1$ as in an unweighted FFT,
$h_i$ and $h_{n-i}$ cannot be separated as $\hat{h}_i = \hat{h}_{l-i}$. 
With $a=-1$, they cannot be separated either, as 
$\hat{h}_i = - \hat{h}_{l-i}$. 
With, e.g., $a=\sqrt{-1}$, we have
\begin{eqnarray*}
a (h_i + ah_{l-i}) - (ah_i + h_{l-i}) & = & (a^2-1) h_{l-i} \\
 & = & -2 h_{l-i}
\end{eqnarray*}
so we can separate $h_i$ and $h_{l-i}$. Hence after processing
\begin{eqnarray*}
\hat{h}_{l-i} & := & -(a \hat{h}_i - \hat{h}_{l-i})/2 \\
\hat{h}_i & := & \hat{h}_i - a\hat{h}_{l-i}
\end{eqnarray*}
for $l-d_h \leq i < l/2$, we have the desired coefficients $h_i$ for 
$0 \leq i \leq d_h$ in $\hat{h}_i$.


% Vecdiv(a,b)={if(length(a)!=length(b),error("Vecdiv: vectors of unequal length"));vector(length(a),i,a[i]/b[i])}
% Vecmul(a,b)={if(length(a)!=length(b),error("Vecdiv: vectors of unequal length"));vector(length(a),i,a[i]*b[i])}
% Vecdiv(Vecrev(Polrev([f0,f1*w,f2*w^2,f3*w^3,0*w^4,f3*w^-3,f2*w^-2,f1*w^-1])^2 % (x^8-1)), [1,w^1,w^-6,w^-5,w^-4,w^-3,w^-2,w^-1])

% (Vecrev((f0 + f1*(x+1/x) + f2*(x^2+1/x^2) + f3*(x^3+1/x^3) + f4*(x^4+1/x^4))^2 % (x^10-w^10)) - Vecdiv(Vecrev(Polrev(Vecmul([f0,f1,f2,f3,f4,0,f4,f3,f2,f1],[1,w,w^2,w^3,w^4,w^5,w^-4,w^-3,w^-2,w^-1]))^2 % (x^10-1)), [1,w^1,w^2,w^3,w^4,w^5,w^6,w^7,w^8,w^9]))

% (Vecrev(f(x)^2 % (x^10-w^10)) - Vecdiv(Vecrev(Polrev(Vecrev(f(w*x) % (x^10-1)))^2 % (x^10-1)), [1,w^1,w^2,w^3,w^4,w^5,w^6,w^7,w^8,w^9]))

To obtain the desired product $\bmod{(x^l-1/a)}$, we can choose a  
constant $w^l=1/a$ and compute the product
$h(wx) \bmod{(x^l-1)} = 
(f(wx) \bmod{(x^l-1)} \cdot g(wx) \bmod{(x^l-1)}) \bmod {(x^l-1)}$. We have
\begin{eqnarray*}
 h(wx) \bmod{(x^l-1)} & = & \sum_{i=0}^{d_h} (w^i h_i x^i + w^{-i} h_i x^{l-i}) \\
 & = & \sum_{i=0}^{d_h} w^i h_i x^i + \sum_{i=0}^{d_h} w^{-i} h_i x^{l-i} \\
 & = & \sum_{i=0}^{d_h} w^i h_i x^i + \sum_{i=l-d_h}^l w^{i-l} h_{l-i} x^i \\
 & = & \sum_{i=0}^{l-d_h-1} w^i h_i x^i + \sum_{i=l-d_h}^l (w^i h_i + w^{i-l} h_{l-i}) x^i \\
 & = & \sum_{i=0}^{l-d_h-1} w^i h_i x^i + \sum_{i=l-d_h}^l w^i (h_i + a h_{l-i}) x^i
\end{eqnarray*}
so dividing the $i$-th coefficient by $w^i$ yields $\hat{h}_i$ as desired.

For example, with two degree $6$ reciprocal Laurent polynomials and a length 
$8$ convolution, the coefficient vectors of $f(wx) \bmod{(x^8-1)}$ and 
$g(wx) \bmod{(x^8-1)}$ 
are $(f_0, w f_1, w^2 f_2, w^3 f_3, 0, w^{-3} f_3, w^{-2} f_2, w^{-1} f_1)$ 
and \\ 
$(g_0, w g_1, w^2 g_2, w^3 g_3, 0, w^{-3} g_3, w^{-2} g_2, w^{-1} g_1)$, 
respectively.
Their product $\bmod{(x^8-1)}$ has coefficient vector
\begin{eqnarray*}
&& (f_0 g_0 + 2(f_1 g_1 + f_2 g_2 + f_3 g_3), \\
&& (f_1 g_0 + (f_0 + f_2) g_1 + (f_1 + f_3) g_2 + f_2 g_3) w, \\
&& ((f_2 g_0 + (f_1 + f_3) g_1 + f_0 g_2 + f_1 g_3) w^8 + f_3 g_3)/w^6, \\
&& ((f_3 g_0 + f_2 g_1 + f_1 g_2 + f_0 g_3) w^8 + f_3 g_2 + f_2 g_3)/w^5, \\
&& ((f_3 g_1 + f_2 g_2 + f_1 g_3) w^8 + f_3 g_1 + f_2 g_2 + f_1 g_3)/w^4, \\
 && ((f_3 g_2 + f_2 g_3) w^8 + f_3 g_0 + f_2 g_1 + f_1 g_2 + f_0 g_3)/w^3, \\
&& (f_3 g_3 w^8 + f_2 g_0 + (f_1 + f_3) g_1 + f_0 g_2 + f_1 g_3)/w^2, \\
&& (f_1 g_0 + (f_0 + f_2) g_1 + (f_1 + f_3) g_2 + f_2 g_3)/w).
\end{eqnarray*}
With 
\begin{eqnarray*}
&& h_0 = f_0 g_0 + 2 f_1 g_1 + 2 f_2 g_2 + 2 f_3 g_3 \\
&& h_1 = f_1 g_0 + (f_0 + f_2) g_1 + (f_1 + f_3) g_2 + f_2 g_3 \\
&& h_2 = f_2 g_0 + (f_1 + f_3) g_1 + f_0 g_2 + f_1 g_3 \\
&& h_3 = f_3 g_0 + f_2 g_1 + f_1 g_2 + f_0 g_3 \\
&& h_4 = f_3 g_1 + f_2 g_2 + f_1 g_3 \\
&& h_5 = f_3 g_2 + f_2 g_3 \\
&& h_6 = f_3 g_3 \\
\end{eqnarray*}
this vector is equal to
$(h_0, h_1 w, h_2 w^2 + h_6 w^{-6}, h_3 w^3 + h_5 w^{-5}, h_4 w^4 + h_4 w^{-4},
  h_5 w^5 + h_3 w^{-3}, h_6 w^6 + h_2 w^{-2}, h_1 w^{-1})$ 
and after dividing the $i$-th coefficient by $w^i$ is equal to
$(h_0, h_1, h_2 + h_6 w^{-8}, h_3 + h_5 w^{-8}, h_4 + h_4 w^{-8},
  h_5 + h_3 w^{-8}, h_6 + h_2 w^{-8}, h_1 w^{-8}) = (h_0, h_1, h_2 + a h_6, 
h_3 + a h_5, h_4 + a h_4, h_5 + a h_3, h_6 + a h_2, a h_1)$.
For $0 \leq i < l - d_h = 2$, the coefficients $h_i$ can be read directly. 
For $l - d_h \leq i \leq l/2$, coefficients $i$ and $l-i$ overlap and must
be separated as shown above.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{With DWT (rewrite)}

Let $Q(X) = q_0 + \sum_{j=1}^{d_q} q_j (X^j + X^{-j})$ be an RLP of degree
$2d_q \leq 2l-2$ and likewise $R(X)$ an RLP of degreee $2d_r \leq 2l-2$.
To obtain the product RLP $S(X) = Q(x) R(x)$ of degree $2d_s = 2(d_q+d_r)$,
we can perform a weighted
convolution product by computing $\tilde{S}(wx) = Q(wX) R(wX) \bmod{(X^l-1)}$,
$w^l = \sqrt{-1}$, and separating the wrapped-around coefficients in 
$\tilde{S}(X)$ to obtain the coefficients of $S(X)$. Since $w$ is a $4l$-th
root of unity we need to choose NTT primes $p_j \equiv 1 \pmod{4l}$.

With $\tilde{S}(wX) = \sum_{i=0}^{l-1} \tilde{s}_j x^j = S(wX) \bmod (X^l-1)$ 
we have
\begin{eqnarray*}
  \tilde{S}(wX) 
  & = & \sum_{j=0}^{d_s} (w^j s_j x^j + w^{-j} s_j x^{l-i}) \\
  & = & \sum_{j=0}^{d_s} w^j s_j x^j + \sum_{j=0}^{d_s} w^{-j} s_j x^{l-i} \\
  & = & \sum_{j=0}^{d_s} w^j s_j x^j + \sum_{j=l-d_s}^l w^{i-l} s_{l-i} x^i \\
  & = & \sum_{j=0}^{l - d_s - 1} w^j s_j x^j + 
        \sum_{j=l-d_s}^l w^j (s_j + w^l s_{l-j}) x^j
\end{eqnarray*}
Hence for $0 \leq i < l - d_s$ we can read $s_i$ directly off $\tilde{s}_i$ 
while the remaining coefficients first need to be separated by using
$w^l \tilde{s}_j - \tilde{s}_{l-j} = w^l (s_j + w^l s_{l-j}) - 
(w^l s_j + s_{l-j}) =  (w^{2l}-1) s_{l-j} = -2 s_{l-j}.$

These ideas yield the algorithm shown in figure~\ref{DWTNTT}. Since we need 
only squaring of RLPs in section \ref{}, we show only the squaring algorithm
here.

\begin{figure}[ht]
\begin{center}
\begin{tabular}{l}
\hline  
  Input:  RLP $Q(X) = \sum_{j=0}^{d_q} q_j (x^j + x^{-j})$ of degree $2d_q$ 
    in standard basis \\
  Output: RLP $S(X) = \sum_{j=0}^{d_s} s_j (x^j + x^{-j}) = Q(X)^2$ of degree 
    $2d_s = 4d_q$ \\ in standard basis \\
  Auxiliary storage: A length $l > 2d_s$ NTT array $M$ with \\
  separate memory for vectors mod each $p_j$\\
\hline
  For each prime $p_j$ \\
\qquad  Compute $w$ with $w^{2l} \equiv -1 \pmod{p_j}$  \\
\qquad  For $0 \leq i \leq d_q$ store $w^i q_i \pmod{p_j}$ in $M_i$ \\
\qquad  For $1 \leq i \leq d_q$ store $w^{-i} q_i \pmod{p_j}$ in $M_{l-i}$ \\
\qquad Perform forward NTT modulo $p_j$ on $M$, square elementwise and \\
\qquad \quad perform inverse 
NTT \\
\qquad For $1 \leq i \leq d_s$ set $M_i := w^{-i} M_i \pmod{p_j}$ \\
\qquad For $l - d_s < i \leq l/2$ \\
\qquad \qquad Set $M_{l-i} := -(w^l M_i - M_{l-i})/2 \pmod{p_j}$ \\
\qquad \qquad If $i < l/2$ set $M_i := M_i - w^l M_{l-i} \pmod{p_j}$ \\
 For $0 \leq i \leq d_s$ perform CRT on $M_i$ residues to obtain $s_i$, store
  in output \\
\hline
\end{tabular}
\end{center}
\caption{NTT based squaring algorithm for Reciprocal Laurent polynomials}
\label{DWTNTT}
\end{figure}

\end{document}
