This repository has been archived on 2022-06-22. You can view files and clone it, but cannot push or open issues or pull requests.
usaco-guide/content/7_Advanced/FFT-ext.mdx
2020-06-26 17:12:04 -04:00

305 lines
No EOL
7.1 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: {
general: [
new Problem("Plat", "Tree Depth", "974", "Hard", false, ["genfunc"], ""),
new Problem("Plat", "Exercise", "1045", "Hard", false, ["permutation"], "genfuncs not required but possibly helpful"),
],
}
};
<info-block title="Pro Tip">
Does not appear on USACO.
</info-block>
## More Complex Operations
- [cp-algo - Operations on Polynomials & Series](https://cp-algorithms.com/algebra/polynomial.html)
- [My Implementations](https://github.com/bqi343/USACO/tree/master/Implementations/content/numerical/Polynomials)
## yosupo
https://judge.yosupo.jp/problem/partition_function
https://judge.yosupo.jp/problem/bernoulli_number
## Counting
<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. :|