Skip to content

Operations on polynomials and series

Problems in competitive programming, especially the ones involving enumeration some kind, are often solved by reducing the problem to computing something on polynomials and formal power series.

This includes concepts such as polynomial multiplication, interpolation, and more complicated ones, such as polynomial logarithms and exponents. In this article, a brief overview of such operations and common approaches to them is presented.

Basic Notion and facts

In this section, we focus more on the definitions and "intuitive" properties of various polynomial operations. The technical details of their implementation and complexities will be covered in later sections.

Polynomial multiplication

Definition

Univariate polynomial is an expression of form A(x)=a0+a1x++anxn.

The values a0,,an are polynomial coefficients, typically taken from some set of numbers or number-like structures. In this article, we assume that the coefficients are taken from some field, meaning that operations of addition, subtraction, multiplication and division are well-defined for them (except for division by 0) and they generally behave in a similar way to real numbers.

Typical example of such field is the field of remainders modulo prime number p.

For simplicity we will drop the term univariate, as this is the only kind of polynomials we consider in this article. We will also write A instead of A(x) wherever possible, which will be understandable from the context. It is assumed that either an0 or A(x)=0.

Definition

The product of two polynomials is defined by expanding it as an arithmetic expression:

A(x)B(x)=(i=0naixi)(j=0mbjxj)=i,jaibjxi+j=k=0n+mckxk=C(x).

The sequence c0,c1,,cn+m of the coefficients of C(x) is called the convolution of a0,,an and b0,,bm.

Definition

The degree of a polynomial A with an0 is defined as degA=n.

For consistency, degree of A(x)=0 is defined as degA=.

In this notion, degAB=degA+degB for any polynomials A and B.

Convolutions are the basis of solving many enumerative problems.

Example

You have n objects of the first kind and m objects of the second kind.

Objects of first kind are valued a1,,an, and objects of the second kind are valued b1,,bm.

You pick a single object of the first kind and a single object of the second kind. How many ways are there to get the total value k?

Solution

Consider the product (xa1++xan)(xb1++xbm). If you expand it, each monomial will correspond to the pair (ai,bj) and contribute to the coefficient near xai+bj. In other words, the answer is the coefficient near xk in the product.

Example

You throw a 6-sided die n times and sum up the results from all throws. What is the probability of getting sum of k?

Solution

The answer is the number of outcomes having the sum k, divided by the total number of outcomes, which is 6n.

What is the number of outcomes having the sum k? For n=1, it may be represented by a polynomial A(x)=x1+x2++x6.

For n=2, using the same approach as in the example above, we conclude that it is represented by the polynomial (x1+x2++x6)2.

That being said, the answer to the problem is the k-th coefficient of (x1+x2++x6)n, divided by 6n.

The coefficient near xk in the polynomial A(x) is denoted shortly as [xk]A.

Formal power series

Definition

A formal power series is an infinite sum A(x)=a0+a1x+a2x2+, considered regardless of its convergence properties.

In other words, when we consider e.g. a sum 1+12+14+18+=2, we imply that it converges to 2 when the number of summands approach infinity. However, formal series are only considered in terms of sequences that make them.

Definition

The product of formal power series A(x) and B(x), is also defined by expanding it as an arithmetic expression:

A(x)B(x)=(i=0aixi)(j=0bjxj)=i,jaibjxi+j=k=0ckxk=C(x),

where the coefficients c0,c1, are define as finite sums

ck=i=0kaibki.

The sequence c0,c1, is also called a convolution of a0,a1, and b0,b1,, generalizing the concept to infinite sequences.

Thus, polynomials may be considered formal power series, but with finite number of coefficients.

Formal power series play a crucial role in enumerative combinatorics, where they're studied as generating functions for various sequences. Detailed explanation of generating functions and the intuition behind them will, unfortunately, be out of scope for this article, therefore the curious reader is referenced e.g. here for details about their combinatorial meaning.

However, we will very briefly mention that if A(x) and B(x) are generating functions for sequences that enumerate some objects by number of "atoms" in them (e.g. trees by the number of vertices), then the product A(x)B(x) enumerates objects that can be described as pairs of objects of kinds A and B, enumerates by the total number of "atoms" in the pair.

