【NOI2019模拟2019.6.27】B （生成函数+整数划分dp|多项式exp）

Description：

---

(1<=n,k<=1e5,mod~1e9+7)

题解：

---

考虑最经典的排列dp，每次插入第(i)大的数，那么可以增加的逆序对个数是(0-i-1)。

不难得到生成函数：

(Ans=prod_{i=0}^{n-1}(sum_{j=0}^ix^j)[x^k])

(=prod_{i=1}^{n}{1-x^iover 1-x}[x^k])

分母是一个经典的生成函数：

({1over 1-x}^n=(sum_{i>=0}x^i)^n=sum_{i>=0}C_{i+n-1}^{n-1})

那么问题就变为了求：

(prod_{i=1}^{n}{1-x^i})的前k项。

考虑利用整数划分dp，相当于把k划分成若干不同且<=n的数和，系数是((-1)^{数的个数})。

不难得出dp：

设(f[i][j])表示已经划分了i个数，和为j的所有方案系数和。

转移(f[i][j]=f[i][j-i]-f[i-1][j-i]+f[i-1][j-(n+1)])

由于(i<=sqrt {2k})，所以复杂度是(O(ksqrt k))。

另一种多项式exp的做法：

不妨对这个式子进行ln最后再exp回去。

我们知道有：

(ln(1+x))

$=int ~ln(1+x)' $

(=int~{1over 1+x})

(=int ~ sum_{i>=0}(-1)^ix^i)

(=sum_{i>=1}{(-1)^{i-1}x^iover i})

(Ans=exp(sum_{i=1}^n(ln(1-x^i)-ln(1-x)))[x^k])

(ln(1-x^i)=-sum_{j>=1}{x^{ij} over j})

所以暴力展开只有调和级数个有用项。

(ln(1-x))同理暴力展开后乘上n即可。

复杂度(O(n~log~n))，但是要写MTT，所以跑得巨慢，又难写。

Code：

#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, B = y; i <= B; i ++)
#define ff(i, x, y) for(int i = x, B = y; i <  B; i ++)
#define fd(i, x, y) for(int i = x, B = y; i >= B; i --)
#define ll long long
#define pp printf
#define hh pp("
")
using namespace std;

const int mo = 1e9 + 7;

ll ksm(ll x, ll y) {
	ll s = 1;
	for(; y; y /= 2, x = x * x % mo)
		if(y & 1) s = s * x % mo;
	return s;
}

const int N = 1e5 + 5;

int n, k, m;
ll fac[N * 2], nf[N * 2];
ll f[450][N];
ll g[N];
void calc(int n) {
	fac[0] = 1; fo(i, 1, n) fac[i] = fac[i - 1] * i % mo;
	nf[n] = ksm(fac[n], mo - 2); fd(i, n, 1) nf[i - 1] = nf[i] * i % mo;
}
ll C(int n, int m) {
	return fac[n] * nf[n - m] % mo * nf[m] % mo;
}

