Hi codeforces! I want to show that how to implement the algorithm from noshi91 who showed that we could modify the Bostan--Mori's algorithm to make use in FPS composition and compositional inverse.

noshi91's blog: FPS の合成と逆関数、冪乗の係数列挙 Θ(n (log(n))^2) in Japanese.

UPD: my submission. relative code from line 874 to 923.

Compositional Inverse

The key idea is consider the bivariate formal power series

Let $$$f(0)=0$$$, we have

Apply Bostan--Mori's algorithm here:

so

The degree of $$$y$$$ doubled after each time of iteration, but the degree of $$$x$$$ halved.

Combined with Lagrange Inversion Formula, we have the following algorithm:

Composition

Given $$$f=\sum_{k\geq 0}f_kx^k\in\mathbb{C}\left\lbrack\left\lbrack x\right\rbrack\right\rbrack,g\in x\mathbb{C}\left\lbrack\left\lbrack x\right\rbrack\right\rbrack$$$ then

Consider again

Our goal is computing

Apply Bostan--Mori's algorithm here:

We want to compute $$$\dfrac{P(y)}{Q(x,y)Q(-x,y)}$$$ first, and we only need to track finitly many coefficients of $$$y^j$$$ for $$$j\leq 0$$$.

and the composition algorithm is just

One can try with composition_of_formal_power_series_large at Library Checker.

Reference

• Alin Bostan, Ryuhei Mori. A Simple and Fast Algorithm for Computing the N-th Term of a Linearly Recurrent Sequence.
• noshi91. FPS の合成と逆関数、冪乗の係数列挙 Θ(n (log(n))^2).

Hi, Codeforces! I want to introduce a simple way to understand the subset convolution and a simple optimization from scratch.

Actually, I do not know this trick is well-known or just out-of-date.

If there are any mistakes, I am happy to fix them (if I am able to) and thanks for pointing out.

We know that subset convolution is equivalent to truncated multivariate polynomial multiplication like

simply written

we know that we could compute

by using the FFT-like way (Computing $$$A(\mathbf{x})\bmod{(\lbrace x_1,x_1-1\rbrace \times \cdots \times\lbrace x_n,x_n-1\rbrace)}$$$ and using Chinese remainder algorithm to restore the result. This is called the Zeta transform and the Moebius transform respectively). But the coefficients are somehow “mixed”, we can not tell the coefficients of $$$\mathbf{x}^2$$$ and $$$\mathbf{x}$$$. An idea is that we could introduce a new variable $$$t$$$ and computing $$$A(\mathbf{x})B(\mathbf{x})\bmod{\left(\mathbf{x}^2-\mathbf{x}t\right)}$$$, simply let $$$t\gets 0$$$ to get the result of subset convolution. Here is the main idea of the whole algorithm we may familar with:

First, we define $$$s$$$ for each term of $$$A(\mathbf{x}),B(\mathbf{x})$$$:

and

then we compute

For the $$$O(n^2)$$$ times of computation of form $$$A_i(\mathbf{x})B_j(\mathbf{x})\bmod{(\mathbf{x}^2-\mathbf{x}t)}$$$, we have an invariant that $$$s(y)+\deg_t(y)=i+j$$$ for any term $$$y$$$ of $$$A_i(\mathbf{x})B_j(\mathbf{x})\bmod{(\mathbf{x}^2-\mathbf{x}t)}$$$. Even we only do the computation of $$$A_i(\mathbf{x})B_j(\mathbf{x})\bmod{(\mathbf{x}^2-\mathbf{x})}$$$, we could restore the result of $$$A_i(\mathbf{x})B_j(\mathbf{x})\bmod{(\mathbf{x}^2-\mathbf{x}t)}$$$, and we do not use these coefficients that with variable $$$t$$$.

Another observation is that $$$\deg_t(A_i(\mathbf{x})B_j(\mathbf{x})\bmod{(\mathbf{x}^2-\mathbf{x}t)})\leq (i+j)/2$$$. For example, if we have $$$\left(\sum_{i+j=4}A_i(\mathbf{x})B_j(\mathbf{x})\right)+\left(\sum_{i+j=9}A_i(\mathbf{x})B_j(\mathbf{x})\right)\bmod{(\mathbf{x}^2-\mathbf{x}t)}$$$, it is still possible for us to tell apart. With this simple trick, we could reduce the times of Moebius transform for about a half.

This is a submission for the optimized code. I did not implement it carefully.

upd: submission reinplemented.

A simple way to understand the transposed algo of multi-eval

A simple way to understand the transposed algorithm of multi-point evaluation of polynomials in Tellegen’s Principle into Practice.

