Output looks like Codeforces titles.

Codeforces rating -> title

The title and the fact that the input is represented by $$$c$$$ is a hint.

What well-known function has a fixed point at $$$-40$$$?

Celsius -> Fahrenheit

Sample 2 has output $$$28$$$, which is $$$\binom{8}{2}$$$.

$$$28$$$ is the number of nonempty substrings of sample $$$2$$$.

Number of distinct nonempty substrings.

We are looking for a word from the grids.

Exactly one word appears in all the grids with answer YES and does not appear in all the grids with answer NO.

Output YES iff "vapid" can be found in the grid (vertically, horizontally or diagonally).

The title is asking you to factor the sample outputs and spot a pattern.

Look for factors coming from well-known sequences.

The answer is Fib(n)^2*Catalan(n). Originally I planned to use Fib(n)*Catalan(n) but unfortunately it appears in OEIS for some reason.

The queries look like subrectangle queries. Play around with the samples. # contributes to the answer somehow. Answer is always even.

Compute the "perimeter" of the # region in the query subrectangle.

The input makes no sense and seems unrelated to the title. The title seems to be our only clue?

Notorious coincidence.

Google “Coloring Trees codeforces” (or something similar) will lead you to https://mirror.codeforces.com/problemset/problem/711/C

Input $$$1$$$ is $$$\frac{1}{6} \pmod{998244353}$$$.

We can guess that the output is a rational number mod $$$998244353$$$.

By brute forcing all fractions with numerator and denominator $$$\le 3000$$$ (sth that can run in $$$\tilde{O}(n^2)$$$) tells us that the answer to sample $$$2$$$ is $$$\frac{691}{2730}$$$.

Google $$$\frac{691}{2730}$$$.

Output the absolute value of the $$$(2n-2)$$$-th Bernoulli number (using the indexing from Wikipedia). For $$$n \le 15$$$ you can just hard code it. For $$$n \le 2 \cdot 10^5$$$, you can compute it using its EGF (Exponential Generating Function) in $$$O(n\log n)$$$.

While Eulerian paths and the input format might suggest this is a graph problem, the ? in the title suggests that we should think twice (deja vu from season 1?).

// is very weird. Why is it in the title?

It's a geometry problem.

Output seems to contain two pairs of triangles.

You can reinterpret the title as "Parallel Euler Lines".

Find two pairs of triangles from the set of points with parallel Euler lines.

You are playing a two-player game on a $$$n$$$ by $$$n$$$ grid where on each turn, you can color a $$$1 \times k$$$ or $$$k \times 1$$$ subrectangle. The game ends when all cells are colored, and you lose if you can’t make a move. From here, it is a normal CP problem.

Symmetric Strategy

You can always win as P1 when $$$n$$$ or $$$m$$$ is odd by starting with the center row/column and using a symmetric strategy. Similarly, you win as P2 if $$$n, m$$$ are even.

What does $$$a,b,c$$$ stand for in math?

Side lengths of triangle

Given $$$X$$$, find a triangle with integer side lengths with area $$$\sqrt{X}$$$.

The score here only depends on $$$|s|$$$. $$$|s| = 81$$$ for full points. You get higher partials if the length is a square, otherwise the partial credit roughly indicates how close you are to the nearest square number below. Trying strings of different lengths will eventually lead you to this.

You need to submit a string $$$s$$$ of perfect square length to get nonzero points in B or C. This suggests that we need to use the perfect square condition somehow.

Since the length is a square number, we can arrange the digits in a square grid (for length=81, 9x9 grid).

This criterion counts the number of rows and columns which contain all distinct digits and gives a higher score if you have more such rows/columns. It also gives a higher score for larger grids. To get full score, you need a $$$9 \times 9$$$ Latin square.

This criteria is a bit more complicated. Having a lot of composite numbers will tend to give you a higher score.

First, read each row and column as a (9-digit) number. Your score depends on the average number of (not necessarily distinct) prime factors is (higher if more prime factors). To get a full score, you need an average of $$$\ge \frac{100}{9}$$$.

First, we need to do some basic counting. For a set of vertices $$$S$$$, let $$$t(S)$$$ denote the minimal subtree that contains all vertices of $$$S$$$.

Fix $$$k$$$. It is clear that we have to somehow double count the sum of $$$f(S)$$$. If we look at a vertex $$$v$$$, it is not that easy to count the number of sets $$$S$$$ of size $$$k$$$ such that $$$t(S)$$$ contains $$$v$$$. However, if we look at an edge $$$e$$$, then it is easy to see that $$$t(S)$$$ contains $$$e$$$ if and only if all elements of $$$S$$$ is in the same connected component formed by deleting the edge $$$e$$$ from the tree. In other words, if the edge $$$e$$$ splits the tree into two components of size (number of vertices) $$$a$$$ and $$$n-a$$$ respectively, then the number of $$$S$$$ with $$$|S|=k$$$ and $$$e \in t(S)$$$ is exactly $$$\binom{n}{k} - \binom{a}{k} - \binom{n-a}{k}$$$. Thus, the sum of $$$f(S)$$$ over all subsets $$$S$$$ is just $$$\binom{n}{k}$$$ (since number of vertices = number of edges + 1) plus the sum of $$$\binom{n}{k} - \binom{a}{k} - \binom{n-a}{k}$$$ over all edges $$$e$$$.

We can find the value of $$$a$$$ for each edge $$$e$$$ by simple dfs. Suppose we have computed the value of $$$a$$$ for each edge $$$e$$$, say $$$a_{1},a_{2},...,a_{n-1}$$$. If we can compute $$$\displaystyle\sum_{i=1}^{n-1}\binom{a_{i}}{k}$$$ fast for each $$$1 \le k \le n$$$, then we are done, so we will focus on this task.

Let $$$c_{i}$$$ denote the number of times $$$i$$$ appear in the array $$$a_{1},a_{2},...,a_{n-1}$$$. Hence, we need to find $$$ans_{k} = \displaystyle\sum_{i=0}^{n}c_{i}\binom{i}{k}$$$.

It is almost customary to try and write binomial coefficients in terms of factorials and look for simplifications.

$$$ans_{k} = \displaystyle\sum_{i=0}^{n}c_{i}\binom{i}{k} = \displaystyle\sum_{i=0}^{n}c_{i}\frac{i!}{k!(i-k)!} = \frac{1}{k!}\displaystyle\sum_{i=0}^{n}c_{i}(i!) \cdot \frac{1}{(i-k)!}$$$

Obviously we only need to know how to compute the sum $$$\displaystyle\sum_{i=0}^{n}c_{i}(i!) \cdot \frac{1}{(i-k)!}$$$ fast for all $$$k$$$. If you have studied about generating functions, you see that our sum looks very much like the convolution of two sequences. If we can write our summand in terms of $$$f(k)g(i-k)$$$, then we can define $$$F(x) = \displaystyle\sum_{i \ge 0}f(i)x^{i}$$$ and $$$G(x) = \displaystyle\sum_{i \ge 0}g(i)x^{i}$$$ and read off $$$\displaystyle\sum_{i \ge 0}f(i)g(k-i)$$$ from the coefficients of $$$F(x)G(x)$$$.

Clearly, we can set $$$f(i)=c_{i}(i!)$$$. We want to set $$$g(i)=\frac{1}{(-i)!}$$$. However, you might worry that $$$(-i)!$$$ is not defined for $$$i > 0$$$. It's ok, since we can always shift our sequence can set $$$g(M+i)=\frac{1}{(-i)!}$$$ instead for some large integer $$$M$$$. Then, we have $$$g(i)=\frac{1}{(M-i)!}$$$, which we can now safely compute. Instead of reading off the coefficient of $$$x^{k}$$$ from $$$F(x)G(x)$$$, we read off the coefficient of $$$x^{k+M}$$$.

This problem was created $$$6$$$ years ago and has a 3100 rating on CF, but if you know about generating functions in 2020 (in the era where polynomial operations template is available) it is almost trivial.

Firstly, the obvious step is to let $$$f_{s}$$$ denote the answer for $$$s$$$ and $$$F(x)$$$ as the OGF of $$$f$$$. Let us derive a recurrence for $$$f_{s}$$$.

Clearly, $$$f_{0} = 0$$$. For $$$s>0$$$, we iterate over all possible weights for the root, then iterate over all possible weights of the left subtree, giving the recurrence $$$f_{s} = \displaystyle\sum_{c \in C}\displaystyle\sum_{i \ge 0}f_{i}f_{s-c-i}$$$ (for convenience, $$$f_{i}=0$$$ for $$$i<0$$$).

Speaking in generating functions, this is just the equation

$$$\displaystyle\sum_{n \ge 1}f_{n}x^{n} = \displaystyle\sum_{n \ge 1}\displaystyle\sum_{c \in C}x^{c}\displaystyle\sum_{i \ge 0}f_{i}x^{i}f_{s-c-i}x^{s-c-i}$$$.

This motivates us to define $$$F(x) = \displaystyle\sum_{n \ge 0}f_{n}x^{n}$$$ and $$$C(x) = \displaystyle\sum_{c \in C}x^{c}$$$. Thus, $$$F(x) - 1 = C(x)F(x)^2$$$ (you should be able to deduce this directly from the equation above).

Our goal is to find $$$F(x)$$$, so we solve for $$$F(x)$$$ using the quadratic formula (remember what we did to obtain the OGF of Catalan numbers?) to obtain $$$F(x) = \frac{1 \pm \sqrt{1-4C(x)}}{2C(x)}$$$.

Similar to the case with Catalan numbers, we choose the negative sign since otherwise as $$$x \rightarrow 0$$$, the numerator goes to $$$2$$$ while the denominator goes to $$$0$$$, but we should converge to a finite limit since $$$F(0)=1$$$. Thus, $$$F(x) = \frac{1 - \sqrt{1-4C(x)}}{2C(x)}$$$.

In the language of generating functions, we are pretty much done, but there is one minor implementation detail here. $$$C(x)$$$ has constant term $$$0$$$, and thus we cannot take the reciprocal directly. However, we just need to "rationalize" the numerator and rewrite our function as

$$$F(x) = \frac{1 - \sqrt{1-4C(x)}}{2C(x)} = \frac{(1 - \sqrt{1-4C(x)})(1 + \sqrt{1-4C(x)})}{2C(x)(1 + \sqrt{1-4C(x)})} = \frac{4C(x)}{2C(x)(1 + \sqrt{1-4C(x)})} = \frac{2}{1 + \sqrt{1-4C(x)}}$$$. You can verify that the constant term of the denominator is nonzero, so we can calculate the first $$$m$$$ terms of the series with usual polynomial operations and solve the problem in $$$O(m\log m)$$$.

To be honest, it is unfair to say that the problem is easy because the main difficulty was to compute the square root and reciprocal of a power series. However, in the present times, these algorithms can be easily found online and in online rounds you can just use templates to perform these operations and thus I would consider this problem to be straightforward.

There is a simple $$$O\left(\frac{n^2}{k}\right)$$$ dp solution (try to find it), but that won't be sufficient for this problem. We need to use generating functions. Let $$$f(n)$$$ denote the answer ($$$k$$$ is fixed).

For a permutation $$$p = p_1, p_2, ..., p_n$$$, define the descent set as the set of integers $$$1 \le i \le n-1$$$ such that $$$p_{i}>p_{i+1}$$$. We are looking for the number of length $$$n$$$ permutations with descent set equal to $$$S = \{k,2k,3k,4k,...\}$$$.

The idea is that counting "exact" conditions is hard, but counting "at least" conditions is easier. Let $$$f(S)$$$ denote the number of permutations with descent set exactly equal to $$$S$$$ and $$$g(S)$$$ denote the number of permutations with descent set that is a subset of $$$S$$$. Clearly, we have $$$g(S) = \displaystyle\sum_{T \subseteq S}f(T)$$$. By the Principle of Inclusion-Exclusion, we have $$$f(S) = \displaystyle\sum_{T \subseteq S}(-1)^{|S| - |T|}g(T)$$$ (see here, this is a commonly used form of PIE).

$$$g(S)$$$ is much easier to count. Suppose $$$S = \{s_1,s_2,...,s_m\}$$$ where $$$s_{i}<s_{i+1}$$$ for all $$$1 \le i \le k-1$$$. This means that $$$p_{j}<p_{j+1}$$$ must hold for all $$$j$$$ that is not in $$$S$$$. A better way to visualize this is that we divide our permutation into several increasing blocks, the first block has size $$$s_{1}$$$, the second block has size $$$s_{2}-s_{1}$$$ and so on. The last block has size $$$n-s_{m}$$$. The only restrictions is that the elements in each block must be in increasing order. Hence, it is clear that the probability that a random permutation satisfies this condition is just $$$\frac{1}{s_{1}!(s_{2}-s_{1})!(s_{3}-s_{2})!...(n-s_{m})!}$$$ (multiply the probability that each block is ordered correctly). Hence, $$$g(S) = \frac{n!}{s_{1}!(s_{2}-s_{1})!(s_{3}-s_{2})!...(n-s_{m})!}$$$.

Let's substitute this back to our equation. For simplicity, let $$$D = {k,2k,3k,...,} \cap [n-1]$$$. Our problem reduces to finding $$$f(D)$$$.

Any subset $$$T = \{s_1,s_2,...,s_{m-1}\}$$$ (elements sorted in increasing order) of $$$D$$$ can be described by a sequence of positive integers $$$b_{1},b_{2},...,b_{m}$$$ where $$$b_{i} = s_{i}-s_{i-1}$$$ ($$$s_{0}=0$$$, $$$s_{m}=n$$$ for simplicity) denote the gap between consecutive elements of $$$T$$$. For example, when $$$T = \{3,9,12\}$$$ and $$$n=13$$$, we can describe it with the sequence $$$3,6,3,1$$$. Note that $$$k$$$ divides $$$b_{1},b_{2},...,b_{m-1}$$$ and $$$b_{m}$$$ has the same remainder as $$$n$$$ mod $$$k$$$ (call such sequences good).

This allows us to simplify the formula of $$$g(T)$$$ to $$$\frac{n!}{b_{1}!b_{2}!...b_{m}!}$$$.

Hence,

$$$f(D) = \displaystyle\sum_{T \subseteq D}(-1)^{|D| - |T|}g(T)$$$

$$$= \displaystyle\sum_{\sum b_{i} = n, b_{i} \text{ good}}(-1)^{|D|-(m-1)}\frac{n!}{b_{1}!b_{2}!...b_{m}!}$$$

For simplicity, let $$$n = kq+r$$$ where $$$1 \le r \le k$$$. Since we only care about the answer for $$$k \mid n - r$$$, we look at the EGF on these terms only, i.e. consider $$$F(x) = \displaystyle\sum_{q \ge 0}\frac{f(kq+r)}{q!}x^{q}$$$.

We have

$$$F(x) = \displaystyle\sum_{q \ge 0}\frac{f(kq+r)}{(kq+r)!}x^{kq+r}$$$

$$$= \displaystyle\sum_{q \ge 0}\displaystyle\sum_{m \ge 1}\displaystyle\sum_{\sum b_{i} = kq+r, b_{i} \text{ good}}(-1)^{q-m+1}\frac{(kq+r)!}{b_{1}!b_{2}!...b_{m}!} \cdot \frac{x^{kq+r}}{(kq+r)!}$$$

$$$= \displaystyle\sum_{m \ge 1}\displaystyle\sum_{q \ge 0}\displaystyle\sum_{\sum b_{i} = kq+r, b_{i} \text{ good}}(-1)^{q-m+1}\frac{x^{b_{1}} \cdot x^{b_{2}} \cdot ... \cdot x^{b_{m}}}{b_{1}!b_{2}!...b_{m}!}$$$

$$$= \displaystyle\sum_{m \ge 1}\displaystyle\left(\displaystyle\sum_{i \ge 1}\frac{(-1)^{i-1} \cdot x^{ki}}{(ki)!}\right)^{m-1} \cdot \left(\displaystyle\sum_{i \ge 0}\frac{(-1)^{i} \cdot x^{ki+r}}{(ki+r)!}\right)$$$

Take a moment to digest the last identity (try expanding the brackets). The idea is that we are able to make $$$b_{1},b_{2},...,b_{m}$$$ independent of each other and simplify our sum into the convolution (or power) of several polynomials.

Continuing our simplifications, we have

$$$= \left(\displaystyle\sum_{i \ge 0}\frac{(-1)^{i} \cdot x^{ki+r}}{(ki+r)!}\right) \cdot \displaystyle\sum_{m \ge 1}\displaystyle\left(\displaystyle\sum_{i \ge 1}\frac{(-1)^{i-1} \cdot x^{ki}}{(ki)!}\right)^{m-1}$$$

$$$ = \left(\displaystyle\sum_{i \ge 0}\frac{(-1)^{i} \cdot x^{ki+r}}{(ki+r)!}\right) \cdot \frac{1}{1 - \left(\displaystyle\sum_{i \ge 1}\dfrac{(-1)^{i-1} \cdot x^{ki}}{(ki)!}\right)}$$$

$$$ = \frac{\displaystyle\sum_{i \ge 0}\frac{(-1)^{i} \cdot x^{ki+r}}{(ki+r)!}}{\displaystyle\sum_{i \ge 0}\dfrac{(-1)^{i} \cdot x^{ki}}{(ki)!}}$$$

We can compute the first $$$n$$$ terms of the last expression in $$$O(n\log n)$$$ time, and we're done.

I think exponential generating functions are especially good at handling sums involving multinomial coefficients as you can separate the factorials into different polynomials and reduce it to convolution of EGFs.

