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.


2 Answers 2


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)$.

  • $\begingroup$ I'm not sure if there's missing context or what, but this doesn't seem to hold and I'm not sure it can hold, given that you would expect some form of binomial expansion in the second step $\endgroup$
    – b3m2a1
    Apr 14, 2023 at 1:01
  • $\begingroup$ I'm using the substitution $i = i+d$, and the fact that the indices are modulo $n$ (since $\omega^n = 1$). $\endgroup$ Apr 20, 2023 at 18:42
  • $\begingroup$ To be a bit clearer, maybe you could provide an example, because when I try this approach out using Mathematica to expand the polynomial coefficients in the x+d polynomial the approach as you detailed it the Fourier coefficients differ significantly. Here's a sample polynomial for which it does not work 10 + 2 n + 4 n^2 + 3 n^3 - 2 n^4 $\endgroup$
    – b3m2a1
    Apr 21, 2023 at 17:13
  • $\begingroup$ Unfortunately that's beyond my paygrade. $\endgroup$ Apr 24, 2023 at 15:39
  • $\begingroup$ Okay, but not to put too fine a point on this, I think what you posted is incorrect and I'd love to be proven wrong but I have no idea why someone would accept an answer that simply doesn't work $\endgroup$
    – b3m2a1
    Apr 24, 2023 at 18:29

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_] :=
  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_] :=
   nco = Length@coeffs - 1
     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_] :=
     ArrayPad[c1, {0, # - 1} & /@ Dimensions[c2]], 
     FourierParameters -> {1, -1}
     ArrayPad[c2, {0, # - 1} & /@ Dimensions[c1]], 
     FourierParameters -> {1, -1}
   FourierParameters -> {1, -1}
tayShiftFourierND[coeffs_, shift_] :=
   ncos = Dimensions[coeffs] - 1,
  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[
     ][[Sequence @@ Table[n ;;, {n, Dimensions@coeffs}] ]];

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


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


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.