I learned this from Bernstein's paper Scaled Remainder Trees. If there are any mistakes, please tell me, thanks! (Sorry for my poor English)

Let $$$\mathbb{R}((x))$$$ denote the ring of formal Laurent series, $$$\mathbb{R}\lbrack x\rbrack$$$ denote the ring of polynomial. $$$f(x)\bmod 1:=f_{-1}x^{-1}+f_{-2}x^{-2}+\cdots$$$.

Bernstein showed that division in $$$\mathbb{R}\lbrack x\rbrack$$$ is division in $$$\mathbb{R}((x^{-1}))$$$. Here is a concrete example. Let $$$A(x):=1+x+4x^2+5x^3+x^4+4x^5\in\mathbb{R}\lbrack x\rbrack$$$ and $$$B(x):=1+9x+8x^2+x^3\in\mathbb{R}\lbrack x\rbrack$$$, we want to compute $$$Q(x),R(x)\in\mathbb{R}\lbrack x\rbrack$$$ such that $$$A(x)=B(x)Q(x)+R(x)$$$ and $$$\deg R(x)\lt \deg B(x)$$$ (with $$$\deg 0=-\infty$$$). Let's do this in $$$\mathbb{R}((x^{-1}))$$$. First, we should find the multiplicative inverse of $$$B(x)$$$ which is

and compute

We could find that $$$Q(x)$$$ is exactly $$$A(x)B^{-1}(x)-(A(x)B^{-1}(x)\bmod 1)$$$, and the leading coefficient of $$$R(x)$$$ equals $$$[x^{-1}]A(x)B^{-1}(x)$$$. That's not a coincidence. Consider the simplest case

So we have $$$f(v)=[x^{-1}]f/(x-v)$$$ cause

The problem of multi-point evaluation is as follows: Given polynomial $$$f(x)=\sum_{i=0}^nf_ix^i$$$ and $$$v_0,v_1,\dots ,v_n$$$, print $$$f(v_0),f(v_1),\dots ,f(v_n)$$$.