It is hard to compute when a game ends, and it is also hard to compute the probability that the game ends in exactly $$$k$$$ moves. However, it is relatively easier to compute the probability that the state of switches are all off after exactly $$$k$$$ moves. Let $$$a(k)$$$ denote the probability that all switches are off after exactly $$$k$$$ moves, and $$$A(x)$$$ be the EGF (we'll see why we choose EGF soon) of $$$a$$$. How to find $$$A(x)$$$?

Suppose we fix a sequence $$$a_{1}, a_{2}, ..., a_{n}$$$ such that $$$\sum a_{i} = k$$$ where $$$a_{i}$$$ denotes the number of flips of the $$$i$$$-th switch. The probability of achieving this sequence is $$$\frac{k!}{a_{1}!a_{2}!...a_{n}!}\left(\frac{p_1}{S}\right)^{a_{1}}\left(\frac{p_2}{S}\right)^{a_{2}}...\left(\frac{p_n}{S}\right)^{a_{n}}$$$ (the product of the number of sequences of switch flips such that switch $$$i$$$ is flipped exactly $$$a_{i}$$$ times and the probability the sequence of switch flips occur). Let $$$q_{i} = \frac{p_i}{S}$$$ for convenience. Hence, $$$\frac{a(k)}{k!} = \frac{q_{1}^{a_{1}}}{a_{1}!} \cdot \frac{q_{2}^{a_{2}}}{a_{2}!} \cdot ... \cdot \frac{q_{n}^{a_{n}}}{a_{n}!}$$$.

Let's assume that all switches are off at the initial state for now. Then, the EGF is $$$A(x) = \displaystyle\prod_{i=1}^{n}\left(\frac{(q_{i}x)^{0}}{0!} + \frac{(q_{i}x)^{2}}{2!} + ...\right)$$$, since we require each switch to be flipped an even number of times. In the general case, our EGF is similar, but some switches are required to be flipped an odd number of times instead. Motivated by this special case, we let $$$E_{i}(x) = \displaystyle\sum_{j \text{ even}}\frac{(q_{i}x)^{j}}{j!}$$$ and $$$O_{i}(x) = \displaystyle\sum_{j \text{ odd}}\frac{(q_{i}x)^{j}}{j!}$$$. Then, if $$$s_{i}=1$$$ (the switch is initially on), we choose $$$O_{i}(x)$$$ as the $$$i$$$-th term of our product (in the formula for $$$A(x)$$$), while if $$$s_{i}=0$$$, we choose $$$E_{i}(x)$$$ as the $$$i$$$-th term of our product.

There's an even more compact way to write this "observation". Recall that $$$\frac{e^{x}+e^{-x}}{2} = \cosh x = 1 + \frac{x^2}{2!} + \frac{x^4}{4!} + ...$$$. We can use a similar idea here. To express $$$E_{i}(x)$$$, we can try to add or subtract $$$\exp(q_{i}x)$$$ with $$$\exp(-q_{i}x)$$$ to "filter out" the even or odd power terms (we will see a generalization of this trick called the roots of unity filter later in the post). Verify that $$$E_{i}(x) = \frac{\exp(q_{i}x) + \exp(-q_{i}x)}{2}$$$ and $$$O_{i}(x) = \frac{\exp(q_{i}x) - \exp(-q_{i}x)}{2}$$$. Thus, we can even write this as $$$\frac{\exp(q_{i}x) + (-1)^{s_{i}}\exp(-q_{i}x)}{2}$$$.

To summarize, $$$A(x) = \displaystyle\prod_{i=1}^{n}\frac{[\exp(q_{i}x) + (-1)^{s_{i}}\exp(-q_{i}x)]}{2}$$$.

Ok, so we can find $$$a(k)$$$. Let $$$c(k)$$$ denote the probability that the game ends (all switches are off for the first time) after exactly $$$k$$$ moves. How can we relate $$$c(k)$$$ and $$$a(k)$$$? Here is the trick. Consider any sequence of $$$k$$$ moves resulting in all switches being off. After $$$i$$$ moves (possibly $$$i=k$$$), the switches are all off for the first time and the game ends. For the remaining $$$k-i$$$ moves, we need to flip each switch an even number of times. Thus, if we let $$$b(k)$$$ denote the probability that we flip each switch an even number of times after exactly $$$k$$$ moves, then $$$a(k) = \displaystyle\sum_{i=0}^{k}c(i)b(k-i)$$$. Does this look familiar? Yes, it is just normal convolution (but on OGFs instead of EGFs).

Firstly, let's find the EGF of $$$b(k)$$$. This is just a special case of $$$a(k)$$$ when $$$s_{i}=0$$$, so we have $$$B(x) = \displaystyle\prod_{i=1}^{n}\frac{[\exp(q_{i}x) + \exp(-q_{i}x)]}{2}$$$.

To relate $$$c(k)$$$ with $$$a(k), b(k)$$$, we need to look at the OGFs of $$$a$$$ and $$$b$$$ (call them $$$A_{o}(x)$$$ and $$$B_{o}(x)$$$), since the recurrence we found is related to the convolution of OGFs. Thus, defining $$$C(x)$$$ and $$$C_{o}(x)$$$ analogously, we have $$$A_{o}(x) = C_{o}(x)B_{o}(x)$$$, so $$$C_{o}(x) = \frac{A_{o}(x)}{B_{o}(x)}$$$.

Our answer is $$$\displaystyle\sum_{k=0}^{\infty}kc(k)$$$ (recall the definition of expected value and $$$c(k)$$$). This is just $$$C_{o}'(1)$$$. By Quotient Rule, this is equivalent to finding $$$\frac{A_{o}'(1)B_{o}(1) - A_{o}(1)B_{o}'(1)}{B_{o}(1)^2}$$$.

Let's see if we can find $$$A_{o}(x)$$$ from $$$A(x)$$$. It is hard to extract the coefficients of $$$A(x)$$$ if we are looking at a product of sums like $$$\displaystyle\prod_{i=1}^{n}\frac{[\exp(q_{i}x) + (-1)^{s_{i}}\exp(-q_{i}x)]}{2}$$$. To turn this into a sum, let's expand the brackets! Note that since $$$q_{i} = \frac{p_{i}}{S}$$$, if we expand the whole thing, we will get a sum where each term is of the form $$$c_{i}\exp(\frac{i}{S}x)$$$ for some $$$-S \le i \le S$$$ (since we are multiplying terms of the form $$$\exp(\frac{j}{S}x)$$$ and the sum of $$$p_{i}$$$ is $$$S$$$). In other words, $$$A(x) = \displaystyle\sum_{-S}^{S}a_{i}\exp\left(\frac{i}{S}x\right)$$$ for some constants $$$a_{i}$$$.

How do we expand the terms quickly? We can use dp! Go through the brackets one by one, and maintain $$$dp[i][j]$$$ which is the coefficient of $$$\exp\left(\frac{j}{S}x\right)$$$ after expanding $$$i$$$ brackets. You can do dp transitions in $$$O(1)$$$ to get the final answer in $$$O(nS)$$$ time.

From $$$A(x) = \displaystyle\sum_{i=-S}^{S}a_{i}\exp\left(\frac{i}{S}x\right)$$$, we can find a formula for $$$A_{o}(x)$$$. Indeed, $$$\exp\left(\frac{i}{S}x\right) = \sum_{j \ge 0}\left(\frac{(\frac{i}{S}x)^j}{j!}\right)$$$, so by removing the $$$j!$$$ terms we get a formula for $$$A_{o}(x)$$$. Specifically,

$$$A_{o}(x) = \displaystyle\sum_{i=-S}^{S}a_{i}\sum_{j \ge 0}\left(\frac{i}{S}x\right)^{j} = \displaystyle\sum_{i=-S}^{S}\frac{a_{i}}{1 - \frac{i}{S}x}$$$.

Similarly, we can derive $$$B_{o}(x) = \displaystyle\sum_{i=-S}^{S}\frac{b_{i}}{1 - \frac{i}{S}x}$$$ for some constants $$$b_{i}$$$.

We want to compute $$$A_{o}(1)$$$, $$$A_{o}'(1)$$$ (and similar for $$$B$$$). However, we run into a slight problem of dividing by zero if we try to do it directly, since $$$1 - \frac{S}{S}(1) = 0$$$. However, since $$$C_{o}(x) = \frac{A_{o}(x)}{B_{o}(x)}$$$, we can multiply both $$$A_{o}(x)$$$ and $$$B_{o}(x)$$$ by the troublesome $$$1-x$$$ term. Formally, let $$$E(x) = (1-x)A_{o}(x)$$$ and $$$F(x) = (1-x)B_{o}(x)$$$. Then, we only need to compute $$$E(1), E'(1), F(1)$$$ and $$$F'(1)$$$ (and as we shall see they are computable and easy to compute).

$$$E(1)$$$ is trivial, since $$$(1-x)A_{o}(x) = \displaystyle\sum_{i=-S}^{S}\frac{a_{i}(1-x)}{1 - \frac{i}{S}x} = \displaystyle\sum_{i=-S}^{S-1}\frac{a_{i}(1-x)}{1 - \frac{i}{S}x} + a_{S}$$$. Since all the terms except $$$a_{S}$$$ when we substitute $$$x=1$$$, $$$E(1) = a_{S}$$$.

$$$E'(1)$$$ is not hard either. By Quotient Rule,

$$$E'(x) = \left[\displaystyle\sum_{i=-S}^{S-1}\frac{a_{i}(1-x)}{1 - \frac{i}{S}x}\right]' = \displaystyle\sum_{i=-S}^{S-1}\frac{-a_{i}\left(1 - \frac{i}{S}x\right) - a_{i}(1-x)\left(-\frac{i}{S}\right)}{\left(1 - \frac{i}{S}x\right)^2}$$$. Substituting $$$x=1$$$ gives us $$$\displaystyle\sum_{i=-S}^{S-1}\frac{-a_{i}}{1 - \frac{i}{S}}$$$, which is easy to compute in $$$O(S)$$$ time.

Thus, we have an easy-to-code $$$O(Sn)$$$ time solution.

This is a very standard problem in math olympiads, but let's see how to solve it in CP. Firstly, we need to know the trick behind the problem. Let's look at the case $$$m=2$$$, $$$r=0$$$, i.e. find $$$\displaystyle\sum_{i \equiv 0 \pmod{2}}\binom{n}{i}$$$.

It is well-known that this is just $$$2^{n-1}$$$, but where does it come from. Let us look at the generating function $$$\displaystyle\sum_{i \ge 0}\binom{n}{i}x^{i}$$$. If we plug in $$$x=1$$$, we get the sum of all binomial coefficients. However, we want to "filter" out the terms with even (or odd) power. What values can we easily substitute? Another easy candidate to try is $$$x=-1$$$, which gives us $$$\displaystyle\sum_{i \ge 0}\binom{n}{i}(-1)^{i}$$$. If we write out the terms, we get the equations:

$$$(1+1)^{n} = \binom{n}{0} + \binom{n}{1} + \binom{n}{2} + \binom{n}{3} + ...$$$

$$$(1-1)^{n} = \binom{n}{0} - \binom{n}{1} + \binom{n}{2} - \binom{n}{3} + ...$$$

Notice that the odd-numbered terms are gone when we add up the equations! Thus, adding up the equations and dividing by $$$2$$$, we obtain $$$\binom{n}{0}+\binom{n}{2}+... = 2^{n-1}$$$.

How to generalize the above method? Let's say we want to find $$$\binom{n}{0} + \binom{n}{3} + \binom{n}{6} + ...$$$.

We can split our sum into three parts:

$$$\binom{n}{0}x^{0} + \binom{n}{3}x^{3} + \binom{n}{6}x^{6} + ...$$$

$$$\binom{n}{1}x^{1} + \binom{n}{4}x^{4} + \binom{n}{7}x^{7} + ...$$$

$$$\binom{n}{2}x^{2} + \binom{n}{5}x^{5} + \binom{n}{8}x^{8} + ...$$$

To leave only the sum of binomial coefficients in each group, we need to substitute $$$x$$$ so that $$$x^{3}=1$$$.

Clearly, $$$x=1$$$ works. What other values of $$$x$$$ work?

The values of $$$x$$$ such that $$$x^{3} = 1$$$ are also called the $$$3$$$rd roots of unity. In this case, $$$x = e^{\frac{2\pi i}{3}}, e^{\frac{4\pi i}{3}}$$$ are the other two roots of unity.

Let $$$S(n,i)$$$ denote the sum of $$$\binom{n}{k}$$$ for all $$$k \equiv i \pmod{3}$$$. (so we want to find $$$S(n,0), S(n,1), S(n,2)$$$).

Let $$$\omega = e^{\frac{2\pi i}{3}}$$$, then $$$1, \omega, \omega^{2}$$$ are the 3rd roots of unity. Note that $$$\omega^{3} = 1$$$ by definition. We substitute $$$x = \omega^{i}$$$ for $$$0 \le i \le 2$$$, to get the following system of equations:

$$$S(n,0)+S(n,1)+S(n,2) = (1+1)^{n}$$$

$$$S(n,0)+\omega S(n,1) + \omega^{2}S(n,2) = (1 + \omega)^{n}$$$

$$$S(n,0) + \omega^{2}S(n,1) + \omega^{4}S(n,2) = (1 + \omega^{2})^{n}$$$

How to solve this system of equations? We need an important identity of the roots of unity, which is $$$1 + \omega^{k} + \omega^{2k} + ... + \omega^{(m-1)k} = 0$$$ whenever $$$k$$$ is not divisible by $$$m$$$ (if $$$\omega^{m}=1$$$). This is because by the geometric series formula, we have $$$1 + \omega^{k} + \omega^{2k} + ... + \omega^{(m-1)k} = \frac{1 - \omega^{mk}}{1 - \omega^{k}} = 0$$$ if $$$\omega^{k} \neq 1$$$.

Hence, in this specific case, $$$1 + \omega + \omega^{2} = 0$$$ and $$$1 + \omega^{2} + \omega^{4} = 0$$$.

Summing all three equations gives us:

$$$3S(n,0) = (1+1)^{n} + (1 + \omega)^{n} + (1 + \omega^{2})^{n}$$$

How to obtain $$$S(n,1)$$$ and $$$S(n,2)$$$ easily? Here is a simple trick. Instead of looking at $$$(x+1)^{n}$$$ only, we look at $$$x(x+1)^{n}$$$ and $$$x^{2}(x+1)^{n}$$$ as well and repeat the same process. Now, all coefficients are shifted properly and thus we can take the sum and divide by $$$3$$$ to find $$$S(n,1)$$$ and $$$S(n,2)$$$ as in the previous case.

In summary, you should get something like:

$$$3S(n,0) = (1+1)^{n} + (1 + \omega)^{n} + (1 + \omega^{2})^{n}$$$

$$$3S(n,2) = (1+1)^{n} + \omega(1 + \omega)^{n} + \omega^{2}(1 + \omega^{2})^{n}$$$

$$$3S(n,1) = (1+1)^{n} + \omega^{2}(1 + \omega)^{n} + \omega(1 + \omega^{2})^{n}$$$ (note that $$$\omega^{4} = \omega$$$)

The remaining problem is how to evaluate the values from this formula. We have a problem because $$$\omega$$$ seems to represent a complex number here ($$$\omega = e^{\frac{2\pi i}{3}}$$$).

The idea is that we do not necessarily need to work in the world of complex numbers. What we require of $$$\omega$$$ is for it to satisfy $$$\omega^{3}=1$$$ and $$$\omega^{k} \neq 1$$$ for $$$1 \le k \le 2$$$. Let us compute our answer as a polynomial in $$$\omega$$$, but modulo $$$\omega^{3}-1$$$ (which means that we will have $$$\omega^{3}=1$$$, $$$\omega^{4}=\omega$$$, etc...). Also, obviously the coefficients will be computed modulo $$$MOD$$$.

For example, $$$(1+2\omega+\omega^{2})(1+\omega+3\omega^{2}) = 3\omega^{4} + 7\omega^{3} + 6\omega^{2} + 3\omega + 1 = 3\omega + 7 + 6\omega^{2} + 3\omega + 1 = 6\omega^{2} + 6\omega + 8$$$.

Hence, at any point of time we will always have a polynomial in $$$\omega$$$ with degree $$$\le 2$$$.

To compute something like $$$(1 + \omega)^{n}$$$, we can use the usual divide-and-conquer method for computing large powers. Multiplication of polynomials can be implemented naively.

We can generalize this method for any $$$m$$$. Let $$$\omega$$$ be the $$$m$$$-th root of unity, so $$$\omega^{m}=1$$$. Let $$$S(n,r)$$$ denote the sum of $$$\binom{n}{k}$$$ over all $$$k \equiv r \pmod{m}$$$.

$$$(1+w)^{n} = \displaystyle\sum_{i \ge 0}\binom{n}{i}w^{i}$$$ for any number $$$w$$$. We want to make the coefficients of the form $$$\binom{n}{jm+r}$$$ to match with the powers $$$w^{0}, w^{m}, w^{2m}, ...$$$ because these will help us obtain the sum $$$S(n,r)$$$. So, it is more helpful to consider the polynomial $$$w^{m-r}(1+w)^{n} = \displaystyle\sum_{i \ge 0}\binom{n}{i}w^{i+m-r}$$$.

Substituting $$$w = 1, \omega, \omega^{2}, ..., \omega^{m-1}$$$ and summing up, we get

$$$\displaystyle\sum_{j=0}^{m-1}w^{j(m-r)}(1+w^{j})^{n} = \displaystyle\sum_{j=0}^{m-1}\displaystyle\sum_{i \ge 0}\binom{n}{i}\omega^{j(i+m-r)} = \displaystyle\sum_{i \ge 0}\binom{n}{i}\displaystyle\sum_{j=0}^{m-1}(\omega^{i+m-r})^{j}$$$.

Recall that $$$1 + w + w^{2} + ... + w^{m-1} = 0$$$ whenever $$$w = \omega, \omega^{2}, ..., \omega^{m-1}$$$ (by the geometric series formula) and $$$= m$$$ whenever $$$w = 1$$$. Note that $$$\omega^{i+m-r} = 1$$$ if and only if $$$i \equiv r \pmod{m}$$$. Hence,

$$$\displaystyle\sum_{i \ge 0}\binom{n}{i}\displaystyle\sum_{j=0}^{m-1}(\omega^{i+m-r})^{j} = m \cdot \displaystyle\sum_{i \equiv r \pmod{m}}\binom{n}{i} = mS(n,r)$$$ (note the multiple of $$$m$$$).

It remains to compute $$$\frac{1}{m}\displaystyle\sum_{j=0}^{m-1}w^{j(m-r)}(1+w^{j})^{n}$$$ modulo $$$MOD$$$. The easiest way to do this is to first calculate the polynomial $$$F(x) = x^{m-r}(1+x)^{n}$$$ modulo $$$x^{m}-1$$$. We can do this via binary exponentiation and multiplying polynomials naively in $$$O(m^{2}\log n)$$$ time (remembering to set $$$x^{m} = 1, x^{m+1} = x, ...$$$ after each multiplication).

After we obtained the polynomial, we are actually done. The sum we want to find is $$$\frac{1}{m}[F(1)+F(\omega)+F(\omega^{2})+...+F(\omega^{m-1})]$$$. Letting $$$F(x) = \displaystyle\sum_{i=0}^{m-1}a_{i}x^{i}$$$, we realize that the desired sum is

$$$\frac{1}{m}\displaystyle\sum_{i=0}^{m-1}a_{i}\displaystyle\sum_{j=0}^{m-1}(\omega^{j})^i = a_{0}$$$ since all the terms with $$$i \ge 1$$$ sum up to $$$0$$$. Hence, we just need to output the constant term of $$$F(x)$$$.

Note: This problem is also solvable in $$$O(m\log m\log n)$$$ if $$$MOD$$$ is a FFT-friendly modulo by using FFT to multiply the polynomials during exponentiation.

This problem is mentioned in jqdai0815's blog, but I will explain it in detail here. A funny story is that I came up with this problem and was trying to solve it one day, but then the next day I saw this problem on his post :) Anyway, I think this is a nice application of generating functions which deserves to be mentioned here.

By now, it should be clear what the first step should be. We want to seperate the terms $$$\frac{1}{x_{i}!}$$$ into different polynomials and convolute them, and the most obvious way to do it is to consider the following product:

$$$\left(1 + \frac{x^d}{d!} + \frac{x^{2d}}{2d!} + ... +\right)\left(1 + \frac{x^d}{d!} + \frac{x^{2d}}{2d!} + ... +\right)...\left(1 + \frac{x^d}{d!} + \frac{x^{2d}}{2d!} + ... +\right)$$$ ($$$k$$$ times).

It is easy to see that the coefficient of $$$x^{n}$$$ is precisely the sum of $$$\frac{1}{x_{1}!x_{2}!...x_{k}!}$$$ over all valid tuples of $$$x_{i}$$$.

Let $$$F(x) = 1 + \frac{x^d}{d!} + \frac{x^{2d}}{2d!} + ...$$$. How do we write $$$F(x)$$$ is a simpler way? For the case $$$d=2$$$, we know that this is just $$$\cosh x = \frac{e^{x} + e^{-x}}{2}$$$. The idea is the same as the previous problem: we want to substitute different values of $$$x$$$ into our "original function" (in this case it is $$$e^{x} = 1 + \frac{x}{1!} + \frac{x^2}{2!} + ...$$$) to "filter out" the powers that are not divisible by $$$d$$$.

Let's try to play with roots of unity again. Let $$$\omega$$$ be the $$$m$$$-th root of unity (so $$$\omega^{m}=1$$$ and $$$1 + \omega + \omega^{2} + ... + \omega^{m-1} = 0$$$). We substitute $$$x$$$ as $$$\omega^{i}x$$$ for $$$0 \le i \le d-1$$$ into $$$e^{x}$$$ and see what happens (note that this is where $$$e^{-x}$$$ comes from for $$$d=2$$$).

$$$e^{x} = 1 + \frac{x}{1!} + \frac{x^2}{2!} + ... + \frac{x^{d}}{d!} + ...$$$

$$$e^{\omega x} = 1 + \frac{\omega x}{1!} + \frac{\omega^{2} x^2}{2!} + ... + \frac{\omega^{d} x^{d}}{d!} + ...$$$

...

$$$e^{\omega^{d-1} x} = 1 + \frac{\omega^{d-1} x}{1!} + \frac{\omega^{2(d-1)} x^2}{2!} + ... + \frac{\omega^{d(d-1)} x^{d}}{d!} + ...$$$

If we sum all these equations, by the identity $$$1 + (\omega^{i}) + (\omega^{i})^2 + ... + (\omega^{i})^{d-1} = 0$$$ for $$$i$$$ not divisible by $$$d$$$, we obtain

$$$e^{x} + e^{\omega x} + e^{\omega^{2} x} + ... + e^{\omega^{d-1} x} = d\left(1 + \frac{x^d}{d!} + \frac{x^{2d}}{2d!} + ... +\right)$$$

Back to our problem, our goal is to find the coefficient of $$$x^n$$$ in the following product:

$$$\frac{1}{d^{k}}\left(e^{x} + e^{\omega x} + e^{\omega^{2} x} + ... + e^{\omega^{d-1} x}\right)^{k}$$$.

The next step requires some knowledge on roots of unity. The good thing about $$$d=4,6$$$ is that $$$\phi(4) = \phi(6) = 2$$$, so the $$$d$$$-th cyclotomic polynomial has degree $$$2$$$ (it is the minimal polynomial of the primitive $$$d$$$-th roots of unity). It can be obtained via the expansion of $$$\displaystyle\prod_{\gcd(i,d)=1}(x-\omega^{i})$$$ (you can find more details on the wikipedia page).

For example, for $$$d=4$$$, we have $$$\omega^{2} + 1 = 0$$$ and for $$$d=6$$$, we have $$$\omega^{2} - \omega + 1 = 0$$$. In both cases, we can always represent $$$\omega^{i}$$$ in the form of $$$a\omega + b$$$ by repeatedly reducing the maximal power using this equation.

Thus, if we let $$$\omega^{i} = a_{i} + b_{i}\omega$$$, our goal is to find $[x^{n}]\frac{1}{d^{k}}\left(\displaystyle\sum_{i=0}^{d-1} e^{(a_i+b_i\omega)x}\right)^{k}.

Notice that the terms of the form $$$e^{ax}$$$ and $$$e^{b\omega x}$$$ are in some sense separated. We make the substitution $$$u = e^{x}$$$ and $$$v = e^{\omega x}$$$ to make things simpler:

$$$[x^{n}]\frac{1}{d^{k}}\left(\displaystyle\sum_{i=0}^{d-1} u^{a_i}v^{b_i}\right)^{k}$$$.

Notice that $$$-1 \le a_{i}, b_{i} \le 1$$$ (you can explicitly write out these values from the definition). If we expand this bivariate polynomial in $$$u$$$ and $$$v$$$, we'll get a sum where all the terms are of the form $$$u^{i}v^{j}$$$ if $$$-k \le i,j \le k$$$ (consider the way to choose the terms in the expansion).

To avoid dealing with negative terms, let us multiply by $$$(uv)^{k}$$$. Let $$$F(u,v) = \displaystyle\sum_{i=0}^{d-1}u^{a_{i}+1}v^{b_{i}+1}$$$. We want to find $$$G(u,v) = F(u,v)^{k}$$$. Note that $$$G(u,v) = \displaystyle\sum_{0 \le i,j \le 2k}g_{i,j}u^{i}v^{j}$$$ for some constants $$$c_{i,j}$$$ (define $$$f_{i,j}$$$ similarly).

While we can do something like 2D FFT, it is probably not going to pass the time limit. The next step is a magical trick mentioned in jqdai0815's article. The idea is that we want to find a recurrence on the coefficients $$$g_{i,j}$$$ so that we can use dp to compute them in $$$O(k^2)$$$ time. Let's differentiate both sides of $$$G(u,v) = F(u,v)^{k}$$$ with respect to $$$u$$$ (or $$$v$$$, it doesn't matter). Let $$$g(u,v)$$$ and $$$f(u,v)$$$ denote the partial derivative of $$$G(u,v)$$$ and $$$F(u,v)$$$ with respect to $$$u$$$. Then, by the Chain Rule,

$$$g(u,v) = kF(u,v)^{k-1}f(u,v)$$$

Noting that $$$F(u,v)^{k-1} = \frac{G(u,v)}{F(u,v)}$$$, we have

$$$F(u,v)g(u,v) = kG(u,v)f(u,v)$$$

Comparing the coefficients of $$$u^{i}v^{j}$$$ on both sides, and noting that the coefficient of $$$u^{i}v^{j}$$$ of $$$f(u,v)$$$ is $$$(i+1)f_{i+1,j}$$$, we obtain the recurrence

$$$\displaystyle\sum_{0 \le i_1 \le i, 0 \le j_1 \le j} (i+1-i_{1})g_{i+1-i_{1},j-j_{1}}f_{i_{1},j_{1}} = k\displaystyle\sum_{0 \le i_1 \le i, 0 \le j_1 \le j} (i_{1}+1)f_{i_{1}+1,j_{1}}g_{i-i_{1},j-j_{1}}$$$

The good thing about this equation is that if we find the values of $$$g_{i,j}$$$ in increasing order of $$$i$$$ followed by increasing order of $$$j$$$, the value of $$$g_{i,j}$$$ is only dependent on previous values! A subtle note is that $$$f_{0,0} = 0$$$ but $$$f_{0,1} \neq 0$$$. Thus, if you only leave the summand with $$$i_{1}=0$$$ and $$$j_{1}=1$$$ on the LHS and throw everything to the RHS, you can compute $$$g_{i,j}$$$ with dp. Note that there are only $$$O(1)$$$ nonzero values of $$$f_{i,j}$$$, so you can actually just iterate over $$$i_{1}, j_{1} \le 2$$$ to compute the dp transitions in $$$O(1)$$$.

The base case is $$$g_{i,j}$$$ with $$$i$$$ or $$$j$$$ equal to $$$0$$$, which can be found by binomial theorem after substituting $$$v=0$$$ (I leave this as an exercise).

To summarize, you can find the expansion of $$$[x^{n}]\frac{1}{d^{k}}\left(\displaystyle\sum_{i=0}^{d-1} u^{a_i}v^{b_i}\right)^{k}$$$ in $$$O(k^2)$$$ time.

How to compute the answer after you reduce it to $$$[x^{n}]\frac{1}{d^{k}}\displaystyle\sum_{-k \le i,j \le k}g_{i+k,j+k}u^{i}v^{j}$$$?. Let's look at each individual term $$$[x^{n}]u^{i}v^{j}$$$. It is exactly equal to $$$[x^n]\exp(x(i + j\omega)) = \frac{1}{n!}(i+j\omega)^{n}$$$. The $$$\frac{1}{n!}$$$ cancels off with the numerator of $$$\frac{n!}{x_{1}!x_{2}!...x_{k}!}$$$. Hence, it is sufficient to compute $$$(i+j\omega)^n$$$ for $$$-k \le i,j \le k$$$.

We can try to use the same trick as the previous problem: Maintaining a polynomial of the form $$$a+b\omega$$$ while doing binary exponentiation. However, it turns out to be a bit too slow (at least my implementation of it). Here, we can exploit the fact that we are computing the answer modulo $$$p = 19491001$$$.

You can find that $$$7$$$ is a primitive root of $$$19491001$$$, and $$$p-1$$$ is divisible by both $$$4$$$ and $$$6$$$. Thus, if we take $$$\omega = 7^{\frac{p-1}{d}}$$$, we will preserve the property that $$$\omega^{d} = 1$$$ and $$$\omega^{i} \neq 1$$$ for $$$1 \le i \le d-1$$$ (equivalently, $$$7^{\frac{p-1}{d}}$$$ is a primitive $$$d$$$-th root of unity in $$$\mathbb{Z}_{p}$$$). Thus, the computations can be done in integers which is faster.

In any case, this gives an $$$O(k^2\log n)$$$ solution with relatively small constants.

Recall that for a permutation $$$p_1, p_2, ..., p_n$$$, an inversion is a pair of indices $$$i<j$$$ such that $$$p_i > p_j$$$.

In my post $$$4$$$ years ago, I mentioned how to solve an easier variant of this problem using a doubling trick to perform dp. Now, let's describe a much simpler and faster solution using generating functions.

Firstly, let us rephrase the problem. Suppose we start with the sequence $$$1$$$ and add the elements $$$2, 3, ..., n$$$ to the permutation one by one. Note that by adding $$$2$$$, we can increase the number of inversions by $$$0$$$ or $$$1$$$ in exactly one way each. Similar, among the $$$i$$$ possible ways to add $$$i$$$ into the sequence, there is exactly one way to increase the number of inversions by $$$0$$$ to $$$i-1$$$ each. Thus, our problem is equivalent to finding the number of sequences $$$a_{0},a_{1},...,a_{n-1}$$$ such that $$$a_{i} \le i$$$ for all $$$i$$$ and $$$\sum a_{i} = k$$$.

Let us rephrase this in the language of generating functions. The idea is that we can "pick" any element from $$$0$$$ to $$$i-1$$$ in the $$$i$$$-th "bracket", so it is natural to consider the function

$$$F(x) = (1)(1+x)(1+x+x^2)(1+x+x^2+x^3)...(1+x+x^2+...+x^{n-1})$$$

The coefficient of $$$x^k$$$ in $$$F(x)$$$ gives us the answer.

Unfortunately, this polynomial is too large to compute directly. Let us rewrite it using the geometric series formula.

$$$F(x) = \frac{1-x}{1-x} \cdot \frac{1-x^2}{1-x} \cdot \frac{1-x^3}{1-x} \cdot ... \cdot \frac{1-x^n}{1-x}$$$

$$$= \frac{(1-x)(1-x^2)...(1-x^n)}{(1-x)^n}$$$

$$$= \prod_{i=1}^{n}(1-x^i) \cdot (1-x)^{-n}$$$

We know how to find the coefficient of $$$[x^i]$$$ in $$$(1-x)^{-n}$$$ for $$$0 \le i \le k$$$ (recall that $$$[x^i]$$$ $$$(1-x)^{-n} = \binom{n+i-1}{i}$$$). Hence, it is sufficient to find $$$[x^j]\prod_{i=1}^{n}(1-x^i)$$$ for all $$$0 \le j \le k$$$.

As $$$n \ge k$$$, we actually only need to compute $$$(1-x)(1-x^2)...(1-x^k)$$$ (as larger terms here don't contribute to the coefficients of lower powers of $$$x$$$). We will present an $$$O(k)$$$ solution which can solve the problem even when $$$k \le 10^6$$$.

The trick is that since the larger terms $$$1-x^{k+1}$$$, $$$1-x^{k+2}$$$, etc doesn't matter in our product, why not just add all of them and consider the infinite product $$$\displaystyle\prod_{n \ge 1}(1 - x^{n})$$$. It turns out that this is a well-known generating function, whose series expansion has a simple form given by the Pentagonal number theorem.

The theorem states that $$$\displaystyle\prod_{n \ge 1}(1 - x^{n}) = 1 + \displaystyle\sum_{i \ge 1}(-1)^{i}\left(x^{\frac{i(3i+1)}{2}} + x^{\frac{i(3i-1)}{2}}\right)$$$. There is a nice bijective proof of this (which you can find on Wikipedia or Enumerative Combinatorics Volume 1), but we only need to use the formula here.

From the series expansion, it is obvious how we can extract the coefficients in $$$O(k)$$$ time (actually even in $$$O(\sqrt{k})$$$ time).

Note that since $$$n$$$ is large, to compute $$$\binom{n+i-1}{i}$$$ for all $$$i$$$, you need to write it in the form $$$\frac{(n+i-1)(n+i-2)...(n)}{i!}$$$ and calculate it in increasing order of $$$i$$$ by multiplying a suitable factor each time $$$i$$$ increases.

Overall, we get an $$$O(k)$$$ time solution, which is more than enough for this problem.

Let's get straight to the generating functions. Define $$$P_{n}(x)$$$ as the OGF for $$$d_{n,k}$$$. $$$d_{1,k}$$$ is just the Fibonacci sequence, so we know that $$$P_{1}(x) = \frac{x}{1-x-x^2}$$$.

How to obtain $$$P_{n}(x)$$$? Thanks to our wonderful "Prefix Sum Trick", we know that multiplying $$$P_{1}(x)$$$ by $$$\frac{1}{1-x}$$$ gives us $$$P_{2}(x)$$$ because $$$d_{2,k}$$$ is just the prefix sum of $$$d_{1,k}$$$. Similarly, we have $$$P_{n}(x) = \frac{1}{(1-x)^{n-1}}P_{1}(x) = \frac{1}{(1-x)^{n-1}} \cdot \frac{x}{1-x-x^2}$$$.

However, now we have some problems, because we need to calculate the coefficient of $$$x^m$$$ in this function with $$$m$$$ up to $$$10^{18}$$$. There is no way we can expand this naively and thus we need to do something clever.

The main annoyance is that we are dealing with a product in the denominator. We know how to compute $$$[x^m]\frac{1}{(1-x)^{n-1}}$$$ and $$$[x^m]\frac{x}{1-x-x^2}$$$ fast (the former is just some binomial coefficient while the latter is just the Fibonacci sequence). However, we don't know how to compute their convolution fast. The trick here is to forcefully separate these two functions by partial fractions. Note that the theory of partial fractions tell us that

$$$\frac{x}{(1-x)^{n-1} \cdot (1-x-x^2)} = \frac{A(x)}{(1-x)^{n-1}} + \frac{B(x)}{1 - x - x^2}$$$

where $$$A(x)$$$ is a polynomial of degree $$$\le n-2$$$ and $$$B(x)$$$ is a polynomial of degree $$$\le 1$$$. If we can find $$$A(x)$$$ and $$$B(x)$$$, then our problem will be much easier to solve. However, $$$A(x)$$$ and $$$B(x)$$$ is hard to find explicitly on paper (though you can guess $$$B(x)$$$ by some pattern from small cases). How do we proceed?

Here is a painless way to do it. Firstly, we clear the denominators to obtain the identity

$$$A(x)(1-x-x^2) + B(x)(1-x)^{n-1} = x$$$.

Since this is an identity, it remains true for any value of $$$x$$$ we substitute! What convenient values of $$$x$$$ can we substitute? We want to make either $$$1-x-x^2$$$ or $$$1-x$$$ equal to $$$0$$$ to leave us with only one polynomial to deal with. Substituting $$$x=1$$$ doesn't tell us too much since $$$A(x)$$$ has degree $$$\le n-2$$$ and we can't determine it with only $$$1$$$ point of information. However, what if we let $$$1-x-x^2=0$$$? We can solve the quadratic equation to get two roots $$$a+b\sqrt{5}$$$ and $$$a-b\sqrt{5}$$$ for some constants $$$a,b$$$. Substituting $$$x=a \pm b\sqrt{5}$$$, we have the nice pair of equations

$$$B(a \pm b\sqrt{5})(1 - (a \pm b\sqrt{5}))^{n-1} = a \pm b\sqrt{5}$$$.

Since $$$B(x)$$$ is linear, if we let $$$B(x) = mx+c$$$ we can solve for the coefficients of $$$B$$$ using these $$$2$$$ simultaneous equations! An implementation detail is that since we are dealing with $$$\sqrt{5}$$$ here, it is helpful to store the numbers as a pair $$$(a,b)$$$ which denotes $$$a+b\sqrt{5}$$$ and do arithmetic on these pairs. The value $$$(1 - (a \pm b\sqrt{5}))^{n-1}$$$ can be found via binary exponentiation in $$$O(\log n)$$$ time (or even naive multiplication works here).

In any case, after some work you can find $$$B(x)$$$. How do we find $$$A(x)$$$? Just refer to the identity to obtain $$$A(x) = \frac{x - B(x)(1-x)^{n-1}}{1-x-x^2}$$$, which we can compute in $$$O(n\log n)$$$ time with one FFT (note that you don't even need to compute the reciprocal of a polynomial, as you can just divide using long division since $$$A(x)$$$ is a polynomial).

It remains to compute $$$[x^m]\frac{A(x)}{(1-x)^{n-1}} + \frac{B(x)}{1 - x - x^2}$$$, which is a significantly easier task. For the second term, we can partition it into $$$\frac{Mx}{1-x-x^2}$$$ and $$$\frac{C}{1-x-x^2}$$$ and note that both are generating functions for the Fibonacci numbers and thus we can just compute the coefficient of $$$x^m$$$ as a Fibonacci number in $$$O(\log m)$$$ time (using any divide-and-conquer method you like).

For the first term, we have $$$A(x)(1-x)^{-(n-1)}$$$ and we know how to compute both $$$[x^i]A(x)$$$ and $$$[x^j]$$$ $$$(1-x)^{-(n-1)}$$$ (it is a binomial coefficient) for all $$$i,j$$$. Since $$$A(x)$$$ is of degree $$$\le n-2$$$, we can iterate through all $$$i$$$ and compute the sum of $$$[x^i]A(x) \cdot [x^{m-i}]$$$ $$$(1-x)^{-(n-1)}$$$ (the latter requires large binomial coefficients which should be computed in a similar manner as the previous problem).

Thus, we obtain an $$$O(n\log n + \log m)$$$ solution.

I believe you can generalize this solution to solve other recurrences of a similar form.

Let us derive a simple recurrence first. Let $$$E(i)$$$ denote the expected number of moves to reach $$$0$$$ from $$$i$$$. Clearly, $$$E(0) = 0$$$, and we have $$$E(i) = p_{1}E(i-1) + p_{2}E(i-2) + ... + p_{n}E(i-n) + 1$$$. We use the convention that $$$E(-r) = E(M-r)$$$.

First, we look at the high-level idea of the solution. The idea is that for $$$i \ge n$$$, we can always write $$$E(i)$$$ in terms of $$$c_{0}E(0) + c_{1}E(1) + c_{2}E(2) + ... + c_{n-1}E(n-1) + C$$$ for some constants $$$c_{i}$$$ and $$$C$$$ by repeatedly using the recurrence relation. In particular, $$$E(M+1) = E(1)$$$, $$$E(M+2) = E(2)$$$, ..., $$$E(M+n-1)=E(n-1)$$$ can all be represented in this form in a non-trivial manner using the recurrence (note that the recurrence doesn't hold for multiples of $$$M$$$, but $$$E(M+i)$$$ can still be represented non-trivially in this form for $$$1 \le i \le n-1$$$).

Hence, moving everything unknown to one side and the constants to the other, we have $$$n-1$$$ nontrivial equations of the form $$$c_{i,1}E(1) + c_{i,2}E(2) + ... + c_{i,n-1}E(n-1) = C_{i}$$$ for $$$1 \le i \le n-1$$$. We can solve this system of equations using Gaussian Elimination in $$$O(n^3)$$$ time.

Once we obtained the values of $$$E(1), E(2), ..., E(n-1)$$$, we just need to represent $$$E(k)$$$ in terms of $$$E(0), E(1), .., E(n-1)$$$ and a constant and we are done.

The remaining problem here is how do we get the representation of $$$E(m)$$$ in terms of $$$E(0), E(1), .., E(n-1)$$$ and a constant $$$C$$$ for any $$$m$$$ in an efficient manner.

Let's assume for the time being that $$$E(M), E(2M), E(3M), ...$$$ also satisfy the linear recurrence $$$E(i) = p_{1}E(i-1) + p_{2}E(i-2) + ... + p_{n}E(i-n) + 1$$$. Thus, the values of $$$E(M)$$$, $$$E(M+1)$$$, ... might not be what we want now but we will deal with this issue later.

Now, we use the same trick as in this blog. For a polynomial $$$f(x) = a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^3+...+a_{k}x^{k}$$$ define its valuation $$$val(f)$$$ as $$$a_{0}+a_{1}E(1)+...+a_{k}E(k)$$$. Since $$$E(i) - p_{1}E_(i-1) - p_{2}E(i-2) - ... - p_{n}E(i-n) = 1$$$ for all $$$i \ge n$$$, we have $$$val(x^{i} - p_{1}x^{i-1} - p_{2}x^{i-2} - ... - p_{n}x^{i-n}) = 1$$$ for all $$$i \ge n$$$. Let $$$P(x) = x^{n} - p_{1}x^{n-1} - p_{2}x^{n-2} - ... - p_{0}x^{0}$$$ for convenience. Then, $$$val(x^{k}P(x)) = 1$$$ for all $$$k \ge 0$$$. Since $$$val$$$ is additive, for any polynomial $$$Q(x)$$$, we have $$$val(Q(x)P(x)) = Q(1)$$$ (sum of coefficients, since $$$x^{k}$$$ corresponds to $$$1$$$. If this is unclear, try writing $$$Q(x)$$$ as $$$q_{0}x^{0} + q_{1}x^{1} + ...$$$).

Our goal is to find $$$val(x^{m})$$$ for some integer $$$m$$$. By the division algorithm, we can write $$$x^{m} = P(x)Q(x)+R(x)$$$ where $$$R(x)$$$ is a polynomial of degree $$$\le n-1$$$. Hence, $$$val(x^{m}) = val(P(x)Q(x)+R(x)) = val(P(x)Q(x))+val(R(x)) = Q(1)+val(R(x))$$$. Notice that $$$val(R(x))$$$ is already a linear combination of $$$E(1), E(2), ..., E(n-1)$$$, while $$$Q(1)$$$ is a constant. Thus, if we can somehow find $$$Q(x)$$$ and $$$R(x)$$$, we can represent $$$E(m)$$$ as a linear combination of $$$E(1)$$$, $$$E(2)$$$, ..., $$$E(n-1)$$$ and a constant.

To find $$$Q(1)$$$ and $$$R(x)$$$, we can use a divide-and-conquer algorithm similar to binary exponentiation. Consider the function $$$solve(m)$$$ that returns a pair denoting $$$R(x)$$$ and $$$Q(1)$$$ (a polynomial and a constant). If $$$m$$$ is even, let $$$R_{1}(x)$$$, $$$Q_{1}(1)$$$ be the return value of $$$solve\left(\frac{m}{2}\right)$$$. Let $$$x^{\frac{m}{2}} = P(x)Q_{1}(x) + R_{1}(x)$$$. We have

$$$x^{m} = (P(x)Q_{1}(x) + R_{1}(x))(P(x)Q_{1}(x) + R_{1}(x)) = P(x)^{2}Q_{1}(x)^{2} + P(x)[2R_{1}(x)Q_{1}(x)] + R_{1}(x)^2 = P(x)[P(x)Q_{1}(x)^{2} + 2R_{1}(x)Q_{1}(x)] + R_{1}(x)^2$$$.

Let $$$R_{1}(x)^{2} = Q_{2}(x)P(x) + R_{2}(x)$$$ (just do long division). It follows from the equation above that we can take $$$R(x) = R_{2}(x)$$$ and $$$Q(1) = P(1)Q_{1}(1)^{2} + 2R_{1}(1)Q_{1}(1) + Q_{2}(1)$$$. The case where $$$m$$$ is odd is similar. Thus, we can compute $$$val(x^{m})$$$ in $$$O(n^{2}\log m)$$$ time.

Also, note that if we have computed $$$val(x^{m})$$$, we can compute $$$val(x^{m+1})$$$ in $$$O(n^{2})$$$ time using the same trick. Thus, we can use this method to compute a representation of $$$E(M-n+1)$$$, $$$E(M-n+2)$$$, ..., $$$E(M-1)$$$ in terms of $$$E(1),E(2),...,E(n-1)$$$ and a constant in $$$O(n^{2}\log m + n^{3})$$$ time. The reason we cannot compute $$$E(M+1)$$$, $$$E(M+2)$$$, ..., $$$E(M+n-1)$$$ directly using $$$val$$$ is because the recurrence doesn't hold for $$$E(M)$$$ as stated before. However, once we have the representation for $$$E(M-n+1)$$$, $$$E(M-n+2)$$$, ..., $$$E(M-1)$$$, we can now plug those into the original recurrence $$$E(i) = p_{1}E(i-1) + p_{2}E(i-2) + ... + p_{n}E(i-n) + 1$$$ and obtain the representations for $$$E(M+1)$$$, $$$E(M+2)$$$, ..., $$$E(M+n-1)$$$ in $$$O(n^{2})$$$ time each.

This gives us a $$$O(n^{2}(n + \log m))$$$ solution.

If you have heard of Cayley's Formula, you know that the answer is $$$n^{n-2}$$$.

Let us count the number of rooted trees on $$$n$$$ vertices. Call this number $$$t(n)$$$. If we remove the root from a rooted tree, we get a collection of rooted subtrees. This allows us to get the recurrence

$$$t(n+1) = (n+1)\displaystyle\sum_{k \ge 0}\sum_{i_{1}+i_{2}+...+i_{k}=n, i_{j} \ge 1}\frac{n!}{i_{1}!i_{2}!...i_{k}!}t(i_{1})t(i_{2})...t(i_{k}) \cdot \frac{1}{k!}$$$, as we have $$$n+1$$$ ways to choose the root, $$$\frac{n!}{i_{1}!i_{2}!...i_{k}!}$$$ ways to assign the non-root nodes to a subtree (we divide by $$$k!$$$ because each set of subtrees is counted $$$k!$$$ times for each permutation), and $$$t(i_{1})t(i_{2})...t(i_{k})$$$ is the number of ways to form each rooted subtree.

Rearranging and multiplying by $$$x^{n+1}$$$, we obtain

$$$\frac{t(n+1)}{(n+1)!}x^{n+1} = x\displaystyle\sum_{k \ge 0}\frac{1}{k!} \cdot \sum_{i_{1}+i_{2}+...+i_{k}=n, i_{j} \ge 1}\frac{t(i_{1})x^{i_1}}{i_{1}!} \cdot \frac{t(i_{2})x^{i_2}}{i_{2}!} \cdot ... \cdot \frac{t(i_{k})x^{i_k}}{i_{k}!}$$$.

Hence, letting $$$T(x)$$$ be the EGF of $$$t(n)$$$ (and define $$$t(0)=0$$$ for simplicity), we have

$$$T(x) = \displaystyle\sum_{n \ge 0}\frac{t(n+1)}{(n+1)!}x^{n+1} = x\displaystyle\sum_{k \ge 0}\frac{1}{k!} \cdot \displaystyle\sum_{n \ge 0}\sum_{i_{1}+i_{2}+...+i_{k}=n, i_{j} \ge 1}\frac{t(i_{1})x^{i_1}}{i_{1}!} \cdot \frac{t(i_{2})x^{i_2}}{i_{2}!} \cdot ... \cdot \frac{t(i_{k})x^{i_k}}{i_{k}!}$$$

$$$= x\displaystyle\sum_{k \ge 0} \frac{T(x)^{k}}{k!}$$$ (verify that we get the previous line by expanding this)

$$$= xe^{T(x)}$$$

Hence, we have the functional equation $$$T(x) = xe^{T(x)}$$$. It is not easy to solve this equation directly, however. However, we can see that we have a function in $$$T(x)$$$ which is equal to $$$x$$$, which motivates us to write

$$$T(x)e^{-T(x)} = x$$$

and let $$$F(x) = xe^{-x}$$$, $$$G(x) = T(x)$$$ in Lagrange Inversion Formula, we obtain

$$$[x^{n}]T(x) = \frac{1}{n}[x^{-1}]\frac{1}{(xe^{-x})^n} = \frac{1}{n}[x^{-1}]x^{-n}e^{nx} = \frac{1}{n}[x^{n-1}]e^{nx} = \frac{1}{n} \cdot \frac{n^{n-1}}{(n-1)!} = \frac{n^{n-1}}{n!}$$$.

Thus, $$$t(n) = n^{n-1}$$$.

Finally, to count the number of unrooted labelled trees, simply divide $$$t(n)$$$ by $$$n$$$ as each unrooted tree is counted $$$n$$$ times in $$$t(n)$$$ by the $$$n$$$ choices of root. Hence, the answer is $$$n^{n-2}$$$.

Let's warmup with an easier problem. Suppose we want to count the number of labelled connected graphs on $$$n$$$ vertices. There are different ways to compute this, but let's show a method using EGFs as it will be useful later. Let $$$C(x)$$$ be the EGF of the number of connected labelled graphs.

Connected graphs on $$$n$$$ vertices are hard to count, but labelled graphs on $$$n$$$ vertices are trivial: there are exactly $$$2^{\binom{n}{2}}$$$ labelled graphs on $$$n$$$ vertices since we can either choose each edge or not. Let $$$G(x)$$$ denote the EGF of the number of labelled graphs. The nice thing is that labelled graphs are made up of several connected components, so we can use a similar argument as above (I will skip this step) to obtain $$$G(x) = \exp(C(x))$$$, which gives $$$C(x) = \ln(G(x))$$$. Since we know how to find $$$G(x)$$$, we can find $$$C(x)$$$ in $$$O(n\log n)$$$ time using polynomial operations.

Ok, now let's return to our problem. Let $$$b(n)$$$ denote the number of $$$2$$$-edge connected graphs on $$$n$$$ vertices and $$$B(x)$$$ be its EGF. Our goal is to find $$$B(x)$$$.

The idea is to relate $$$b(n)$$$ with $$$c(n)$$$. Suppose we have a labelled connected graph on $$$n$$$ vertices, say $$$G$$$. Any connected graph $$$G$$$ can be decomposed into a bridge tree, where each vertex is a $$$2$$$-edge connected component and the edges of the tree are the bridges of the graph. Let $$$s$$$ be the size of the $$$2$$$-edge connected component containing vertex $$$1$$$, and fix the $$$2$$$-edge connected component containing vertex $$$1$$$ as the root of the bridge tree. There are $$$\binom{n-1}{s-1}$$$ ways to choose the other elements in the component and $$$b(s)$$$ ways to connect edges within the component. Now, in the bridge tree, let $$$a_{1}, a_{2}, ..., a_{k}$$$ be the total weight of subtrees of the children of the root (we define the weight of a vertex in the bridge tree as the size of the $$$2$$$-edge connected component represented by it and the weight of a subtree as the sum of weights of all vertices in the subtree). Then, there are $$$\frac{(n-s)!}{a_{1}!a_{2}!...a_{k}!}$$$ ways to assign the remaining $$$n-s$$$ vertices to each subtree. Each subtree represented a general connected graph on $$$a_{i}$$$ vertices, so there are $$$c(a_{i})$$$ ways to connect edges in the $$$i$$$-th subtree. Finally, there are $$$sa_{i}$$$ ways to choose the "bridge" between subtree $$$i$$$ and the root, because we need to pick exactly one vertex from the subtree and the root component to connect.

Summing over all tuples $$$(a_1,a_2,...,a_k)$$$ with sum $$$n-s$$$, and dividing by $$$k!$$$ to account for the fact that each set of subtrees is counting $$$k!$$$ times, we obtain the recurrence

$$$c(n) = \displaystyle\sum_{s=1}^{n} b(s) \cdot \binom{n-1}{s-1} \cdot (n-s)! \cdot \displaystyle\sum_{k \ge 0} \frac{1}{k!} \displaystyle\sum_{a_{1}+a_{2}+...+a_{k}=n-s} \displaystyle\prod_{j=1}^{k}\frac{c(a_{j}) \cdot a_{j} \cdot s}{a_{j}!}$$$

$$$= \displaystyle\sum_{s=1}^{n} b(s) \cdot \frac{(n-1)!}{(s-1)!} \cdot \displaystyle\sum_{k \ge 0} \frac{s^{k}}{k!} \displaystyle\sum_{a_{1}+a_{2}+...+a_{k}=n-s} \displaystyle\prod_{j=1}^{k}\frac{c(a_{j}) \cdot a_{j}}{a_{j}!}$$$

Hence,

$$$\frac{nc(n)}{n!} = \displaystyle\sum_{s=1}^{n} \frac{b(s)}{(s-1)!} \cdot \displaystyle\sum_{k \ge 0} \frac{s^{k}}{k!} \displaystyle\sum_{a_{1}+a_{2}+...+a_{k}=n-s} \displaystyle\prod_{j=1}^{k}\frac{c(a_{j}) \cdot a_{j}}{a_{j}!}$$$

Note that we have $$$nc(n)$$$ and $$$a_{j} \cdot c(a_{j})$$$ appearing in the summand. It seems like it is easier to consider the EGF of the sequence $$$nc_{n}$$$, say $$$C_{1}(x)$$$. Note that $$$C_{1}(x) = xC'(x)$$$ by the "multiplication by $$$n$$$" rule. In any case, we work in generating functions to obtain

$$$C_{1}(x) = \displaystyle\sum_{n \ge 0}\frac{nc(n)}{n!}x^{n} = \displaystyle\sum_{n \ge 0}x^{n}\displaystyle\sum_{s=1}^{n} \frac{b(s)}{(s-1)!} \cdot \displaystyle\sum_{k \ge 0} \frac{s^{k}}{k!} \displaystyle\sum_{a_{1}+a_{2}+...+a_{k}=n-s} \displaystyle\prod_{j=1}^{k}\frac{c(a_{j}) \cdot a_{j}}{a_{j}!}$$$

$$$= \displaystyle\sum_{n \ge 0}\displaystyle\sum_{s=1}^{n} \frac{b(s)x^{s}}{(s-1)!} \cdot \displaystyle\sum_{k \ge 0} \frac{s^{k}}{k!} \displaystyle\sum_{a_{1}+a_{2}+...+a_{k}=n-s} \displaystyle\prod_{j=1}^{k}\frac{c(a_{j}) \cdot a_{j} \cdot x^{a_{j}}}{a_{j}!}$$$

Simplifying the interior sum and product by noting its relevance to the coefficients of $$$C_{1}(x)$$$, we obtain

$$$= \displaystyle\sum_{n \ge 0}\displaystyle\sum_{s=1}^{n} \frac{b(s)x^{s}}{(s-1)!} \cdot \displaystyle\sum_{k \ge 0} \frac{s^{k}x^{n-s}}{k!} [x^{n-s}]C_{1}(x)^k$$$

$$$= \displaystyle\sum_{n \ge 0}\displaystyle\sum_{s=1}^{n} \frac{b(s)x^{s}}{(s-1)!} \cdot x^{n-s}[x^{n-s}]\displaystyle\sum_{k \ge 0} \frac{s^{k}}{k!}C_{1}(x)^k$$$

$$$= \displaystyle\sum_{n \ge 0}\displaystyle\sum_{s=1}^{n} \frac{b(s)x^{s}}{(s-1)!} \cdot x^{n-s}[x^{n-s}]\exp(sC_{1}(x))$$$

Swapping the order of summation, we get

$$$= \displaystyle\sum_{s \ge 0}\frac{b(s)x^{s}}{(s-1)!}\displaystyle\sum_{n \ge s}x^{n-s}[x^{n-s}]\exp(sC_{1}(x))$$$

$$$= \displaystyle\sum_{s \ge 0}\frac{b(s)x^{s}}{(s-1)!}\displaystyle\sum_{n \ge 0}x^{n}[x^{n}]\exp(sC_{1}(x))$$$

$$$= \displaystyle\sum_{s \ge 0}\frac{b(s)x^{s}}{(s-1)!}\exp(sC_{1}(x))$$$.

$$$= \displaystyle\sum_{s \ge 0}\frac{sb(s)(x\exp(C_{1}(x))^{s}}{s!}$$$.

Let $$$B_{1}(x) = xB'(x)$$$ be the EGF of $$$sb(s)$$$. Thus, we obtain $$$C_{1}(x) = B_{1}(x\exp(C_{1}(x)))$$$.

We know how to find $$$C_{1}$$$ and our aim now is to find the coefficients of $$$B_{1}$$$.

Let $$$P(x) = C_{1}(x)$$$ and $$$Q(x) = x\exp(C_{1}(x))$$$. We have the relation $$$P(x) = B_{1}(Q(x))$$$ and want to find $$$[x^n]B_{1}(x)$$$. To make $$$B_{1}(x)$$$ appear, substitute $$$x$$$ as $$$Q^{-1}(x)$$$ (the compositional inverse exist because $$$Q(x)$$$ is a power series with a nonzero $$$x$$$ term and no constant term). Thus, $$$B_{1}(x) = P(Q^{-1}(x))$$$. This looks very similar to Lagrange Inversion Formula!

Indeed, we let $$$P = H$$$ and $$$Q^{-1} = G$$$. Then, $$$Q = F$$$, so we have

$$$[x^{n}]B_{1}(x) = [x^{n}]P(Q^{-1}(x)) = \frac{1}{n}[x^{-1}]P'(x)\frac{1}{Q(x)^{n}} = \frac{1}{n}[x^{-1}]C_{1}'(x)\frac{1}{x^{n}\exp(C_{1}(x))^{n}} = \frac{1}{n}[x^{n-1}]\frac{C_{1}'(x)}{\exp(C_{1}(x))^{n}}$$$.

This can be computed via standard polynomial operations in $$$O(n\log n)$$$ time.

The official editorial has a solution using PIE which ends up with an application of Lagrange Inversion Formula very similar to the example just shown. You can try to see the connection between the previous example and the trick used in the official editorial (for the Lagrange Inversion part). Here, I will demonstrate Elegia's solution described here which to me is a more intuitive way to approach the problem (thanks to Elegia for helping me with some parts I didn't understand and showing that this appraoch can lead to a full solution!).

The first step is to reduce the problem into something simpler. If we look at the samples, we see that curiously the sum of all the answers in the output is $$$n \cdot n!$$$, which suggest that there are only $$$n!$$$ good sequences. This cannot be coincidence, and a bijection between permutations and good sequences should exist.

In fact, this is IMO 2002 Shortlist C3. Since the last occurrence of $$$k$$$ must occur after the first occurrence of $$$k-1$$$, we can consider iterating through the good sequence from right to left several times. In the first run, we record the positions of the number $$$1$$$s, from right to left. On the second run, we record the positions of the number $$$2$$$s, from right to left and so on. At the end of the process, $$$n$$$ distinct numbers from $$$1$$$ to $$$n$$$ are recorded. This will be our permutation.

We claim that this is the bijection we seek. The idea is that if we read the values in our final permutation $$$p$$$ from left to right, every time we meet an ascent ($$$p_{i} < p_{i+1}$$$), it indicates that we have finished recording all occurrences of the current number and is now going from right to left again. The condition of good sequences guarantees that every time we record the occurrences of a new number, there will be a new ascent in our permutation. Thus, we can easily recover the good sequence by marking the ascents of the permutation.

This also gives us a nice way to formulate the problem. Let $$$ans(i)$$$ denote the $$$i$$$-th answer for our problem. For any permutation $$$p$$$, divide it into a minimum number of decreasing blocks. Let's say the block sizes are $$$b_{1},b_{2},...,b_{k}$$$. Then, this permutation contributes $$$b_{1}$$$ to $$$ans(1)$$$, $$$b_{2}$$$ to $$$ans(2)$$$ and so on. For example, if $$$p = (5, 3, 2, 7, 1, 4, 6)$$$, then we divide it into blocks $$$[5,3,2]$$$, $$$[7,1]$$$, $$$[4]$$$, $$$[6]$$$. $$$p$$$ increases $$$b_{1}$$$, $$$b_{2}$$$, $$$b_{3}$$$ and $$$b_{4}$$$ by $$$3$$$, $$$2$$$, $$$1$$$, $$$1$$$ respectively.

Now, let's do some double counting. Fix a position $$$1 \le j \le n$$$. Can we calculate how many times this position contributes to $$$ans(k)$$$ over all $$$n!$$$ permutations (for a fixed $$$k$$$)?

Note that for position $$$j$$$ to contribute to the answer, the prefix $$$p_{1}$$$, $$$p_{2}$$$, ..., $$$p_{j}$$$ must contain exactly $$$k-1$$$ ascents. Additionally, the last $$$n-j$$$ elements can be permuted in any manner. Finally, there are $$$\binom{n}{j}$$$ ways to choose which elements occur in the prefix. Thus, if we let $$$A(n,k)$$$ denote the number of permutations of length $$$n$$$ with exactly $$$k-1$$$ ascents, position $$$j$$$ contributes $$$(n-j)!\binom{n}{j}A(j,k) = \frac{n!}{j!}A(j,k)$$$ to the answer.

Summing over all $$$j$$$, we get the nice-looking formula $$$ans(k) = n!\displaystyle\sum_{j=1}^{n}\frac{A(j,k)}{j!}$$$. From here, you can derive a simple $$$O(n^2)$$$ solution, since there is a simple recurrence for $$$A(n,k)$$$. However, we are looking for more here.

For convenience, now we ignore the $$$n!$$$ term and just let $$$ans(k)$$$ be $$$\displaystyle\sum_{j=1}^{n}\frac{A(j,k)}{j!}$$$. We can multiply each answer by $$$n!$$$ at the end.

What's nicer is that $$$A(n,k)$$$ is actually a well-known sequence called the Eulerian numbers! If you use Wikipedia or Enumerative Combinatorics $$$1$$$, you can find the formulas related to generating functions involving $$$A(n,k)$$$.

Here, I use the notation $$$A(n,k)$$$ in Enumerative Combinatorics 1, which differs a bit from Wikipedia $$$k$$$ is shifted and $$$A(0,0)=1$$$. You can either derive or google the following formula (first thing that comes up when you look for generating functions of Eulerian numbers, though it may look slightly different due to variable shifting):

$$$F(x,y) = \displaystyle\sum_{k \ge 0}\displaystyle\sum_{n \ge 0}A(n,k)\frac{x^{n}}{n!}y^{k} = \frac{1-y}{1 - ye^{(1-y)x}}$$$

What we want to find is $$$\displaystyle\sum_{j=1}^{n}\frac{A(j,k)}{j!}$$$, so let us rewrite

$$$F(x,y) = \displaystyle\sum_{k \ge 0}y^{k}\displaystyle\sum_{n \ge 0}A(n,k)\frac{x^{n}}{n!}$$$ and focus on $$$G_{k}(x) = \displaystyle\sum_{n \ge 0}A(n,k)\frac{x^{n}}{n!}$$$ for a fixed $$$k$$$.

Observe that we want the prefix sum of the coefficients of $$$G_{k}(x)$$$. By the prefix sum trick, we get $$$ans(k) = [x^{n}]\frac{1}{1-x} \cdot G_{k}(x)$$$.

Thus, $$$ans(k) = [x^{n}y^{k}]\frac{1}{1-x}F(x,y) = [x^{n}y^{k}]\frac{1-y}{(1-x)(1 - ye^{(1-y)x})}$$$, which is exactly the same expression quoted in Elegia's post.

The hard part is finding this fast. I got to this point on my own but didn't really know how to proceed (I wasn't even sure if it was doable) before reading Elegia's post so huge thanks to him!

Let us focus on the expression $$$(1-y)[x^n]\frac{1}{(1-x)(1 - ye^{(1-y)x})}$$$ and keep in mind that we want to find the answer as a polynomial in $$$y$$$, which we can read the answer from.

The first major concern is that we have $$$e^{(1-y)x}$$$, which is an exponential of a multivariable function. This seems painfully awful to deal with especially if we need to expand it in power series form in the end. A nice trick here is to eliminate this nonsenes by making the substitution $$$z = (1-y)x$$$. Then, $$$x = \frac{z}{1-y}$$$ and our expression becomes:

$$$(1-y)[x^n]\frac{1}{(1-x)(1 - ye^{(1-y)x})} = (1-y)[x^n]\frac{1}{\left(1-\frac{z}{1-y}\right)(1 - ye^{z})} = (1-y)^{n+1}[z^n]\frac{1}{\left(1-\frac{z}{1-y}\right)(1 - ye^{z})}$$$

Let us clear the denominator in $$$1 - \frac{z}{1-y}$$$ by multiplying $$$1-y$$$ to the denominator, so we need to find

$$$(1-y)^{n+2}[z^n]\frac{1}{(1-y-z)(1 - ye^{z})}$$$

Ok, this looks simpler, though we still have products of multivariable functions. Life would be simpler if we can separate the two factors in the denominator. As a matter of fact, we can! Let's write the expression in partial fractions. Let

$$$\frac{1}{(1-y-z)(1 - ye^{z})} = \frac{A}{1 - y - z} + \frac{B}{1 - ye^{z}}$$$. We want to find some simple functions $$$A$$$, $$$B$$$ (hopefully in one variable) so that $$$A(1 - ye^{z}) + B(1 - y - z) = 1$$$.

Note that we have a linear function on $$$y$$$ on the LHS if $$$z$$$ is a treated as a constant. Thus, it motivates us to let $$$A$$$ and $$$B$$$ be functions in $$$z$$$ that annilhates the LHS. We can treat $$$z$$$ as a constant and compare the coefficients of $$$y$$$ and $$$1$$$ to obtain $$$A(z) = \frac{1}{1 - e^{z}(1-z)}$$$, $$$B(z) = \frac{-e^{z}}{1 - e^{z}(1-z)}$$$.

Substituting back, we need to find

$$$(1-y)^{n+2}[z^n]\frac{1}{(1-y-z)(1 - ye^{z})} = (1-y)^{n+2}[z^{n}]\left(\frac{1}{(1 - e^{z}(1-z))(1 - y - z)} + \frac{-e^{z}}{(1 - e^{z}(1-z))(1 - ye^{z})}\right)$$$

It remains to compute $$$[z^{n}]$$$ of each fraction fast. Our main goal is to isolate the variable $$$y$$$ as much as possible, treating $$$z$$$ pretty much like a constant.

For example, let's look at the second fraction. We want to find $$$[z^{n}]\frac{-e^{z}}{(1 - e^{z}(1-z))(1 - ye^{z})}$$$. Let $$$f(z) = \frac{-e^{z}}{1 - e^{z}(1-z)}$$$. Thus, we need to compute $$$[z^{n}]\frac{f(z)}{1 - ye^{z}}$$$.

Similarly, the first fraction can be rewritten as

$$$[z^{n}]\frac{1}{(1-e^{z}(1-z))(1-z-y)} = \frac{\frac{1}{1-z}}{(1-e^{z}(1-z)\left(1 - \frac{y}{1-z}\right)}$$$ (note that we divide by $$$1-z$$$ to make something of the form $$$\frac{1}{1-h(z)y}$$$ to make it similar to our second fraction).

Letting $$$g(z) = \frac{1}{(1-z)(1 - e^{z}(1-z))}$$$, the first fraction reduces to $$$[z^{n}]\frac{g(z)}{1 - \frac{y}{1-z}}$$$.

In both cases, we have to compute something of the form $$$[z^{n}]\frac{F(z)}{1 - G(z)y}$$$ for some functions $$$F$$$ and $$$G$$$. You can see that this is similar to $$$[x^{k}]\frac{1}{1-uf(x)}$$$ in our previous example.

Here, we have $$$G(z) = e^{z}$$$ and $$$\frac{1}{1-z}$$$. There is a subtle detail here which is that $$$[z^{0}]G(z) = 1 \neq 0$$$ in both cases, so $$$G$$$ does not have a compositional inverse and we can't apply what we did just now directly. However, $$$G(z)-1$$$ does have a compositional inverse in both cases, so let $$$A(z) = e^{z}-1$$$ and $$$B(z) = \frac{1}{1-z} - 1$$$. Thus, we need to compute $$$[z^{n}]\frac{f(z)}{1 - (A(z)+1)y}$$$ and $$$[z^{n}]\frac{g(z)}{1 - (B(z)+1)y}$$$.

Now we use the same approach as the previous example. Let $$$A^{-1}(z)$$$ denote the compositional inverse of $$$A(z)$$$ (you can find it explicitly by the definition of inverse). We let $$$H(z) = \frac{f(A^{-1}(z))}{1 - (z+1)y}$$$, $$$G(z) = A(z)$$$ and $$$F(z) = A^{-1}(z)$$$ in Lagrange Inversion Formula. Then,

$$$[z^{n}]\frac{f(z)}{1 - (A(z)+1)y} = [z^{n}]H(G(z)) = \frac{1}{n}[x^{-1}]H'(x)\frac{1}{F(x)^{n}} = \frac{1}{n}[x^{n-1}]\left(\frac{f(A^{-1}(x))}{1 - (x+1)y}\right)'\frac{1}{\left(\frac{A^{-1}(x)}{x}\right)^{n}}$$$. Let $$$C(x) = f(A^{-1}(x))$$$ and $$$D(x) = \frac{1}{\left(\frac{A^{-1}(x)}{x}\right)^{n}}$$$. Take a moment to check that we can still compute the first $$$n$$$ terms of $$$C(x)$$$ and $$$D(x)$$$ using polynomial operations in $$$O(n\log n)$$$ (find $$$A^{-1}(x)$$$ and substitute back into their definitions).

Now, our problem reduces to finding $$$\frac{1}{n}[x^{n-1}]\left(\frac{C(x)}{1-(x+1)y}\right)'D(x)$$$. Using quotient rule, we have (differentiating with respect to $$$x$$$),

$$$\left(\frac{C(x)}{1-(x+1)y}\right)' = \frac{C'(x)[1-(x+1)y] - C(x)(-y)}{(1 - (x+1)y)^2} = \frac{C'(x)}{1 - (x+1)y} + \frac{C(x)y}{(1 - (x+1)y)^2}$$$.

Thus,

$$$\frac{1}{n}[x^{n-1}]\left(\frac{C(x)}{1-(x+1)y}\right)'D(x) = \frac{1}{n}[x^{n-1}]\frac{C'(x)D(x)}{1 - (x+1)y} + [x^{n-1}]\frac{C(x)D(x)y}{(1 - (x+1)y)^2}$$$.

We have two subproblems, finding $$$[x^{n}]\frac{P(x)}{1 - (x+1)y}$$$ and $$$[x^{n}]\frac{P(x)y}{(1-(x+1)y)^2}$$$ for some computable functions $$$P$$$. Both can be dealt with using the geometric series formula. We have

$$$[x^{n}]\frac{P(x)}{1 - (x+1)y} = [x^{n}]P(x)\displaystyle\sum_{i \ge 0}(x+1)^{i}y^{i}$$$

Let $$$P(x) = p_{0} + p_{1}x + p_{2}x^{2} + ...$$$ be the series expansion, then

$$$[x^{n}]P(x)\displaystyle\sum_{i \ge 0}(x+1)^{i}y^{i} = \displaystyle\sum_{i \ge 0}y^{i}[x^{n}]\left(P(x)(x+1)^{i}\right) = \displaystyle\sum_{i \ge 0}y^{i}\displaystyle\sum_{j \ge 0}p_{n-j}\binom{i}{j}$$$.

Hence, we need to compute $$$ans(i) = \displaystyle\sum_{j \ge 0}p_{n-j}\binom{i}{j} = i!\displaystyle\sum_{j \ge 0}\frac{p_{n-j}}{j!(i-j)!}$$$ fast. But this is almost the same thing as what we did in the Atcoder problem mentioned at the very beginning of this article! Define $$$E(x) = \displaystyle\sum_{i \ge 0}\frac{p_{n-i}}{i!}x^{i}$$$ and $$$F(x) = \displaystyle\sum_{i \ge 0}\frac{1}{i!}x^{i}$$$ for some large enough $$$M$$$. Then, our answer is the coefficient of $$$x^{i}$$$ in $$$E(x)F(x)$$$. This part can be solved in $$$O(n\log n)$$$ time.

Finally, we need to compute $$$[x^{n}]\frac{P(x)y}{(1 - (x+1)y)^2}$$$. With a similar expansion,

$$$[x^{n}]\frac{P(x)y}{(1 - (x+1)y)^2} = [x^{n}]P(x)y\displaystyle\sum_{i \ge 0}(i+1)(x+1)^{i}y^{i}$$$

$$$= \displaystyle\sum_{i \ge 0}(i+1)y^{i+1}[x^{n}]\left(P(x)(x+1)^{i}\right)$$$

$$$= \displaystyle\sum_{i \ge 0}(i+1)y^{i+1}[x^{n}]\left(P(x)(x+1)^{i}\right)$$$

$$$= \displaystyle\sum_{i \ge 0}(i+1)y^{i+1}\displaystyle\sum_{j \ge 0}p_{n-j}\binom{i}{j}$$$.

Thus, we can compute this in the same way as before in $$$O(n\log n)$$$ using FFT.

Putting it altogether, we obtain a (albeit complicated) solution in $$$O(n\log n)$$$ time!

I hope this explanation makes Elegia's solution more intuitive to understand. :) Thanks to Elegia for the wonderful solution.

This is unfortunately the most standard problem of the set. Obviously, we can model the friends as vertices and friendships as edges in an undirected graph. The problem basically asks us to answer queries of the form: "For a pair of vertices $$$u, v$$$, find the number of vertices $$$w \neq u, v$$$ such that removing $$$w$$$ from the graph disconnects $$$u$$$ and $$$v$$$".

Removing vertices and disconnecting graphs should remind one of articulation points. The data structure to solve this problem is block-cut tree. Each biconnected component is considered as a block. An articulation point might belong to multiple blocks, but a non-articulation point only belongs to one block. We build a tree where the vertices are blocks and the articulation points. We connect each block to the articulation points it contains. If the graph was already biconnected, the answer is obviously $$$0$$$.

After we build the block-cut tree, let’s see how to answer the queries. Clearly, only removing articulation points have an effect on the connectivity of $$$u$$$ and $$$v$$$. With the block-cut tree, it can be easily seen that the articulation points that disconnects $$$u$$$ and $$$v$$$ are precisely the articulation points (not equal to $$$u$$$ or $$$v$$$) on the path from the block corresponding to $$$u$$$ (or $$$u$$$ itself if $$$u$$$ is an articulation point) and the block corresponding to $$$v$$$ (or $$$v$$$ itself if $$$v$$$ is an articulation point). This can be easily handled using LCA and calculating the number of articulation points on the path to root. The time complexity is $$$O(M + Q\log N)$$$.

Subtask 1 can be solved by brute forcing which vertices to delete and Subtask 2 can be solved with LCA.

Hopefully, this problem brings forth a slightly uncommon topic, the block-cut tree, for more inexperienced participants.

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
 
using namespace std;
using namespace __gnu_pbds;
 
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key
 
typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;

const int N = 200000;
int A[N*2+10];
int at[N*2+10];

struct Tree
{
	struct data
	{
		ll w;
	};
	
	struct node
	{
		int p; //parent
		ll w; //modify for different problems
	};
	
	struct edge
	{
		int v; data dat;
	};
	
	vector<vector<edge> > adj;
	int n;
	
	Tree(int _n)
	{
		adj.resize(_n);
		n = _n;
	}
	
	vi level;
	vi depth;
	vi h;
	vi h2;
	vi euler;
	vi firstocc;
	vector<vi> rmqtable;
	vi subsize;
	vi start; vi en;
	vector<vector<node> > st;
	
	void addedge(int u, int v)
	{
		edge tmp; tmp.v = v;
		adj[u].pb(tmp);
		tmp.v = u;
		adj[v].pb(tmp);
	}
	
	void reset(int _n)
	{
		adj.clear();
		level.clear();
		depth.clear();
		euler.clear();
		rmqtable.clear();
		subsize.clear();
		start.clear();
		en.clear();
		st.clear();
		firstocc.clear();
		adj.resize(_n);
		n = _n;
	}
	
	void dfssub(int u, int p)
	{
		subsize[u] = 1;
		for(int i = 0; i < adj[u].size(); i++)
		{
			int v = adj[u][i].v;
			if(v == p) continue;
			dfssub(v, u);
			subsize[u] += subsize[v];
		}
	}
	
	void calcsub()
	{
		subsize.resize(n);
		dfssub(0, -1);
	}
	
	int timer;
	
	void dfsstartend(int u, int p)
	{
		start[u] = ++timer;
		if(p == -1) h[u] = 0;
		else h[u] = h[p] + 1;
		for(int i = 0; i < adj[u].size(); i++)
		{
			int v = adj[u][i].v;
			if(v == p) continue;
			dfsstartend(v, u);
		}
		en[u] = ++timer;
	}
	
	void calcstartend()
	{
		timer = 0;
		start.resize(n); en.resize(n); h.resize(n);
		dfsstartend(0, -1);
	}
	
	int eulercnt;
	
	void dfseuler(int u, int p)
	{
		euler[eulercnt] = u; eulercnt++;
		if(p == -1) {depth[u] = 0;}
		else {depth[u] = depth[p] + 1;}
		firstocc[u] = eulercnt-1;
		for(int i = 0; i < adj[u].size(); i++)
		{
			int v = adj[u][i].v;
			if(v == p) continue ;
			dfseuler(v, u);
			euler[eulercnt] = u; eulercnt++;
		}
	}
	
	void calceuler()
	{
		eulercnt = 0;
		level.assign(2*n+1, 0);
		euler.assign(2*n+1, 0);
		depth.assign(n, 0);
		firstocc.resize(n);
		dfseuler(0, -1);
	}

	void filllevel()
	{
		int LG = 0;
		while((1<<LG) <= n*2) LG++;
		rmqtable.resize(LG);
		for(int i = 0; i < LG; i++) rmqtable[i].resize(eulercnt);
		for(int i = 0; i < eulercnt; i++)
		{
			level[i] = depth[euler[i]];
		}
		level[eulercnt] = 1000000000;
		for(int j = 0; j < LG; j++)
		{
			for(int i = 0; i < eulercnt; i++)
			{
				rmqtable[j][i] = eulercnt;
				if(i + (1<<j) - 1 < eulercnt)
				{
					if(j == 0)
					{
						rmqtable[j][i] = i;
					}
					else
					{
						if(level[rmqtable[j - 1][i]] < level[rmqtable[j-1][i + (1<<(j-1))]])
						{
							rmqtable[j][i] = rmqtable[j-1][i];
						}
						else
						{
							rmqtable[j][i] = rmqtable[j-1][i + (1<<(j-1))];
						}
					}
				}
			}
		}
	}

	int rmq(int l, int r)
	{
		int k = 31 - __builtin_clz(r-l);
		//cout << l << ' ' << r << ' ' << rmqtable[l][k] << ' ' << rmqtable[r - (1<<k) + 1][k] << endl;
		if(level[rmqtable[k][l]] < level[rmqtable[k][r - (1<<k) + 1]])
		{
			return rmqtable[k][l];
		}
		else
		{
			return rmqtable[k][r - (1<<k) + 1];
		}
	}

	int lcaeuler(int u, int v)
	{
		if(firstocc[u] > firstocc[v]) swap(u, v);
		//cerr << firstocc[u] << ' ' << firstocc[v] << ' ' << rmq(firstocc[u], firstocc[v]) << ' ' << euler[rmq(firstocc[u], firstocc[v])] << endl;
		return euler[rmq(firstocc[u], firstocc[v])];
	}
	
	bool insub(int u, int v) //is u in the subtree of v?
	{
		if(start[v] <= start[u] && en[u] <= en[v]) return true;
		return false;
	}
	
	void dfspar(int u, int p)
	{
		//cerr << u << ' ' << p << '\n';
		st[0][u].p = p;
		if(p == -1) 
		{
			h2[u] = A[u]; h[u]=0;
		}
		else 
		{
			h2[u] = h2[p] + A[u];
			h[u]=h[p]+1;
		}
		//cerr<<"H: "<<u<<' '<<h[u]<<' '<<A[u]<<'\n';
		for(int i = 0; i < adj[u].size(); i++)
		{
			int v = adj[u][i].v;
			if(v == p) continue;
			dfspar(v, u);
		}
	}
	
	int LOG;
	
	void calcpar()
	{
		h.resize(n); h2.resize(n);
		int LG = 0; LOG = 0;
		while((1<<LG) <= n) {LG++; LOG++;}
		st.resize(LG);
		for(int i = 0; i < LG; i++)
		{
			st[i].resize(n);
		}
		dfspar(0, -1);
		//cerr << "HER" << ' ' << LG << endl;
		for(int i = 1; i < LG; i++)
		{
			for(int j = 0; j < n; j++)
			{
				if(st[i-1][j].p == -1) st[i][j].p = -1;
				else st[i][j].p = st[i-1][st[i-1][j].p].p;
			}
		}
	}
	
	int getpar(int u, ll k)
	{
		for(int i = LOG - 1; i >= 0; i--)
		{
			if(k&(1<<i))
			{
				u = st[i][u].p;
			}
		}
		return u;
	}
	
	int lca(int u, int v)
	{
		if(h[u] > h[v]) swap(u, v);
		for(int i = LOG - 1; i >= 0; i--)
		{
			if(st[i][v].p != -1 && h[st[i][v].p] >= h[u])
			{
				v = st[i][v].p;
			}
		}
		if(u == v) return u;
		for(int i = LOG - 1; i >= 0; i--)
		{
			if(st[i][v].p != -1 && st[i][v].p != st[i][u].p)
			{
				u = st[i][u].p;
				v = st[i][v].p;
			}
		}
		return st[0][u].p;
	}

	int distance(int u, int v)
	{
		int lc = lca(u, v);
		return (h[u]+h[v]-2*h[lc]);
	}
};

Tree T(1);
struct graph
{
	int n;
	vector<vector<int>> adj;
 
	graph(int n) : n(n), adj(n) {}
 
	void add_edge(int u, int v)
	{
		adj[u].push_back(v);
		adj[v].push_back(u);
	}
 
	int add_node()
	{
		adj.push_back({});
		return n++;
	}
 
	vector<int>& operator[](int u) { return adj[u]; }
};
 
vector<int> id;

void biconnected_components(graph &adj)
{
	int n = adj.n;
 
	vector<int> num(n), low(n), art(n), stk;
	vector<vector<int>> comps;
 
	function<void(int, int, int&)> dfs = [&](int u, int p, int &t)
	{
		num[u] = low[u] = ++t;
		stk.push_back(u);
 
		for (int v : adj[u]) if (v != p)
		{
			if (!num[v])
			{
				dfs(v, u, t);
				low[u] = min(low[u], low[v]);
 
				if (low[v] >= num[u])
				{
					art[u] = (num[u] > 1 || num[v] > 2);
 
					comps.push_back({u});
					while (comps.back().back() != v)
						comps.back().push_back(stk.back()), stk.pop_back();
				}
			}
			else low[u] = min(low[u], num[v]);
		}
	};
 
	for (int u = 0, t; u < n; ++u)
		if (!num[u]) dfs(u, -1, t = 0);
 
	// build the block cut tree
	graph tree(0);
	id.resize(n);

	for (int u = 0; u < n; ++u)
	{
		if(art[u]) 
		{
			at[u]=art[u];
			id[u]=tree.add_node();
			A[id[u]]=1;
		}
	}

	for (auto &comp : comps)
	{
		int node = tree.add_node();
		for (int u : comp)
		{
			if (!art[u]) id[u] = node;
			else 
			{
				T.addedge(node,id[u]);
			}
		}
	}
 
}

//main part of solution

int query(int u, int v)
{
	int ans = -at[u]-at[v];
	u=id[u]; v=id[v];
	if(u==v) return 0;
	int lc = T.lca(u,v);
	ans+=T.h2[u]+T.h2[v]+A[lc]-2*T.h2[lc];
	return ans;
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	int n,m; cin>>n>>m;
	graph G(n);
	for(int i=0;i<m;i++)
	{
		int u,v; cin>>u>>v; u--; v--;
		G.add_edge(u,v);
	}
	T.reset(2*n+10);
	biconnected_components(G);
	T.calcpar();
	int q; cin>>q;
	for(int z=0;z<q;z++)
	{
		int u,v; cin>>u>>v; u--; v--;
		cout<<query(u,v)<<'\n';
	}
}

#include <bits/stdc++.h>
using namespace std;
 
#define ll long long
#define ar array
 
const int mxN=2e5;
int n, m, bccI, dt, tin[mxN], low[mxN], st[mxN], sth, c[2*mxN], d[2*mxN], anc[2*mxN][19], q;
vector<int> adj1[mxN], adj2[2*mxN];
 
void dfs1(int u=0, int p=-1) {
	tin[u]=low[u]=++dt;
	st[sth++]=u;
	for(int v : adj1[u]) {
		if(v==p)
			continue;
		if(!tin[v]) {
			dfs1(v, u);
			if(low[v]>=tin[u]) {
				adj2[u].push_back(n+bccI);
				do {
					adj2[n+bccI].push_back(st[sth-1]);
				} while(st[--sth]^v);
				++bccI;
			}
			low[u]=min(low[u], low[v]);
		} else
			low[u]=min(low[u], tin[v]);
	}
}
 
void dfs2(int u=0) {
	for(int i=1; i<19; ++i)
		anc[u][i]=anc[anc[u][i-1]][i-1];
	c[u]+=u<n;
	for(int v : adj2[u]) {
		c[v]=c[u];
		d[v]=d[u]+1;
		anc[v][0]=u;
		dfs2(v);
	}
}
 
int lca(int u, int v) {
	if(d[u]<d[v])
		swap(u, v);
	for(int i=18; ~i; --i)
		if(d[u]-(1<<i)>=d[v])
			u=anc[u][i];
	if(u==v)
		return u;
	for(int i=18; ~i; --i) {
		if(anc[u][i]^anc[v][i]) {
			u=anc[u][i];
			v=anc[v][i];
		}
	}
	return anc[u][0];
}
 
int main() {
	ios::sync_with_stdio(0);
	cin.tie(0);
 
	cin >> n >> m;
	for(int i=0, u, v; i<m; ++i) {
		cin >> u >> v, --u, --v;
		adj1[u].push_back(v);
		adj1[v].push_back(u);
	}
	dfs1();
	dfs2();
	cin >> q;
	for(int u, v; q--; ) {
		cin >> u >> v, --u, --v;
		int w=lca(u, v);
		cout << (c[u]+c[v]-c[w]-(w?c[anc[w][0]]:0)-2) << "\n";
	}
}

Subtask 1 can be solved by trying all permutations. From now on, we will refer to the graph as the graph $$$i \rightarrow p_i$$$ for all $$$i$$$.

For subtask 2, note that for $$$N > 2$$$, the first person will always fail and we can reduce our graph into a chain. The idea is that after any set of moves, our graph will be a set of independent chains where the head of the chain might or might not have confessed. Thus, we can compute the answer using $$$dp[n][0], dp[n][1]$$$ denoting the answer for a chain of length $$$n$$$ and whether the head of the chain has confessed. This leads to a direct $$$O(N^2)$$$ solution for this subtask. There are also faster solutions for this subtask but as we will not use it for the main solution, I will leave them as an exercise for the reader. :)

Now, let’s solve the general case. My idea is to use linearity of expectation and analyze the probability that student $$$i$$$ will be accepted by their crush. For simplicity I will ignore the existence of 2-cycles here, but the solution idea remains the same.

Fix a student $$$i$$$. We look at the chain of crushes $$$i, p_i, p_{p_i}, …$$$ until someone occurs in the list twice. For convenience, label them as $$$x_0=i, x_1, …, x_c, x_{c+1}, …, x_{d}, x_{d+1}=x_c$$$ (where $$$c$$$ can possibly be $$$0$$$).

Let’s look at the conditions necessary for student $$$i$$$ to be accepted when they confess. Firstly, $$$p_i$$$ must be to the left of $$$i$$$ in the permutation $$$a$$$. However, that itself is insufficient, for two reasons. Firstly, let $$$y$$$ be another person with a crush on $$$p_i$$$. Then, if $$$y$$$ appears between $$$p_i$$$ and $$$i$$$ in the permutation, then if $$$i$$$ would have succeeded, $$$y$$$ would succeed first, a contradiction. Thus, let $$$S_{x}$$$ be the set of students with a crush on student $$$x$$$. Then, all elements in $$$S_{p_{i}} \setminus \{i\}$$$ must not lie between $$$p_i$$$ and $$$i$$$ in the permutation.

Are these conditions sufficient? No! It might also happen that student $$$p_i=x_1$$$ has already confessed successfully and is thus unavailable. Now we ask ourselves, when does this situation occur? This looks like the same problem we were trying to solve! However, there are a few additional conditions which make it not as simple: $$$x_0$$$ must occur to the right of $$$x_1$$$ in the permutation, and all elements of $$$S_{x_{1}}$$$ must not lie strictly between $$$x_0$$$ and $$$x_1$$$ in the permutation.

To compute the number of permutations where the additional conditions are satisfied and $$$x_1$$$ confessed successfully, we follow essentially the same process. $$$x_2$$$ must be to the left of $$$x_1$$$ in the permutation, and all elements of $$$S_{x_{2}}$$$ must not be between $$$x_1$$$ and $$$x_2$$$. However, we need to subtract the number of ways when $$$x_2$$$ confessed successfully, and thus we need to repeat the same process again but with $$$x_2$$$.

In general, we have to solve a problem of the following form: Count the number of permutations where $$$x_a, x_{a-1}, …, x_{0}$$$ appear from left to right in this order and any element of $$$S_{x_{i}} \setminus \{x_{i-1}\} $$$ for $$$i \ge 1$$$ does not appear between $$$x_{i}$$$ and $$$x_{i-1}$$$ in the permutation.

Note that $$$S_{x_{i}} \setminus \{x_{i-1}\}$$$ are all disjoint and does not contain elements in our chain $$$x_0, x_1, …, x_{d}$$$, with the sole exception of possibly $$$i = c$$$ and $$$x_{d} \in S_{x_{c}}$$$. However, it turns out that we can handle this sole exception easily as it is already near the end of our chain. Hence, only the sizes of $$$S_{x_{i}}$$$ are important to us and the actual elements can be ignored.

To solve our subproblem, we will use the idea of PIE (Principle of Inclusion-Exclusion). Here’s the general sketch:

First, we start with the sequence $$$x_a, x_{a-1}, …, x_{0}$$$. Now, we want to insert $$$b_{0}$$$ 0s, $$$b_{1}$$$ 1s, …, $$$b_{a-1}$$$ (a-1)s into the sequence such that there is no $$$i$$$ between $$$x_{i}$$$ and $$$x_{i+1}$$$ for all $$$i$$$ (the values $$$b_{i}$$$ correspond to the size of the set $$$S_{x_{i+1}}$$$ subtracted by $$$1$$$). We call a value $$$i$$$ between the elements $$$x_{i}$$$ and $$$x_{i+1}$$$ as a violation. For any final sequence $$$F$$$, we will count the number of pairs $$$(F, V)$$$ where $$$V$$$ is a subset of violations in $$$F$$$.

It turns out that we can compute the number of pairs $$$(F, V)$$$ using some $$$dp[i][j]$$$, where $$$i$$$ denotes that we have considered the numbers from $$$0$$$ to $$$i$$$, and now we are trying to add violations with the value $$$i+1$$$. $$$j$$$ denotes the total number of violations we have added till now. With some binomial coefficients, we can iterate through all values from $$$0$$$ to $$$b_{i+1}$$$ denoting the number of violations between $$$x_{i+1}$$$ and $$$x_{i+2}$$$ we add to $$$V$$$, and compute the dp in $$$O(N^{3})$$$ time total.

Naively doing this for every $$$a$$$ will give a solution that takes $$$O(N^{5})$$$ total (since we have to do it separately for each $$$x_{0} = i$$$), but it is easy to see that we only have to do our dp once (or actually twice to handle the edge case of the final person in the chain as mentioned above) to compute the answer for all $$$a$$$. A simple implementation of the ideas above gives an $$$O(N^4)$$$ solution, which is sufficient to pass.

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
 
using namespace std;
using namespace __gnu_pbds;
 
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key
 
typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;

const int N = 200;
const int MOD = (1e9 + 7);
int add(int a, int b)
{
	a+=b;
	while(a>=MOD) a-=MOD;
	return a;
}
void radd(int &a, int b)
{
	a=add(a,b); 
}
int mult(int a, int b)
{
	return (a*1LL*b)%MOD;
}
void rmult(int &a, int b)
{
	a=mult(a,b);
}
int modpow(int a, int b)
{
	int r=1;
	while(b)
	{
		if(b&1) r=mult(r,a);
		a=mult(a,a);
		b>>=1;
	}
	return r;
}
int inverse(int a)
{
	return modpow(a,MOD-2);
}

int ncr[N+10][N+10];
int fact[N+10];
int inv[N+10];
int ifact[N+10];

int choose(int n, int r)
{
	if(r>n||r<0) return 0;
	if(r==0||r==n) return 1;
	if(ncr[n][r]!=-1) return ncr[n][r];
	return (ncr[n][r]=add(choose(n-1,r),choose(n-1,r-1)));
}

int dp[N+10][N+10];
int solve_fast(vi a)
{
	int n=a.size();
	vi indeg(n,0);
	for(int i=0;i<n;i++) indeg[a[i]]++;
	int ans=0;
	for(int u=0;u<n;u++)
	{
		if(a[a[u]]==u)
		{
			ans=add(ans,inv[2]); //auto couple
			continue;
		}
		//check for the path+cycle combo
		bitset<N+10> visited;
		visited.reset();
		int cur=u;
		vi path;
		path.pb(cur);
		while(1)
		{
			visited[cur]=1;
			cur=a[cur];
			path.pb(cur);
			if(visited[cur]) break;
		}
		int id=-1;
		for(int i=0;i<path.size();i++)
		{
			if(path[i]==path.back())
			{
				id=i; break;
			}
		}
		assert(id!=-1);
		//probability that cur will succeed?
		vi b;
		int res=0; //count # of valid permutations
		memset(dp,0,sizeof(dp));
		dp[0][0]=1;
		int bsum=0;
		for(int i=0;i+1<int(path.size())-1;i++) 
		{
			int v=path[i+1]; //person u success
			if(id+2==int(path.size())-1) //cycle of size 2
			{
				if(i==id-1) break; //sure fail
			}
			//last person
			if(i+2==int(path.size())-1&&id>0) 
			{
				b.clear(); bsum=0;
				for(int j=0;j<=i;j++)
				{
					int v=path[j+1];
					if(v==path.back()) b.pb(indeg[v]-2);
					else b.pb(indeg[v]-1);
					bsum+=b.back();
				}
				memset(dp,0,sizeof(dp));
				dp[0][0]=1;
				for(int i2=0;i2<=i;i2++)
				{
					for(int j=0;j<=n;j++)
					{
						if(dp[i2][j]==0) continue;
						int v=dp[i2][j];
						for(int k=0;k<=b[i2];k++)
						{
							radd(dp[i2+1][j+k],mult(v,mult(fact[k],choose(b[i2],k))));
						}
					}
				}				
			}
			else
			{
				b.pb(indeg[v]-1);
				bsum+=indeg[v]-1;
				for(int j=0;j<=n;j++)
				{
					if(dp[i][j]==0) continue;
					int v=dp[i][j];
					for(int k=0;k<=b[i];k++)
					{
						radd(dp[i+1][j+k],mult(v,mult(fact[k],choose(b[i],k))));
					}
				}
			}
			int curans=0;
			for(int j=0;j<=n;j++)
			{
				if(j%2==0) radd(curans,mult(dp[i+1][j],mult(fact[bsum-j],choose(bsum+i+2,bsum-j))));
				else radd(curans,MOD-mult(dp[i+1][j],mult(fact[bsum-j],choose(bsum+i+2,bsum-j))));
			}
			curans=mult(curans,mult(choose(n,bsum+i+2),fact[n-(bsum+i+2)]));
			if(i%2==0) radd(res,curans);
			else radd(res,MOD-curans);
		}
		res=mult(res,ifact[n]);
		ans=add(ans,res);
	}
	return ans;
}

vi read()
{
	int n; cin>>n;
	vi a(n);
	for(int i=0;i<n;i++) 
	{
		cin>>a[i]; a[i]--;
	}
	return a;
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	fact[0]=ifact[0]=1;
	for(int i=1;i<=N+1;i++) 
	{
		fact[i]=mult(fact[i-1],i);
		ifact[i]=inverse(fact[i]);
		inv[i]=inverse(i);
	}
	memset(ncr,-1,sizeof(ncr));
	vi a = read();
	cout<<solve_fast(a)<<'\n';
}

For subtask $$$1$$$, brute forcing all possible paths work. For subtask $$$2$$$, we can do a simple $$$dp[n][x][y]$$$ denoting the number of valid paths from $$$(x, y)$$$ in $$$O(N^3)$$$. The real question is how to exploit the small value of $$$D$$$ to solve the problem is sub-cubic time.

This turns out to be a simple exercise on PIE (Principle of Inclusion-Exclusion). We call a point $$$(x, y)$$$ bad if $$$|x| + |y| = D$$$. Note that a path is valid if and only if it does not pass through a bad point. Then, a random path might go through bad points several times. For an arbitrary path (not necessarily good), let $$$\{b_{1}, b_{2}, …, b_{k}\}$$$ be the multiset of all bad points it passes through, at times $$$t_{1}, t_{2}, …, t_{k}$$$ respectively. We will count the pairs $$$(P, B)$$$, where $$$P$$$ is any path and $$$B$$$ is a subset of $$${(b_{1}, t_{1}), (b_{2}, t_{2}), …, (b_{k}, t_{k})}$$$. Let $$$dp[i][j][k]$$$ denote the number of pairs $$$(P, B)$$$ such that $$$P$$$ is a path of length $$$i$$$ that ends at bad point $$$j$$$ and the size of $$$B$$$ is $$$k$$$. To compute this dp, we can iterate over all smaller $$$i’$$$ and all other $$$j’$$$ and use the value of $$$dp[i’][j’][k-1]$$$ to update $$$dp[i][j][k]$$$ (special handling for the case $$$k = 1$$$). However, there is a small problem. We need to compute the number of ways to go from some $$$(a, b)$$$ to $$$(c, d)$$$ in exactly $$$M$$$ moves fast (without any other conditions).

Fortunately, this is doable in $$$O(1)$$$ with binomial coefficients, as we claim that the answer is $$$\displaystyle\binom{M}{\frac{M+c+d-a-b}{2}} \displaystyle\binom{M}{\frac{M+c-d-a+b}{2}}$$$. The basic idea is to biject each path $$$P$$$ of length $$$M$$$ from $$$(a, b)$$$ to $$$(c, d)$$$ into a pair of up-right paths of length $$$M$$$ $$$(P_1, P_2)$$$. For each right move in $$$P$$$, biject it to up in $$$P_1$$$ and $$$P_2$$$. For each left move in $$$P$$$, biject it to right in $$$P_1$$$ and $$$P_2$$$. For each up move in $$$P$$$, biject it to up in $$$P_1$$$ and right in $$$P_2$$$. For each down move in $$$P$$$, biject it to right in $$$P_1$$$ and up in $$$P_2$$$. Then, it is easy to see that $$$P_1$$$ has $$$c-a+d-b$$$ more up moves than right moves and $$$P_2$$$ has $$$c-a+b-d$$$ more up moves than right moves. It can be readily seen that this map is a bijection. Counting the up-right paths gives us the claimed formula.

Armed with the formula, we are able to do our dp transitions in $$$O(ND)$$$ time. However, we still have $$$O(N^2D)$$$ states in total. We can easily reduce it to $$$O(ND)$$$ states by noticing that only the parity of $$$k$$$ matters, since we are doing inclusion-exclusion. This gives a $$$O(N^2D^2)$$$ dp solution.

Note that we can actually omit the parameter $$$k$$$ in our dp altogether by reversing signs in our calculations appropriately (see code), but the existence of $$$k$$$ makes the explanation easier.

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
 
using namespace std;
using namespace __gnu_pbds;
 
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key
 
typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;

const int N = 2011;
const int MOD = (1e9 + 7);
int add(int a, int b)
{
	a+=b;
	while(a>=MOD) a-=MOD;
	return a;
}
void radd(int &a, int b)
{
	a=add(a,b); 
}
int mult(int a, int b)
{
	return (a*1LL*b)%MOD;
}
void rmult(int &a, int b)
{
	a=mult(a,b);
}
int modpow(int a, int b)
{
	int r=1;
	while(b)
	{
		if(b&1) r=mult(r,a);
		a=mult(a,a);
		b>>=1;
	}
	return r;
}
int inverse(int a)
{
	return modpow(a,MOD-2);
}

int ncr[N+10][N+10];
int fact[N+10];
int inv[N+10];
int ifact[N+10];

int choose(int n, int r)
{
	if(r>n||r<0) return 0;
	if(r==0||r==n) return 1;
	if(ncr[n][r]!=-1) return ncr[n][r];
	return (ncr[n][r]=add(choose(n-1,r),choose(n-1,r-1)));
}

int x,y,n,d; 

void read()
{
	cin>>x>>y>>n>>d;
}

int dx[4] = {1,-1,0,0};
int dy[4] = {0,0,1,-1};

int ways(int a, int b, int c, int d, int n)
{
	if((n+a+b+c+d)%2==0)
	{
		return mult(choose(n,(n+c+d-a-b)/2),choose(n,(n+c-d-a+b)/2));
	}
	else return 0;
}

vector<ii> badpt;
const int D = 5;
int dp[N+10][4*D+10];
int solve_n2()
{
	int mand = abs(x)+abs(y);
	if(mand<=d) return 0;
	if(mand-n>d)
	{
		return modpow(4,n);
	}
	badpt.clear();
	for(int i=-d;i<=d;i++)
	{
		for(int j=-d;j<=d;j++)
		{
			if(abs(i)+abs(j)==d) badpt.pb({i,j});
		}
	}
	int badsiz=badpt.size();
	assert(badsiz<4*D+10);
	memset(dp,0,sizeof(dp));
	for(int i=1;i<=n;i++)
	{
		for(int j=0;j<badsiz;j++)
		{
			dp[i][j]=(MOD-ways(x,y,badpt[j].fi,badpt[j].se,i))%MOD;
		}
	}
	for(int i=2;i<=n;i++)
	{
		for(int j=0;j<badsiz;j++)
		{
			for(int k=1;k<i;k++)
			{
				for(int l=0;l<badsiz;l++)
				{
					if(dp[k][l]==0) continue;
					radd(dp[i][j],MOD-mult(dp[k][l],ways(badpt[l].fi,badpt[l].se,badpt[j].fi,badpt[j].se,i-k)));
				}
			}
		}
	}
	int ans=modpow(4,n);
	for(int i=1;i<=n;i++)
	{
		for(int j=0;j<badsiz;j++)
		{
			radd(ans,mult(dp[i][j],modpow(4,n-i)));
		}
	}
	return ans;
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	memset(ncr,-1,sizeof(ncr));
	read();
	cout<<solve_n2()<<'\n';
}

Basically, the problem asks you to find the number of positive integers $$$T$$$ such that $$$(2k+1)T$$$ is in set $$$A$$$ (comprised of $$$X$$$ disjoint intervals) for all $$$k \ge 0$$$ while $$$2kT$$$ is in set $$$B$$$ (comprised of $$$Y$$$ disjoint intervals) for all $$$k \ge 1$$$.

For Subtask 1, we can literally try all possible values of $$$T$$$, since it takes $$$O(\frac{N}{T})$$$ time to check if $$$T$$$ works, and thus the total complexity is $$$O(N\log N)$$$ since $$$\frac{N}{1} + \frac{N}{2} + … + \frac{N}{N} \approx N\log N$$$.

How to solve the general case? First, we show an $$$O(X+Y)$$$ time solution to test a fixed value of $$$T$$$. The key idea is to look at the complement. Let $$$A’, B’$$$ be the complements of $$$A$$$ and $$$B$$$ respectively (which also consists of $$$O(X)$$$ and $$$O(Y)$$$ intervals). A value $$$T$$$ is invalid iff some interval in $$$A$$$ contains a number of the form $$$(2k+1)T$$$ or some interval in $$$B$$$ contains a number of the form $$$2kT$$$. Checking these conditions for an interval takes $$$O(1)$$$ time, so we can test a fixed value of $$$T$$$ in $$$O(X+Y)$$$ time.

However, we have $$$O(N)$$$ values of $$$T$$$, so our $$$O(N(X+Y))$$$ time solution is still too slow. The key idea here is that if $$$T$$$ is large, the value $$$k$$$ above will be of order $$$\frac{N}{T}$$$, which is small. Thus, we can use a square root decomposition idea. For $$$T \le \sqrt{N}$$$, we use the naive $$$O(X+Y)$$$ solution to test it. Now, we want to find all the bad values of $$$T$$$ larger than $$$\sqrt{N}$$$ .

Let’s fix an interval from $$$A’$$$, say $$$[l, r]$$$, and see which values of $$$T$$$ it will invalidate. A value $$$T$$$ fails iff there exist $$$k \ge 0$$$ such that $$$l \le (2k+1)T \le r$$$, or $$$\frac{l}{2k+1} \le T \le \frac{r}{2k+1}$$$. Sine $$$k \le \sqrt{N}$$$ when $$$T \ge \sqrt{N}$$$, we can loop through all the values of $$$k \le \sqrt{N}$$$ and add the interval of $$$T$$$ it invalidates to our set of candidate bad values. Thus, for each of the $$$X+Y$$$ intervals, we add $$$O(\sqrt{N})$$$ intervals of $$$T$$$. Finally, we sort these intervals and extract their union to find the number of distinct $$$T$$$ larger than $$$\sqrt{N}$$$ which are bad.

The solution works in $$$O((X+Y)\sqrt{N}\log({N(X+Y)}))$$$ time.

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
 
using namespace std;
using namespace __gnu_pbds;
 
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key
 
typedef long long ll;
typedef pair<int,int> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;

vector<ii> clean(vector<ii> vec)
{
	sort(vec.begin(),vec.end());
	vector<ii> V;
	int R=-1; int L=-1;
	for(int i=0;i<vec.size();i++)
	{
		int l=vec[i].fi; int r=vec[i].se;
		if(l>r) continue;
		if(l>R+1)
		{
			if(R>=0) V.pb({L,R});
			L=l;R=r; continue;
		}
		R=max(R,r);
	}
	if(R>=0) V.pb({L,R});
	return V;
}

vector<ii> complement(vector<ii> vec, int n)
{
	if(vec.empty()) return {{1,n}};
	vector<ii> S;
	S.pb({1,vec[0].fi-1});
	for(int i=0;i+1<vec.size();i++)
	{
		S.pb({vec[i].se+1,vec[i+1].fi-1});
	}
	S.pb({vec.back().se+1,n});
	S=clean(S);
	return S;
}

int a[1111];
int b[1111];
vector<ii> A,B;
int n; 

void read()
{
	cin>>n; A.clear(); B.clear();
	int s1; cin>>s1;
	for(int i=0;i<s1;i++)
	{
		int l,r; cin>>l>>r;
		A.pb({l,r});
	}
	int s2; cin>>s2;
	for(int i=0;i<s2;i++)
	{
		int l,r; cin>>l>>r;
		B.pb({l,r});
	}
}

bool test(int T)
{
	for(ii x:A)
	{
		int l=x.fi; int r=x.se;
		l+=T; r+=T;
		if((l-1)/(2*T)!=(r/(2*T))) return false;
	}
	for(ii x:B)
	{
		int l=x.fi; int r=x.se;
		if((l-1)/(2*T)!=(r/(2*T))) return false;
	}
	return true;
}

const int C = 35000; //C*C>10^9
int solve_fast()
{
	A = complement(A,n);
	B = complement(B,n);
	int ans=0;
	for(int i=1;i<=min(C,n);i++)
	{
		if(test(i)) ans++;
	}
	if(n<=C) return ans;
	//only looking for numbers larger than C
	vector<ii> bad;
	for(ii x:A)
	{
		int l=x.fi; int r=x.se;
		for(int i=1;i<=C;i+=2)
		{
			int L = (l+i-1)/i;
			int R = r/i;
			L=max(L,C+1);
			if(L<=R) bad.pb({L,R});
		}
	}
	for(ii x:B)
	{
		int l=x.fi; int r=x.se;
		for(int i=2;i<=C;i+=2)
		{
			int L = (l+i-1)/i;
			int R = r/i;
			L=max(L,C+1);
			if(L<=R) bad.pb({L,R});
		}
	}
	bad = clean(bad);
	ans+=n-C;
	for(ii x:bad)
	{
		ans-=(x.se-x.fi+1);
	}
	return ans;
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	read();
	int ans = solve_fast();
	cout<<ans<<'\n';
}

There can be different approaches for this problem, which might work for various ranges of $$$X$$$, but I will mainly demonstrate my approach that works for all $$$X$$$ up to $$$7995843$$$.

Firstly, the theoretical upper bound of $$$X$$$ is $$$N^3$$$ for an $$$N \times N$$$ grid by counting the number of $$$1 \times k$$$ and $$$k \times 1$$$ strips. In our problem, $$$N \le 200$$$ and our constraint is $$$X \le 7995051 = N^3 - 4949$$$, which suggests that we need a tight solution.

Let’s fix a grid size, say $$$N \times M$$$ and see what we can do with this grid. Obviously, there are $$$NM$$$ $$$1 \times 1$$$ subrectangles, and we will subtract it from $$$X$$$, thereby ignoring $$$1 \times 1$$$ subrectangles from now on.

We should look at the 1-dimensional case first. Note that we can divide a 1D strip into chains of monotonic segments (that may overlap at a square). This enables us to do a simple dp to generate all values of $$$X’$$$ (the number of monotonic $$$1 \times k$$$ subrectangles of the 1D strip with $$$k \ge 2$$$) a 1D strip of size $$$N$$$ can generate. If you run the simple dp, you will find that you can actually generate most possible values of $$$X’$$$ up to the theoretical maximum (except for some values close to the maximum). Let $$$S_{N}$$$ denote the set of values (# of nontrivial horizontal substrips) you can generate using a 1D strip of size $$$N$$$. Then, our observation (by running our dp) is that $$$S_{N}$$$ contains all nonnegative integers less than $$$T$$$ for some $$$T$$$ close to $$$\frac{N(N-1)}{2}$$$.

Inspired by this observation, let’s suppose we don’t have to care about vertical strips first, and we have a target value $$$H$$$ of the number of nontrivial (area > 1) horizontal strips. We want to search for $$$N$$$ values in $$$S_{M}$$$ that sum up to $$$H$$$. Since our set $$$S_{M}$$$ contains almost all possible values, intuitively if we greedily choose the largest value in $$$S_{M}$$$ not exceeding our current target value of $$$H$$$ every time, then after $$$M$$$ iterations we will most likely get our answer. In fact, this greedy approach works very well and we will use this idea in our full solution.

To achieve all possible values of $$$X \le N^3 - 4949$$$ however, it is clear that we need to take vertical strips into account. The only hurdle is how to combine vertical and horizontal strips, since our task is somewhat easy if we can consider them independently. It turns out that there is a neat way to overcome this issue (there may also be many other ways but here is what I used): We add a suitable multiple of 1000 to each row, to adjust the pattern of the vertical strips. We also ensure that no two consecutive rows get added by the same multiple of 1000. Thus, all vertical strips will have the same pattern, and our monotonic chains cannot be separated by equal elements in between. This turns out to be not an issue, however, since we already have enough degrees of freedom. Let $$$S_{N}’$$$ denote the set of possible values generated by a $$$1 \times N$$$ strip where any two adjacent elements are not equal. Then, we want $$$X = NM + H + MV$$$, where $$$H \in S_{M}+S_{M}+...+S_{M}$$$ ($$$N$$$ times) and $$$V \in S_{N}’$$$. For a fixed pair $$$(N, M)$$$, we can greedily search from the largest valid value of $$$V$$$ (that still keeps the sum $$$\le X$$$), and greedily search for a valid combination that forms $$$H$$$ as described in the previous paragraph.

For $$$N = M = 200$$$, this approach can generate all sufficiently large $$$X \le N^3 - 4949$$$ (say $$$> 700000$$$), but what should we do for smaller values of $$$X$$$. Just test different sets of values $$$(N, M)$$$ depending on your range of $$$X$$$ (e.g. $$$(30, 30), (50, 50), (100, 100)$$$ and $$$N, M \le 20$$$ for small $$$X$$$). A neater solution would be to choose the grid size as $$$N \times 200$$$ where $$$N$$$ is minimal such that the theoretical upper bound of the grid size is at least $$$X$$$. There are a lot of flexibilities in choosing the grid size for smaller values of $$$X$$$, so any decent approach should work.

How to prove that all cases $$$X \le 7995001$$$ have a valid solution? Just test all such cases with a program >_<. The subtasks are chosen to (hopefully) reflect various levels of progress in this problem.

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
 
using namespace std;
using namespace __gnu_pbds;
 
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key
 
typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;

const int N = 200;
const int S = (N*(N-1))/2+1;
int dp[N+10][S+10];
int ty[N+10][S+10];
int dp2[N+10][S+10];
int a[N+10][N+10];

int c2(int x)
{
	return (x*(x-1))/2;
}

vi solve_dp1(int n, int S) //return a sequence of length n that works
{
	assert(dp[n][S]!=-1);
	int cur=n;
	vector<ii> sub;
	while(cur>0)
	{
		int las = dp[cur][S];
		if(las>0) sub.pb({las,ty[cur][S]});
		if(ty[cur][S])
		{
			S-=c2(cur-las+1);
		}
		else
		{
			S-=c2(cur-las);
		}
		cur=las;
	}
	reverse(sub.begin(),sub.end());
	vi vec;
	int ptr=0; int sgn=1;
	for(int i=1;i<n;i++)
	{
		if(ptr<sub.size()&&i==sub[ptr].fi) 
		{
			if(sub[ptr].se==0)
			{
				vec.pb(0); sgn*=-1; ptr++; continue;
			}
			else
			{
				sgn*=-1; ptr++;
			}
		}
		vec.pb(sgn);
	}
	vi V; int mn=0;
	int s=0; V.pb(s);
	for(int i=0;i<n-1;i++)
	{
		if(vec[i]==0) V.pb(s);
		else if(vec[i]==-1) V.pb(--s);
		else V.pb(++s);
		mn=min(mn,s);
	}
	for(int i=0;i<n;i++) V[i]-=mn;
	return V;
}

vi solve_dp2(int n, int S) //return a sequence of length n that works
{
	assert(dp2[n][S]!=-1);
	int cur=n;
	vi sub;
	while(cur>0)
	{
		int las = dp2[cur][S];
		if(las>0) S-=c2(cur-las+1);
		else S-=c2(cur-las);
		if(las>0) sub.pb(las); 
		cur=las;
	}
	reverse(sub.begin(),sub.end());
	vi vec;
	int ptr=0; int sgn=1;
	for(int i=1;i<n;i++)
	{
		if(ptr<sub.size()&&i==sub[ptr]) 
		{
			sgn*=-1; ptr++;
		}
		vec.pb(sgn);
	}
	vi V; int mn=0;
	int s=0; V.pb(s);
	for(int i=0;i<n-1;i++)
	{
		if(vec[i]==0) V.pb(s);
		else if(vec[i]==-1) V.pb(--s);
		else V.pb(++s);
		mn=min(mn,s);
	}
	for(int i=0;i<n;i++) V[i]-=mn;
	return V;
}

bool test(int s, int n, int m) //n*m + (S_m)^n + m*S''
{
	int T = s-n*m;
	if(T<0) return false;
	if(T==0) 
	{
		for(int i=0;i<n;i++)
		{
			for(int j=0;j<m;j++)
			{
				a[i][j]=1;
			}
		}
		return true;
	}
	int mx = T/m + 2;
	for(int i=mx;i>=n-1;i--)
	{
		if(dp2[n][i]==-1) continue;
		int T2 = T-m*i;
		if(T2<0) continue;
		vi V;
		for(int j=0;j<n;j++)
		{
			if(T2==0) break;
			for(int k=min(T2,c2(m));k>=0;k--)
			{
				if(dp[m][k]>=0)
				{
					T2-=k; 
					V.pb(k); 
					break;
				}
			}
		}
		if(T2==0) 
		{
			while(int(V.size())<n) V.pb(0);
			//the n rows follow these pattern
			vi rowcode = solve_dp2(n,i);
			const int C = 1000;
			for(int i=0;i<n;i++)
			{
				vi colcode = solve_dp1(m,V[i]);
				for(int j=0;j<m;j++)
				{
					a[i][j]=rowcode[i]*C+colcode[j];
				}
			}
			return true;
		}
	}
	return false;
}

void output(int i, int j)
{
	cout<<i<<' '<<j<<'\n';
	for(int r=0;r<i;r++)
	{
		for(int c=0;c<j;c++)
		{
			cout<<a[r][c];
			if(c+1<j) cout<<' ';
		}
		cout<<'\n';
	}
}

void solve(int x)
{
	if(x<=2000)
	{
		for(int i=1;i<=20;i++)
		{
			for(int j=1;j<=20;j++)
			{
				if(test(x,i,j))
				{
					output(i,j);
					return ;
				}
			}
		}
		assert(0);
	}
	if(x<=25000)
	{
		if(test(x,30,30)) 
		{
			output(30,30);
			return ;
		}
		if(test(x,50,50)) 
		{
			output(50,50);
			return ;
		}
		assert(0);
	}
	if(x<=700000)
	{
		if(test(x,100,100)) 
		{
			output(100,100); return ;
		}
		assert(0);
	}
	if(test(x,200,200)) 
	{
		output(200,200); return ;
	}
	assert(0);
}

void precompute()
{
	int n = N;
	memset(dp,-1,sizeof(dp));
	memset(dp2,-1,sizeof(dp2));
	dp[0][0]=1; dp2[0][0]=1;
	for(int i=1;i<=n;i++)
	{
		for(int j=0;j<i;j++)
		{
			for(int k=0;k<S;k++)
			{
				if(dp[j][k]==-1) continue;
				if(j>0) //last guy is still the same
				{
					dp[i][k+c2(i-j+1)]=j;
					ty[i][k+c2(i-j+1)]=1;
				}
				//last guy is ignored
				dp[i][k+c2(i-j)]=j;
				ty[i][k+c2(i-j)]=0;
			}
		}
	}
	for(int i=1;i<=n;i++)
	{
		for(int j=0;j<i;j++)
		{
			for(int k=0;k<S;k++)
			{
				if(dp2[j][k]==-1) continue;
				if(j>0) //last guy is still the same
				{
					dp2[i][k+c2(i-j+1)]=j;
				}			
				else
				{
					dp2[i][k+c2(i-j)]=j;
				}
			}
		}
	}
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	precompute();
	int t; cin>>t;
	while(t--) 
	{
		int x; cin>>x; 
		solve(x);
	}
}

#include <bits/stdc++.h>
using namespace std;
 
#define ll long long
#define ar array
 
const int M=19900;
bitset<M+1> dp[201];
int t, x;
 
void pa(int i, int o) {
	int k=200;
	while(k) {
		int j=1;
		while(!dp[k-j][i-j*(j-1)/2])
			++j;
		k-=j;
		i-=j*(j-1)/2;
		while(j--)
			cout << o++ << " ";
		--o;
	}
	cout << "\n";
}
 
void solve() {
	cin >> x;
	if(x<200) {
		cout << "1 " << x << "\n";
		for(int i=0; i<x; ++i)
			cout << 0 << " \n"[i==x-1];
		return;
	}
	int n=1;
	while(n*(n+1)/2*200+n*200*199/2-4949<x)
		++n;
	cout << n << " 200\n";
	x-=n*(n+1)/2*200;
	for(int i=0; i<n; ++i) {
		int j=min(M, x);
		while(!dp[200][j])
			--j;
		pa(j, i*200);
		x-=j;
	}
}
 
int main() {
	ios::sync_with_stdio(0);
	cin.tie(0);
 
	dp[0][0]=1;
	for(int i=1; i<=200; ++i)
		for(int j=1; j<=i; ++j)
			dp[i]|=dp[i-j]<<j*(j-1)/2;
	cin >> t;
	while(t--)
		solve();
}

Call a position in the string bad if there exists a way for Nibutani to place one of the letters ‘L’, ‘O’, ‘V’, ‘E’ at that position such that a new substring of the form “LOVE” is formed. The intuitive idea is to always play at bad positions if possible (Player 1 should use it to form LOVE while Player 2 should block it). Indeed, this immediately gives the lower bound of ceil(# of bad positions/2) as the number of additional LOVE strings that can be formed by Player 1, since Player 2 can only block at most one bad position per turn. Hence, as Player 1, your strategy is simple: Place a character at a bad position which will form the substring LOVE while at least one bad position remains.

Now, we prove that Player 2 can always prevent Player 1 from forming more LOVE substrings. The key idea is that Player 2 can always place a character at any position (whether it is bad or not) without forming a new substring LOVE or creating a new bad position. You can prove this using casework, noting that there are always 4 possible choices to place a new character but less than 4 characters are effectively banned. Hence, your program can just iterate over all characters and check naively which character works.

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
 
using namespace std;
using namespace __gnu_pbds;
 
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key
 
typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;

const int N = 200000;
string s; 
string love = "LOVE";
map<char,int> ma;
int n;
set<int> bad; //set of bad positions
set<int> emp;
const int BAD = int(1e9)+2;

void check(int x) //check if the substring starting from x contributes to a bad position
{
	if(x+int(love.length())-1>=n) return ;
	int cnt=0; int id=-1;
	for(int i=0;i<love.length();i++)
	{
		if(s[x+i]=='?')
		{
			cnt++; id=x+i;
		}
		else if(ma[s[x+i]]!=i) return ;
	}
	if(cnt==1) bad.insert(id);
}

int check2(int x)
{
	if(x+int(love.length())-1>=n) return int(1e9);
	int cnt=0; int id=-1;
	for(int i=0;i<love.length();i++)
	{
		if(s[x+i]=='?')
		{
			cnt++; id=x+i;
		}
		else if(ma[s[x+i]]!=i) return int(1e9);
	}
	if(cnt==0) return BAD;
	return cnt;
}

void fill(int x) //position x was just filled, now update
{
	for(int i=-3;i<=0;i++)
	{
		if(x+i>=0) check(x+i);
	}
}

void out(int x)
{
	cout<<x+1<<' '<<s[x]<<'\n'; fflush(stdout);
}

void place(int pl, int x)
{
	if(pl==1)
	{
		int mncnt=int(1e9)+10; int stpos=-1;
		for(int i=max(0,x-3);i<=x;i++)
		{
			int res = check2(i);
			if(res<mncnt)
			{
				mncnt=res; stpos=i;
			}
		}
		s[x]=love[x-stpos];
	}
	else
	{
		//make sure I don't create any extra problems (e.g. LO??VE)
		for(int i=0;i<4;i++)
		{
			s[x]=love[i];
			if(x-i>=0&&(check2(x-i)==1||check2(x-i)==BAD))
			{
				continue;
			}
			break;
		}
	}
	out(x);
	fill(x);
}

void play(int pl)
{
	if(bad.empty())
	{
		int x=(*emp.begin()); emp.erase(x);
		place(pl,x); return ;
	}
	int x=(*bad.begin());
	emp.erase(x); bad.erase(x); place(pl,x); return ;
}

void getmove()
{
	int x; cin>>x; x--; 
	char c; cin>>c;
	s[x]=c; emp.erase(x); bad.erase(x);
	fill(x);
}

int main()
{
	ma['L']=0; ma['O']=1; ma['V']=2; ma['E']=3;
	cin>>s; n=s.length();
	int cnt=0;
	for(int i=0;i<n;i++) 
	{
		if(s[i]=='?')
		{
			cnt++; emp.insert(i);
		}
	}
	for(int i=0;i+3<n;i++)
	{
		check(i);
	}
	int pl; cin>>pl; //pl = 1 or 2
	pl%=2;
	int curp = 1;
	for(int i=0;i<cnt;i++)
	{
		if(curp!=pl)
		{
			getmove();
		}
		else
		{
			play(pl);
		}
		curp^=1;
	}
}

Let $$$B_{i} = b_{1} + b_{2} + … + b_{i}, D_{i} = (a_1 - b_1) + (a_2 - b_2) + … + (a_i - b_i)$$$.

Observe that the desired sum can be rewritten as $$$\displaystyle\sum_{i = l}^{r}\max(D_{i} - D_{l-1}, 0) + (B_{i} - B_{l-1}) = \displaystyle\sum_{i = l}^{r}\max(D_{i}, D_{l-1}) + (B_{i} - B_{l-1} - D_{l-1})$$$. Note that the $$$B_{i} - B_{l-1} - D_{l-1}$$$ part can be computed in $$$O(1)$$$ with a prefix sum, so the hard part is computing $$$\displaystyle\sum_{i = l}^{r}\max(D_{i}, D_{l-1})$$$.

Basically, our problem reduces to answering the following queries: Given $$$r, V$$$, find $$$\displaystyle\sum_{i=1}^{r}\max(D_{i}, V)$$$(since the case of general $$$l$$$ can be reduced to the difference between two queries of this type). This subproblem can be solved in many ways. For example, we can abuse the fact that the given queries are offline, and maintain some Fenwick trees on the sorted order of $$$D_{i}$$$, and process queries in increasing order of $$$r$$$. The offline solution works in $$$O((N+Q)\log N)$$$.

Subtask 1 can be solved with pure brute force. There exists a $$$O(N\log{N} + Q\log^{2}N)$$$ solution using binary search on segment tree of sorted vectors which should solve up to Subtask 4.

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
 
using namespace std;
using namespace __gnu_pbds;
 
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key
 
typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;

const int N = 500000;
const int Q = 500000;

ll a[N+10];
ll b[N+10];
int n,q; 
ll ans[Q+10];
ii queries[Q+10];
ll d[N+10];
ll bsum[N+10];
ll wbsum[N+10];

void read()
{
	scanf("%d %d",&n,&q);
	for(int i=0;i<n;i++)
	{
		scanf("%lld",a+i);
	}
	for(int i=0;i<n;i++)
	{
		scanf("%lld",b+i);
	}
	for(int i=0;i<q;i++)
	{
		int l,r; scanf("%d %d",&l,&r);
		l--; r--;
		queries[i]={l,r};
	}
}

ll sumb(int l, int r)
{
	if(l==0) return bsum[r];
	return bsum[r]-bsum[l-1];
}

ll sumwb(int l, int r)
{
	if(l==0) return wbsum[r]-(n-r-1)*1LL*sumb(l,r);
	return wbsum[r]-wbsum[l-1]-(n-r-1)*1LL*sumb(l,r);
}

struct Fenwick
{
	vector<ll> t;
    Fenwick(int n)
    {
        t.assign(n+1,0);
    }
    void reset(int n)
    {
		t.assign(n+1, 0);
	}
    void update(int p, ll v)
    {
        for (; p < (int)t.size(); p += (p&(-p))) t[p] += v;
    }
    ll query(int r) //finds [1, r] sum
    {                     
        ll sum = 0;
        for (; r; r -= (r&(-r))) sum += t[r];
        return sum;
    }
    ll query(int l, int r) //finds [l, r] sum
    {
		if(l == 0) return query(r);
		return query(r) - query(l-1);
	}
};

ll dsum[N+11];
vector<pair<pair<ll,int> ,int> > qs[N+11];

int main()
{
	read();
	for(int i=0;i<n;i++) d[i]=a[i]-b[i];
	vector<ii> dsrt;
	for(int i=0;i<n;i++)
	{
		bsum[i]=b[i];
		dsum[i]=a[i]-b[i];
		wbsum[i]=(n-i)*1LL*b[i];
		if(i>0) 
		{
			wbsum[i]+=wbsum[i-1];
			bsum[i]+=bsum[i-1];
			dsum[i]+=dsum[i-1];
		}
		dsrt.pb({dsum[i],i});
	}
	dsrt.pb({0,-1});
	for(int i=0;i<q;i++)
	{
		int l=queries[i].fi; int r=queries[i].se;
		ans[i]+=sumwb(l,r);
		if(l>0) ans[i]-=(r-l+1)*1LL*dsum[l-1];
		qs[l-1].pb({{r,i},1});
		if(l>0) qs[l-1].pb({{l-1,i},-1});
	}
	Fenwick f1(n+10); Fenwick f2(n+10);
	sort(dsrt.rbegin(),dsrt.rend());
	for(ii x:dsrt)
	{
		ll D = x.fi; int id = x.se;
		for(auto query:qs[id])
		{
			int sgn = query.se;
			int til = query.fi.fi;
			int lab = query.fi.se;
			int cnt = til+1;
			int cntlarger = f1.query(til+1);
			ans[lab]+=sgn*((cnt-cntlarger)*1LL*D+f2.query(til+1));
		}
		if(id>=0) 
		{
			f1.update(id+1,1);
			f2.update(id+1,D);
		}
	}
	for(int i=0;i<q;i++)
	{
		cout<<ans[i]<<'\n';
	}
}

Author : zscoder

This is a simple implementation problem. We iterate through all banknotes one by one and check if Oleg can take each of them. If a banknote is at position x, then Oleg can take it if and only if b < x < c. This can be checked in O(1) time. Thus, the total complexity is O(n). Note that the information on the starting position of Oleg is useless here.

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<ll> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	int a, b, c; cin>>a>>b>>c;
	int n; cin>>n;
	int ans=0;
	for(int i=0;i<n;i++)
	{
		int x; cin>>x;
		if(x>b&&x<c) ans++;
	}
	cout<<ans<<'\n';
}

Author : zscoder

Let's find the value of x_i explicitly. Suppose we make the i-th cut and distance x_i from the apex. Then, the ratio of similitude of the isosceles triangle with apex equal to the apex of the carrot and the base equal to the i-th cut and the whole carrot is . Since the area of this smaller isosceles triangle is the sum of areas of the first i pieces, which is of the whole carrot. Thus, , which is equivalent to . Thus, and we can find each x_i in O(1) time. The total complexity is O(n).

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<ll> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	ll n, h; cin>>n>>h;
	for(int i=1;i<=n-1;i++)
	{
		cout<<fixed<<setprecision(12)<<sqrt(ld(i)/ld(n))*ld(h);
		if(i<n-1) cout<<' ';
	}	
	cout<<'\n';
}

Author : zscoder

First, it is clear that Oleg will place letters and Igor will place letters. Next, it is clear that Oleg and Igor will both choose their smallest and biggest letters respectively to place in the final string. Thus, we now consider that Oleg places his smallest letters and Igor places his largest letters.

Consider the following greedy strategy. When it's Oleg's turn, he will replace the frontmost question mark with his smallest letter. When it's Igor's turn, he will replace the frontmost question mark with his largest letter. At first glance, you might think that this works. However, there's another case that we haven't considered.

Suppose Oleg has the letters {x, y, z} and Igor has the letters {a, b, c}. According to our previous strategy, Oleg will place x as the first letter. However, that's not optimal. He can place his letters at the back and force Igor to place the first letter. The reason is because the largest letter of Igor is not larger than the smallest letter of Oleg. Thus, it is beneficial for Oleg to place his letters at the back and force Igor to place his letters in front.

So, what exactly will the final string look like? We'll look at the moves one by one. If at some point Oleg's smallest letter is still strictly smaller than Igor's largest letter, then both player must put their smallest (largest if it's Igor) letter as the frontmost letter. Why? Suppose not, then on the next turn the other player will occupy that spot with their best (smallest if Oleg, largest if Igor) letter, and the resulting string will be worse for the current player. This proves that greedy is correct in this case.

Now, what if Oleg's smallest letter is not smaller than Igor's largest letter. In this case, both players will want to force the other player to place their own letter at the beginning of the string. It can be proven that in this case, each person will place their current worst (largest if Oleg, smallest if Igor) letter at the back of the string in the optimal strategy. Thus, we can calculate the final string starting from this point and after that reverse this part and combine it with the first part of the string where both players greedily place their best letters in the beginning.

Time Complexity : O(n)

Many people failed on pretest 6 initially because they didn't consider the second case.

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<ll> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	string ansl;
	string ansr;
	string s,t;
	cin>>s>>t;
	sort(s.begin(),s.end());
	sort(t.begin(),t.end());
	reverse(t.begin(),t.end());
	int n = s.length();
	deque<char> a,b;
	for(int i=0;i<(n+1)/2;i++)
	{
		a.pb(s[i]);
	}
	for(int i=0;i<n/2;i++)
	{
		b.pb(t[i]);
	}
	bool mode=0;
	for(int i=0;i<n;i++)
	{
		if(i&1)
		{
			if(!a.empty()&&a[0]>=b[0])
			{
				mode=1;
			}
			if(mode)
			{
				ansr+=b.back();
				b.pop_back();
			}
			else
			{
				ansl+=b[0];
				b.pop_front();
			}
		}
		else //P1's turn
		{
			if(!b.empty()&&a[0]>=b[0])
			{
				mode=1;
			}
			if(mode)
			{
				ansr+=a.back();
				a.pop_back();
			}
			else
			{
				ansl+=a[0];
				a.pop_front();
			}
		}
	}
	reverse(ansr.begin(),ansr.end());
	ansl+=ansr;
	cout<<ansl<<'\n';
}

Author : AnonymousBunny

Add each vertex to its own adjacency list. Now, we claim that if it is possible to label the cities to satisfy the problem conditions, then it is possible to do so so that for every two cities with the same adjacency list, they're labelled with the same number.

Indeed, if they have the same adjacency list, they must be neighbours. Thus, the difference between their labels is at most 1. Suppose we label the first vertex u with number i and the second vertex v with the number i + 1. Note that since their adjacency lists are equal, a vertex x is a neighbour of u iff it's a neighbour of v. Thus, u and v can't have neighbours with labels i - 1 or i + 2, or else it will contradict the condition. Thus, all neighbours of u and v have labels i or i + 1. Thus, we can safely change the label of the second vertex v to i and the conditions will still hold.

Thus, we can sort the set of adjacency lists of each vertex, and then group the vertices with the same adjacency list together. Suppose there are k such groups. For simplicity, we can create a new graph where each group represent a vertex of the new graph. Connect two groups i and j if and only if there exist some vertex in group i that connects to a vertex in group j. Note that the graph will have at most O(m) edges. Now, if a vertex has degree  ≥ 3, we can't assign a number to that vertex properly, as one of its neighbours will not have a label which have a difference  ≤ 1 from it. Thus, all vertices in the new graph must have degree  ≤ 2. Since it's connected, it must be either a cycle or a path. However, it can be easily seen that there is no labelling if it's a cycle. Thus, it must be a path. Now, we can just assign the labels to the graph from one end of the path to the other end by the numbers 1 to k. Finally, the label of a vertex is simply the label of its group.

This solution can be implemented in time.

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<int,int> ii;
typedef vector<ll> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int N = 300000;

pair<vi,int> adj[N+2];
int ans[N+2];
int lab[N+2];
int lab2[N+2];
set<int> adj2[N+2];
vector<ii> edges;

struct DSU
{
	int S;
	
	struct node
	{
		int p; ll sum;
	};
	vector<node> dsu;
	
	DSU(int n)
	{
		S = n;
		for(int i = 0; i < n; i++)
		{
			node tmp;
			tmp.p = i; tmp.sum = 0;
			dsu.pb(tmp);
		}
	}
	
	void reset(int n)
	{
		dsu.clear();
		S = n;
		for(int i = 0; i < n; i++)
		{
			node tmp;
			tmp.p = i; tmp.sum = 0;
			dsu.pb(tmp);
		}
	}
	
	int rt(int u)
	{
		if(dsu[u].p == u) return u;
		dsu[u].p = rt(dsu[u].p);
		return dsu[u].p;
	}
	
	void merge(int u, int v)
	{
		u = rt(u); v = rt(v);
		if(u == v) return ;
		if(rand()&1) swap(u, v);
		dsu[v].p = u;
		dsu[u].sum += dsu[v].sum;
	}
	
	bool sameset(int u, int v)
	{
		if(rt(u) == rt(v)) return true;
		return false;
	}
	
	ll getstat(int u)
	{
		return dsu[rt(u)].sum;
	}
};

deque<int> chain;

void dfs(int u, int p, bool type)
{
	if(type) chain.pb(u);
	else chain.push_front(u);
	int c=0;
	for(sit it = adj2[u].begin(); it != adj2[u].end(); it++)
	{
		int v = (*it);
		if(v==p) continue;
		if(p!=-1)
		{
			dfs(v,u,type);
		}
		else
		{
			dfs(v,u,c);
			c++;
		}
	}
}

int main()
{
	//ios_base::sync_with_stdio(0); cin.tie(0);
	int n, m; scanf("%d %d", &n, &m);
	for(int i = 0; i < m; i++)
	{
		int u, v; scanf("%d %d", &u, &v);
		u--; v--;
		adj[u].fi.pb(v);
		adj[v].fi.pb(u);
		edges.pb(mp(u,v));
	}
	for(int i = 0; i < n; i++)
	{
		adj[i].fi.pb(i);
		adj[i].se = i;
		sort(adj[i].fi.begin(),adj[i].fi.end());
	}
	sort(adj,adj+n);
	int cnt = 1;
	for(int i = 0; i < n; i++)
	{
		if(i==0) 
		{
			lab[adj[i].se] = cnt;
		}
		else
		{
			if(adj[i].fi==adj[i-1].fi)
			{
				lab[adj[i].se]=cnt;
			}
			else
			{
				lab[adj[i].se]=++cnt;
			}
		}
	}
	if(cnt==1)
	{
		printf("YES\n");
		for(int i = 0; i < n; i++)
		{
			printf("%d ",lab[i]);
		}
		printf("\n");
		return 0;
	}
	DSU dsu(cnt+1);
	for(int i = 0; i < m; i++)
	{
		int u = edges[i].fi; int v = edges[i].se;
		if(lab[u]!=lab[v])
		{
			adj2[lab[u]].insert(lab[v]);
			adj2[lab[v]].insert(lab[u]);
			dsu.merge(lab[u],lab[v]);
		}
	}
	bool pos = 1;
	for(int i = 1; i <= cnt; i++)
	{
		if(dsu.rt(i)!=dsu.rt(1))
		{
			pos=0;
			break;
		}
	}
	if(!pos)
	{
		printf("NO\n");
		return 0;
	}
	int d1 = 0;
	for(int i = 1; i <= cnt; i++)
	{
		if(adj2[i].size()>2)
		{
			printf("NO\n");
			return 0;
		}
		if(adj2[i].size()==1) d1++;
		else assert(adj2[i].size()==2);
	}
	if(d1==2)
	{
		printf("YES\n");
		dfs(1,-1,0);
		for(int i = 0; i < chain.size(); i++)
		{
			lab2[chain[i]] = i+1;
		}
		for(int i = 0; i < n; i++)
		{
			printf("%d ",lab2[lab[i]]);
		}
		printf("\n");
	}
	else
	{
		printf("NO\n");
		return 0;
	}
}

Author : zscoder

First, we solve the problem when no one has any extra turns.

Suppose we're binary searching the answer. Let all the numbers  ≥ x be equal to 1 and all the numbers  < x be equal to 0. Both players can remove one number from one end of the row. The goal of the first player is to let the remaining number be 1 and the goal of the second player is to leave 0 in the end. If the first player can win, this means that the answer is at least x. Thus, we first try to solve this simpler problem.

We claim that the first player wins if and only if :

1. n is even and one of the two middle numbers is 1.

2. n is odd, the middle digit is 1 and at least one of the digits beside the middle digit is 1 (unless n = 1, for which first players wins when the only carrot is labelled 1)

Indeed, once we deduce this, we can easily prove this by induction on n. The proof is just doing casework and considering all possible moves.

Once we have this fact, we realize we don't actually have to binary search the answer. If n is even, the answer is while if n ≥ 3 is odd, the answer is . (If n = 1 then the answer is obviously a₁.)

Now, we have to take extra moves into account. Fortunately, it's not very difficult. Having k extra moves just means that Player 1 can choose to start the game in any subsegment of length n - k. Thus, we just have to compute the maximum answer for all subsegments of length n - k for all 0 ≤ k ≤ n - 1. With the formula above, you can find all the answers in O(n) time or even time if you use sparse table for range maximum query.

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<ll> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

int a[300001];
int ans[300001];
int b[300001];

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	int n; cin>>n;
	int mx=0;
	for(int i=0;i<n;i++) 
	{
		cin>>a[i];
		mx=max(mx,a[i]);
	}
	for(int i=0;i<n-2;i++)
	{
		b[i]=min(a[i+1],max(a[i],a[i+2]));
	}
	ans[n-1]=mx;
	int odd=n;int even=n;
	if(n&1) even=n-1;
	else odd=n-1;
	mx=0;
	for(int i=even;i>=2;i-=2)
	{
		int l = (i-1)/2; int r=n-i/2;
		assert(l<=r);
		mx=max(mx,max(a[l],a[r]));
		if(i==even)
		{
			assert(r-l<=2);
			if(r-l==2)
			{
				mx=max(mx,a[l+1]);
			}
		}
		ans[n-i]=mx;
	}
	mx=0;
	for(int i=odd;i>=3;i-=2)
	{
		int l = i/2-1; int r=n-2-i/2;
		assert(l<=r);
		if(i==odd) assert(r-l<=1);
		mx=max(mx,max(b[l],b[r]));
		ans[n-i]=mx;
	}
	for(int i=0;i<n;i++)
	{
		cout<<ans[i];
		if(i<n-1) cout<<' ';
	}
	cout<<'\n';
}

Author : hloya_ygrt

We use a segment tree to solve this problem. For each node, it is sufficient to store two arrays : sum[i], denoting the total contribution of the digit i in the current segment (if a digit is in the tens digit then it contributes 10 to the sum and etc...), and also nxt[i], what all the digits i in the current segment are changed to.

Maintaining these arrays is quite straightforward with lazy propogation. When we push an update down a node, we need to update the nxt array of the children. First, we change st[id].nxt[u] to v, where the current update is to change all digits u to v. Then, we change st[id * 2].nxt[i] to st[id].nxt[st[id * 2].nxt[i]], where st[id] is the current node and st[id * 2] is one of the children nodes. (Do the same for the right children). You can see the code if you need more details. Finally, update the sum array of the current segment.

The total complexity of the code is , which is fast enough.

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;
const int N = 100000;
int a[N+1][10];

struct node
{
	int nxt[10];
	ll sum[10];
};

node st[N*4+1];

void combine(int id)
{
	for(int i = 0; i < 10; i++)
	{
		st[id].sum[i]=st[id*2].sum[i]+st[id*2+1].sum[i];
	}
}

void build(int id, int l, int r)
{
	if(r-l<2)
	{
		for(int i = 0; i < 10; i++) st[id].sum[i]=a[l][i];
		for(int i = 0; i < 10; i++)
		{
			st[id].nxt[i] = i;
		}
		return ;
	}
	for(int i = 0; i < 10; i++)
	{
		st[id].nxt[i] = i;
	}
	int mid=(l+r)>>1;
	build(id*2,l,mid);
	build(id*2+1,mid,r);
	combine(id);
}

int nxt1[10];
int nxt2[10];
ll sum[10];

void push(int id, int l, int r)
{
	memset(sum,0,sizeof(sum));
	if(r-l>=2)
	{
		for(int i = 0; i < 10; i++)
		{
			nxt1[i] = st[id].nxt[st[id*2].nxt[i]];
			nxt2[i] = st[id].nxt[st[id*2+1].nxt[i]];
		}
		for(int i=0;i<10;i++)
		{
			st[id*2].nxt[i]=nxt1[i];
			st[id*2+1].nxt[i]=nxt2[i];
		}
	}
	for(int i=0;i<10;i++)
	{
		sum[st[id].nxt[i]]+=st[id].sum[i];
	}
	for(int i=0;i<10;i++) 
	{
		st[id].sum[i]=sum[i];
		st[id].nxt[i]=i;
	}
}

void update(int id, int l, int r, int ql, int qr, int u, int v)
{
	push(id,l,r);
	if(ql>=r||l>=qr) return ;
	if(ql<=l&&r<=qr)
	{
		st[id].nxt[u]=v;
		push(id,l,r);
		return ;
	}
	int mid=(l+r)>>1;
	update(id*2,l,mid,ql,qr,u,v);
	update(id*2+1,mid,r,ql,qr,u,v);
	combine(id);
}

ll query(int id, int l, int r, int ql, int qr)
{
	push(id,l,r);
	if(ql>=r||l>=qr) return 0;
	if(ql<=l&&r<=qr)
	{
		ll sum=0;
		for(int i=1;i<10;i++)
		{
			sum+=ll(i)*st[id].sum[i];
		}
		return sum;
	}
	int mid=(l+r)>>1;
	return (query(id*2,l,mid,ql,qr)+query(id*2+1,mid,r,ql,qr));
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	int n, q; cin>>n>>q;
	for(int i = 0; i < n; i++)
	{
		int x; cin>>x;
		int cur=1;
		for(int j = 0; j < 9; j++)
		{
			a[i][x%10] += cur;
			x/=10;
			cur*=10;
			if(x==0) break;
		}
	}
	build(1,0,n);
	for(int i = 0; i < q; i++)
	{
		int type;
		cin>>type;
		if(type==1)
		{
			int l,r,u,v;
			cin>>l>>r>>u>>v;
			l--; r--;
			update(1,0,n,l,r+1,u,v);
		}
		else
		{
			int l,r; cin>>l>>r;
			l--; r--;
			ll sum = query(1,0,n,l,r+1);
			cout<<sum<<'\n';
		}
	}
}

Author : zscoder

First, we solve the problem when there're no question marks, i.e. we find a way to calculate the number of good pairs of strings fast for a constant pair of strings A and B.

Call a pair of strings (S, T) where |S| ≤ |T| coprime if S = T or S is a prefix of T, and if T = S + X, then (X, S) is also coprime. (S, T) where |S| > |T| is coprime iff (T, S) is coprime.

If A = B, then all possible strings work. Thus, we assume A ≠ B from now on. We remove the longest common prefix of A and B. Thus, we can assume A[0] ≠ B[0]. Thus, either S is a prefix of T or T is a prefix of S. WLOG, S is a prefix of T. Let T = S + X. Now, A and B consists of only S and X. Using this, we can prove by induction on |S| + |T| that S and T must be coprime.

One important property of coprime strings is that S + T = T + S holds. (again induction works here)

Now, since the strings S and T needs to be coprime, we have S + T = T + S. This allows us to swap any neighbouring Ss and Ts (or 'A's and 'B's) in A and B, as the resulting strings will still be equal. Thus, swapping repeatedly allows us to sort the strings A and B. (the 'A's appear in front and 'B's appear at the back) Let x_A, x_B, y_A, y_B denote the number of As and Bs in the first string and second string respectively. If (x_A, x_B) > (y_A, y_B), then the answer is 0. We'll handle the case (x_A, x_B) = (y_A, y_B) later. Now, assume x_A > y_A, x_B < y_B. Thus, we have to solve the equation (x_A - y_A) copies of S = (y_B - x_B) copies of T.

Now, let x = x_A - y_A, y = y_B - x_B. If x = y, then the solution is S = T. Otherwise, assume x > y. Then, |S| < |T|. So, by comparing, we again have T = S + X, for some nonempty binary string X. Note that S and X must be coprime too, so we can sort the second string as well. We cancel off the Ss on both sides to get (x - y)S = yX. Thus, this means that if (S, T) is a solution for (x, y), then (S, X) is a solution for (x - y, y). Note that repeating this process will eventually lead us to (1, 1). (this process is similar to Euclidean Algorithm)

The answer for (1, 1) is the number of solutions to S = T. Let's denote the solution here as X. Doing some backtracking, we realize that the answer for (x, y) is equal to (X....X (y times), X...X (x times)). Note that we still have the condition |S|, |T| ≤ N, so we can translate this to an appropriate condition on the length of X and the answer is simply the number of binary strings of length not exceeding the maximum possible length of X.

The only case that remains is that (x_A, x_B) = (y_A, y_B). In this case, any pair of coprime strings S and T will work. Thus, our task reduces to calculating the number of coprime pair of strings with length not exceeding N.

We claim that the number of coprime pair of strings (S, T) with |S| = p, |T| = q is .

If p = q the claim is obviously true. Otherwise, we can induct on p + q agin. If q > p, we can write T = S + X and then the number of coprime pairs of (S, T) is equal to the number of coprime pairs of (S, X), which by induction is equal to . This proves the claim.

Thus, we just need to compute the sum of for all 1 ≤ p, q ≤ N.

Indeed, since N ≤ 3·10⁵, it is enough to count the number of pairs (p, q) with gcd = g for all g.

However, this is quite easy. Let cnt[i] denote the number of pairs (p, q) such that p and q are both divisible by i. Let ans[i] denote the number of pairs (p, q) with gcd = i. Then, . Thus, this can be computed in .

Now, we need to find out how to calculate the sum of all these values on two strings X and Y with question marks. Handle the case when the two strings become equal separately.

Let's first make a summary of the number of good pairs of strings for constant strings A and B. In fact, note that the formulaes above only depends on (d_A, d_B), the difference between the number of As in A and B, and the difference between the number of Bs in A and B (note that d_A, d_B can be negative)

If d_A = d_B = 0, then the answer is the sum of for all 1 ≤ p, q ≤ N, which as we have just saw can be precomputed in time.

Otherwise, if d_A, d_B ≥ 0 or d_A, d_B ≤ 0, then there are no good pair of strings.

Finally, in other cases, let p = |d_A|, q = |d_B|. Then, the answer is .

This also means that we can compute the answer if we know d_A and d_B very fast. (worst case is )

Now, suppose in the strings X and Y, we have a and b question marks respectively. Additionally, suppose the current difference between the number of As and Bs of these strings is (p, q).

If we choose x and y of the question marks from X and Y to be replaced with As, then the difference between As and Bs in the strings become (p + x - y, q + (a - b) - (x - y)). Let's denote q as q + a - b for simplicity. Thus, the difference is now written as (p + (x - y), q - (x - y)). The values of x and y can be any integer in the range [0, a] and [0, b] respectively. Suppose for all  - b ≤ d ≤ a, we know how many ways to assign the question marks have x - y = d. Then, we can iterate through all the ds one by one and compute the answer fast for each d.

Thus, the final hurdle is to calculate the number of ways to obtain x - y = d for all possible d so that 0 ≤ x ≤ a, 0 ≤ y ≤ b. This is just the sum of for all 0 ≤ x ≤ a. However, this is equal to , as the number of ways to choose b + d objects from a + b objects is the same as the sum of the product of the number of ways to choose x objects from the first a objects and the number of ways to choose b + d - x objects from the first b objects for all 0 ≤ a ≤ x. Thus, this value can be computed in O(1) with precomputed factorials and inverse factorials (or you can maintain this value when we iterate through all d).

Finally, don't forget to take care of the cases where it is possible for both strings to be equal.

The time complexity of the solution is .

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<ll> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int MOD = 1e9 + 7;

struct NumberTheory
{
	vector<ll> primes;
	vector<bool> prime;
	vector<ll> totient;
	vector<ll> sumdiv;
	vector<ll> bigdiv;
	void Sieve(ll n)
	{
		prime.assign(n+1, 1);
		prime[1] = false;
		for(ll i = 2; i <= n; i++)
		{
			if(prime[i])
			{
				primes.pb(i);
				for(ll j = i*2; j <= n; j += i)
				{
					prime[j] = false;
				}
			}
		}
	}
	
	ll phi(ll x)
	{
		map<ll,ll> pf;
		ll num = 1; ll num2 = x;
		for(ll i = 0; primes[i]*primes[i] <= x; i++)
		{
			if(x%primes[i]==0)
			{
				num2/=primes[i];
				num*=(primes[i]-1);
			}
			while(x%primes[i]==0)
			{
				x/=primes[i];
				pf[primes[i]]++;
			}
		}
		if(x>1)
		{
			pf[x]++; num2/=x; num*=(x-1);
		}
		x = 1;
		num*=num2;
		return num;
	}
	
	bool isprime(ll x)
	{
		if(x==1) return false;
		for(ll i = 0; primes[i]*primes[i] <= x; i++)
		{
			if(x%primes[i]==0) return false;
		}
		return true;
	}

	void SievePhi(ll n)
	{
		totient.resize(n+1);
		for (int i = 1; i <= n; ++i) totient[i] = i;
		for (int i = 2; i <= n; ++i)
		{
			if (totient[i] == i)
			{
				for (int j = i; j <= n; j += i)
				{
					totient[j] -= totient[j] / i;
				}
			}
		}
	}
	
	void SieveSumDiv(ll n)
	{
		sumdiv.resize(n+1);
		for(int i = 1; i <= n; ++i)
		{
			for(int j = i; j <= n; j += i)
			{
				sumdiv[j] += i;
			}
		}
	}
	
	ll getPhi(ll n)
	{
		return totient[n];
	}
	
	ll getSumDiv(ll n)
	{
		return sumdiv[n];
	}
	
	ll modpow(ll a, ll b, ll mod)
	{
		ll r = 1;
		if(b < 0) b += mod*100000LL;
		while(b)
		{
			if(b&1) r = (r*a)%mod;
			a = (a*a)%mod;
			b>>=1;
		}
		return r;
	}
	
	ll inv(ll a, ll mod)
	{
		return modpow(a, mod - 2, mod);
	}
	
	ll invgeneral(ll a, ll mod)
	{
		ll ph = phi(mod);
		ph--;
		return modpow(a, ph, mod);
	}
	
	void getpf(vector<ii>& pf, ll n)
	{
		for(ll i = 0; primes[i]*primes[i] <= n; i++)
		{
			int cnt = 0;
			while(n%primes[i]==0)
			{
				n/=primes[i]; cnt++;
			}
			if(cnt>0) pf.pb(ii(primes[i], cnt));
		}
		if(n>1)
		{
			pf.pb(ii(n, 1));
		}
	}

	//ll op;
	void getDiv(vector<ll>& div, vector<ii>& pf, ll n, int i)
	{
		//op++;
		ll x, k;
		if(i >= pf.size()) return ;
		x = n;
		for(k = 0; k <= pf[i].se; k++)
		{
			if(i==int(pf.size())-1) div.pb(x);
			getDiv(div, pf, x, i + 1);
			x *= pf[i].fi;
		}
	}
};

NumberTheory nt;

ll modpow(ll a, ll b)
{
	ll r = 1;
	while(b)
	{
		if(b&1) r=(r*a)%MOD;
		a=(a*a)%MOD;
		b>>=1;
	}
	return r;
}

ll inv(ll a)
{
	return modpow(a,MOD-2);
}

ll n;
ll cnt[300001];
ll mob[300001];

ll mobius(ll x)
{
	int cc = 0;
	for(int i=0;nt.primes[i]*nt.primes[i]<=x;i++)
	{
		int z=0;
		while(x%nt.primes[i]==0)
		{
			z++;
			x/=nt.primes[i];
		}
		if(z>=2) return 0;
		if(z>0) cc++;
	}
	if(x>1) cc++;
	if(cc&1) return -1;
	else return 1;
}

ll solve(ll x, ll y)
{
	if(x==0&&y==0)
	{
		for(int i=1;i<=n;i++)
		{
			cnt[i]=ll(n/i)*ll(n/i);
		}
		for(int i=1;i<=n;i++)
		{
			for(int j=2*i;j<=n;j+=i)
			{
				cnt[i]+=mob[j/i]*cnt[j];
			}
		}	
		ll ans = 0;
		ll cur = 2;
		for(int i=1;i<=n;i++)
		{
			cnt[i]%=MOD;
			if(cnt[i]<0) cnt[i]+=MOD;
			//cerr<<i<<' '<<cnt[i]<<'\n';
			ans=(ans+(cur*cnt[i])%MOD)%MOD;
			if(ans<0) ans+=MOD;
			cur=(cur*2)%MOD;
			if(cur<0) cur+=MOD;
		}
		return ans;
	}
	else if(x>=0&&y>=0)
	{
		return 0;
	}
	else if(x<=0&&y<=0)
	{
		return 0;
	}
	else
	{
		x=abs(x); y=abs(y);
		ll g = __gcd(x,y);
		x/=g; y/=g;
		ll k = n/max(x,y);
		ll ans = modpow(2,k+1)+MOD-2;
		while(ans>=MOD) ans-=MOD;
		return ans;
	}
}

ll fact[600001];
ll ifact[600001];
ll inverse[600001];

ll choose(ll n, ll r)
{
	if(r==0) return 1;
	ll ans = fact[n];
	ans=(ans*ifact[r])%MOD;
	ans=(ans*ifact[n-r])%MOD;
	if(ans<0) ans+=MOD;
	return ans;
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	string s, t;
	cin>>s>>t;
	cin>>n;
	fact[0]=1; ifact[0]=1;
	for(int i=1;i<=600000;i++)
	{
		fact[i]=(fact[i-1]*i)%MOD;
		if(fact[i]<0) fact[i]+=MOD;
		ifact[i]=inv(fact[i]);
		inverse[i]=inv(i);
	}
	nt.Sieve(300001);
	for(int i=2;i<=n;i++)
	{
		mob[i]=mobius(i);
	}
	ll sa, sb, sc; //sa = # of As in s, sb = # of Bs in s, sc = # of ?s in s
	ll ta, tb, tc;
	sa=sb=sc=ta=tb=tc=0;
	ll same = 1; //number of ways to fill in ?s such that |S| = |T|
	if(s.length()!=t.length()) same=0;
	else
	{
		for(int i=0;i<s.length();i++)
		{
			if(s[i]=='?'&&t[i]=='?') same=(same*2)%MOD;
			else if(s[i]=='?'||t[i]=='?')
			{
				
			}
			else if(s[i]==t[i])
			{
				
			}
			else
			{
				same=0;
			}
		}
	}
	for(int i=0;i<s.length();i++)
	{
		if(s[i]=='A') sa++;
		else if(s[i]=='B') sb++;
		else sc++;
	}
	for(int i=0;i<t.length();i++)
	{
		if(t[i]=='A') ta++;
		else if(t[i]=='B') tb++;
		else tc++;
	}
	ll ans = 0;
	ll c = 1;
	int cntt=0;
	for(ll i = sa - ta - tc; i <= sa - ta + sc; i++)
	{
		if(i==0)
		{
			ll cc = (c-same)%MOD;
			if(cc<0) cc+=MOD;
			ans=(ans+(cc*solve(i,sa+sb+sc-ta-tb-tc-i))%MOD)%MOD;
			ll tmp = modpow(2,n+1)+MOD-2;
			while(tmp>=MOD) tmp-=MOD;
			tmp=(tmp*tmp)%MOD;
			ans=(ans+(same*tmp)%MOD)%MOD;
		}
		else
		{
			ans=(ans+(c*solve(i,sa+sb+sc-ta-tb-tc-i))%MOD)%MOD;
		}
		if(ans<0) ans+=MOD;
		c=(c*inverse[cntt+1])%MOD;
		c=(c*(sc+tc-cntt))%MOD;
		if(c<0) c+=MOD;
		cntt++;
	}
	cout<<ans<<'\n';
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

vector<pair<char,int> > vec;
const int N = 100000;

void gen()
{
	for(int i = 1; i <= N*10; i++)
	{
		int x=i;
		string r;
		while(x)
		{
			r+=char('0'+x%10);
			x/=10;
		}
		for(int j = int(r.size())-1;j>=0;j--)
		{
			vec.pb(mp(r[j],i));
		}	
	}
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	gen();
	string s; cin>>s;
	int ptr=0; int p2 = 0;
	while(p2<s.length()&&ptr<vec.size())
	{
		if(s[p2]==vec[ptr].fi) 
		{
			p2++;
			ptr++;
			continue;
		}
		ptr++;
	}
	ptr--; //remember this
	cout<<vec[ptr].se;
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;
const int N = 100011;
const int T = 301;
ii dp[N][T];
ii a[N];

ii maxi(ii a, ii b)
{
	if(a.fi!=b.fi)
	{
		if(a.fi>b.fi) return a;
		else return b;
	}
	else
	{
		if(a.se<b.se) return a;
		else return b;
	}
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	int n; cin>>n;
	int t = 300;
	for(int i = 1; i <= n; i++)
	{
		cin>>a[i].fi>>a[i].se;
	}
	sort(a+1,a+n+1);
	for(int i = 0; i <= n; i++)
	{
		for(int j = 0; j <= t; j++)
		{
			dp[i][j]=mp(-1,-1);
		}
	}
	dp[0][0] = mp(0,0);
	ii maxans = mp(0,0);
	for(int i = 1; i <= n; i++)
	{
		for(int j = 0; j <= 300; j++)
		{
			dp[i][j] = dp[i-1][j];
			if(j-a[i].fi>=0&&dp[i-1][j-a[i].fi].fi!=-1)
			{
				dp[i][j] = maxi(dp[i][j], mp(dp[i-1][j-a[i].fi].fi+1, dp[i-1][j-a[i].fi].se+j+20LL*a[i].se));
			}
			if(i==n) maxans = maxi(maxans, dp[i][j]);
		}
	}
	cout<<maxans.fi<<' '<<maxans.se<<'\n';
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

int flip(int l, int r)
{
	cout<<l<<' '<<r<<'\n';
	fflush(stdout);
	int cnt; 
	cin>>cnt;
	return cnt;
}

int main()
{
	int n; cin>>n;
	int c1 = flip(1,1);
	if(c1==n)
	{
		cout<<-1<<'\n'; fflush(stdout);
		return 0;
	}
	int c2 = flip(1,1);
	bool pre = 0; //does previous need a flip?
	if(c2<c1)
	{
		pre=1;
	}
	if(c2==n)
	{
		cout<<-1<<'\n'; fflush(stdout);
		return 0;
	}
	int cnt = c2;
	for(int i = 2; i <= n - 1; i++)
	{
		int prevcnt = cnt;
		int newcnt;
		if(pre)
		{
			newcnt = flip(i-1,i);
			if(newcnt==prevcnt)
			{
				pre=1;
			}
			else
			{
				pre=0;
			}
		}
		else
		{
			newcnt = flip(i,i);
			if(newcnt<prevcnt)
			{
				pre=1;
			}
			else
			{
				pre=0;
			}
		}
		cnt=newcnt;
		if(newcnt==n)
		{
			cout<<-1<<'\n'; fflush(stdout);
			return 0;
		}
	}
	if(cnt==n)
	{
		cout<<-1<<'\n'; fflush(stdout);
		return 0;
	}
	if(cnt==n-2)
	{
		flip(n-1,n);
		cout<<-1<<'\n'; fflush(stdout);
		return 0;
	}
	assert(cnt==n-1);
	if(pre)
	{
		flip(n-1,n-1);
		cout<<-1<<'\n'; fflush(stdout);
	}
	else 
	{
		flip(n,n);
		cout<<-1<<'\n'; fflush(stdout);
	}
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

bool good(int n) //returns whether n can be represented as ax + by
{
	cout<<"? "<<n<<'\n';
	fflush(stdout);
	int x; cin>>x;
	return x;
}

const int N = 500000;
int main()
{
	int a, b;
	a=20;
	int t;
	cin>>t;
	while(t--)
	{
		a=20;
		for(int i = 2; i <= 19; i++)
		{
			if(good(i))
			{
				a=i;
				break;
			}
		}
		int lo = 1; int hi = N/a+2;
		while(lo<=hi)
		{
			int mid = (lo+hi)>>1;
			int pos = -1;
			for(int i = 1; i <= a - 1; i++)
			{
				if(good(mid*a+i))
				{
					pos=mid*a+i;
					break;
				}
			}
			if(pos==-1)
			{
				lo=mid+1;
			}
			else
			{
				b=pos;
				hi=mid-1;
			}
		}
		cout<<"! "<<a<<' '<<b<<'\n';
		fflush(stdout);
		//int x; cin>>x;
	}
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

bool good(int n) //returns whether n can be represented as ax + by
{
	cout<<"? "<<n<<'\n';
	fflush(stdout);
	int x; cin>>x;
	return x;
}

struct NumberTheory
{
	vector<ll> primes;
	vector<bool> prime;
	vector<ll> totient;
	vector<ll> sumdiv;
	vector<ll> bigdiv;
	void Sieve(ll n)
	{
		prime.assign(n+1, 1);
		prime[1] = false;
		for(ll i = 2; i <= n; i++)
		{
			if(prime[i])
			{
				primes.pb(i);
				for(ll j = i*2; j <= n; j += i)
				{
					prime[j] = false;
				}
			}
		}
	}
	
	ll phi(ll x)
	{
		map<ll,ll> pf;
		ll num = 1; ll num2 = x;
		for(ll i = 0; primes[i]*primes[i] <= x; i++)
		{
			if(x%primes[i]==0)
			{
				num2/=primes[i];
				num*=(primes[i]-1);
			}
			while(x%primes[i]==0)
			{
				x/=primes[i];
				pf[primes[i]]++;
			}
		}
		if(x>1)
		{
			pf[x]++; num2/=x; num*=(x-1);
		}
		x = 1;
		num*=num2;
		return num;
	}
	
	bool isprime(ll x)
	{
		if(x==1) return false;
		for(ll i = 0; primes[i]*primes[i] <= x; i++)
		{
			if(x%primes[i]==0) return false;
		}
		return true;
	}

	void SievePhi(ll n)
	{
		totient.resize(n+1);
		for (int i = 1; i <= n; ++i) totient[i] = i;
		for (int i = 2; i <= n; ++i)
		{
			if (totient[i] == i)
			{
				for (int j = i; j <= n; j += i)
				{
					totient[j] -= totient[j] / i;
				}
			}
		}
	}
	
	void SieveSumDiv(ll n)
	{
		sumdiv.resize(n+1);
		for(int i = 1; i <= n; ++i)
		{
			for(int j = i; j <= n; j += i)
			{
				sumdiv[j] += i;
			}
		}
	}
	
	ll getPhi(ll n)
	{
		return totient[n];
	}
	
	ll getSumDiv(ll n)
	{
		return sumdiv[n];
	}
	
	ll modpow(ll a, ll b, ll mod)
	{
		ll r = 1;
		if(b < 0) b += mod*100000LL;
		while(b)
		{
			if(b&1) r = (r*a)%mod;
			a = (a*a)%mod;
			b>>=1;
		}
		return r;
	}
	
	ll inv(ll a, ll mod)
	{
		return modpow(a, mod - 2, mod);
	}
	
	ll invgeneral(ll a, ll mod)
	{
		ll ph = phi(mod);
		ph--;
		return modpow(a, ph, mod);
	}
	
	void getpf(vector<ii>& pf, ll n)
	{
		for(ll i = 0; primes[i]*primes[i] <= n; i++)
		{
			int cnt = 0;
			while(n%primes[i]==0)
			{
				n/=primes[i]; cnt++;
			}
			pf.pb(ii(primes[i], cnt));
		}
		if(n>1)
		{
			pf.pb(ii(n, 1));
		}
	}

	//ll op;
	void getDiv(vector<ll>& div, vector<ii>& pf, ll n, int i)
	{
		//op++;
		ll x, k;
		if(i >= pf.size()) return ;
		x = n;
		for(k = 0; k <= pf[i].se; k++)
		{
			if(i==int(pf.size())-1) div.pb(x);
			getDiv(div, pf, x, i + 1);
			x *= pf[i].fi;
		}
	}
};

NumberTheory nt;

const int N = 500000;
int main()
{
	nt.Sieve(N+1);
	int a, b;
	a=20;
	int t; cin>>t;
	while(t--)
	{
		a=20;
		for(int i = 2; i <= 19; i++)
		{
			if(good(i))
			{
				a=i;
				break;
			}
		}
		int MX = N*a-1;
		int lo = 1; int hi = N-1;
		int ans = 0;
		while(lo<=hi)
		{
			int mid=(lo+hi)>>1;
			if(good(MX-mid*a))
			{
				ans=mid;
				lo=mid+1;
			}
			else
			{
				hi=mid-1;
			}
		}
		int big = MX - ans*a; //this thing is a multiple of b
		vector<ii> ppf;
		nt.getpf(ppf, big);
		for(int i = 0; i < ppf.size(); i++)
		{
			int p = ppf[i].fi;
			while(big%p==0)
			{
				if(!good(big/p)) break;
				big/=p;
			}
		}
		b=big;
		cout<<"! "<<a<<' '<<b<<'\n';
		fflush(stdout);
		//int x; cin>>x;
	}
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int N = 1000011;

bool iscycle[N];
int visited[N];
ll ans[N];
int n;
vi adj[N];
int a[N];
vi cyc;
ll subsize[N];
vector<vector<ll> > vec; //nodes sorted by height
ll subsizesum;

void dfs2(int u)
{
	cyc.pb(u);
	iscycle[u]=1;
	visited[u] = 3;
	if(visited[a[u]] == 3) return ;
	dfs2(a[u]);
}

void findcyc(int u)
{
	visited[u] = 2;
	if(visited[a[u]] == 0)
	{
		findcyc(a[u]);
	}
	else if(visited[a[u]] == 1)
	{
		visited[u] = 1;
		return ;
	}
	else
	{
		dfs2(u);
	}
	visited[u] = 1;
}

void upd(int l, int r, ll s)
{
	ans[r+1]-=s;
	ans[l]+=s;
}

int ptr;
void dfs(int u, int d)
{
	if(vec[ptr].size()<=d) vec[ptr].pb(0);
	vec[ptr][d]++;
	subsize[u]=1;
	for(int i = 0; i < adj[u].size(); i++)
	{
		int v = adj[u][i];
		dfs(v,d+1);
		subsize[u]+=subsize[v];
	}
	if(d>0) subsizesum+=subsize[u];
}

int main()
{
	//ios_base::sync_with_stdio(0); cin.tie(0);
	scanf("%d", &n);
	for(int i = 0; i < n; i++)
	{
		//cin>>a[i];
		scanf("%d", a + i);
		a[i]--;
	}
	findcyc(0);
	memset(visited,0,sizeof(visited));
	for(int i = 0; i < n; i++)
	{
		if(!(iscycle[a[i]]&&iscycle[i]))
		{
			adj[a[i]].pb(i);
		}
	}
	vec.resize(int(cyc.size()));
	int c = cyc.size();
	for(int i = 0; i < cyc.size(); i++)
	{
		dfs(cyc[i],0);
		ptr++;
	}
	upd(c,c,ll(n-1)*ll(n-c) - subsizesum);
	for(int i = 0; i < c; i++)
	{
		for(int j = 0; j < vec[i].size(); j++)
		{
			//j is the height
			upd(j+1,j+c,vec[i][j]);
		}
	}
	upd(c,c,-c); //subtract original cycle edges
	ll cur=0;
	for(int i = 1; i <= n; i++)
	{
		cur+=ans[i];
		printf("%I64d ", cur);
	}
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<int,int> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int N = 1e5 + 3;
ll a[N];
int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	int n, c; cin >> n >> c;
	for(int i = 0; i < n; i++)
	{
		cin >> a[i];
	}
	//sort(a, a + n);
	int cnt = 0;
	for(int i = 0; i < n; i++)
	{
		if(i == 0) cnt++;
		else
		{
			if(a[i] - a[i - 1] <= c) cnt++;
			else cnt = 1;
		}
	}
	cout << cnt;
	return 0;
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<int,int> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int N = 10000;
int cnt[27];
string s; int n;

bool valid()
{
	for(int i = 0; i < 26; i++)
	{
		if(cnt[i] >= 2) return false;
	}
	return true;
}

void fillall()
{
	for(int i = 0; i < n; i++)
	{
		if(s[i] == '?') s[i] = 'A';
	}
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	cin >> s;
	n = s.length();
	if(n < 26) {cout << -1; return 0;}
	for(int i = 25; i < n; i++)
	{
		memset(cnt, 0, sizeof(cnt));
		for(int j = i; j >= i - 25; j--)
		{
			cnt[s[j]-'A']++;
		}
		if(valid())
		{
			//cout << "GG " << i << '\n';
			int cur = 0;
			while(cnt[cur]>0) cur++;
			for(int j = i - 25; j <= i; j++)
			{
				if(s[j] == '?')
				{
					s[j] = cur + 'A';
					cur++;
					while(cnt[cur]>0) cur++;
				}
			}
			fillall();
			cout << s;
			return 0;
		}
	}
	cout << -1;
	return 0;
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<int,int> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int N = 10000;
int cnt[27];
string s; int n;

bool valid()
{
	for(int i = 0; i < 26; i++)
	{
		if(cnt[i] >= 2) return false;
	}
	return true;
}

void fillall()
{
	for(int i = 0; i < n; i++)
	{
		if(s[i] == '?') s[i] = 'A';
	}
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	cin >> s;
	n = s.length();
	if(n < 26) {cout << -1; return 0;}
	for(int i = 0; i < 26; i++) cnt[s[i]-'A']++;
	if(valid())
	{
		int cur = 0;
		while(cnt[cur]>0) cur++;
		for(int i = 0; i < 26; i++)
		{
			if(s[i] == '?')
			{
				s[i] = cur + 'A';
				cur++;
				while(cnt[cur]>0) cur++;
			}
		}
		fillall();
		cout << s;
		return 0;
	}
	for(int i = 26; i < n; i++)
	{
		cnt[s[i]-'A']++; cnt[s[i-26]-'A']--;
		if(valid())
		{
			//cout << "GG " << i << '\n';
			int cur = 0;
			while(cnt[cur]>0) cur++;
			for(int j = i - 25; j <= i; j++)
			{
				if(s[j] == '?')
				{
					s[j] = cur + 'A';
					cur++;
					while(cnt[cur]>0) cur++;
				}
			}
			fillall();
			cout << s;
			return 0;
		}
	}
	cout << -1;
	return 0;
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<int,int> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int N = 50000;
int cnt[27];
string s; int n;
int counter;

bool valid()
{
    //cout << counter << endl;
	return (counter == 26);
}

void fillall()
{
	for(int i = 0; i < n; i++)
	{
		if(s[i] == '?') s[i] = 'A';
	}
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	cin >> s;
	n = s.length();
	if(n < 26) {cout << -1; return 0;}
	counter = 0;
	for(int i = 0; i < 26; i++)
	{
		if(s[i] == '?')
		{
			counter++; continue;
		}
		cnt[s[i]-'A']++;
		if(cnt[s[i]-'A'] == 1) counter++;
	}
	if(valid())
	{
		int cur = 0;
		while(cnt[cur]>0) cur++;
		for(int i = 0; i < 26; i++)
		{
			if(s[i] == '?')
			{
				s[i] = cur + 'A';
				cur++;
				while(cnt[cur]>0) cur++;
			}
		}
		fillall();
		cout << s;
		return 0;
	}
	for(int i = 26; i < n; i++)
	{
		if(s[i] != '?') {cnt[s[i]-'A']++; if(cnt[s[i]-'A']==1) counter++;}
		if(s[i-26] != '?') {cnt[s[i-26]-'A']--; if(cnt[s[i-26]-'A']==0) counter--;}
		if(s[i-26] == '?') counter--;
		if(s[i] == '?') counter++;
		if(valid())
		{
			int cur = 0;
			while(cnt[cur]>0) cur++;
			for(int j = i - 25; j <= i; j++)
			{
				if(s[j] == '?')
				{
					s[j] = cur + 'A';
					cur++;
					while(cnt[cur]>0) cur++;
				}
			}
			fillall();
			cout << s;
			return 0;
		}
	}
	cout << -1;
	return 0;
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	ll n; cin >> n;
	for(ll i = 1; i <= n; i++)
	{
		if(i == 1) cout << 2 << '\n';
		else cout << i*(i+1)*(i+1)-(i-1) << '\n';
	}
	return 0;
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int N = 1001;
const int M = 10001;
const ll INF = ll(1e18);

vector<ii> adj[N];
vector<int> adj2[N];
int L[M]; int R[M];
ll d1[N];
ll d2[N];
int par[N];
ll dist[N][N];

int n, m, l, s, e;

void dijkstra()
{
	d1[s] = 0;
	priority_queue<ii, vector<ii>, greater<ii> > pq;
	pq.push(ii(0, s));
	int cnt = 0;
	while(!pq.empty())
	{
		cnt++;
		int u = pq.top().se; ll d = pq.top().fi; pq.pop();
		for(int i = 0; i < adj[u].size(); i++)
		{
			int v = adj[u][i].fi; ll w = adj[u][i].se;
			if(d + w < d1[v])
			{
				d1[v] = d + w;
				pq.push(ii(d1[v], v));
			}
		}
	}
	cerr << "DIJKSTRA OPERATIONS : " << cnt << '\n';
}

bool dijkstra2()
{
	d2[s] = 0;
	priority_queue<ii, vector<ii>, greater<ii> > pq;
	pq.push(ii(0, s));
	while(!pq.empty())
	{
		int u = pq.top().se; ll d = pq.top().fi; pq.pop();
		for(int i = 0; i < adj2[u].size(); i++)
		{
			int v = adj2[u][i]; ll w = dist[u][v];
			if(d + abs(w) < d2[v])
			{
				d2[v] = d + abs(w);
				par[v] = u;
				pq.push(ii(d2[v], v));
			}
		}
	}
	if(d2[e] > l) return false;
	int cur = e;
	while(cur != s)
	{
		int p = par[cur];
		if(dist[p][cur] < 0)
		{
			dist[p][cur] = -2; dist[cur][p] = -2;
		}
		cur = par[cur];
	}
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < n; j++)
		{
			if(dist[i][j] == -1)
			{
				dist[i][j] = INF;
			}
		}
	}
	cur = e;
	while(cur != s)
	{
		int p = par[cur];
		if(dist[p][cur] < 0)
		{
			dist[p][cur] = -1; dist[cur][p] = -1;
		}
		cur = par[cur];
	}
	return true;
}

void print()
{
	cout << "YES\n";
	for(int i = 0; i < m; i++)
	{
		int u = L[i]; int v = R[i];
		ll d = dist[u][v];
		if(d < 0) d = -d;
		cout << u << ' ' << v << ' ' << d << '\n';
	}
}

bool relax()
{
	for(int i = 0; i < n; i++) d2[i] = INF;
	memset(par, -1, sizeof(par)); //shouldn't be neccesary
	d2[s] = 0;
	priority_queue<ii, vector<ii>, greater<ii> > pq;
	pq.push(ii(0, s));
	while(!pq.empty())
	{
		int u = pq.top().se; ll d = pq.top().fi; pq.pop();
		for(int i = 0; i < adj2[u].size(); i++)
		{
			int v = adj2[u][i]; ll w = dist[u][v];
			if(d + abs(w) < d2[v])
			{
				d2[v] = d + abs(w);
				par[v] = u;
				pq.push(ii(d2[v], v));
			}
		}
	}
	//cerr << d2[e] << '\n';
	if(d2[e] == l) return true;
	int cur = e; bool meet = false;
	while(cur != s)
	{
		int p = par[cur];
		if(!meet && dist[p][cur] < 0)
		{
			ll d = abs(dist[p][cur]);
			dist[p][cur] = d + l - d2[e];
			dist[cur][p] = d + l - d2[e];
			//cerr << d << endl;
			meet = true;
		}
		if(meet) break;
		cur = par[cur];
	}
	return false;
}

int dist1[N][N];

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	cin >> n >> m >> l >> s >> e;
	memset(dist1, -1, sizeof(dist1));
	for(int i = 0; i < m; i++)
	{
		int u, v, w; cin >> u >> v >> w;
		if(w > 0) 
		{
			adj[u].pb(ii(v, w)); adj[v].pb(ii(u, w));
			adj2[u].pb(v); adj2[v].pb(u);
			dist1[u][v] = w; dist1[v][u] = w;
			dist[u][v] = w; dist[v][u] = w;
		}
		else
		{
			adj2[u].pb(v); adj2[v].pb(u);
			dist[u][v] = -1; dist[v][u] = -1;
		}
		L[i] = u; R[i] = v;
	}	
	for(int i = 0; i < n; i++)
	{
		d1[i] = INF; d2[i] = INF;
	}
	dijkstra(); //cerr << d1[e] << '\n';
	if(d1[e] < l)
	{
		cout << "NO\n";
		return 0;
	}
	if(d1[e] == l)
	{
		cout << "YES\n";
		for(int i = 0; i < m; i++)
		{
			int u = L[i]; int v = R[i];
			ll d = dist1[u][v];
			if(d <= 0) d = INF;
			cout << u << ' ' << v << ' ' << d << '\n';
		}
		return 0;
	}
	bool tmp = dijkstra2();
	if(!tmp)
	{
		cout << "NO\n";
		return 0;
	}
	//cerr << d2[e] << '\n';
	int cnt = 0;
	while(!relax()) 
	{
		relax(); cnt++;
	}
	cerr << "RELAXATIONS DONE : " << cnt << '\n';
	print();
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<ll,ll> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int N = 1001;
const int M = 10001;
const ll INF = ll(1e18);

int n, m, l, s, e;

struct edge
{
	int to; ll w; int label;
	edge(int _to, int _w, int _label){to = _to, w = _w, label = _label;}
};

int edgecnt = -1;
vector<edge> adj[N];
ll dist[N];
set<ii> used;

ll dijk(int p, ll val)
{
	for(int i = 0; i < n; i++) dist[i] = INF;
	dist[s] = 0;
	priority_queue<ii, vector<ii>, greater<ii> > pq;
	pq.push(ii(0, s));
	while(!pq.empty())
	{
		int u = pq.top().se; ll d = pq.top().fi; pq.pop();
		for(int i = 0; i < adj[u].size(); i++)
		{
			edge tmp = adj[u][i];
			int v = tmp.to; ll w = tmp.w; int lab = tmp.label;
			if(lab >= 0)
			{
				if(lab < p) w = 1;
				else if(lab == p) w = val;
				else w = ll(1e14);
			}
			if(d + w < dist[v])
			{
				dist[v] = d + w;
				pq.push(ii(dist[v], v));
			}
		}
	}
	return dist[e];
}

void setw(int p, ll val)
{
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < adj[i].size(); j++)
		{
			int lab = adj[i][j].label;
			if(lab >= 0)
			{
				if(lab < p)
				{
					adj[i][j].w = 1;
				}
				else if(lab == p)
				{
					adj[i][j].w = val;
				}
				else
				{
					adj[i][j].w = INF;
				}
			}
		}
	}
}

void print()
{
	cout << "YES\n";
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < adj[i].size(); j++)
		{
			edge tmp = adj[i][j];
			int v = tmp.to; ll w = tmp.w; 
			if(used.find(ii(i, v)) == used.end())
			{
				cout << i << ' ' << v << ' ' << w << '\n';
				used.insert(ii(i, v)); used.insert(ii(v, i));
			}
		}
	}
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	cin >> n >> m >> l >> s >> e;
	for(int i = 0; i < m; i++)
	{
		int u, v, w;
		cin >> u >> v >> w;
		int lab = -1;
		if(w == 0) 
		{
			lab = ++edgecnt;
		}
		adj[u].pb(edge(v, w, lab));
		adj[v].pb(edge(u, w, lab));
	}
	ll x = dijk(edgecnt, 1); ll y = dijk(-1, 1);
	if(!(x <= l && l <= y))
	{
		cout << "NO\n";
		return 0;
	}
	ll lo = -1; ll hi = edgecnt;
	ll mid, ans;
	while(lo <= hi)
	{
		mid = (lo+hi)/2;
		if(dijk(mid, 1) <= l)
		{
			ans = mid;
			hi = mid - 1;
		}
		else
		{
			lo = mid + 1;
		}
	}
	//now [0..ans] as 1 will give <= L whereas [0..ans - 1] as 1 will give > L
	if(ans == -1)
	{
		setw(-1, 0);
		print();
		return 0;
	}
	lo = 1; hi = INF;
	int ans2 = 0;
	while(lo <= hi)
	{
		mid = (lo+hi)>>1;
		if(dijk(ans, mid) <= l)
		{
			ans2 = mid;
			lo = mid + 1;
		}
		else
		{
			hi = mid - 1;
		}
	}
	//cerr << ans << ' ' << ans2 << '\n';
	setw(ans, ans2);
	print();
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<int,int> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int N = 1e5 + 1;
const int MAX = 1e9;
int MOD, n;

bool isprime[100001];
vector<ll> primes;
vector<ii> adj[N];
int subsize[N];
bool visited[N];
int treesize;
vi clrlist;
ll up[N];
ll down[N];
int h[N];
int PHI;
int dppart[N];

ll mult(ll a, ll b)
{
	return (a*b)%MOD;
}

ll add(ll a, ll b)
{
	return (a+b+MOD)%MOD;
}

ll modpow(ll a, ll b)
{
	ll r = 1;
	while(b)
	{
		if(b&1) r=(r*a)%MOD;
		a=(a*a)%MOD;
		b>>=1;
	}
	return r;
}

void Sieve(int n)
{
	memset(isprime, 1, sizeof(isprime));
	isprime[1] = false;
	for(int i = 2; i <= n; i++)
	{
		if(isprime[i])
		{
			primes.pb(i);
			for(int j = 2*i; j <= n; j += i)
			{
				isprime[j] = false;
			}
		}
	}
}

int phi(int n)
{
	ll num = 1; ll num2 = n;
	for(ll i = 0; primes[i]*primes[i] <= n; i++)
	{
		if(n%primes[i]==0)
		{
			num2/=primes[i];
			num*=(primes[i]-1);
		}
		while(n%primes[i]==0)
		{
			n/=primes[i];
		}
	}
	if(n>1)
	{
		num2/=n; num*=(n-1);
	}
	n = 1;
	num*=num2;
	return num;
}

ll inv(ll a)
{
	return modpow(a, PHI-1);
}

void dfs(int u, int par)
{
	if(par == -1) clrlist.clear();
	subsize[u] = 1; clrlist.pb(u);
	for(int i = 0; i < adj[u].size(); i++)
	{
		int v = adj[u][i].fi;
		if(visited[v]) continue;
		if(v == par) continue;
		dfs(v, u);
		subsize[u] += subsize[v];
	}	
	if(par == -1) treesize = subsize[u];
}

int centroid(int u, int par)
{
	for(int i = 0; i < adj[u].size(); i++)
	{
		int v = adj[u][i].fi;
		if(visited[v]) continue;
		if(v == par) continue;
		if(subsize[v]*2 > treesize) return centroid(v, u);
	}
	return u;
}

int parts = 0;
void fill(int u, int p, int cent)
{
	if(p == cent)
	{
		dppart[u] = parts;
		parts++;
	}
	else if(p != -1)
	{
		dppart[u] = dppart[p];
	}
	for(int i = 0; i < adj[u].size(); i++)
	{
		int v = adj[u][i].fi; int w = adj[u][i].se;
		if(v == p || visited[v]) continue;
		down[v] = add(mult(down[u], 10), w);
		up[v] = add(up[u], mult(modpow(10, h[u]), w));
		h[v] = h[u] + 1;
		fill(v, u, cent);
		//cout << v << ' ' << u << ' ' << up[v] << ' ' << up[u] << '\n';
	}
}

ll solve(int cent)
{
	for(int i = 0; i < clrlist.size(); i++)
	{
		up[clrlist[i]] = 0; down[clrlist[i]] = 0; h[clrlist[i]] = 0;
	}
	parts = 0;
	fill(cent, -1, cent); parts--;
	dppart[cent] = -1; 
	map<ll,ll> tot; //only count up
	vector<map<ll,ll> > vec; //only count up, but in specific subtree
	vec.resize(parts+1);
	tot[0]++;
	for(int i = 0; i < clrlist.size(); i++)
	{
		int u = clrlist[i];
		//cout << u << ' ' << up[u] << ' ' << down[u] << '\n';
		if(u == cent) continue;
		tot[up[u]]++;
		vec[dppart[u]][up[u]]++;
	}
	ll ans = 0;
	for(int i = 0; i < clrlist.size(); i++)
	{
		int u = clrlist[i];
		int ht = h[u];
		int pt = dppart[u];
		if(u == cent)
		{
			ans += (tot[0] - 1); //exclude cent as the vertex
		}
		else
		{
			ll val = ((-down[u])%MOD+MOD)%MOD;
			val = mult(val, inv(modpow(10, ht)));
			ans += (tot[val] - vec[pt][val]);
		}
	}
	return ans;
}

ll compsolve(int u)
{
	dfs(u, -1);
	int cent = centroid(u, -1);
	ll ans = solve(cent);
	//cout << u << ' ' << cent << ' ' << ans << '\n';
	visited[cent] = true;
	for(int i = 0; i < adj[cent].size(); i++)
	{
		int v = adj[cent][i].fi;
		if(!visited[v]) ans += compsolve(v);
	}
	return ans;
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	cin >> n >> MOD;
	if(MOD == 1)
	{
		cout << ll(n)*ll(n - 1);
		return 0;
	}
	Sieve(100000); PHI = phi(MOD);
	for(int i = 0; i < n - 1; i++) //tree is 0-indexed
	{
		int u, v, w; cin >> u >> v >> w;
		adj[u].pb(ii(v, w)); adj[v].pb(ii(u, w));
	}
	cout << compsolve(0) << '\n';
	return 0;
}

Our AC solution uses base 6 instead of binary. Write t in base 6. Note that t has at most 24 digits in base 6, so to add a new digit we can increase the dimensions by 2 and the number of deleted edges can be up to 12 per digit. We'll construct such a way. This method is explained in the picture below. The key is to first construct a grid which has 0 in it, then find a way to get 6x + i for all 0 ≤ i ≤ 5 from x by maintaining a wall of 1s around the squares. This method uses a 50 × 50(2·24 + 2 = 50) grid and at most 24·12 + 2 = 290 edge deletions and will get AC.

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<int,int> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int INF = 1e9 + 7;
const int MOD = 1e9 + 7;

typedef pair<ii,ii> move;

set<move> ans;
int curx, cury;

bool isvalid(move x)
{
	if(x.fi.fi > 0 && x.se.fi > 0 && x.fi.se > 0 && x.se.se > 0 && x.fi.fi <= curx && x.fi.se <= cury && x.se.fi <= curx && x.se.se <= cury) return true;
	return false;
}

void edge(int x1, int y1, int x2, int y2)
{
	ans.insert(mp(mp(x1, y1), mp(x2, y2)));
}

void add(int bit)
{
	int x = curx;
	edge(x,x+2,x,x+3);
	edge(x+1,x+2,x+1,x+3);
	edge(x+2,x,x+3,x);
	edge(x+2,x+1,x+3,x+1);
	edge(x-2,x+3,x-1,x+3);
	edge(x,x+4,x+1,x+4);
	edge(x+3,x-2,x+3,x-1);
	edge(x+4,x,x+4,x+1);
	edge(x-1,x+1,x,x+1);
	if(bit%3==0) edge(x-1,x+2,x,x+2);
	if(bit%3!=2) edge(x+2,x-1,x+2,x);
	if(bit<3) edge(x+1,x-1,x+1,x);
	curx += 2; cury += 2;
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	ll t; cin >> t;
	vector<int> digits;
	while(t)
	{
		digits.pb(t%6);
		t/=6;
	}
	reverse(digits.begin(), digits.end());
	edge(1, 2, 2, 2);
	edge(2, 1, 2, 2);
	curx = 2; cury = 2;
	for(int i = 0; i < digits.size(); i++)
	{
		add(digits[i]);
	}
	cout << curx << ' ' << cury << '\n';
	vector<move> clr;
	for(set<move>::iterator it = ans.begin(); it != ans.end(); it++)
	{
		if(!isvalid(*it))
		{
			clr.pb(*it);
		}
	}
	for(int i = 0; i < clr.size(); i++)
	{
		ans.erase(clr[i]);
	}
	cout << ans.size() << '\n';
	for(set<move>::iterator it = ans.begin(); it != ans.end(); it++)
	{
		move tmp = (*it);
		cout << tmp.fi.fi << ' ' << tmp.fi.se << ' ' << tmp.se.fi << ' ' << tmp.se.se << '\n';
	}
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<int,int> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int N = 251;
const int MOD = 998244353;
ll inv2;
ll prt;
ll iprt;

ll dpncr[N][N];
ll fact[N];
ll inverse[N];
ll g[N][N];
ll sumg[N][N][N];

vector<ii> perm;
int A, B, C, D;

ll modpow(ll a, ll b)
{
	ll r = 1;
	while(b)
	{
		if(b&1) r = (r*a)%MOD;
		a = (a*a)%MOD;
		b>>=1;
	}
	return r;
}

ll inv(ll a)
{
	return modpow(a, MOD - 2);
}

ll choose(int n, int m)
{
	if(m < 0) return 0;
	if(n < m) return 0;
	if(m == 0) return 1;
	if(n == m) return 1;
	if(dpncr[n][m] != -1) return dpncr[n][m];
	dpncr[n][m] = choose(n - 1, m - 1) + choose(n - 1, m);
	dpncr[n][m] += MOD; dpncr[n][m] %= MOD;
	return dpncr[n][m];
}

void computefact()
{
	fact[0] = 1;
	for(ll i = 1; i < N; i++)
	{
		fact[i] = (fact[i - 1]*i)%MOD;
	}
	for(ll i = 1; i < N; i++)
	{
		inverse[i] = modpow(i, MOD - 2);
	}
}

void print(vector<ii>& vec)
{
	for(int i = 0; i < vec.size(); i++)
	{
		cout << vec[i].fi << ' ' << vec[i].se << endl;
	}
	cout << "------------------------------------------------" << endl;
}

void printans(vector<ll>& vec)
{
	for(int i = 0; i < vec.size(); i++)
	{
		cout << vec[i] << ' ';
	}
	cout << endl;
}

void printansi(vector<int>& vec)
{
	for(int i = 0; i < vec.size(); i++)
	{
		cout << vec[i] << ' ';
	}
	cout << endl;
}

void calcpos(vector<ii>& pos)
{
	pos.resize(perm.size());
	for(int i = 0; i < perm.size(); i++)
	{
		pos[i] = ii(-1, -1);
	}
	for(int i = 1; i < perm.size(); i++)
	{
		if(perm[i].fi > 0)
		{
			pos[perm[i].fi].fi = i;
		}
		if(perm[i].se > 0)
		{
			pos[perm[i].se].se = i;
		}
	}
}

int reduce()
{
	int n = perm.size() - 1;
	vector<ii> pos;
	int cnt = 0;
	for(int i = 1; i <= n; i++) //Do a reduction
	{
		calcpos(pos);
		//print(pos);
		if(pos[i].fi > 0 && pos[i].se > 0)
		{
			if(pos[i].fi == pos[i].se) 
			{
				cnt++;
				ii tmp1 = perm[pos[i].fi];
				for(vector<ii>::iterator it = perm.begin(); it != perm.end(); it++)
				{
					if((*it) == tmp1)
					{
						perm.erase(it); break;
					}
				}
				continue;
			}
			int p1 = pos[i].se; int l = perm[p1].fi; ii tmp1 = perm[p1];
			int p2 = pos[i].fi; int r = perm[p2].se; ii tmp2 = perm[p2];
			for(vector<ii>::iterator it = perm.begin(); it != perm.end(); it++)
			{
				if((*it) == tmp1)
				{
					perm.erase(it); break;
				}
			}
			for(vector<ii>::iterator it = perm.begin(); it != perm.end(); it++)
			{
				if((*it) == tmp2)
				{
					perm.erase(it); break;
				}
			}
			perm.pb(ii(l, r));
		}
	}
	//count A, B, C, D
	for(int i = 1; i < perm.size(); i++)
	{
		if(perm[i].fi > 0 && perm[i].se > 0)
		{
			assert(perm[i].fi != perm[i].se);
			C++;
		}
		else if(perm[i].fi > 0)
		{
			A++;
		}
		else if(perm[i].se > 0)
		{
			B++;
		}
		else
		{
			D++;
		}
	}
	return cnt;
}

ll mult(ll a, ll b)
{
	ll r = (a*b)%MOD;
	r = (r+MOD)%MOD;
	return r;
}

ll add(ll a, ll b)
{
	ll r = ((a+b)%MOD+MOD)%MOD;
	return r;
}

ll F(ll a, ll b, ll c, ll d)
{
	ll ans = 1;
	if(d == 0)
	{
		if(a == 0 && b == 0 && c == 0) return 1;
		else return 0;
	}
	ans = mult(ans, fact[a]);
	ans = mult(ans, fact[b]);
	ans = mult(ans, fact[c]);
	ans = mult(ans, choose(a+d-1, d-1));
	ans = mult(ans, choose(b+d-1, d-1));
	ans = mult(ans, choose(d, C));
	return ans;
}

const int LG = 9;
const int root_pw = (1<<LG);
void fft (vector<int> & a, bool invert) 
{
	int n = (int) a.size();
 
	for (int i=1, j=0; i<n; ++i) {
		int bit = n >> 1;
		for (; j>=bit; bit>>=1)
			j -= bit;
		j += bit;
		if (i < j)
			swap (a[i], a[j]);
	}
 
	for (int len=2; len<=n; len<<=1) {
		int wlen = invert ? iprt : prt;
		for (int i=len; i<root_pw; i<<=1)
			wlen = int((wlen*1LL*wlen)%MOD);
		for (int i=0; i<n; i+=len) {
			int w = 1;
			for (int j=0; j<len/2; ++j) {
				int u = a[i+j]; int v = int((a[i+j+len/2]*1LL*w)%MOD);
				a[i+j] = u+v < MOD ? u+v : u+v-MOD;
				a[i+j+len/2] = u-v >= 0 ? u-v : u-v+MOD;
				w = int (w * 1LL * wlen % MOD);
			}
		}
	}
	if (invert) {
		ll nrev = inv(n);
		for (int i=0; i<n; ++i)
			a[i] = int((a[i]*1LL*nrev)%MOD);
	}
}

void multiply(vector<int>& a, vector<int>& b, vector<int>& res)
{
	vector<int> fa(a.begin(), a.end()), fb(b.begin(), b.end());
	int n = 1;
	while(n < max(a.size(), b.size())) n <<= 1;
	fa.resize(n); fb.resize(n);
	//cerr << "A : "; printansi(fa); cerr << "B : "; printansi(fb);
	fft(fa, 0); fft(fb, 0);
	//cerr << "INVERT ONCE : ";
	//printans(fa); printans(fb);
	//fft(fa, 1); fft(fb, 1); cerr << "INVERT BACK A: "; printans(fa); cerr << "INVERT BACK B: "; printans(fb);
	res.resize(n);
	for(int i = 0; i < n; i++) res[i] = int((fa[i]*1LL*fb[i])%MOD);
	//printans(fa);
	fft(res, 1);
	//cerr << "CONVOLUTION : "; printansi(res);
	
}

void computeg(int n)
{
	g[0][0] = 1;
	g[1][1] = 1;
	for(int i = 2; i <= n; i++)
	{
		for(int j = 1; j <= i; j++)
		{
			for(int k = 1; k <= i; k++)
			{
				g[i][j] = add(g[i][j], mult(g[i-k][j-1], mult(choose(i, k), fact[k-1])));
			}
			//cerr << g[i][j] << '\n';
			g[i][j] = mult(g[i][j], inverse[j]);
			//cerr << "G : " << i << ' ' << j << ' ' << g[i][j] << '\n';
		}
	}
	for(int i = 0; i <= A; i++)
	{
		for(int j = 0; j <= B; j++)
		{
			vector<int> gi; vector<int> gj;
			gi.resize(n+1); gj.resize(n+1);
			for(int k = 0; k <= n; k++)
			{
				gi[k] = int(g[i][k]);
				gj[k] = int(g[j][k]);
				//cerr << gi[k] << ' ' << gj[k] << '\n';
			}
			vector<int> res;
			multiply(gi, gj, res);
			for(int k = 0; k <= n; k++)
			{
				sumg[i][j][k] = res[k];
			}
		}
	}
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	int n; cin >> n; perm.resize(n+1);
	for(int i = 1; i <= n; i++)
	{
		cin >> perm[i].fi;
	}
	for(int i = 1; i <= n; i++)
	{
		cin >> perm[i].se;
	}
	ll tmpmult = 7*17;
	tmpmult = mult(tmpmult, modpow(2, 23 - LG)); 
	prt = modpow(3, tmpmult);
	inv2 = inv(2);
	iprt = inv(prt);
	A = 0; B = 0; C = 0; D = 0;
	memset(dpncr, -1, sizeof(dpncr));
	memset(sumg, 0, sizeof(sumg));
	memset(g, 0, sizeof(g));
	memset(fact, 0, sizeof(fact));
	memset(inverse, 0, sizeof(inverse));
	int cycles = reduce();
	computefact(); 
	computeg(n);
	//cerr << h(1, 0, 0) << endl;
	
	vector<ll> ans; ans.assign(n+1, 0);
	
	if(D - C < 0)
	{
		for(int i = 0; i < n; i++)
		{
			cout << ans[i] << ' ';
		}
		cout << endl;
		return 0;
	}

	for(int i = 0; i <= A; i++)
	{
		for(int j = 0; j <= B; j++)
		{
			ll coef = 1;
			coef = mult(coef, F(A-i,B-j,C,D));
			if(A > 0) coef = mult(coef, choose(A, i));
			if(B > 0) coef = mult(coef, choose(B, j));
			vector<int> gi; vector<int> gj;
			gi.resize(n+1); gj.resize(n+1);
			for(int k = 0; k <= n - cycles; k++)
			{
				gi[k] = int(g[D][k]);
				gj[k] = int(sumg[i][j][k]);
			}
			vector<int> res;
			multiply(gi, gj, res);
			//cout << gi.size() << ' ' << gj.size() << ' ' << res.size() << ' ' << ans.size() << ' ' << n << ' ' << cycles << '\n';
			for(int k = 0; k <= n - cycles; k++)
			{
				int moves = n - (k + cycles);
				ans[moves] = add(ans[moves], mult(res[k], coef));
			}
		}
	}
	
	for(int i = 0; i < ans.size(); i++)
	{
		ans[i] = mult(ans[i], fact[D-C]);
	}
	
	for(int i = 0; i < n; i++)
	{
		cout << ans[i] << ' ';
	}
	cout << endl;
	
	return 0;
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<int,int> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int N = 251;
const int MOD = 998244353;
ll inv2;
ll prt;
ll iprt;

ll dpncr[N][N];
ll fact[N];
ll inverse[N];
ll g[N][N];
ll sumg[N][N][N];

vector<ii> perm;
int A, B, C, D;

ll modpow(ll a, ll b)
{
	ll r = 1;
	while(b)
	{
		if(b&1) r = (r*a)%MOD;
		a = (a*a)%MOD;
		b>>=1;
	}
	return r;
}

ll inv(ll a)
{
	return modpow(a, MOD - 2);
}

ll choose(int n, int m)
{
	if(m < 0) return 0;
	if(n < m) return 0;
	if(m == 0) return 1;
	if(n == m) return 1;
	if(dpncr[n][m] != -1) return dpncr[n][m];
	dpncr[n][m] = choose(n - 1, m - 1) + choose(n - 1, m);
	dpncr[n][m] += MOD; dpncr[n][m] %= MOD;
	return dpncr[n][m];
}

void computefact()
{
	fact[0] = 1;
	for(ll i = 1; i < N; i++)
	{
		fact[i] = (fact[i - 1]*i)%MOD;
	}
	for(ll i = 1; i < N; i++)
	{
		inverse[i] = modpow(i, MOD - 2);
	}
}

void print(vector<ii>& vec)
{
	for(int i = 0; i < vec.size(); i++)
	{
		cout << vec[i].fi << ' ' << vec[i].se << endl;
	}
	cout << "------------------------------------------------" << endl;
}

void printans(vector<ll>& vec)
{
	for(int i = 0; i < vec.size(); i++)
	{
		cout << vec[i] << ' ';
	}
	cout << endl;
}

void printansi(vector<int>& vec)
{
	for(int i = 0; i < vec.size(); i++)
	{
		cout << vec[i] << ' ';
	}
	cout << endl;
}

void calcpos(vector<ii>& pos)
{
	pos.resize(perm.size());
	for(int i = 0; i < perm.size(); i++)
	{
		pos[i] = ii(-1, -1);
	}
	for(int i = 1; i < perm.size(); i++)
	{
		if(perm[i].fi > 0)
		{
			pos[perm[i].fi].fi = i;
		}
		if(perm[i].se > 0)
		{
			pos[perm[i].se].se = i;
		}
	}
}

int reduce()
{
	int n = perm.size() - 1;
	vector<ii> pos;
	int cnt = 0;
	for(int i = 1; i <= n; i++) //Do a reduction
	{
		calcpos(pos);
		//print(pos);
		if(pos[i].fi > 0 && pos[i].se > 0)
		{
			if(pos[i].fi == pos[i].se) 
			{
				cnt++;
				ii tmp1 = perm[pos[i].fi];
				for(vector<ii>::iterator it = perm.begin(); it != perm.end(); it++)
				{
					if((*it) == tmp1)
					{
						perm.erase(it); break;
					}
				}
				continue;
			}
			int p1 = pos[i].se; int l = perm[p1].fi; ii tmp1 = perm[p1];
			int p2 = pos[i].fi; int r = perm[p2].se; ii tmp2 = perm[p2];
			for(vector<ii>::iterator it = perm.begin(); it != perm.end(); it++)
			{
				if((*it) == tmp1)
				{
					perm.erase(it); break;
				}
			}
			for(vector<ii>::iterator it = perm.begin(); it != perm.end(); it++)
			{
				if((*it) == tmp2)
				{
					perm.erase(it); break;
				}
			}
			perm.pb(ii(l, r));
		}
	}
	//count A, B, C, D
	for(int i = 1; i < perm.size(); i++)
	{
		if(perm[i].fi > 0 && perm[i].se > 0)
		{
			assert(perm[i].fi != perm[i].se);
			C++;
		}
		else if(perm[i].fi > 0)
		{
			A++;
		}
		else if(perm[i].se > 0)
		{
			B++;
		}
		else
		{
			D++;
		}
	}
	return cnt;
}

ll mult(ll a, ll b)
{
	ll r = (a*b)%MOD;
	r = (r+MOD)%MOD;
	return r;
}

ll add(ll a, ll b)
{
	ll r = ((a+b)%MOD+MOD)%MOD;
	return r;
}

ll F(ll a, ll b, ll c, ll d)
{
	ll ans = 1;
	if(d == 0)
	{
		if(a == 0 && b == 0 && c == 0) return 1;
		else return 0;
	}
	ans = mult(ans, fact[a]);
	ans = mult(ans, fact[b]);
	ans = mult(ans, fact[c]);
	ans = mult(ans, choose(a+d-1, d-1));
	ans = mult(ans, choose(b+d-1, d-1));
	ans = mult(ans, choose(d, C));
	return ans;
}

ll buffer[20001], bufferpos, siz = 1024;
const int LG = 4;

void multiply(int size, ll a[], ll b[], ll r[])
{
	if(size <= (1<<LG))
	{
		for(int i = 0; i < size*2; i++) r[i] = 0;
		for(int i = 0; i < size; i++)
		{
			if(a[i])
			{
				for(int j = 0; j < size; j++)
				{
					r[i+j] += a[i]*b[j];
					r[i+j] %= MOD;
				}
			}
		}
		for(int i = 0; i < size*2; i++)
		{
			r[i] %= MOD;
		}
		return ;
	}
	int s = size/2;
	multiply(s, a, b, r);
	multiply(s, a+s, b+s, r+size);
	ll *a2 = buffer+bufferpos; bufferpos += s;
	ll *b2 = buffer+bufferpos; bufferpos += s;
	ll *r2 = buffer+bufferpos; bufferpos += size;
	for(int i = 0; i < s; i++)
	{
		a2[i] = a[i] + a[i+s];
		if(a2[i]>=MOD) a2[i]-=MOD;
	}
	for(int i = 0; i < s; i++)
	{
		b2[i] = b[i] + b[i+s];
		if(b2[i]>=MOD) b2[i]-=MOD;
	}
	multiply(s, a2, b2, r2);
	for(int i = 0; i < size; i++)
	{
		r2[i] -= (r[i] + r[i+size]);
	}
	for(int i = 0; i < size; i++)
	{
		r[i+s] += r2[i];
		r[i+s]%=MOD;
		if(r[i+s]<0) r[i+s]+=MOD;
	}
	bufferpos -= (s+s+size);
}
ll gi[N+5]; ll gj[N+5];

void computeg(int n)
{
	g[0][0] = 1;
	g[1][1] = 1;
	for(int i = 2; i <= n; i++)
	{
		for(int j = 1; j <= i; j++)
		{
			for(int k = 1; k <= i; k++)
			{
				g[i][j] = add(g[i][j], mult(g[i-k][j-1], mult(choose(i, k), fact[k-1])));
			}
			g[i][j] = mult(g[i][j], inverse[j]);
		}
	}
	for(int i = 0; i <= A; i++)
	{
		for(int j = 0; j <= B; j++)
		{
			siz = 512;
			while(siz/2 >= n+1) siz>>=1;
			for(int k = 0; k < siz; k++)
			{
				if(k <= n)
				{
					gi[k] = g[i][k];
					gj[k] = g[j][k];
				}
				else
				{
					gi[k] = gj[k] = 0;
				}
			}
			ll *res = buffer+bufferpos;
			bufferpos+=2*siz;
			multiply(siz,gi,gj,res);
			for(int k = 0; k <= n; k++)
			{
				sumg[i][j][k] = res[k];
			}
			bufferpos-=2*siz;
		}
	}
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	int n; cin >> n; perm.resize(n+1);
	
	for(int i = 1; i <= n; i++)
	{
		cin >> perm[i].fi;
	}
	for(int i = 1; i <= n; i++)
	{
		cin >> perm[i].se;
	}
	
	A = 0; B = 0; C = 0; D = 0;
	memset(dpncr, -1, sizeof(dpncr));
	memset(sumg, 0, sizeof(sumg));
	memset(g, 0, sizeof(g));
	memset(fact, 0, sizeof(fact));
	memset(inverse, 0, sizeof(inverse));
	int cycles = reduce();
	computefact(); computeg(n);
	vector<ll> ans; ans.assign(n+1, 0);
	
	if(D - C < 0)
	{
		for(int i = 0; i < n; i++)
		{
			cout << ans[i] << ' ';
		}
		cout << endl;
		return 0;
	}
	
	for(int i = 0; i <= A; i++)
	{
		for(int j = 0; j <= B; j++)
		{
			ll coef = 1;
			coef = mult(coef, F(A-i,B-j,C,D));
			if(A > 0) coef = mult(coef, choose(A, i));
			if(B > 0) coef = mult(coef, choose(B, j));
			siz = 512;
			while(siz/2 >= n-cycles+1) siz>>=1;
			for(int k = 0; k < siz; k++)
			{
				if(k <= n-cycles)
				{
					gi[k] = g[D][k];
					gj[k] = sumg[i][j][k];
				}
				else
				{
					gi[k] = gj[k] = 0;
				}
			}
			ll *res = buffer+bufferpos;
			bufferpos+=2*siz;
			multiply(siz,gi,gj,res);
			for(int k = 0; k <= n - cycles; k++)
			{
				int moves = n - (k + cycles);
				ans[moves] = add(ans[moves], mult(res[k], coef));
			}
			bufferpos-=2*siz;
		}
	}
	
	for(int i = 0; i < ans.size(); i++)
	{
		ans[i] = mult(ans[i], fact[D-C]);
	}
	
	for(int i = 0; i < n; i++)
	{
		cout << ans[i] << ' ';
	}
	cout << endl;
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<int,int> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int INF = 1e9 + 7;
const int MOD = 1e9 + 7;
const int N = 1000;
char bus[N][5];

void printbus(int n)
{
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < 5; j++)
		{
			cout << bus[i][j];
		}
		cout << '\n';
	}
}

void yes()
{
	cout << "YES" << '\n';
}

void no()
{
	cout << "NO" << '\n';
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	int n; cin >> n;
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < 5; j++)
		{
			cin >> bus[i][j];
		}
	}
	bool possible = false;
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < 2; j++)
		{
			if(bus[i][j*3] == 'O' && bus[i][j*3+1] == 'O')
			{
				bus[i][j*3] = '+';
				bus[i][j*3+1] = '+';
				possible = true;
				break;
			}
		}
		if(possible) break;
	}
	if(!possible)
	{
		no();
		return 0;
	}
	else
	{
		yes();
		printbus(n);
		return 0;
	}
	return 0;
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<int,int> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int INF = 1e9 + 7;
const int MOD = 1e9 + 7;
const int LG = 20;

ll a[1001][1001];
ll r[1001]; //row sum
ll c[1001]; //column sum
int n;

void no()
{
	cout << -1 << '\n';
}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	cin >> n;
	int x, y; ll diagonal1 = 0; ll diagonal2 = 0;
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < n; j++)
		{
			cin >> a[i][j];
			if(a[i][j] == 0)
			{
				x = i; y = j;
			}
			else
			{
				r[i] += a[i][j];
				c[j] += a[i][j];
				if(i == j)
				{
					diagonal1 += a[i][j];
				}
				if(i + j == n - 1)
				{
					diagonal2 += a[i][j];
				}
			}
		}
	}
	if(n == 1)
	{
		cout << 1 << '\n';
		return 0;
	}
	ll commonsum = r[0];
	if(x == 0) commonsum = r[1];
	//cout << commonsum << '\n';
	ll rowsum = -1; ll colsum = -1; ll d1sum = -1; ll d2sum = -1;
	for(int i = 0; i < n; i++)
	{
		if(i != x)
		{
			if(r[i] != commonsum)
			{
				no();
				return 0;
			}
		}
		else
		{
			rowsum = r[i];
		}
	}
	for(int i = 0; i < n; i++)
	{
		if(i != y)
		{
			if(c[i] != commonsum)
			{
				no(); return 0;
			}
		}
		else
		{
			colsum = c[i];
		}
	}
	bool isdiagonal1 = false; bool isdiagonal2 = false;
	if(x == y) isdiagonal1 = true;
	if(x + y == n - 1) isdiagonal2 = true;
	if(!isdiagonal1)
	{
		if(diagonal1 != commonsum)
		{
			no();
			return 0;
		}
	}
	else
	{
		d1sum = diagonal1;
	}
	if(!isdiagonal2)
	{
		if(diagonal2 != commonsum)
		{
			no();
			return 0;
		}
	}
	else
	{
		d2sum = diagonal2;
	}
	if(rowsum == colsum)
	{
		if(isdiagonal1 && d1sum != rowsum)
		{
			no();
			return 0;
		}
		if(isdiagonal2 && d2sum != rowsum)
		{
			no();
			return 0;
		}
		ll value = commonsum - rowsum;
		if(value > 0)
		{
			cout << value << '\n';
			return 0;
		}
		else
		{
			no();
			return 0;
		}
	}
	else
	{
		no();
		return 0;
	}
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
 
using namespace std;
using namespace __gnu_pbds;
 
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key
 
typedef long long ll;
typedef pair<int,int> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int N = 101;
const int MOD = 1e9 + 7;
const ll INF = ll(1e18);

ll dp[N][N][N];
int c[N];
ll cost[N][N];

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	int n, m, k; cin >> n >> m >> k;
	for(int i = 1; i <= n; i++)
	{
		cin >> c[i];
	}
	for(int i = 0; i <= n; i++)
	{
		for(int j = 0; j <= k; j++)
		{
			for(int a = 0; a <= m; a++)
			{
				dp[i][j][a] = INF;
			}
		}
	}
	for(int i = 1; i <= n; i++)
	{
		for(int j = 1; j <= m; j++)
		{
			cin >> cost[i][j];
		}
	}
	if(c[1] == 0)
	{
		for(int i = 1; i <= m; i++)
		{
			dp[1][1][i] = cost[1][i];
		}
	}
	else
	{
		dp[1][1][c[1]] = 0;
	}
	for(int i = 2; i <= n; i++)
	{
		for(int j = 1; j <= k; j++)
		{
			if(c[i] == 0)
			{
				for(int a = 1; a <= m; a++)
				{
					dp[i][j][a] = min(dp[i][j][a], dp[i-1][j][a] + cost[i][a]);
					for(int b = 1; b <= m; b++)
					{
						if(b != a) dp[i][j][a] = min(dp[i][j][a], dp[i-1][j-1][b] + cost[i][a]);
					}
				}
			}
			else
			{
				dp[i][j][c[i]] = min(dp[i][j][c[i]], dp[i-1][j][c[i]]);
				for(int b = 1; b <= m; b++)
				{
					if(b != c[i]) dp[i][j][c[i]] = min(dp[i][j][c[i]], dp[i-1][j-1][b]);
				}
				//cout << i << ' ' << j << ' ' << c[i] << ' ' << dp[i][j][c[i]] << '\n';
			}
		}
	}
	ll ans = INF;
	for(int i = 1; i <= m; i++)
	{
		ans = min(ans, dp[n][k][i]);
	}
	if(ans >= INF) ans = -1;
	cout << ans;
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
 
using namespace std;
using namespace __gnu_pbds;
 
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key
 
typedef long long ll;
typedef pair<int,int> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int N = 101;
const int MOD = 1e9 + 7;
const ll INF = ll(1e18);

ll dp[N][N][N];
int c[N];
ll cost[N][N];
ll idx[N][N];
ll m1[N][N];
ll m2[N][N];

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0);
	int n, m, k; cin >> n >> m >> k;
	for(int i = 1; i <= n; i++)
	{
		cin >> c[i];
	}
	for(int i = 0; i <= n; i++)
	{
		for(int j = 0; j <= k; j++)
		{
			m1[i][j] = INF; m2[i][j] = INF; idx[i][j] = -1;
			for(int a = 0; a <= m; a++)
			{
				dp[i][j][a] = INF;
			}
		}
	}
	for(int i = 1; i <= n; i++)
	{
		for(int j = 1; j <= m; j++)
		{
			cin >> cost[i][j];
		}
	}
	if(c[1] == 0)
	{
		for(int i = 1; i <= m; i++)
		{
			dp[1][1][i] = cost[1][i];
			if(dp[1][1][i] <= m1[1][1])
			{
				if(dp[1][1][i] == m1[1][1])
				{
					idx[1][1] = -2;
				}
				else
				{
					idx[1][1] = i;
				}
				m2[1][1] = m1[1][1];
				m1[1][1] = dp[1][1][i];
			}
			else if(dp[1][1][i] <= m2[1][1])
			{
				m2[1][1] = dp[1][1][i];
			}
		}
	}
	else
	{
		dp[1][1][c[1]] = 0;
		m1[1][1] = 0; idx[1][1] = c[1];
	}
	for(int i = 2; i <= n; i++)
	{
		for(int j = 1; j <= k; j++)
		{
			if(c[i] == 0)
			{
				for(int a = 1; a <= m; a++)
				{
					dp[i][j][a] = min(dp[i][j][a], dp[i-1][j][a] + cost[i][a]);
					ll tmp = INF;
					if(a == idx[i-1][j-1])
					{
						tmp = m2[i-1][j-1];
					}
					else
					{
						tmp = m1[i-1][j-1];
					}
				    dp[i][j][a] = min(dp[i][j][a], tmp + cost[i][a]);
				}
			}
			else
			{
				dp[i][j][c[i]] = min(dp[i][j][c[i]], dp[i-1][j][c[i]]);
				for(int b = 1; b <= m; b++)
				{
					if(b != c[i]) dp[i][j][c[i]] = min(dp[i][j][c[i]], dp[i-1][j-1][b]);
				}
				//cout << i << ' ' << j << ' ' << c[i] << ' ' << dp[i][j][c[i]] << '\n';
			}
			for(int a = 1; a <= m; a++)
			{
				if(dp[i][j][a] <= m1[i][j])
				{
					if(dp[i][j][a] == m1[i][j])
					{
						idx[i][j] = -2;
					}
					else
					{
						idx[i][j] = a;
					}
					m2[i][j] = m1[i][j];
					m1[i][j] = dp[i][j][a];
				}
				else if(dp[i][j][a] <= m2[i][j])
				{
					m2[i][j] = dp[i][j][a];
				}
			}
		}
	}
	ll ans = INF;
	for(int i = 1; i <= m; i++)
	{
		ans = min(ans, dp[n][k][i]);
	}
	if(ans >= INF) ans = -1;
	cout << ans;
}

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