Example

Let A(x)=i=02ixi enumerate packs of stones, each stone colored in one of 2 colors (so, there are 2i such packs of size i) and B(x)=j=03jxj enumerate packs of stones, each stone colored in one of 3 colors. Then C(x)=A(x)B(x)=k=0ckxk would enumerate objects that may be described as "two packs of stones, first pack only of stones of type A, second pack only of stones of type B, with total number of stones being k" for ck.

In a similar way, there is an intuitive meaning to some other functions over formal power series.

Long polynomial division

Similar to integers, it is possible to define long division on polynomials.

Definition

For any polynomials A and B0, one may represent A as

A=DB+R, degR<degB,

where R is called the remainder of A modulo B and D is called the quotient.

Denoting degA=n and degB=m, naive way to do it is to use long division, during which you multiply B by the monomial anbmxnm and subtract it from A, until the degree of A is smaller than that of B. What remains of A in the end will be the remainder (hence the name), and the polynomials with which you multiplied B in the process, summed together, form the quotient.

Definition

If A and B have the same remainder modulo C, they're said to be equivalent modulo C, which is denoted as

AB(modC).

Polynomial long division is useful because of its many important properties:

  • A is a multiple of B if and only if A0(modB).

  • It implies that AB(modC) if and only if AB is a multiple of C.

  • In particular, AB(modCD) implies AB(modC).

  • For any linear polynomial xr it holds that A(x)A(r)(modxr).

  • It implies that A is a multiple of xr if and only if A(r)=0.

  • For modulo being xk, it holds that Aa0+a1x++ak1xk1(modxk).

Note that long division can't be properly defined for formal power series. Instead, for any A(x) such that a00, it is possible to define an inverse formal power series A1(x), such that A(x)A1(x)=1. This fact, in turn, can be used to compute the result of long division for polynomials.

Basic implementation

Here you can find the basic implementation of polynomial algebra.

It supports all trivial operations and some other useful methods. The main class is poly<T> for polynomials with coefficients of type T.

All arithmetic operation +, -, *, % and / are supported, % and / standing for remainder and quotient in Euclidean division.

There is also the class modular<m> for performing arithmetic operations on remainders modulo a prime number m.

Other useful functions:

  • deriv(): computes the derivative P(x) of P(x).
  • integr(): computes the indefinite integral Q(x)=P(x) of P(x) such that Q(0)=0.
  • inv(size_t n): calculate the first n coefficients of P1(x) in O(nlogn).
  • log(size_t n): calculate the first n coefficients of lnP(x) in O(nlogn).
  • exp(size_t n): calculate the first n coefficients of expP(x) in O(nlogn).
  • pow(size_t k, size_t n): calculate the first n coefficients for Pk(x) in O(nlognk).
  • deg(): returns the degree of P(x).
  • lead(): returns the coefficient of xdegP(x).
  • resultant(poly<T> a, poly<T> b): computes the resultant of a and b in O(|a||b|).
  • bpow(T x, size_t n): computes xn.
  • bpow(T x, size_t n, T m): computes xn(modm).
  • chirpz(T z, size_t n): computes P(1),P(z),P(z2),,P(zn1) in O(nlogn).
  • vector<T> eval(vector<T> x): evaluates P(x1),,P(xn) in O(nlog2n).
  • poly<T> inter(vector<T> x, vector<T> y): interpolates a polynomial by a set of pairs P(xi)=yi in O(nlog2n).
  • And some more, feel free to explore the code!

Arithmetic

Multiplication

The very core operation is the multiplication of two polynomials. That is, given the polynomials A and B:

A=a0+a1x++anxn
B=b0+b1x++bmxm

You have to compute polynomial C=AB, which is defined as

C=i=0nj=0maibjxi+j=c0+c1x++cn+mxn+m.

It can be computed in O(nlogn) via the Fast Fourier transform and almost all methods here will use it as subroutine.

Inverse series

If A(0)0 there always exists an infinite formal power series A1(x)=q0+q1x+q2x2+ such that A1A=1. It often proves useful to compute first k coefficients of A1 (that is, to compute it modulo xk). There are two major ways to calculate it.