The traditional algorithm for multi-point evaluation is that we make an almost balanced binary (sub)product tree with leaves are polynomials $$$p_i(x)=x-v_i$$$ and non-leaf nodes are the product of its children. Then we compute $$$f\bmod (p_1\cdots p_n)$$$ and then compute $$$f\bmod (p_1\cdots p_{n/2})$$$ and $$$f\bmod (p_{n/2+1}\cdots p_n)$$$ cause $$$(f\bmod (p_1\cdots p_n))\bmod (p_1\cdots p_{n/2})=f\bmod (p_1\cdots p_{n/2})$$$ (suppose $$$n$$$ is even) and then ... The transposed algorithm is that we compute $$$f/(p_1\cdots p_n)$$$ and $$$f/(p_1\cdots p_n)\cdot (p_{n/2+1}\cdots p_{n})=f/(p_1\cdots p_{n/2})$$$, finally we will have $$$f/p_i$$$ at the leaf of the tree, the coefficient of $$$x^{-1}$$$ is exactly what we want. The only remaining problem is how many terms should we track when doing the computation. It's clear that we don't care about the coefficients of $$$x^0,x^1,\dots$$$ in $$$f/(p_1\cdots p_n)$$$ cause they do nothing with the coefficient of $$$x^{-1},x^{-2},\dots$$$ when $$$f/(p_1\cdots p_n)$$$ was multiplied by a polynomial like $$$g_0+g_1x+g_2x^2+\cdots$$$. Another observation is that for $$$f/P\bmod 1=a_{-1}x^{-1}+a_{-2}x^{-2}+\cdots +a_{-\deg P}x^{\deg P}+\cdots$$$, only $$$a_{-1},\dots ,a_{-\deg P}$$$ are enough. Consider a very simple case that we have $$$f/(p_1p_2)\bmod 1=c_{-1}x^{-1}+c_{-2}x^{-2}+c_{-3}x^{-3}+\cdots$$$ and want to compute $$$f/p_1\bmod 1$$$ via multiplying $$$f/(p_1p_2)\bmod 1$$$ by $$$p_2$$$, the coefficients of $$$x^{-3},x^{-4},\dots$$$ have no effect on the coefficient of $$$x^{-1}$$$ (it's more complicated for the integer case). So we could almost halve the size of problem after one multiplication.

About how to compute $$$(p_1\cdots p_n)^{-1}$$$ in $$$\mathbb{R}((x^{-1}))$$$, we could simply use the method that we compute the multiplicative inverse of formal power series $$$x^{\deg(p_1\cdots p_n)}p_1(x^{-1})\cdots p_n(x^{-1})$$$ in $$$\mathbb{R}\lbrack \lbrack x\rbrack \rbrack$$$.

Simple C++(17) implementation for univariate formal power series computation in competitive programming.

a.hpp with main function.

Idea from fjzzq2002's blog (in Chinese). My code is about 2 times slower than fjzzq2002's, but mine is shorter and simpler.

I am glad if you would tell me about bugs or optimizations.

There are several mistakes in the comments of the code ... sorry.

Usage

GF for Catalan numbers A000108.

#include <iostream>
#include "a.hpp"
int main() {
  using lib::Mint;
  using lib::FPS;
  FPS<Mint<998244353>> A;
  auto X = FPS<Mint<998244353>>::idm();
  A.reset().set(X / (1 - A));
  for (int i = 0; i <= 10; ++i) std::cout << A[i] << ' ';
}

GF for alkyl A000598.

#include <iostream>
#include "a.hpp"
int main() {
  using lib::Mint;
  using lib::FPS;
  FPS<Mint<998244353>> A;
  auto X = FPS<Mint<998244353>>::idm();
  A.reset().set((A.pow(3) / 6 + A * A(X ^ 2) / 2 + A(X ^ 3) / 3) * X + 1);
  for (int i = 0; i <= 10; ++i) std::cout << A[i] << ' ';
}

GF for unlabeled trees A000081.

#include <iostream>
#include "a.hpp"
int main() {
  using lib::Mint;
  using lib::FPS;
  FPS<Mint<998244353>> A;
  auto X = FPS<Mint<998244353>>::idm();
  A.reset().set(A.Exp() * X);
  auto B = A - (A * A - A(X ^ 2)) / 2;
  for (int i = 0; i <= 10; ++i) std::cout << B[i] << ' ';
}

Apply Euler transform several times A144036.

#include <iostream>
#include "a.hpp"
int main() {
  using lib::Mint;
  using lib::FPS;
  FPS<Mint<998244353>> A;
  auto X = FPS<Mint<998244353>>::idm();
  A.reset().set((((A.Exp() - 1).Exp() - 1).Exp() - 1).Exp() * X);
  for (int i = 0; i <= 10; ++i) std::cout << A[i] << ' ';
}

License

CC0 1.0.

Reference

• Daniel J. Bernstein. The tangent FFT.
• J. van der Hoeven. Relax, but don’t be too lazy. JSC, 34:479–542, 2002.
• van der Hoeven, J., 2003a. New algorithms for relaxed multiplication. Tech. Rep. 2003–44, Universit´e Paris-Sud, Orsay, France.
• Philippe Flajolet and Robert Sedgewick. Analytic Combinatorics.

Suppose $$$p$$$ is an odd prime and $$$a$$$ is a quadratic residue modulo $$$p$$$. Cipolla's algorithm shows that $$$b:=x^{(p+1)/2}\bmod{(x^2-tx+a)}$$$ such that $$$b^2=a$$$ for some irreducible polynomial $$$x^2-tx+a\in\mathbb{F}_p[x]$$$.

I don't know how to prove this properly. If there are any mistakes or typos, please let me know, thanks! Let $$$\alpha$$$ be a zero of $$$x^2-tx+a$$$, $$$\mathbb{F}_p(\alpha):=\lbrace a_0+a_1\alpha :a_0,a_1\in\mathbb{F}_p\rbrace$$$, and let $$$\beta$$$ be a zero of $$$x^2-(t^2-4a)$$$, $$$\mathbb{F}_p(\beta):=\lbrace a_0+a_1\beta :a_0,a_1\in\mathbb{F}_p\rbrace$$$. We may easily find homomorphisms between $$$\mathbb{F}_p(\alpha)$$$ and $$$\mathbb{F}_p(\beta)$$$. So I think we just need to show that $$$((t+\beta )/2)^{p+1}=a$$$.

A lot of blogs show that $$$(t+\beta)^p=t-\beta$$$ according to binomial theorem and Fermat's little theorem, so

Bostan and Mori's paper shows that the computation of $$$x^n$$$ modulo a monic polynomial sometimes is equivalent to the computation of one term of a rational function $$$P(x)/Q(x)$$$ where $$$P,Q$$$ are both polynomials and $$$\deg(P(x))\lt \deg(Q(x))$$$.

We have

and

This may save some multiplications.

submission(Library Checker)

Actually, the initial value of $$$k_1$$$ does no matter to the result. So we could mainly focus on the computation of one term of the inverse of formal power series. I don't know if we could do better.