using namespace std;
using namespace __gnu_pbds;

#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define fbo find_by_order
#define ook order_of_key

typedef long long ll;
typedef pair<int,int> ii;
typedef vector<int> vi;
typedef long double ld; 
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> pbds;
typedef set<int>::iterator sit;
typedef map<int,int>::iterator mit;
typedef vector<int>::iterator vit;

const int INF = 1e9 + 7;
const int MOD = 1e9 + 7;
const int N = 1e6 + 3;

int a[N];
int visited[N];
ll ans;
vector<int> cycles;
ll dp[N];
int cyclecnt;

void dfs2(int u)
{
	cycles[cyclecnt]++;
	visited[u] = 3;
	if(visited[a[u]] == 3) return ;
	dfs2(a[u]);
}

void dfs(int u)
{
	visited[u] = 2;
	if(visited[a[u]] == 0)
	{
		dfs(a[u]);
	}
	else if(visited[a[u]] == 1)
	{
		visited[u] = 1;
		return ;
	}
	else
	{
		cycles.pb(0);
		dfs2(u);
		cyclecnt++;
	}
	visited[u] = 1;
}

int main()
{
	//ios_base::sync_with_stdio(0); cin.tie(0);
	int n; scanf("%d", &n);
	for(int i = 1; i <= n; i++)
	{
		scanf("%d", a + i);
	}
	dp[0] = 1;
	for(int i = 1; i <= n; i++)
	{
		dp[i] = (dp[i-1]*2LL)%MOD;
	}
	ans = 1;
	memset(visited, 0, sizeof(visited));
	for(int i = 1; i <= n; i++)
	{
		if(visited[i] == 0)
		{
			dfs(i);
		}
	}
	ll cnt = n;
	for(int i = 0; i < cycles.size(); i++)
	{
		cnt -= cycles[i];
		ans = (ans*(dp[cycles[i]]-2+MOD))%MOD;
	}
	ans = (ans*dp[cnt])%MOD;
	if(ans < 0) ans += MOD;
	int ans2 = ans;
	printf("%d\n", ans2);
	return 0;
}

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;
typedef vector<int> vi;

