---
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. :|