313 lines
7.3 KiB
Text
313 lines
7.3 KiB
Text
---
|
|
id: fft-ext
|
|
title: "More Complex Operations Using FFT"
|
|
author: Benjamin Qi
|
|
prerequisites:
|
|
description: "?"
|
|
frequency: 0
|
|
---
|
|
|
|
import { Problem } from "../models";
|
|
|
|
export const metadata = {
|
|
problems: {
|
|
ys: [
|
|
new Problem("YS", "Partition Function", "partition_function", "?", false, [], ""),
|
|
new Problem("YS", "Bernoulli Number", "bernoulli_number", "?", false, [], ""),
|
|
],
|
|
general: [
|
|
new Problem("Plat", "Tree Depth", "974", "Hard", false, ["genfunc"], ""),
|
|
new Problem("Plat", "Exercise", "1045", "Hard", false, ["permutation"], "genfuncs not required but possibly helpful"),
|
|
],
|
|
}
|
|
};
|
|
|
|
## More Complex Operations
|
|
|
|
<resources>
|
|
<resource source="cp-algo" title="Operations on Polynomials & Series" url="algebra/polynomial.html"> </resource>
|
|
</resources>
|
|
|
|
### Implementations
|
|
|
|
<resources>
|
|
<resource source="Benq" title="Polys" url="https://github.com/bqi343/USACO/tree/master/Implementations/content/numerical/Polynomials"> </resource>
|
|
</resources>
|
|
|
|
### Problems
|
|
|
|
<problems-list problems={metadata.problems.ys} />
|
|
|
|
## Counting
|
|
|
|
No advanced knowledge about generating functions is required for either of these problems.
|
|
|
|
<problems-list problems={metadata.problems.general} />
|
|
|
|
- [zscoder GenFunc Pt 1](https://codeforces.com/blog/entry/77468)
|
|
- [zscoder GenFunc Pt 2](https://codeforces.com/blog/entry/77551)
|
|
|
|
## [F2: Slime & Sequences](https://codeforces.com/contest/1349/problem/F2)
|
|
|
|
Providing some code to make the explanation for the last problem from zscoder's GenFunc Pt 2 clearer. It's quite easy to lose some factor of $y$ or $z$ somewhere during the calculation and end up with gibberish.
|
|
|
|
Of course, this will be split into several parts to make it easier to follow.
|
|
|
|
### Part 0: Includes
|
|
|
|
<warning-block>400+ lines of polynomial template</warning-block>
|
|
|
|
(link?)
|
|
|
|
### Part 1: What do we want to calculate?
|
|
|
|
The problem reduces to finding $n!$ times the coefficients of $y^1,y^2,\ldots,y^n$ in
|
|
$$
|
|
(1-y)[x^n]\frac{1}{(1-x)(1-ye^{(1-y)x})}.
|
|
$$
|
|
Since we don't care about powers of $y$ greater than $n$, we can expand:
|
|
$$
|
|
\frac{1}{(1-ye^{(1-y)x})}=\sum_{t=0}^ny^te^{t(1-y)x}
|
|
$$
|
|
and
|
|
$$
|
|
e^{t(1-y)x}=\frac{\sum_{i=0}^{\infty}x^it^i(1-y)^i}{i!}.
|
|
$$
|
|
We can naively implement this as follows. This gets MLE on F1, but at least it gives the correct results.
|
|
|
|
|
|
<spoiler title="Slow Solution">
|
|
|
|
```cpp
|
|
int n;
|
|
poly py[5005]; // powers of 1-y
|
|
|
|
vector<poly> po(int t) { // y^t * e^{t*(1-y)*x}
|
|
vector<poly> res(n+1);
|
|
F0R(i,n+1) res[i] = RSZ(shift(pow(mi(t),i)*py[i]*ifac[i],t),n+1);
|
|
return res;
|
|
}
|
|
|
|
poly brute() {
|
|
py[0] = {1}; FOR(i,1,n+1) py[i] = poly{1,-1}*py[i-1];
|
|
vector<poly> ans(n+1); // store polynomial in y for each degree of x
|
|
F0R(i,n+1) {
|
|
auto a = po(i); // dbg("OOPS",i,a);
|
|
F0R(j,n+1) ans[j] += a[j];
|
|
}
|
|
poly res; F0R(i,n+1) res += ans[i];
|
|
res = RSZ(res*poly{1,-1},n+1);
|
|
FOR(i,1,n+1) res[i] *= fac[n];
|
|
return vmi(1+all(res));
|
|
}
|
|
|
|
int main() {
|
|
re(n); genFac(n+10);
|
|
poly res = brute();
|
|
FOR(i,1,n+1) {
|
|
res[i] *= fac[n];
|
|
pr(res[i],' ');
|
|
}
|
|
ps();
|
|
}
|
|
|
|
```
|
|
|
|
</spoiler>
|
|
|
|
### Part 2: Setting up a function to compute the same thing more quickly
|
|
|
|
Now we'll make a function `smart()` that does the same thing as `brute()` (but quickly). We need to compute
|
|
|
|
$$
|
|
(1-y)^{n+2}[z^n]\frac{1}{(1-y-z)(1-ye^z)}
|
|
$$
|
|
|
|
It suffices to compute
|
|
$$
|
|
[z^n]\frac{1}{(1-y-z)(1-ye^z)}=[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).
|
|
$$
|
|
|
|
We'll let `frac1()` compute the polynomial in $y$ corresponding to
|
|
|
|
$$
|
|
[z^n]n\cdot \frac{1}{1-e^z(1-z))(1-y-z)}
|
|
$$
|
|
|
|
and `frac2()` compute the polynomial in $y$ corresponding to
|
|
|
|
$$
|
|
[z^n]n\cdot \frac{e^z}{(1-e^z(1-z))(1-ye^z)}.
|
|
$$
|
|
|
|
The reasons for the factors of $n$ will become apparent later on.
|
|
|
|
```cpp
|
|
poly smart() {
|
|
poly res = frac1()-frac2(); res /= n;
|
|
poly py; // (1-y)^{n+2}
|
|
F0R(i,n+3) {
|
|
py.pb(comb(n+2,i));
|
|
if (i&1) py.bk *= -1;
|
|
}
|
|
return RSZ(conv(res,py),n+1); // multiply res by (1-y)^{n+1}
|
|
}
|
|
```
|
|
|
|
### Part 3: Finding $[x^n]\frac{P(x)}{1-(x+1)y}$
|
|
|
|
This is equal to
|
|
$$
|
|
[x^n]P(x)\sum_{i\ge 0}(x+1)^iy^i.
|
|
$$
|
|
Letting $Q(x)=\sum_{i\le n}p_ix^{n-i}$, we have
|
|
$$
|
|
[x^ny^k]P(x)\sum_{i\ge 0}(x+1)^iy^i=\sum_{i=0}^{\infty}q_i\binom{k}{i}.
|
|
$$
|
|
|
|
```cpp
|
|
poly coef1(poly p, int n) {
|
|
p.rsz(n+1); reverse(all(p));
|
|
poly difs; F0R(i,n+1) difs.pb(ifac[i]);
|
|
F0R(i,n+1) p[i] *= ifac[i];
|
|
p = conv(p,difs);
|
|
F0R(i,n+1) p[i] *= fac[i];
|
|
return RSZ(p,n+1);
|
|
}
|
|
```
|
|
|
|
### Part 4: Finding $[x^n]\frac{P(x)}{(1-(x+1)y)^2}$
|
|
|
|
This is equal to
|
|
$$
|
|
[x^n]P(x)\sum_{i\ge 0}(i+1)(x+1)^iy^i.
|
|
$$
|
|
We get the same expression as above except the coefficient of $y^i$ is multiplied by $i+1$.
|
|
|
|
```cpp
|
|
poly coef2(poly p, int n) {
|
|
poly q = coef1(p,n);
|
|
F0R(i,n+1) q[i] *= i+1;
|
|
return q;
|
|
}
|
|
```
|
|
|
|
### Part 5: First Fraction
|
|
|
|
Letting $g(z)=\frac{1}{(1-z)(1-e^z(1-z))}$, we want to compute $[z^n]\frac{g(z)}{1-\frac{y}{1-z}}$.
|
|
|
|
Let
|
|
$$
|
|
B(z)=\frac{1}{1-z}-1, B^{-1}(z)=\frac{z}{z+1}.
|
|
$$
|
|
Then
|
|
$$
|
|
C(x)=g(B^{-1}(x))=\frac{1}{\frac{1}{1+x}\cdot (1-\frac{e^{x/(1+x)}}{1+x})}=\frac{(1+x)^2}{1+x-e^{x/(x+1)}}.
|
|
$$
|
|
$$
|
|
D(x)=\left(\frac{1}{1+x}\right)^{-n}=(1+x)^n
|
|
$$
|
|
|
|
We want $\frac{1}{n}$ times the following:
|
|
$$
|
|
[x^{n-1}]\left(\frac{C'(x)D(x)}{1-(x+1)y}+\frac{C(x)D(x)y}{(1-(x+1)y)^2}\right).
|
|
$$
|
|
Note that $x^2$ divides the denominator of $C(x)$ so we need to be careful about how we manipulate it.
|
|
|
|
```cpp
|
|
poly deriv(poly a, int b) { // 0-th element of vector corresponds to x^{b}
|
|
assert(b < 0 && a[0] != 0);
|
|
poly ans; F0R(i,sz(a)) ans.pb((i+b)*a[i]);
|
|
return ans;
|
|
}
|
|
|
|
// [z^n]1/(1-e^z(1-z))/(1-y-z)
|
|
poly frac1() {
|
|
poly C, D;
|
|
{
|
|
poly ex{0}; FOR(i,1,n+5) ex.pb(i&1?1:-1);
|
|
poly dem = poly{1,1}-exp(ex,n+5);
|
|
assert(dem[0] == 0 && dem[1] == 0 && dem[2] != 0);
|
|
dem = poly(2+all(dem));
|
|
dem = inv(dem,n+5);
|
|
C = poly{1,1}*poly{1,1}*dem;
|
|
}
|
|
{
|
|
F0R(i,n+1) D.pb(comb(n,i));
|
|
}
|
|
poly C2 = deriv(C,-2);
|
|
poly X = conv(C2,D); // lowest deg term is x^{-3}
|
|
poly Y = conv(C,D); // lowest deg term is x^{-2}
|
|
return coef1(X,n+2)+shift(coef2(Y,n+1),1);
|
|
}
|
|
```
|
|
|
|
### Part 6: Second Fraction
|
|
|
|
Let $f(z)=\frac{e^z}{1-e^z(1-z)}$. We want to compute $[z^n]\frac{f(z)}{1-ye^z}$.
|
|
|
|
Let
|
|
$$
|
|
A(z)=e^z-1, A^{-1}(z)=\ln (z+1).
|
|
$$
|
|
Then by Lagrange Inversion (??)
|
|
$$
|
|
[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}.
|
|
$$
|
|
Ok so let
|
|
$$
|
|
C(x)=f(A^{-1}(x))=\frac{1+x}{1-(1+x)(1-\ln(x+1))}
|
|
$$
|
|
and
|
|
$$
|
|
D(x)=\frac{1}{\left(\frac{A^{-1}(x)}{x}\right)^n}=\left(\frac{\ln(x+1)}{x}\right)^{-n}.
|
|
$$
|
|
and we want $\frac{1}{n}$ times the following:
|
|
$$
|
|
[x^{n-1}]\left(\frac{C'(x)D(x)}{1-(x+1)y}+\frac{C(x)D(x)y}{(1-(x+1)y)^2}\right).
|
|
$$
|
|
Note that $x^2$ divides the denominator of $C(x)$ so we need to be careful about how we manipulate it.
|
|
|
|
```cpp
|
|
poly po(poly a, int b) {
|
|
assert(a[0] == 1);
|
|
poly x = log(a,n+5);
|
|
return exp(b*x,n+5);
|
|
}
|
|
|
|
// [z^n]e^z/(1-e^z(1-z))(1-ye^z)
|
|
poly frac2() {
|
|
poly LN; // x-x^2/2+x^3/3-x^4/4+...
|
|
LN.pb(0);
|
|
FOR(i,1,n+5) {
|
|
LN.pb(invs[i]);
|
|
if (i%2 == 0) LN.bk *= -1;
|
|
}
|
|
poly C, D;
|
|
{ // calculating C*x^2
|
|
poly oops = poly{1}-LN;
|
|
poly dem = poly{1}-poly{1,1}*oops;
|
|
assert(dem[0] == 0 && dem[1] == 0 && dem[2] != 0);
|
|
dem = poly(2+all(dem));
|
|
dem = inv(dem,n+5);
|
|
C = poly{1,1}*dem;
|
|
}
|
|
{
|
|
poly oops = LN; oops.erase(begin(oops));
|
|
D = po(oops,-n);
|
|
}
|
|
poly C2 = deriv(C,-2);
|
|
poly X = conv(C2,D); // x^{-3}
|
|
poly Y = conv(C,D); // x^{-2}
|
|
return coef1(X,n+2)+shift(coef2(Y,n+1),1);
|
|
}
|
|
```
|
|
|
|
### Part 7: Summary
|
|
|
|
All in all, this solution (barely) passes the time limit. :| |