\documentclass{article}
\usepackage{amsmath}
\begin{document}

\title{Polynomial multipoint evaluation along a geometric progression}
\author{A. Kruppa}
\maketitle

Let
\begin{displaymath}
F(x) = \sum_{j=0}^d f_j x^j
\end{displaymath}

To evaluate $F$ on $n$ points in geometric progression, i.e. 
$F(1), F(q), F(q^2), \ldots ,F(q^{n-1})$, we can do a convolution of the 
polynomials
\begin{eqnarray*}
B(x) & = & \sum_{j=0}^{l - 1} b_j x^j \\
     &   & \textrm{with } b_j = q^{b(j)}\\
C(x) & = & \sum_{j=0}^d c_j x^j \\
     &   & \textrm{with } c_j = f_j q^{c(j)}
\end{eqnarray*}
where $l = n + d$ and get the product
\begin{displaymath}
P(x) = B(x)C(x) = \sum_{i=0}^{l+d-1} p_i x^i
\end{displaymath}
with
\begin{eqnarray*}
p_i & = & \sum_{\substack{0\leq j\leq d,\\0\leq i-j < l}} b_{i-j} c_{j} \\
    & = & \sum_{\substack{0\leq j\leq d,\\0\leq i-j < l}} f_j q^{b(i-j)} q^{c(j)}.
\end{eqnarray*}
If $d \leq i < l$ and $0\leq j\leq d$, $0 \leq i-j < l$, 
so the condition on $j$ suffices for this range of $i$:
\begin{eqnarray*}
p_i & = & \sum_{0\leq j\leq d} f_j q^{b(i-j) + c(j)}, \textrm{ for } d \leq i < l.
\end{eqnarray*}

Since we only want the coefficients $p_i$ for $d\leq i < l$, it is
permissible to compute $B(x)C(x) \% (x^l-1)$, i.e. by a length $l$ cyclic 
convolution product. This way, the coefficients
$p_i$, $l \leq i < l + d$ will overlap with the coefficients $p_i$, 
$0 \leq i < d$, but we aren't interested in any of them.


\subsection{Simple progression}
Now we would like 
\begin{displaymath}
q^{h(i)} p_{i+d} =  \sum_{j=0}^{d} f_j q^{ji}
\end{displaymath}
so that $q^{h(i)} p_{i+d} = F(q^i)$, for $0 \leq i < n$.
%
Hence we can equate the exponents
\begin{displaymath}
  b(i+d-j) + c(j) + h(i) = ij.
\end{displaymath}
If $b(x)$ is a polynomial, it must be of degree at least $2$ to produce the 
$ij$ term. Let $b(x)=-(x-d)^2/2$. Then 
\begin{eqnarray*}
  -i^2/2 + ji - j^2/2 + c(j) + h(i) & = & ij \\
  -i^2/2 - j^2/2 + c(j) + h(i) & = & 0
\end{eqnarray*}
so with $c(x) = h(x) = x^2/2$, the equality is satisfied.

The values of $f(x)=-(x-d)^2/2$ are symmetric around $x=d$, so only $d+1$ of 
the $b_j = q^{f(j)}, 0 \leq j \leq d+n$ need to be computed. For example, if
$n = d$, we
can start at $j=d$ and work up to $j=2d$. Computing $q^{f(j)}$ for 
successive $j$ is done by $f(j+1) - f(j) = -j -1/2 + d$,
so $b_{j + 1} = b_j \cdot q^{-j -1/2 + d}$ which for $d \leq j \leq 2d$
requires the sequence $b'_{j'} = q^{-j' -1/2}$ for $0 \leq j' \leq d$.

The values $q^{g(j)}, 0 \leq j \leq d,$ for the $c_j$ sequence can be computed 
for sucessive values of $j$ as well, but it is profitable to do it in reverse 
order, starting at $j = d$. Then we have 
$q^{g(j - 1)} = q^{g(j)} \cdot q^{-j + 1/2}$, and the $q^{-j + 1/2}$ values 
are identical to those computed for the $c'_{j'}$ sequence when $j = j' + 1$,
so we can reuse all the $c'_{j'}$ except for $j' = d$ and only need to compute
$q^{1/2}$ afresh. This way, computing all $c_j$ for $0 \leq j \leq d$
costs $d - 1 + o(d)$ multiplications for the $q^{g(j)}$ values, plus $d$ for 
the multiplication with $f_i$, for a total of $2d + o(d)$.