Divide and conquer

This algorithm was mentioned in Schönhage's article and is inspired by Graeffe's method. It is known that for B(x)=A(x)A(x) it holds that B(x)=B(x), that is, B(x) is an even polynomial. It means that it only has non-zero coefficients with even numbers and can be represented as B(x)=T(x2). Thus, we can do the following transition:

A1(x)1A(x)A(x)A(x)A(x)A(x)T(x2)(modxk)

Note that T(x) can be computed with a single multiplication, after which we're only interested in the first half of coefficients of its inverse series. This effectively reduces the initial problem of computing A1(modxk) to computing T1(modxk/2).

The complexity of this method can be estimated as

T(n)=T(n/2)+O(nlogn)=O(nlogn).

Sieveking–Kung algorithm

The generic process described here is known as Hensel lifting, as it follows from Hensel's lemma. We'll cover it in more detail further below, but for now let's focus on ad hoc solution. "Lifting" part here means that we start with the approximation B0=q0=a01, which is A1(modx) and then iteratively lift from modxa to modx2a.

Let BkA1(modxa). The next approximation needs to follow the equation ABk+11(modx2a) and may be represented as Bk+1=Bk+xaC. From this follows the equation

A(Bk+xaC)1(modx2a).

Let ABk1+xaD(modx2a), then the equation above implies

xa(D+AC)0(modx2a)DAC(modxa)CBkD(modxa).

From this, one can obtain the final formula, which is

xaCBkxaDBk(1ABk)(modx2a)Bk+1Bk(2ABk)(modx2a)

Thus starting with B0a01(modx) we will compute the sequence Bk such that ABk1(modx2k) with the complexity

T(n)=T(n/2)+O(nlogn)=O(nlogn).

The algorithm here might seem a bit more complicated than the first one, but it has a very solid and practical reasoning behind it, as well as a great generalization potential if looked from a different perspective, which would be explained further below.

Euclidean division

Consider two polynomials A(x) and B(x) of degrees n and m. As it was said earlier you can rewrite A(x) as

A(x)=B(x)D(x)+R(x),degR<degB.

Let nm, it would imply that degD=nm and the leading nm+1 coefficients of A don't influence R. It means that you can recover D(x) from the largest nm+1 coefficients of A(x) and B(x) if you consider it as a system of equations.

The system of linear equations we're talking about can be written in the following form:

[anam+1am]=[bm00bm0bm1bm][dnmd1d0]

From the looks of it, we can conclude that with the introduction of reversed polynomials

AR(x)=xnA(x1)=an+an1x++a0xn
BR(x)=xmB(x1)=bm+bm1x++b0xm
DR(x)=xnmD(x1)=dnm+dnm1x++d0xnm

the system may be rewritten as

AR(x)BR(x)DR(x)(modxnm+1).

From this you can unambiguously recover all coefficients of D(x):

DR(x)AR(x)(BR(x))1(modxnm+1)

And from this, in turn, you can recover R(x) as R(x)=A(x)B(x)D(x).

Note that the matrix above is a so-called triangular Toeplitz matrix and, as we see here, solving system of linear equations with arbitrary Toeplitz matrix is, in fact, equivalent to polynomial inversion. Moreover, inverse matrix of it would also be triangular Toeplitz matrix and its entries, in terms used above, are the coefficients of (BR(x))1(modxnm+1).

Calculating functions of polynomial

Newton's method

Let's generalize the Sieveking–Kung algorithm. Consider equation F(P)=0 where P(x) should be a polynomial and F(x) is some polynomial-valued function defined as

F(x)=i=0αi(xβ)i,

where β is some constant. It can be proven that if we introduce a new formal variable y, we can express F(x) as

F(x)=F(y)+(xy)F(y)+(xy)2G(x,y),

where F(x) is the derivative formal power series defined as

F(x)=i=0(i+1)αi+1(xβ)i,

and G(x,y) is some formal power series of x and y. With this result we can find the solution iteratively.

Let F(Qk)0(modxa). We need to find Qk+1Qk+xaC(modx2a) such that F(Qk+1)0(modx2a).

Substituting x=Qk+1 and y=Qk in the formula above, we get

