Efficient algorithm to translate/shift polynomials

Question

I have a polynomial $P(x)$, and given some constant $d$, I need to find the polynomial $P(x+d)$. For example, if $P(x)=x^2$ and $d=1$, then the result would be $P(x+1)=(x+1)^2=x^2+2x+1$ (with the coefficients stored in an array/vector).

I know the algorithm with $O(n^2)$ time complexity, which is based on Horner's method: $$a_0+(x+d)(a_1+(x+d)(a_2+...+(x+d)a_n))$$ However, it is not efficient enough for my needs. Is there an algorithm that is more efficient for solving this problem?

P.S. The coefficients and $d$ will always be integers, if that helps.

Answer 1

Use FFT. Suppose that $\deg P \leq n$ and we use FFT on $n$ points (usually $n$ would be the smallest power of 2 which is at least $\deg P$). Let $\omega$ be an $n$th root of unity. The Fourier coefficients of $Q(x) = P(x+d)$ are $$ \hat{Q}(m) = \sum_{i=0}^{n-1} Q(i) \omega^{mi} = \sum_{i=0}^{n-1} P(i+d) \omega^{mi} = \sum_{i=0}^{n-1} P(i) \omega^{m(i-d)} = \omega^{-md} \hat{P}(m). $$ This shows how to compute the Fourier transform of $Q$ from that of $P$ in linear time. Since FFT takes time $O(n\log n)$ (for suitable $n$), the overall running time is $O(n\log n)$. Choosing $n$ as indicated above, this becomes $O(\deg P \log \deg P)$.

Answer 2

Here's another FFT-based approach that I can confirm works. This builds off of Algorithm F in Fast algorithms for Taylor shifts and certain difference equations and has a natural extension to multivariate polynomials and shifts

Assuming our polynomial is of degree $n$, given the polynomial coefficient vector ${p_k}$ and a shift $a$, we'll construct the intermediate polynomial coefficient vectors

$$ \begin{align} u_k &= (n-k)!\ p_{(n-k)} \\ v_k &= a^k / k! \end{align} $$

By treating these as polynomial coefficients and taking their direct convolution, $g_k$ (which an be done in $O(n \log{n})$ time via FFT with appropriate zero-padding) we get our final polynomial coefficients as

$$ q_k = g_{(n-k)} / k! $$

Here's Mathematica code for that

fourierConvolve[c1_, c2_] :=
 InverseFourier[
  Fourier[PadRight[c1, Length@c1 + Length@c2 - 1], FourierParameters -> {1, -1}] *
    Fourier[PadRight[c2, Length@c1 + Length@c2 - 1], FourierParameters -> {1, -1}],
  FourierParameters -> {1, -1}
  ]

tayShiftFourier[coeffs_, a_] :=
 With[{
   nco = Length@coeffs - 1
   },
  fourierConvolve[
     Reverse[coeffs*Factorial[Range[0, nco]]], 
     a^Range[0, nco]/Factorial[Range[0, nco]]
     ][[nco + 1 ;; 1 ;; -1]]/Factorial[Range[0, nco]]
  ]

As noted, this generalizes easily to multivariate polynomials. Considering that one could apply this approach to each dimension recursively, we pretty naturally find that we can replace our coefficient sequence $P$ with a tensor of coefficients, our factorials can be replaced with the outer-product tensor of factorials in each dimension, and the powers of the shifts

Here's the N-dimensional analog of the prior code, although I did some small modifications in orders of operations to avoid having to reverse my coefficient tensors

Clear[tayShiftFourierND, fourierConvolveND];
fourierConvolveND[c1_, c2_] :=
  InverseFourier[
   Fourier[
     ArrayPad[c1, {0, # - 1} & /@ Dimensions[c2]], 
     FourierParameters -> {1, -1}
     ]*
    Fourier[
     ArrayPad[c2, {0, # - 1} & /@ Dimensions[c1]], 
     FourierParameters -> {1, -1}
     ],
   FourierParameters -> {1, -1}
   ];
tayShiftFourierND[coeffs_, shift_] :=
 Block[
  {
   ncos = Dimensions[coeffs] - 1,
   facTensor,
   shiftTensor,
   revFacTensor,
   convolvedCoeffs
   },
  facTensor = Outer[Times, ##] & @@ Table[Factorial[Range[n, 0, -1]], {n, ncos}];
  shiftTensor = 
   Outer[Times, ##] & @@ MapThread[#^Range[#2, 0, -1] &, {shift, ncos}];
  revFacTensor = Map[Reverse, facTensor, Range[0, Length@ncos - 1]];
  convolvedCoeffs = fourierConvolveND[
      coeffs*revFacTensor, 
     shiftTensor/facTensor
     ][[Sequence @@ Table[n ;;, {n, Dimensions@coeffs}] ]];
  convolvedCoeffs/revFacTensor
  ]

and here's some proof that it works

coeffs3D = BlockRandom@RandomReal[{-10, 10}, {5, 3, 4}];
poly3D = Total@
   Flatten[coeffs3D*Outer[Times, x^Range[0, 4], y^Range[0, 2], z^Range[0, 3]]];

tayShiftFourierND[coeffs3D, {1, 2, 3}] // Chop // Flatten // Norm

4976.12

poly3D /. {x -> x + 1, y -> y + 2, z -> z + 3} // Expand // 
    CoefficientList[#, {x, y, z}] & // Chop // Flatten // Norm

4976.12

where I just used the norm of the flattened tensor of polynomial coefficients after shifting as effectively a hash