For the P-1 factoring algorithm, the $q^{h(i)}$ values need not be computed 
as $q \perp N$ ($\perp$ meaning coprime) and so 
$q^{h(i)} \perp N$, and $\gcd(F(q^i) q^{h(i)}, N) = \gcd(F(q^i), N)$.
Hence the total cost for the multipoint evaluation is $4d + o(d)$ 
multiplications and a convolution of length $l$.


\subsection{More general progression}
To evaluate $F(q^{\alpha}), F(q^{\alpha+\beta}), ..., 
F(q^{\alpha+(n-1)\beta})$, we would like instead
\begin{displaymath}
q^{h(i)} g_{i+d} = \sum_{j=0}^{d} f_j q^{j(\alpha+i\beta)}
\end{displaymath}
and thus
\begin{displaymath}
  f(i+d-j) + g(j) + h(i) = \beta ij + \alpha j
\end{displaymath}
Setting $f(x)=-\beta (x-d)^2/2$, we get
\begin{eqnarray*}
-\beta i^2/2 + \beta ij - \beta j^2/2 + g(j) + h(i) & = & \beta ij + \alpha j\\
-\beta i^2/2 - \beta j^2/2 + g(j) + h(i)& = & \alpha j
\end{eqnarray*}
%
which admits a solution with $g(x) = \beta x^2/2 + \alpha x$ and 
$h(x) = \beta x^2/2$. The $f(x)$ values are still symmetric around $x=d$ but
due to the $\alpha$ term in $g(x - 1) - g(x) = -\beta x + \beta/2 - \alpha$,
we cannot reuse the $b'_{j'}$ sequence as before and computing the $c_j$ values
requires $3$ multiplies each, for a total cost for the multipoint evaluation
of $5d + o(d)$ multiplications and a convolution of length $l$.


\subsection{Progression in Montgomery's fast P-1 stage 2}
Montgomery wants to evaluate, in his notation, $f(X)$ with 
$X=b_1^{2k_2 + (2m+1)P}$, for $m_1 \leq m \leq m_2$. In our notation, 
this means $\alpha = 2k_2 + (2m_1+1)P$ and $\beta = 2P$. Further, he requires 
that the sequences $b_j$ and $c_j$ are symmetric so that by suitable shifting 
by $s, t$ we have $b_{j+s} = b_{-(j+s)}$ and $c_{j+t} = c_{-(j+t)}$. 
Let $l = n + d$.

He sets $f(x)=\alpha (n-1+d/2-x)+\beta(n-1+d/2-x)^2/2$ and 
$g(x)=-\beta(x-d/2)^2/2$ and $h(x) = -\beta x^2/2$.
This way, 
\begin{eqnarray*}
g_{2d-1-i} & = & \sum_{0\leq j\leq d} b_{i-j} c_j \\
           & = & \sum_{0\leq j\leq d} q^{f(i-j)} q^{g(j)} f_j
\end{eqnarray*}


\subsection{What about exponents that are powers}
If we want $F(q^{i^S})$, $S>1$, we'll need $f(x)$, $g(x)$ and $h(x)$ so that
\begin{displaymath}
  f(i+d-j) + g(j) = h(i) + i^{S}j
\end{displaymath}
In order to get the $i^{S}j$ term, $f(x)$ must have degree at least $S+1$, 
which for $S>1$ produces at least one term in $i^k j^l$, $k,l>0$ other
than the $i^{S}j$ we want, and those extra terms cannot be absorbed by $g(j)$ 
and $h(i)$. It may work if we convolute more than two sequences, 
but this will require a greater convolution length.

\subsection{What if $F(x)$ is symmetric?}
Let $F(x) = \sum_{j=1}^{d} f_j (x^j + x^{-j})$ of degree $2d$.
We would like
\begin{eqnarray*}
g_{i+d} & = & q^{h(i)} F(q^i) \\
 & = & h_i \sum_{j=1}^{d} f_j \left(q^{ij} + q^{-ij} \right) 
\end{eqnarray*}
where $h_i$ may depend on $q$, $d$ and $i$, but not $j$.

I am hopelessly stuck here. I can't find $b_j$ and $c_j$ so that 
$c_{i-j}b_j = h_i f_j \left(q^{ij} + q^{-ij} \right)$.
\end{document}