整数划分 题目分析以及衍生出来的一系列做法
整数划分 题目分析以及衍生出来的一系列做法
推荐看背包计数问题的多项式优化。
题目概述
将 \(n\) 分为若干个不同整数的和，问有多少种方案。
分析题目性质
想一想，注意到是不同整数的和，也就是说我们分成的 \(k\) 部满足：
\[k(k+1)\leq 2n
\]
约等一下，\(k\leq \sqrt{2n}+1\)，即 \(k_{\max}=\sqrt{2n}+1.\)
因此可以考虑 \(\mathcal{O}(nk)\) 算法。
思路
动态规划小优化
设 \(f_{i,j}\) 表示数字 \(i\) 被分为 \(j\) 个部分的方案数。
考虑怎么转移。
为了使复杂度不过高，我们可以观察 \(1\) 的贡献。
若前面的方案中有 \(1\)，我们强制使前面的方案每个数 \(+1\) 变成当前方案再加上 \(1\) 的方案，即 \(f_{i-j,j-1}.\)
第二种就是使前面的方案每个数 \(+1\) 但是我不加 \(1\) 的方案，即 \(f_{i-j,j}.\)
显然有：
\[f_{i,j}=f_{i-j,j}+f_{i-j,j-1}
\]
于是这道题就做完了。
#include <cstdio>
#include <cstring>
#include <iostream>
#include <cmath>
#define int long long
#define N 50005
using namespace std;
const int mod = 1e9 + 7;
int p[N][405];
signed main() {
int n, k;
while (~scanf("%lld", &n)) {// 分成 k 个部分。
memset(p, 0, sizeof(p));
p[0][0] = 1;
k = sqrt(2 * n) + 1;// 本题可以分成 sqrt2n
int i;
for (i = 1;i <= n;i ++)
for (int j = 1;j <= k;j ++)
if (i - j >= 0)
p[i][j] = (p[i - j][j] + p[i - j][j - 1]) % mod;
int res = 0;
for (int i = 1;i <= k;i ++) res = (res + p[n][i]) % mod;
printf("%lld\n", res);
}
return 0;
}
但细想一下，本题还可以衍生出许多做法，甚至可以达到 \(\mathcal{O}(n\log n).\)
生成函数思路
类似背包计数
\(dp\) 这么做似乎没什么前途，我们考虑生成函数（或者说使母函数）。
本质上是求：
\[[x^n]\prod_{k=1}^\infty (1+x^k)
\]
我们发现这不类似于 01 背包计数问题，有兴趣的可以参考一下引用，其为更加一般的形式。
其的生成函数为：
\[[x^m]\prod_{k=1}^n(1+x^{a_k}) \]我们还是可以想到之前的一道题，下面会讲。
事实上，这个式子与 \(\prod_{k=1}^n(x+a_k)\) 很像，这个式子不难用分治 NTT 解决，但是对于上述的式子，其次数过高，不好处理（底层会炸）。
注意到有：
\[xy=\exp\left(\ln\left(xy\right)\right)=\exp\left(\ln x+\ln y\right)
\]
其中 \(\exp y=e^y,\ln y=\log_ey.\)
所以原式不难化为：
\[\prod_{k=1}^\infty (1+x^k)=\exp\left(\sum_{k=1}^\infty\ln\left(1+x^k\right)\right)
\]
根据麦克劳森公式，或者泰勒公式（可以考虑用倒数和泰勒级数证明）：
\[\ln(1+x)=x-\frac{x^2}2+\frac{x^3}3-\dots+(-1)^{(n+1)}\frac{x^n}n=\sum_{k=1}^n(-1)^{k+1}\frac{x^k}k
\]
可以将 \(x\) 替换成 \(x^k\)：
\[\ln(1+x^k)=\sum_{p=1}^n(-1)^{p+1}\frac{x^{pk}}k
\]
在总量为 \(n\) 的条件下（背包是 \(m\)），这个数 \(k\) 的式子中有用的项只有 \(\frac{n}k\)，其余项为 \(0.\)
记录每个数有几个，那么每个数最多只会计算一次，而 \(\sum\frac{n}k\) 的复杂度是调和级数 \(n\log n\) 的，因此，我们用 \(\mathcal{O}(n\log n)\) 的时间复杂度解决了这个问题。
纯数学推导
欲求：
\[[x^n]\prod_{k=1}^\infty (1+x^k)=[x^n]\prod_{k=1}^\infty\frac{1-x^{2k}}{1-x^k}
\]
上述等式成立的证明如下：
\[\begin{align} (1+x^k) & = \sum_{r\geq 0}x^{rk}-\sum_{r\geq 2}x^{rk}\\ & = \sum_{r\geq 0}\left(x^{rk}-x^{(r+2)k}\right)\\ & = \sum_{r\geq 0}x^{rk}\left(1-x^{2k}\right)\\ & = \frac{1-x^{2k}}{1-x^k} \end{align} \]套个 \(\prod\) 就行，于是由上可知证毕。
根据欧拉五边形数定理：
\[\prod_{k=1}^\infty\left(1-x^n\right)=\sum_{k\in \Z}\left(-1\right)^kx^{k\left(3k-1\right)\div2}
\]
于是原式 \(S\) 就有：
\[S=[x^n]\frac{\sum_{k\in\Z}\left(-1\right)^kx^{k(3k-1)}}{\sum_{k\in\Z}\left(-1\right)^kx^{k(3k-1)\div2}}
\]
于是可以对分母求逆，暴力是 \(\mathcal{O}(n^\frac{3}{2})\) 的，考虑用 NTT 优化有 \(\mathcal{O}(n\log n).\)
ai 的代码如下：
#include <bits/stdc++.h>
using namespace std;
const int MOD = 1000000007;
// A function to add two mod numbers.
inline int add_mod(int a, int b) {
int res = a + b;
if(res >= MOD) res -= MOD;
return res;
}
// A function to subtract two mod numbers.
inline int sub_mod(int a, int b) {
int res = a - b;
if(res < 0) res += MOD;
return res;
}
// A function to multiply two mod numbers.
inline int mul_mod(int a, int b) {
return (int)((1LL * a * b) % MOD);
}
// The main function.
int main(){
ios::sync_with_stdio(false);
cin.tie(nullptr);
int n;
cin >> n;
// We need two power series up to degree n:
// Denom: D(x) = ∏_{k>=1}(1-x^k) = sum_{m in Z} (-1)^m x^{g(m)}
// where for m>=0: g(m) = (3*m*m - m) / 2 and for m>=1: the negative terms contribute with exponent g(-m) = (3*m*m + m) / 2.
vector<int> D(n+1, 0);
D[0] = 1; // m = 0 term
// For m positive, add contribution for exponent1 = (3*m^2 - m)/2
for (int m = 1;; m++){
int exp1 = (3LL * m * m - m) / 2;
if(exp1 > n) break;
// (-1)^m: if m is odd, value = -1 mod MOD; if even, value = 1.
int coef = (m % 2 == 1) ? MOD - 1 : 1;
D[exp1] = add_mod(D[exp1], coef);
}
// For m positive, add contribution for exponent2 = (3*m^2 + m)/2
for (int m = 1;; m++){
int exp2 = (3LL * m * m + m) / 2;
if(exp2 > n) break;
int coef = (m % 2 == 1) ? MOD - 1 : 1;
D[exp2] = add_mod(D[exp2], coef);
}
// Next, compute the numerator series:
// N(x) = ∏_{k>=1}(1-x^(2*k))
// By replacing x-> x^2 in Euler’s theorem we have:
// N(x) = sum_{m in Z} (-1)^m x^{m(3*m-1)}
// So for m = 0,1,2,..., and also for negative m.
vector<int> Npoly(n+1, 0);
Npoly[0] = 1; // m = 0 term
// for m>=1, add term with exponent1 = m*(3*m - 1)
for (int m = 1;; m++){
int exp1 = m * (3LL * m - 1);
if(exp1 > n) break;
int coef = (m % 2 == 1) ? MOD - 1 : 1;
Npoly[exp1] = add_mod(Npoly[exp1], coef);
}
// for m>=1, add term with exponent2 = m*(3*m + 1)
for (int m = 1;; m++){
int exp2 = m * (3LL * m + 1);
if(exp2 > n) break;
int coef = (m % 2 == 1) ? MOD - 1 : 1;
Npoly[exp2] = add_mod(Npoly[exp2], coef);
}
// Now, let H(x) = Y(x) be the inverse series of D(x) modulo x^(n+1).
// That is, we want H(x) such that D(x)*H(x) = 1 (mod x^(n+1)).
vector<int> H(n+1, 0);
H[0] = 1; // because D(0)=1.
// Instead of summing over all j from 1 to i, we take advantage that D is sparse.
// Precompute all indices j (>= 1) for which D[j] != 0.
vector<int> nonzero;
for (int j = 1; j <= n; j++) {
if (D[j] != 0) nonzero.push_back(j);
}
// Inversion: for each i >= 1, H[i] = - sum_{j in nonzero, j<=i} D[j]*H[i-j] (mod MOD)
for (int i = 1; i <= n; i++){
long long cur = 0;
for (int j : nonzero) {
if(j > i) break;
cur = (cur + 1LL * D[j] * H[i-j]) % MOD;
}
// H[i] = -cur mod MOD (mod normalization)
H[i] = (cur % MOD + MOD) % MOD; // cur is computed mod MOD; then take negative.
H[i] = (MOD - H[i]) % MOD;
}
// Now the final generating function for our answer is:
// F(x) = N(x) * H(x)
// And the answer is the coefficient of x^n in F(x).
long long ans = 0;
// Since Npoly is sparse, loop over indices where Npoly[i] is nonzero.
for (int i = 0; i <= n; i++){
if(Npoly[i] != 0){
int j = n - i;
if(j >= 0 && j <= n){
ans = (ans + 1LL * Npoly[i] * H[j]) % MOD;
}
}
}
cout << (ans % MOD + MOD) % MOD << "\n";
return 0;
}
NTT 优化过后的代码：
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MOD = 998244353; // NTT-friendly prime; note: original problem stated mod=1e9+7,
// but for NTT we use mod 998244353. For mod 1e9+7 one must use multi-mod approaches.
// Modular exponentiation
int modexp(int a, int b) {
int res = 1;
while(b) {
if(b & 1) res = (ll)res * a % MOD;
a = (ll)a * a % MOD;
b >>= 1;
}
return res;
}
// NTT implementation
// Based on iterative in-place Cooley-Tukey NTT.
// Root for MOD=998244353 is 3.
void ntt(vector<int>& a, bool invert) {
int n = 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 = modexp(3, (MOD - 1) / len);
if (invert) {
wlen = modexp(wlen, MOD - 2);
}
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)((ll)a[i+j+len/2] * 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)((ll)w * wlen % MOD);
}
}
}
if (invert) {
int inv_n = modexp(n, MOD - 2);
for (int & x : a)
x = (int)((ll)x * inv_n % MOD);
}
}
// Multiply two polynomials using NTT
vector<int> poly_mul(const vector<int>& a, const vector<int>& b) {
int n = 1;
while(n < (int)(a.size() + b.size() - 1))
n <<= 1;
vector<int> fa(n, 0), fb(n, 0);
for (size_t i = 0; i < a.size(); i++) fa[i] = a[i];
for (size_t i = 0; i < b.size(); i++) fb[i] = b[i];
ntt(fa, false); ntt(fb, false);
for (int i = 0; i < n; i++) {
fa[i] = (int)((ll)fa[i] * fb[i] % MOD);
}
ntt(fa, true);
// Remove trailing zeros
while(!fa.empty() && fa.back() == 0) fa.pop_back();
return fa;
}
// Polynomial inversion using doubling.
// Given polynomial f with f[0] = 1,
// computes g such that f*g = 1 mod x^n.
vector<int> poly_inv(const vector<int>& f, int n) {
vector<int> g(1, modexp(f[0], MOD - 2)); // g[0] = inverse of f[0]
while((int)g.size() < n) {
int m = g.size() << 1;
vector<int> fcut(f.begin(), f.begin() + min((int)f.size(), m));
vector<int> gg = poly_mul(g, g);
vector<int> prod = poly_mul(gg, fcut);
g.resize(m);
for (int i = 0; i < m; i++) {
int tmp = 0;
if(i < (int)g.size()) tmp = g[i];
// g = g + g - prod (mod x^m)
g[i] = ((2LL * g[i] % MOD - (i < (int)prod.size() ? prod[i] : 0)) % MOD + MOD) % MOD;
}
g.resize(m);
}
g.resize(n);
return g;
}
// Main
// We build polynomials D(x) and N(x) as described:
// D(x) = sum_{k in Z} (-1)^k * x^(g(k)), where for k>0 we add terms for exp = k*(3k-1)/2 and k*(3k+1)/2.
// N(x) = sum_{k in Z} (-1)^k * x^(f(k)), where for k>0 we add terms for exp = k*(3k-1) and k*(3k+1).
int main(){
ios::sync_with_stdio(false);
cin.tie(nullptr);
int n;
cin >> n;
// Build D(x) with degree up to n.
vector<int> D(n+1, 0);
// k = 0 term
D[0] = 1;
// For positive k: add term for exponent1 = (3k^2 - k)/2
for (int k = 1; ; k++){
ll exp1 = (ll)k * (3LL * k - 1) / 2;
if(exp1 > n) break;
int coef = (k & 1) ? MOD - 1 : 1; // (-1)^k mod MOD
D[exp1] = (D[exp1] + coef) % MOD;
}
// For positive k: add term for exponent2 = (3k^2 + k)/2
for (int k = 1; ; k++){
ll exp2 = (ll)k * (3LL * k + 1) / 2;
if(exp2 > n) break;
int coef = (k & 1) ? MOD - 1 : 1;
D[exp2] = (D[exp2] + coef) % MOD;
}
// Build numerator polynomial N(x) with degree up to n.
vector<int> Npoly(n+1, 0);
Npoly[0] = 1;
for (int k = 1; ; k++){
ll exp1 = (ll)k * (3LL * k - 1);
if(exp1 > n) break;
int coef = (k & 1) ? MOD - 1 : 1;
Npoly[exp1] = (Npoly[exp1] + coef) % MOD;
}
for (int k = 1; ; k++){
ll exp2 = (ll)k * (3LL * k + 1);
if(exp2 > n) break;
int coef = (k & 1) ? MOD - 1 : 1;
Npoly[exp2] = (Npoly[exp2] + coef) % MOD;
}
// Invert D(x): compute H(x) such that D(x)*H(x) = 1 mod x^(n+1).
// Since D is given up to degree n, we need H of length n+1.
vector<int> invD = poly_inv(D, n+1);
// The final generating function F(x) = N(x) * H(x).
vector<int> F = poly_mul(Npoly, invD);
int ans = 0;
if(n < (int)F.size())
ans = F[n] % MOD;
else
ans = 0;
cout << ans % MOD << "\n";
return 0;
}
资料推荐
求多项式 exp
专题·多项式基本操作【including 多项式逆元，多项式ln，牛顿迭代，多项式exp，多项式快速幂
浙公网安备 33010602011771号