int main() {
	freopen("b.in", "r", stdin);
	freopen("b.out", "w", stdout);
	calc(200000);
	scanf("%d %d", &n, &k);
	m = sqrt(2 * k);
	f[0][0] = 1;
	fo(i, 1, m) {
		fo(j, i, k) {
			f[i][j] = f[i][j - i] - f[i - 1][j - i];
			if(j >= n + 1) f[i][j] += f[i - 1][j - (n + 1)];
			f[i][j] %= mo;
		}
	}
	ll ans = 0;
	fo(i, 0, k) {
		fo(j, 0, m) g[i] += f[j][i];
		g[i] %= mo;
		ans += g[i] * fac[n - 1 + (k - i)] % mo * nf[k - i] % mo;
	}
	ans = (ans % mo * nf[n - 1] % mo + mo) % mo;
	pp("%lld
", ans);
}

#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, B = y; i <= B; i ++)
#define ff(i, x, y) for(int i = x, B = y; i <  B; i ++)
#define fd(i, x, y) for(int i = x, B = y; i >= B; i --)
#define ll long long
#define pp printf
#define hh pp("
")
#define db double
using namespace std;

const int mo = 1e9 + 7;

typedef vector<ll> V;
#define si size()
#define pb push_back

namespace ntt {
	const db pi = acos(-1);
	struct P {
		db x, y;
		P(db _x = 0, db _y = 0) { x = _x, y = _y;}
	};
	P operator + (P a, P b) { return P(a.x + b.x, a.y + b.y);}
	P operator - (P a, P b) { return P(a.x - b.x, a.y - b.y);}
	P operator * (P a, P b) { return P(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}
	const int nm = 1 << 18;
	int r[nm];
	P w[nm], c0[nm], c1[nm], c2[nm], c3[nm];

	void dft(P *a, int n) {
		ff(i, 0, n) {
			r[i] = r[i / 2] / 2 + (i & 1) * (n / 2);
			if(i < r[i]) swap(a[i], a[r[i]]);
		} P b;
		for(int i = 1; i < n; i *= 2) for(int j = 0; j < n; j += 2 * i)
			ff(k, 0, i) b = a[i + j + k] * w[i + k], a[i + j + k] = a[j + k] - b, a[j + k] = a[j + k] + b;
	}
	void rev(P *a, int n) {
		reverse(a + 1, a + n);
		ff(i, 0, n) a[i].x /= n, a[i].y /= n;
	}
	P conj(P a) { return P(a.x, -a.y);}
	void fft(V &a, V b) {
		#define qz(x) ((ll) round(x))
		int n = a.si;
//		ff(i, 0, n) c0[i] = P(a[i], 0), c1[i] = P(b[i], 0);
//		dft(c0, n); dft(c1, n);
//		ff(i, 0, n) c0[i] = c0[i] * c1[i];
//		dft(c0, n); rev(c0, n);
//		ff(i, 0, n) a[i] = qz(c0[i].x) % mo;
//		return;
		ff(i, 0, n) {
			c0[i] = P(a[i] & 32767, a[i] >> 15);
			c1[i] = P(b[i] & 32767, b[i] >> 15);
		}
		dft(c0, n); dft(c1, n);
		ff(i, 0, n) {
			int j = (n - i) & (n - 1);
			P k, d0, d1, d2, d3;
			k = conj(c0[j]);
			d0 = (k + c0[i]) * P(0.5, 0);
			d1 = (k - c0[i]) * P(0, 0.5);
			k = conj(c1[j]);
			d2 = (k + c1[i]) * P(0.5, 0);
			d3 = (k - c1[i]) * P(0, 0.5);
			c2[i] = d0 * d2 + d1 * d3 * P(0, 1);
			c3[i] = d1 * d2 + d0 * d3;
		}
		dft(c2, n); dft(c3, n);
		rev(c2, n); rev(c3, n);
		ff(i, 0, n) {
			a[i] = qz(c2[i].x) + (qz(c2[i].y) % mo << 30) + (qz(c3[i].x) % mo << 15);
			a[i] %= mo;
		}
	}
	void build() {
		for(int i = 1; i < nm; i *= 2) ff(j, 0, i)
			ntt :: w[i + j] = P(cos(pi * j / i), sin(pi * j / i));	
	}
}

ll ksm(ll x, ll y) {
	ll s = 1;
	for(; y; y /= 2, x = x * x % mo)
		if(y & 1) s = s * x % mo;
	return s;
}

V operator + (V a, V b) {
	a.resize(max(a.si, b.si));
	ff(i, 0, a.si) a[i] = (a[i] + b[i]) % mo;
	return a;
}
V operator - (V a, V b) {
	a.resize(max(a.si, b.si));
	ff(i, 0, a.si) a[i] = (a[i] - b[i] + mo) % mo;
	return a;
}
V operator * (V a, int x) {
	ff(i, 0, a.si) a[i] = a[i] * x % mo;
	return a;
}
V operator * (V a, V b) {
	int sa = a.si + b.si - 1, n = 1;
	for(; n <= sa; n *= 2);
	a.resize(n); b.resize(n);
	ntt :: fft(a, b);
	a.resize(sa);
	return a;
}

V qni(V a) {
	int n0 = 1; for(; n0 < a.si; n0 *= 2);
	V b; b.resize(1); b[0] = ksm(a[0], mo - 2);
	for(int n = 2; n <= n0; n *= 2) {
		V c = a; c.resize(n);
		c = c * b; c.resize(n); c = c * b; c.resize(n);
		b = b * 2 - c;
	}
	b.resize(a.si);
	return b;
}

V qd(V a) {
	ff(i, 1, a.si) a[i - 1] = a[i] * i % mo;
	a.resize(a.si - 1);
	return a;
}
V jf(V a) {
	a.pb(0);
	fd(i, a.si - 1, 1) a[i] = a[i - 1] * ksm(i, mo - 2) % mo;
	a[0] = 0;
	return a;
}
V ln(V a) {
	int sa = a.si;
	a = jf(qni(a) * qd(a));
	a.resize(sa);
	return a;
}
V exp(V a) {
	int n0 = 1; for(; n0 < a.si; n0 *= 2);
	V b; b.resize(1); b[0] = 1;
	for(int n = 2; n <= n0; n *= 2) {
		V c = b; c.resize(n);
		V d = a; d.resize(n);
		c = d - ln(c); c[0] ++;
		b = b * c; b.resize(n);
	}
	b.resize(a.si);
	return b;
}

V a;

const int N = 1e5 + 5;

int n, k;

ll ni[N];

int main() {
	freopen("b.in", "r", stdin);
	freopen("b.out", "w", stdout);
	ntt :: build();
	n = 1e5;
	fo(i, 1, n) ni[i] = ksm(i, mo - 2);
	scanf("%d %d", &n, &k);
	a.clear(); a.resize(k + 1);
	fo(j, 1, k) a[j] = ni[j] % mo * (n - 1) % mo;
	fo(i, 2, n) {
		fo(j, 1, k / i) a[i * j] -= ni[j];
	}
	fo(j, 1, k) a[j] %= mo;
	a = exp(a);
	pp("%lld
", (a[k] + mo) % mo);
}