【模板】集合幂级数操作(逐点牛顿迭代法)
正常写法:【模板】子集卷积 - caijianhong - 博客园
由伟大的 Elegia 提出的 Point-wise Newton Iteration(逐点牛顿迭代法):Optimal Algorithm on Polynomial Composite Set Power Series - Codeforces
其它参考资料:feat(math/poly/sps.md): 增加集合幂级数 by hly1204 · Pull Request #5438 · OI-wiki/OI-wiki
以下是集合幂级数的 \(\log\) 的代码实现:注意此写法常数较大!
void fwt(mint vec[], int len, int op) {/*{{{*/if (op == +1) {for (int i = 1; i < len; i <<= 1) {for (int S = i; S < len; S = (S + 1) | i) vec[S] += vec[S ^ i];}} else {for (int i = 1; i < len; i <<= 1) {for (int S = i; S < len; S = (S + 1) | i) vec[S] -= vec[S ^ i];}}
}/*}}}*/
void subset_conv(mint lhs[], mint rhs[], mint res[], int len) {/*{{{*/int n = bitctz(len);vector<vector<mint>> f(n + 1, vector<mint>(1 << n));vector<vector<mint>> g(n + 1, vector<mint>(1 << n));for (int i = 0; i < 1 << n; i++) f[popc(i)][i] = lhs[i];for (int i = 0; i < 1 << n; i++) g[popc(i)][i] = rhs[i];for (int i = 0; i <= n; i++) fwt(f[i].data(), 1 << n, +1);for (int i = 0; i <= n; i++) fwt(g[i].data(), 1 << n, +1);vector<vector<mint>> h(n + 1, vector<mint>(1 << n));for (int i = 0; i <= n; i++) {for (int j = 0; j <= i; j++) {for (int k = 0; k < 1 << n; k++) h[i][k] += f[j][k] * g[i - j][k];}}for (int i = 0; i <= n; i++) fwt(h[i].data(), 1 << n, -1);for (int i = 0; i < 1 << n; i++) res[i] = h[popc(i)][i];
}/*}}}*/
void subset_conv3(mint lhs[], mint rhs[], mint mhs[], mint res[], int len) {/*{{{*/int n = bitctz(len);vector<vector<mint>> f(n + 1, vector<mint>(1 << n));vector<vector<mint>> g(n + 1, vector<mint>(1 << n));vector<vector<mint>> m(n + 1, vector<mint>(1 << n));for (int i = 0; i < 1 << n; i++) f[popc(i)][i] = lhs[i];for (int i = 0; i < 1 << n; i++) g[popc(i)][i] = rhs[i];for (int i = 0; i < 1 << n; i++) m[popc(i)][i] = mhs[i];for (int i = 0; i <= n; i++) fwt(f[i].data(), 1 << n, +1), fwt(g[i].data(), 1 << n, +1);for (int i = 0; i <= n; i++) fwt(m[i].data(), 1 << n, +1);vector<vector<mint>> h(n + 1, vector<mint>(1 << n));for (int i = 0; i <= n; i++) {for (int j = 0; j <= i; j++) {for (int k = 0; k < 1 << n; k++) h[i][k] += f[j][k] * g[i - j][k];}}swap(h, f);for (int i = 0; i <= n; i++) memset(h[i].data(), 0, sizeof(mint) << n);for (int i = 0; i <= n; i++) {for (int j = 0; j <= i; j++) {for (int k = 0; k < 1 << n; k++) h[i][k] += f[j][k] * m[i - j][k];}}for (int i = 0; i <= n; i++) fwt(h[i].data(), 1 << n, -1);for (int i = 0; i < 1 << n; i++) res[i] = h[popc(i)][i];
}/*}}}*/
void sps_inv(mint vec[], mint res[], int len) {int n = bitctz(len);res[0] = 1 / vec[0];for (int i = 0; i < n; i++) {// subset_conv(vec + (1 << i), res, res + (1 << i), 1 << i);// subset_conv(res + (1 << i), res, res + (1 << i), 1 << i);subset_conv3(vec + (1 << i), res, res, res + (1 << i), 1 << i);for (int j = 1 << i; j < 2 << i; j++) res[j] *= -1;}
}
void sps_log(mint vec[], mint res[], int len) {int n = bitctz(len);assert(vec[0] == 1);res[0] = 0;vector<mint> inv(len);sps_inv(vec, inv.data(), len);for (int i = 0; i < n; i++) {subset_conv(vec + (1 << i), inv.data(), res + (1 << i), 1 << i);}
}