const int MOD = 1e6 + 3;

ll power(ll base, ll exp)
{
	ll ans = 1;
    while(exp)
    {
		if(exp&1) ans = (ans*base)%MOD;
		base = (base*base)%MOD;
		exp>>=1;
	}
    return ans;
}

int main()
{
	ios_base::sync_with_stdio(false); cin.tie(0);
	ll n, k;
	cin >> n >> k;
	if(n <= 63 && k > (1LL<<n))
	{
		cout << 1 << " " << 1;
		return 0;
	}
	ll v2 = 0;
	int digits = __builtin_popcountll(k - 1);
	v2 = k - 1 - digits;
	ll ntmp = n % (MOD - 1);
	if(ntmp < 0) ntmp += (MOD - 1);
	ll ktmp = k % (MOD - 1);
	if(ktmp < 0) ktmp += (MOD - 1);
	ll v2tmp = v2 % (MOD - 1);
	if(v2tmp < 0) v2tmp += (MOD - 1);
	ll exponent = ntmp*(ktmp - 1) - v2tmp;
	exponent %= (MOD - 1);
	if(exponent < 0) exponent += MOD - 1;
	ll denom = power(2, exponent);
	ll numpart = 0;
	if(k - 1 >= MOD)
	{
		numpart = 0;
	}
	else
	{
		ll prod = 1;
		ll ntmp2 = power(2, ntmp);
		prod = power(2, v2tmp);
		prod = power(prod, MOD - 2);
		if(prod < 0) prod += MOD;
		for(ll y = 1; y <= k - 1; y++)
		{
			prod = (prod * (ntmp2 - y))%MOD;
		}
		numpart = prod;
	}
	ll num = (denom - numpart)%MOD;
	num %= MOD; denom %= MOD;
	if(num < 0) num += MOD;
	if(denom < 0) denom += MOD;
	cout << num << " " << denom;
	return 0;
}