F(Qk+1)F(Qk)+(Qk+1Qk)F(Qk)+(Qk+1Qk)2G(x,y)(modx)2a.

Since Qk+1Qk0(modxa), it also holds that (Qk+1Qk)20(modx2a), thus

0F(Qk+1)F(Qk)+(Qk+1Qk)F(Qk)(modx2a).

The last formula gives us the value of Qk+1:

Qk+1=QkF(Qk)F(Qk)(modx2a)

Thus, knowing how to invert polynomials and how to compute F(Qk), we can find n coefficients of P with the complexity

T(n)=T(n/2)+f(n),

where f(n) is the time needed to compute F(Qk) and F(Qk)1 which is usually O(nlogn).

The iterative rule above is known in numerical analysis as Newton's method.

Hensel's lemma

As was mentioned earlier, formally and generically this result is known as Hensel's lemma and it may in fact used in even broader sense when we work with a series of nested rings. In this particular case we worked with a sequence of polynomial remainders modulo x, x2, x3 and so on.

Another example where Hensel's lifting might be helpful are so-called p-adic numbers where we, in fact, work with the sequence of integer remainders modulo p, p2, p3 and so on. For example, Newton's method can be used to find all possible automorphic numbers (numbers that end on itself when squared) with a given number base. The problem is left as an exercise to the reader. You might consider this problem to check if your solution works for 10-based numbers.

Logarithm

For the function lnP(x) it's known that:

(lnP(x))=P(x)P(x)

Thus we can calculate n coefficients of lnP(x) in O(nlogn).

Inverse series

Turns out, we can get the formula for A1 using Newton's method. For this we take the equation A=Q1, thus:

F(Q)=Q1A
F(Q)=Q2
Qk+1Qk(2AQk)(modx2k+1)

Exponent

Let's learn to calculate eP(x)=Q(x). It should hold that lnQ=P, thus:

F(Q)=lnQP
F(Q)=Q1
Qk+1Qk(1+PlnQk)(modx2k+1)

k-th power

Now we need to calculate Pk(x)=Q. This may be done via the following formula:

Q=exp[klnP(x)]

Note though, that you can calculate the logarithms and the exponents correctly only if you can find some initial Q0.

To find it, you should calculate the logarithm or the exponent of the constant coefficient of the polynomial.

But the only reasonable way to do it is if P(0)=1 for Q=lnP so Q(0)=0 and if P(0)=0 for Q=eP so Q(0)=1.

Thus you can use formula above only if P(0)=1. Otherwise if P(x)=αxtT(x) where T(0)=1 you can write that:

Pk(x)=αkxktexp[klnT(x)]

Note that you also can calculate some k-th root of a polynomial if you can calculate αk, for example for α=1.

Evaluation and Interpolation

Chirp-z Transform

For the particular case when you need to evaluate a polynomial in the points xr=z2r you can do the following:

A(z2r)=k=0nakz2kr

Let's substitute 2kr=r2+k2(rk)2. Then this sum rewrites as:

A(z2r)=zr2k=0n(akzk2)z(rk)2

Which is up to the factor zr2 equal to the convolution of the sequences uk=akzk2 and vk=zk2.

Note that uk has indexes from 0 to n here and vk has indexes from n to m where m is the maximum power of z which you need.

Now if you need to evaluate a polynomial in the points xr=z2r+1 you can reduce it to the previous task by the transformation akakzk.

It gives us an O(nlogn) algorithm when you need to compute values in powers of z, thus you may compute the DFT for non-powers of two.

Another observation is that kr=(k+r2)(k2)(r2). Then we have

A(zr)=z(r2)k=0n(akz(k2))z(k+r2)

The coefficient of xn+r of the product of the polynomials A0(x)=k=0nankz(nk2)xk and A1(x)=k0z(k2)xk equals z(r2)A(zr). You can use the formula z(k+12)=z(k2)+k to calculate the coefficients of A0(x) and A1(x).

Multi-point Evaluation

