diff --git a/math/power_projection_of_set_power_series/checker.cpp b/math/power_projection_of_set_power_series/checker.cpp new file mode 100644 index 000000000..6a66d5330 --- /dev/null +++ b/math/power_projection_of_set_power_series/checker.cpp @@ -0,0 +1,62 @@ +// https://github.com/MikeMirzayanov/testlib/blob/master/checkers/wcmp.cpp + +// The MIT License (MIT) + +// Copyright (c) 2015 Mike Mirzayanov + +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: + +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. + +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. + +#include "testlib.h" + +using namespace std; + +int main(int argc, char * argv[]) +{ + setName("compare sequences of tokens"); + registerTestlibCmd(argc, argv); + + int n = 0; + string j, p; + + while (!ans.seekEof() && !ouf.seekEof()) + { + n++; + + ans.readWordTo(j); + ouf.readWordTo(p); + + if (j != p) + quitf(_wa, "%d%s words differ - expected: '%s', found: '%s'", n, englishEnding(n).c_str(), compress(j).c_str(), compress(p).c_str()); + } + + if (ans.seekEof() && ouf.seekEof()) + { + if (n == 1) + quitf(_ok, "\"%s\"", compress(j).c_str()); + else + quitf(_ok, "%d tokens", n); + } + else + { + if (ans.seekEof()) + quitf(_wa, "Participant output contains extra tokens"); + else + quitf(_wa, "Unexpected EOF in the participants output"); + } +} diff --git a/math/power_projection_of_set_power_series/gen/example_00.in b/math/power_projection_of_set_power_series/gen/example_00.in new file mode 100644 index 000000000..7f712fe8c --- /dev/null +++ b/math/power_projection_of_set_power_series/gen/example_00.in @@ -0,0 +1,3 @@ +2 10 +1 2 3 4 +5 6 7 8 diff --git a/math/power_projection_of_set_power_series/gen/max_random.cpp b/math/power_projection_of_set_power_series/gen/max_random.cpp new file mode 100644 index 000000000..eef486825 --- /dev/null +++ b/math/power_projection_of_set_power_series/gen/max_random.cpp @@ -0,0 +1,35 @@ +#include +#include + +#include "../params.h" +#include "random.h" + +using namespace std; + +void out(Random& gen, int N, int M) { + printf("%d %d\n", N, M); + + for (int i = 0; i < (1 << N); ++i) { + if (i) printf(" "); + int x = gen.uniform(0, MOD - 1); + printf("%d", x); + } + printf("\n"); + + for (int i = 0; i < (1 << N); ++i) { + if (i) printf(" "); + int x = gen.uniform(0, MOD - 1); + printf("%d", x); + } + printf("\n"); +} + +int main(int, char* argv[]) { + long long seed = atoll(argv[1]); + auto gen = Random(seed); + + const int M = M_MAX; + const int N = N_MAX; + + out(gen, N, M); +} diff --git a/math/power_projection_of_set_power_series/gen/random.cpp b/math/power_projection_of_set_power_series/gen/random.cpp new file mode 100644 index 000000000..91fe48765 --- /dev/null +++ b/math/power_projection_of_set_power_series/gen/random.cpp @@ -0,0 +1,35 @@ +#include +#include + +#include "../params.h" +#include "random.h" + +using namespace std; + +void out(Random& gen, int N, int M) { + printf("%d %d\n", N, M); + + for (int i = 0; i < (1 << N); ++i) { + if (i) printf(" "); + int x = gen.uniform(0, MOD - 1); + printf("%d", x); + } + printf("\n"); + + for (int i = 0; i < (1 << N); ++i) { + if (i) printf(" "); + int x = gen.uniform(0, MOD - 1); + printf("%d", x); + } + printf("\n"); +} + +int main(int, char* argv[]) { + long long seed = atoll(argv[1]); + auto gen = Random(seed); + + const int N = gen.uniform(0, N_MAX); + const int M = gen.uniform(0, M_MAX); + + out(gen, N, M); +} diff --git a/math/power_projection_of_set_power_series/gen/small.cpp b/math/power_projection_of_set_power_series/gen/small.cpp new file mode 100644 index 000000000..064ef6ec4 --- /dev/null +++ b/math/power_projection_of_set_power_series/gen/small.cpp @@ -0,0 +1,37 @@ +#include +#include + +#include "../params.h" +#include "random.h" + +using namespace std; + +void out(Random& gen, int N, int M) { + printf("%d %d\n", N, M); + + for (int i = 0; i < (1 << N); ++i) { + if (i) printf(" "); + int x = gen.uniform(0, MOD - 1); + printf("%d", x); + } + printf("\n"); + + for (int i = 0; i < (1 << N); ++i) { + if (i) printf(" "); + int x = gen.uniform(0, MOD - 1); + printf("%d", x); + } + printf("\n"); +} + +int main(int, char* argv[]) { + long long seed = atoll(argv[1]); + auto gen = Random(seed); + + int ns[] = {0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2}; + int ms[] = {0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 4}; + + int N = ns[seed % 12]; + int M = ms[seed % 12]; + out(gen, N, M); +} diff --git a/math/power_projection_of_set_power_series/hash.json b/math/power_projection_of_set_power_series/hash.json new file mode 100644 index 000000000..8413d4237 --- /dev/null +++ b/math/power_projection_of_set_power_series/hash.json @@ -0,0 +1,54 @@ +{ + "example_00.in": "4c9e5af5790c67342ffddaccdaf84c198dfd391b17343f6f89f04e6ce583eb60", + "example_00.out": "c08aeadb5fe3d77c79b80a3cc2faf1e27af915ab9514852b30a098f57864ea2e", + "max_random_00.in": "ccb066924563134b6bed70fefec0f5e89aabae0aca2a044554ef739c085d9cc0", + "max_random_00.out": "5473c83160b1e231a649872597e42de11f518c0cb2d0226f635838521ef490c4", + "max_random_01.in": "ee754c125f69d7f9bbeebd5f6f2777817783d301793c75435ed7b89b2c520cf9", + "max_random_01.out": "de82efb2d8da462b38edce5b54ba7f5e7463fee426132b85880099d7780cb7e8", + "max_random_02.in": "0f509a2bf123f97a9e841c5494ffa2632a8a4d1c3e9216000838a4e26e2b113b", + "max_random_02.out": "fe92965f22d79fa3b99a514ba91097e6aa8f94e788f07809e09e0f7583428e14", + "random_00.in": "a787f024d631265db6427e2997a3e93846ffd47a33a4aaab83dfa6151459bc96", + "random_00.out": "aa091d3aa02e0f15e7957ffd93e98a9461f2e57a9f8b9fbdc2aa103d0ca24ae7", + "random_01.in": "e701371aa017cf4e20d070870d15539fa46b3c665b101eb8784dbcd96cc22fcd", + "random_01.out": "fc9b6808a9db35210e59d0053186afc0b4a79e894c348386a6000e52de4a8222", + "random_02.in": "693851f7a5db6e1aa8576b17db1872c937108efba0bf688b7db364b945ea53e5", + "random_02.out": "12f42aa658856e322cdc9aa4bd73fa2d0a2815dca4c4f05817c6b79f10146f03", + "random_03.in": "50281d4cd0d60b86ed12c50da7696579978b416769041600020fc297412e910c", + "random_03.out": "bbe7f78a4c354dab9efad6258de7333ea1a3bad00c574dd175f581e1bfcb0b9e", + "random_04.in": "e5ef03ecfb4df40b72bb531d4914fcb29f6d01b3d21c6c02d108a117357c7f17", + "random_04.out": "4339c34605f7a3f4a3df4f8b8522e1c1e169ed18e1ea1eee8af057ba5dfa12f7", + "random_05.in": "d5c82e4be32947dd949db87ce92f787dfa5f8574f6c41c0847bea30f97580e08", + "random_05.out": "c45ff481e94a9bab72ad2c83875d17e7a6065700bb193a859d89d3da133d62d4", + "random_06.in": "336bc08bf897008900b40375a1e7c7fa53d59858ef007d8adeda0861bae15c25", + "random_06.out": "f10156b9a664785023b3487942a2487012d085fc3e684cd0dad38b144703cfba", + "random_07.in": "4f56b10e7039781a6313204fe6804a3b9da18090814f7c329fab10ccadab619c", + "random_07.out": "f69ec8cbabd6101abc28199f57bc53b3e95c17fa142c650d934d7db26b8a1ed4", + "random_08.in": "51bfbec7121b5d80ead7ac5009cd3deeae6c0becab91d9ac1d4c87bf24b42a23", + "random_08.out": "735dd5e6b45409714e3447a13c8c1714e50e277c0f48591357b044b419a14e38", + "random_09.in": "ff95191a217733ec20a3b8b055183018579654734ce63e33787eb3f1ca269ca3", + "random_09.out": "c3a895a4128ca672d336794cabd65da2d4f3f0e85c3622935b946b7c4545e251", + "small_00.in": "be684ebbae295e4abd595990db480324e702a12b53af4cc0337212594dc95383", + "small_00.out": "01ba4719c80b6fe911b091a7c05124b64eeece964e09c058ef8f9805daca546b", + "small_01.in": "d92b4601f9054702a2a1e218a4eea9082c1beadfebd47273f8eda2af594a142e", + "small_01.out": "fd809f85df04159b2c8f002c0147443121246de5e93a79271c22399ec7b3b5d0", + "small_02.in": "706c71af6b74b093abd660725c8f4ae4c35b23b8cb023a444f58b3284934029f", + "small_02.out": "9e9758896b8c700fcb5179e212ef903b30e90396b1e462a7b7c8ca44bd9a3c95", + "small_03.in": "80c28bf717e86686248f2ef8c9b4904564717e368cccfbb2cbcf8ad0b02c68c2", + "small_03.out": "01ba4719c80b6fe911b091a7c05124b64eeece964e09c058ef8f9805daca546b", + "small_04.in": "0999cf8d76eaae173b61aa11c5778e809632789d5b65dc6265093748e08d4b4f", + "small_04.out": "2f315045b26261cafaf670e615d32ecab85d0170d779f3cda6b123cadff84a72", + "small_05.in": "3688942481b7ee54de65e8dd75aa580bccfe57b0a63a01863453ceba654bf7cc", + "small_05.out": "5909990d32741c4ec5cabfb63ca4033d49f50ea752f40d55b68fd93cc286dd4f", + "small_06.in": "ab119679055039e73c76b7947928e4213245c5e1640aec2ae2507dd8b830b2b7", + "small_06.out": "ef3e029ec2d736374a7379afa75f24e8bc192f29e27635746b9899ddff9e3210", + "small_07.in": "d73d357e17992e20124a8fce74c2558ea09c2cf6732ed3fa19ea467a36bc1b62", + "small_07.out": "01ba4719c80b6fe911b091a7c05124b64eeece964e09c058ef8f9805daca546b", + "small_08.in": "0bf8ce1104f9aa2c7375ce5e986f9e2fe09a539f10003d813954c6e05b12298a", + "small_08.out": "27a70f6dc3ff41e5b400df71a202c6c56c28972c8aa23483c7ae68913055ef2c", + "small_09.in": "b916cb15bc2649fbf2210789a2efda82997439513bc1e7cb0febc4e1445b1d6d", + "small_09.out": "1c71ff2998f352bb4376b99d761a6ea0ed39642922b37ecb0fd03afe2423d043", + "small_10.in": "2af25dc20b1c52fa92ef14f7d43c46f91d281e20194a062b8e65a73fc1f0ab5e", + "small_10.out": "3de910e53bca95b9183eb67cf0e8193d6290f14b647bc4dc1691e469a7ba32a1", + "small_11.in": "99f8ae47fd88089937d1c99b04de173ab4f85bcd2a18e7f7ec25e1a7a39c66a7", + "small_11.out": "5a88787a0b88fc64f9780c796bea746e01cc8d45fbe9089d5feb71360b4f314c" +} \ No newline at end of file diff --git a/math/power_projection_of_set_power_series/info.toml b/math/power_projection_of_set_power_series/info.toml new file mode 100644 index 000000000..27447a45d --- /dev/null +++ b/math/power_projection_of_set_power_series/info.toml @@ -0,0 +1,26 @@ +title = 'Power Projection of Set Power Series' +timelimit = 10.0 +forum = "https://github.com/yosupo06/library-checker-problems/issues/975" + +[[tests]] + name = "example.in" + number = 1 +[[tests]] + name = "small.cpp" + number = 12 +[[tests]] + name = "random.cpp" + number = 10 +[[tests]] + name = "max_random.cpp" + number = 3 + +[[solutions]] + name = "naive.cpp" + wrong = false + allow_tle = true + +[params] + M_MAX = 100000 + N_MAX = 20 + MOD = 998244353 diff --git a/math/power_projection_of_set_power_series/sol/correct.cpp b/math/power_projection_of_set_power_series/sol/correct.cpp new file mode 100644 index 000000000..9fa28941f --- /dev/null +++ b/math/power_projection_of_set_power_series/sol/correct.cpp @@ -0,0 +1,287 @@ +#include +#include +#include +#include +#include + +using namespace std; + +using ll = long long; +using u32 = unsigned int; +using u64 = unsigned long long; + +template +using vc = vector; + +#define FOR1(a) for (ll _ = 0; _ < ll(a); ++_) +#define FOR2(i, a) for (ll i = 0; i < ll(a); ++i) +#define FOR3(i, a, b) for (ll i = a; i < ll(b); ++i) +#define FOR4(i, a, b, c) for (ll i = a; i < ll(b); i += (c)) +#define FOR1_R(a) for (ll i = (a)-1; i >= ll(0); --i) +#define FOR2_R(i, a) for (ll i = (a)-1; i >= ll(0); --i) +#define FOR3_R(i, a, b) for (ll i = (b)-1; i >= ll(a); --i) +#define overload4(a, b, c, d, e, ...) e +#define overload3(a, b, c, d, ...) d +#define FOR(...) overload4(__VA_ARGS__, FOR4, FOR3, FOR2, FOR1)(__VA_ARGS__) +#define FOR_R(...) overload3(__VA_ARGS__, FOR3_R, FOR2_R, FOR1_R)(__VA_ARGS__) + +#define FOR_subset(t, s) \ + for (ll t = (s); t >= 0; t = (t == 0 ? -1 : (t - 1) & (s))) +#define all(x) x.begin(), x.end() +#define len(x) ll(x.size()) +#define elif else if + +#define eb emplace_back +#define mp make_pair +#define mt make_tuple +#define fi first +#define se second + +int popcnt(int x) { return __builtin_popcount(x); } +// (0, 1, 2, 3, 4) -> (-1, 0, 1, 1, 2) +int topbit(int x) { return (x == 0 ? -1 : 31 - __builtin_clz(x)); } + +template +vc> ranked_zeta(const vc &f) { + int n = topbit(len(f)); + assert(n <= LIM); + assert(len(f) == 1 << n); + vc> Rf(1 << n); + for (int s = 0; s < (1 << n); ++s) Rf[s][popcnt(s)] = f[s]; + for (int i = 0; i < n; ++i) { + int w = 1 << i; + for (int p = 0; p < (1 << n); p += 2 * w) { + for (int s = p; s < p + w; ++s) { + int t = s | 1 << i; + for (int d = 0; d <= n; ++d) Rf[t][d] += Rf[s][d]; + } + } + } + return Rf; +} + +template +vc ranked_mobius(vc> &Rf) { + int n = topbit(len(Rf)); + assert(len(Rf) == 1 << n); + for (int i = 0; i < n; ++i) { + int w = 1 << i; + for (int p = 0; p < (1 << n); p += 2 * w) { + for (int s = p; s < p + w; ++s) { + int t = s | 1 << i; + for (int d = 0; d <= n; ++d) Rf[t][d] -= Rf[s][d]; + } + } + } + vc f(1 << n); + for (int s = 0; s < (1 << n); ++s) f[s] = Rf[s][popcnt(s)]; + return f; +} + +template +vc subset_convolution(const vc &A, const vc &B) { + auto RA = ranked_zeta(A); + auto RB = ranked_zeta(B); + int n = topbit(len(RA)); + FOR(s, len(RA)) { + auto &f = RA[s], &g = RB[s]; + FOR_R(d, n + 1) { + T x = 0; + FOR(i, d + 1) x += f[i] * g[d - i]; + f[d] = x; + } + } + return ranked_mobius(RA); +} + +template +mint inv(int n) { + static const int mod = mint::get_mod(); + static vector dat = {0, 1}; + assert(0 <= n); + if (n >= mod) n %= mod; + while (len(dat) <= n) { + int k = len(dat); + int q = (mod + k - 1) / k; + dat.eb(dat[k * q - mod] * mint::raw(q)); + } + return dat[n]; +} + +template +mint fact(int n) { + static const int mod = mint::get_mod(); + assert(0 <= n && n < mod); + static vector dat = {1, 1}; + while (len(dat) <= n) dat.eb(dat[len(dat) - 1] * mint::raw(len(dat))); + return dat[n]; +} + +template +mint fact_inv(int n) { + static vector dat = {1, 1}; + if (n < 0) return mint(0); + while (len(dat) <= n) dat.eb(dat[len(dat) - 1] * inv(len(dat))); + return dat[n]; +} + +template +struct modint { + static constexpr u32 umod = u32(mod); + static_assert(umod < u32(1) << 31); + u32 val; + + static modint raw(u32 v) { + modint x; + x.val = v; + return x; + } + constexpr modint() : val(0) {} + constexpr modint(u32 x) : val(x % umod) {} + constexpr modint(u64 x) : val(x % umod) {} + constexpr modint(int x) : val((x %= mod) < 0 ? x + mod : x){}; + constexpr modint(ll x) : val((x %= mod) < 0 ? x + mod : x){}; + bool operator<(const modint &other) const { return val < other.val; } + modint &operator+=(const modint &p) { + if ((val += p.val) >= umod) val -= umod; + return *this; + } + modint &operator-=(const modint &p) { + if ((val += umod - p.val) >= umod) val -= umod; + return *this; + } + modint &operator*=(const modint &p) { + val = u64(val) * p.val % umod; + return *this; + } + modint &operator/=(const modint &p) { + *this *= p.inverse(); + return *this; + } + modint operator-() const { return modint::raw(val ? mod - val : u32(0)); } + modint operator+(const modint &p) const { return modint(*this) += p; } + modint operator-(const modint &p) const { return modint(*this) -= p; } + modint operator*(const modint &p) const { return modint(*this) *= p; } + modint operator/(const modint &p) const { return modint(*this) /= p; } + bool operator==(const modint &p) const { return val == p.val; } + bool operator!=(const modint &p) const { return val != p.val; } + modint inverse() const { + int a = val, b = mod, u = 1, v = 0, t; + while (b > 0) { + t = a / b; + swap(a -= t * b, b), swap(u -= t * v, v); + } + return modint(u); + } + modint pow(ll n) const { + assert(n >= 0); + modint ret(1), mul(val); + while (n > 0) { + if (n & 1) ret *= mul; + mul *= mul; + n >>= 1; + } + return ret; + } + static constexpr int get_mod() { return mod; } +}; + +template +vc convolution_naive(const vc &a, const vc &b) { + int n = int(a.size()), m = int(b.size()); + if (n > m) return convolution_naive(b, a); + if (n == 0) return {}; + vector ans(n + m - 1); + FOR(i, n) FOR(j, m) ans[i + j] += a[i] * b[j]; + return ans; +} + +// for fixed sps s, consider linear map F:a->b = subset-conv(a,s) +// given x, calculate transpose(F)(x) +template +vc transposed_subset_convolution(vc s, vc x) { + /* + sum_{j}x_jb_j = sum_{i subset j}x_ja_is_{j-i} = sum_{i}y_ia_i. + y_i = sum_{j supset i}x_js_{j-i} + (rev y)_i = sum_{j subset i}(rev x)_js_{i-j} + y = rev(conv(rev x), s) + */ + reverse(all(x)); + x = subset_convolution(x, s); + reverse(all(x)); + return x; +} + +// for fixed sps s s.t. s[0] == 0. +// consider linear map F:f->t=f(s) for egf f. +// given x, calcuate transpose(F)(x) +// equivalent: calculate sum_i x_i(s^k/k!)_i for k=0,1,...,N +template +vc transposed_sps_composition_egf(vc &s, vc x) { + const int N = topbit(len(s)); + assert(len(s) == (1 << N) && len(x) == (1 << N) && s[0] == mint(0)); + vc y(N + 1); + y[0] = x[0]; + auto &dp = x; + FOR(i, N) { + vc newdp(1 << (N - 1 - i)); + FOR(j, N - i) { + vc a = {s.begin() + (1 << j), s.begin() + (2 << j)}; + vc b = {dp.begin() + (1 << j), dp.begin() + (2 << j)}; + b = transposed_subset_convolution(a, b); + FOR(k, len(b)) newdp[k] += b[k]; + } + swap(dp, newdp); + y[1 + i] = dp[0]; + } + return y; +} + +// for fixed sps s s.t. s[0] == 0. +// consider linear map F:f->t=f(s) for polynomial f. +// given x, calcuate transpose(F)(x) +// equivalent: calculate sum_i x_i(s^k/k!)_i for k=0,1,...,M-1 +template +vc transposed_sps_composition_poly(vc s, vc x, int M) { + const int N = topbit(len(s)); + assert(len(s) == (1 << N) && len(x) == (1 << N)); + mint c = s[0]; + s[0] -= c; + x = transposed_sps_composition_egf(s, x); + vc g(M); + mint pow = 1; + FOR(i, M) { g[i] = pow * fact_inv(i), pow *= c; } + x = convolution_naive(x, g); + x.resize(M); + FOR(i, M) x[i] *= fact(i); + return x; +} + +using mint = modint<998244353>; + +void solve() { + int N, M; + scanf("%d %d", &N, &M); + vc s(1 << N), w(1 << N); + FOR(i, 1 << N) { + int x; + scanf("%d", &x); + s[i] = x; + } + FOR(i, 1 << N) { + int x; + scanf("%d", &x); + w[i] = x; + } + vc ANS = transposed_sps_composition_poly(s, w, M); + FOR(i, M) { + if (i) printf(" "); + printf("%d", ANS[i].val); + } + printf("\n"); +} + +signed main() { + solve(); + return 0; +} diff --git a/math/power_projection_of_set_power_series/sol/naive.cpp b/math/power_projection_of_set_power_series/sol/naive.cpp new file mode 100644 index 000000000..2649f5e94 --- /dev/null +++ b/math/power_projection_of_set_power_series/sol/naive.cpp @@ -0,0 +1,188 @@ +#include +#include +#include +#include +#include + +using namespace std; + +using ll = long long; +using u32 = unsigned int; +using u64 = unsigned long long; + +template +using vc = vector; + +#define FOR1(a) for (ll _ = 0; _ < ll(a); ++_) +#define FOR2(i, a) for (ll i = 0; i < ll(a); ++i) +#define FOR3(i, a, b) for (ll i = a; i < ll(b); ++i) +#define FOR1_R(a) for (ll i = (a)-1; i >= ll(0); --i) +#define FOR2_R(i, a) for (ll i = (a)-1; i >= ll(0); --i) +#define FOR3_R(i, a, b) for (ll i = (b)-1; i >= ll(a); --i) +#define overload4(a, b, c, d, e, ...) e +#define overload3(a, b, c, d, ...) d +#define FOR(...) overload4(__VA_ARGS__, FOR4, FOR3, FOR2, FOR1)(__VA_ARGS__) +#define FOR_R(...) overload3(__VA_ARGS__, FOR3_R, FOR2_R, FOR1_R)(__VA_ARGS__) + +#define all(x) x.begin(), x.end() +#define len(x) ll(x.size()) +#define elif else if + +#define eb emplace_back + +int popcnt(int x) { return __builtin_popcount(x); } +// (0, 1, 2, 3, 4) -> (-1, 0, 1, 1, 2) +int topbit(int x) { return (x == 0 ? -1 : 31 - __builtin_clz(x)); } + +template +vc> ranked_zeta(const vc &f) { + int n = topbit(len(f)); + assert(n <= LIM); + assert(len(f) == 1 << n); + vc> Rf(1 << n); + for (int s = 0; s < (1 << n); ++s) Rf[s][popcnt(s)] = f[s]; + for (int i = 0; i < n; ++i) { + int w = 1 << i; + for (int p = 0; p < (1 << n); p += 2 * w) { + for (int s = p; s < p + w; ++s) { + int t = s | 1 << i; + for (int d = 0; d <= n; ++d) Rf[t][d] += Rf[s][d]; + } + } + } + return Rf; +} + +template +vc ranked_mobius(vc> &Rf) { + int n = topbit(len(Rf)); + assert(len(Rf) == 1 << n); + for (int i = 0; i < n; ++i) { + int w = 1 << i; + for (int p = 0; p < (1 << n); p += 2 * w) { + for (int s = p; s < p + w; ++s) { + int t = s | 1 << i; + for (int d = 0; d <= n; ++d) Rf[t][d] -= Rf[s][d]; + } + } + } + vc f(1 << n); + for (int s = 0; s < (1 << n); ++s) f[s] = Rf[s][popcnt(s)]; + return f; +} + +template +vc subset_convolution(const vc &A, const vc &B) { + auto RA = ranked_zeta(A); + auto RB = ranked_zeta(B); + int n = topbit(len(RA)); + FOR(s, len(RA)) { + auto &f = RA[s], &g = RB[s]; + FOR_R(d, n + 1) { + T x = 0; + FOR(i, d + 1) x += f[i] * g[d - i]; + f[d] = x; + } + } + return ranked_mobius(RA); +} + +template +struct modint { + static constexpr u32 umod = u32(mod); + static_assert(umod < u32(1) << 31); + u32 val; + + static modint raw(u32 v) { + modint x; + x.val = v; + return x; + } + constexpr modint() : val(0) {} + constexpr modint(u32 x) : val(x % umod) {} + constexpr modint(u64 x) : val(x % umod) {} + constexpr modint(int x) : val((x %= mod) < 0 ? x + mod : x){}; + constexpr modint(ll x) : val((x %= mod) < 0 ? x + mod : x){}; + bool operator<(const modint &other) const { return val < other.val; } + modint &operator+=(const modint &p) { + if ((val += p.val) >= umod) val -= umod; + return *this; + } + modint &operator-=(const modint &p) { + if ((val += umod - p.val) >= umod) val -= umod; + return *this; + } + modint &operator*=(const modint &p) { + val = u64(val) * p.val % umod; + return *this; + } + modint &operator/=(const modint &p) { + *this *= p.inverse(); + return *this; + } + modint operator-() const { return modint::raw(val ? mod - val : u32(0)); } + modint operator+(const modint &p) const { return modint(*this) += p; } + modint operator-(const modint &p) const { return modint(*this) -= p; } + modint operator*(const modint &p) const { return modint(*this) *= p; } + modint operator/(const modint &p) const { return modint(*this) /= p; } + bool operator==(const modint &p) const { return val == p.val; } + bool operator!=(const modint &p) const { return val != p.val; } + modint inverse() const { + int a = val, b = mod, u = 1, v = 0, t; + while (b > 0) { + t = a / b; + swap(a -= t * b, b), swap(u -= t * v, v); + } + return modint(u); + } + modint pow(ll n) const { + assert(n >= 0); + modint ret(1), mul(val); + while (n > 0) { + if (n & 1) ret *= mul; + mul *= mul; + n >>= 1; + } + return ret; + } + static constexpr int get_mod() { return mod; } +}; + +using mint = modint<998244353>; + +void solve() { + int N, M; + scanf("%d %d", &N, &M); + + vc a(1 << N), c(1 << N); + FOR(i, 1 << N) { + int x; + scanf("%d", &x); + a[i] = x; + } + FOR(i, 1 << N) { + int x; + scanf("%d", &x); + c[i] = x; + } + + vc ANS(M); + vc f(1 << N); + f[0] = 1; + + FOR(i, M) { + FOR(j, 1 << N) ANS[i] += f[j] * c[j]; + f = subset_convolution(f, a); + } + + FOR(i, M) { + if (i) printf(" "); + printf("%d", ANS[i].val); + } + printf("\n"); +} + +signed main() { + solve(); + return 0; +} diff --git a/math/power_projection_of_set_power_series/task.md b/math/power_projection_of_set_power_series/task.md new file mode 100644 index 000000000..681797f08 --- /dev/null +++ b/math/power_projection_of_set_power_series/task.md @@ -0,0 +1,39 @@ +## @{keyword.statement} + +@{lang.en} + +Given an integer $M$, an $N$-variable polynomial $\displaystyle s(x _ 0, x _ 1, \dots, x _ {N - 1}) = \sum _ {i = 0} ^ {2 ^ N - 1} a _ i \prod _ {k = 0} ^ {N - 1} x _ k ^ {i _ k} \in \mathbb{F} _ {@{param.MOD}} \lbrack x _ 0, x _ 1, \dots, x _ {N - 1} \rbrack$ and a sequence $(w _ 0, \ldots, w _ {2 ^ N - 1})$ of elements of $\mathbb{F} _ {@{param.MOD}}$. Here, $i_k$ represents the number in the $2 ^ k$'s place when $i$ is expressed in binary. + +For each $m=0,1,\ldots,M-1$, compute $\sum _ {i=0} ^ {2 ^ N - 1} w_ib_i$, where $b_i$ is defined by $s(x _ 0, x _ 1, \dots, x _ {N - 1})^m \equiv \sum _ {i = 0} ^ {2 ^ N - 1} b _ i \prod _ {k = 0} ^ {N - 1} x _ k ^ {i _ k} \pmod{(x _ 0 ^ 2, x _ 1 ^ 2, \dots, x _ {N - 1} ^ 2)}$. + +@{lang.ja} + +整数 $M$ および,$N$ 変数多項式 $\displaystyle s(x _ 0, x _ 1, \dots, x _ {N - 1}) = \sum _ {i = 0} ^ {2 ^ N - 1} a _ i \prod _ {k = 0} ^ {N - 1} x _ k ^ {i _ k} \in \mathbb{F} _ {@{param.MOD}} \lbrack x _ 0, x _ 1, \dots, x _ {N - 1} \rbrack$ および $\mathbb{F} _ {@{param.MOD}}$ の元の列 $(w _ 0, \ldots, w _ {2 ^ N - 1})$ が与えられます. ただし, $i _ k$ は $i$ を $2$ 進法で表した時の $2 ^ k$ の位の値を表します. + +各 $m=0,1,\ldots,M-1$ に対して $s(x _ 0, x _ 1, \dots, x _ {N - 1})^m \equiv \sum _ {i = 0} ^ {2 ^ N - 1} b _ i \prod _ {k = 0} ^ {N - 1} x _ k ^ {i _ k} \pmod{(x _ 0 ^ 2, x _ 1 ^ 2, \dots, x _ {N - 1} ^ 2)}$ により $b _ i$ を定めたときの $\sum _ {i=0} ^ {2 ^ N - 1} w _ ib _ i$ を求めてください. + +@{lang.end} + +## @{keyword.constraints} + +- $0 \leq N \leq @{param.N_MAX}$ +- $0 \leq M \leq @{param.M_MAX}$ +- $0 \leq a_i, w_i \lt @{param.MOD}$ + +## @{keyword.input} + +``` +$N$ $M$ +$a _ 0$ $a _ 1$ $\cdots$ $a _ {2^N-1}$ +$w _ 0$ $w _ 1$ $\cdots$ $w _ {2^N-1}$ +``` + +## @{keyword.output} + +``` +$\mathrm{ans}_0$ $\mathrm{ans}_1$ $\cdots$ $\mathrm{ans}_{M-1}$ +``` + +## @{keyword.sample} + +@{example.example_00} diff --git a/math/power_projection_of_set_power_series/verifier.cpp b/math/power_projection_of_set_power_series/verifier.cpp new file mode 100644 index 000000000..e9f321933 --- /dev/null +++ b/math/power_projection_of_set_power_series/verifier.cpp @@ -0,0 +1,28 @@ +#include "testlib.h" + +#include "params.h" + +int main() { + registerValidation(); + + const int N = inf.readInt(0, N_MAX); + inf.readChar(' '); + inf.readInt(0, M_MAX); + inf.readChar('\n'); + + for (int i = 0; i < 1 << N; i++) { + if (i) inf.readChar(' '); + inf.readInt(0, MOD - 1); + } + inf.readChar('\n'); + + for (int i = 0; i < 1 << N; i++) { + if (i) inf.readChar(' '); + inf.readInt(0, MOD - 1); + } + inf.readChar('\n'); + + inf.readEof(); + + return 0; +}