Assume you need to calculate A(x1),,A(xn). As mentioned earlier, A(x)A(xi)(modxxi). Thus you may do the following:

  1. Compute a segment tree such that in the segment [l,r) stands the product Pl,r(x)=(xxl)(xxl+1)(xxr1).
  2. Starting with l=1 and r=n+1 at the root node. Let m=(l+r)/2. Let's move down to [l,m) with the polynomial A(x)(modPl,m(x)).
  3. This will recursively compute A(xl),,A(xm1), now do the same for [m,r) with A(x)(modPm,r(x)).
  4. Concatenate the results from the first and second recursive call and return them.

The whole procedure will run in O(nlog2n).

Interpolation

There's a direct formula by Lagrange to interpolate a polynomial, given set of pairs (xi,yi):

A(x)=i=1nyijixxjxixj

Computing it directly is a hard thing but turns out, we may compute it in O(nlog2n) with a divide and conquer approach:

Consider P(x)=(xx1)(xxn). To know the coefficients of the denominators in A(x) we should compute products like:

Pi=ji(xixj)

But if you consider the derivative P(x) you'll find out that P(xi)=Pi. Thus you can compute Pi's via evaluation in O(nlog2n).

Now consider the recursive algorithm done on same segment tree as in the multi-point evaluation. It starts in the leaves with the value yiPi in each leaf.

When we return from the recursion we should merge the results from the left and the right vertices as Al,r=Al,mPm,r+Pl,mAm,r.

In this way when you return back to the root you'll have exactly A(x) in it. The total procedure also works in O(nlog2n).

GCD and Resultants

Assume you're given polynomials A(x)=a0+a1x++anxn and B(x)=b0+b1x++bmxm.

Let λ0,,λn be the roots of A(x) and let μ0,,μm be the roots of B(x) counted with their multiplicities.

You want to know if A(x) and B(x) have any roots in common. There are two interconnected ways to do that.

Euclidean algorithm

Well, we already have an article about it. For an arbitrary domain you can write the Euclidean algorithm as easy as:

template<typename T>
T gcd(const T &a, const T &b) {
    return b == T(0) ? a : gcd(b, a % b);
}

It can be proven that for polynomials A(x) and B(x) it will work in O(nm).

Resultant

Let's calculate the product A(μ0)A(μm). It will be equal to zero if and only if some μi is the root of A(x).

For symmetry we can also multiply it with bmn and rewrite the whole product in the following form:

R(A,B)=bmnj=0mA(μj)=bmnamni=0nj=0m(μjλi)=(1)mnanmi=0nB(λi)

The value defined above is called the resultant of the polynomials A(x) and B(x). From the definition you may find the following properties:

  1. R(A,B)=(1)nmR(B,A).
  2. R(A,B)=anmbmn when n=0 or m=0.
  3. If bm=1 then R(ACB,B)=R(A,B) for an arbitrary polynomial C(x) and n,m1.
  4. From this follows R(A,B)=bmdeg(A)deg(ACB)R(ACB,B) for arbitrary A(x), B(x), C(x).

Miraculously it means that resultant of two polynomials is actually always from the same ring as their coefficients!

Also these properties allow us to calculate the resultant alongside the Euclidean algorithm, which works in O(nm).

template<typename T>
T resultant(poly<T> a, poly<T> b) {
    if(b.is_zero()) {
        return 0;
    } else if(b.deg() == 0) {
        return bpow(b.lead(), a.deg());
    } else {
        int pw = a.deg();
        a %= b;
        pw -= a.deg();
        base mul = bpow(b.lead(), pw) * base((b.deg() & a.deg() & 1) ? -1 : 1);
        base ans = resultant(b, a);
        return ans * mul;
    }
}

Half-GCD algorithm

There is a way to calculate the GCD and resultants in O(nlog2n).

The procedure to do so implements a 2×2 linear transform which maps a pair of polynomials a(x), b(x) into another pair c(x),d(x) such that degd(x)dega(x)2. If you're careful enough, you can compute the half-GCD of any pair of polynomials with at most 2 recursive calls to the polynomials which are at least 2 times smaller.

The specific details of the algorithm are somewhat tedious to explain, however you can find its implementation in the library, as half_gcd function.

After half-GCD is implemented, you can repeatedly apply it to polynomials until you're reduced to the pair of gcd(a,b) and 0.

Problems