library
This documentation is automatically generated by online-judge-tools/verification-helper
Multivarate Convolution Cyclic
(Convolution/multivariatecyclic.hpp)
- View this file on GitHub
- Last update: 2024-04-26 03:18:17+09:00
- Include:
#include "Convolution/multivariatecyclic.hpp"
Depends on
- Multipoint Evaluation on Geometric Sequence
(FPS/multievalgeom.hpp)
- Fast Division
(Math/fastdiv.hpp)
- Miller-Rabin
(Math/miller.hpp)
- Pollard-Rho
(Math/pollard.hpp)
- Primitive Function
(Math/primitive.hpp)
- Random
(Utility/random.hpp)
Verified with
Code
#pragma once
#include "Math/primitive.hpp"
#include "FPS/multievalgeom.hpp"
template<typename T>vector<T> MultivariateCyclic
(vector<T> f,vector<T> g,vector<int>& a){
int MO=T::get_mod();
int k=a.size(),n=1;
for(auto& x:a)n*=x;
T pr=getPrimitiveRoot(MO),ipr=T(pr).inv();
int offset=1;
rep(x,0,k){
assert((MO-1)%a[x]==0);
T w=pr.pow((MO-1)/a[x]);
rep(i,0,n)if(i%(offset*a[x])<offset){
vector<T> na(a[x]),nb(a[x]);
rep(j,0,a[x]){
na[j]=f[i+offset*j];
nb[j]=g[i+offset*j];
}
na=MultievalGeomSeq(na,T(1),w,a[x]);
nb=MultievalGeomSeq(nb,T(1),w,a[x]);
rep(j,0,a[x]){
f[i+offset*j]=na[j];
g[i+offset*j]=nb[j];
}
}
offset*=a[x];
}
rep(i,0,n)f[i]*=g[i];
offset=1;
rep(x,0,k){
T iw=ipr.pow((MO-1)/a[x]);
rep(i,0,n)if(i%(offset*a[x])<offset){
vector<T> na(a[x]);
rep(j,0,a[x])na[j]=f[i+offset*j];
na=MultievalGeomSeq(na,T(1),iw,a[x]);
rep(j,0,a[x])f[i+offset*j]=na[j];
}
offset*=a[x];
}
T ninv=T(n).inv();
rep(i,0,n)f[i]*=ninv;
return f;
}
/**
* @brief Multivarate Convolution Cyclic
*/#line 2 "Math/fastdiv.hpp"
struct FastDiv{
using u64=uint64_t;
using u128=__uint128_t;
constexpr FastDiv():m(),s(),x(){}
constexpr FastDiv(int _m)
:m(_m),s(__lg(m-1)),x(((u128(1)<<(s+64))+m-1)/m){}
constexpr int get(){return m;}
constexpr friend u64 operator/(u64 n,const FastDiv& d){
return (u128(n)*d.x>>d.s)>>64;
}
constexpr friend int operator%(u64 n,const FastDiv& d){
return n-n/d*d.m;
}
constexpr pair<u64,int> divmod(u64 n)const{
u64 q=n/(*this);
return {q,n-q*m};
}
int m,s; u64 x;
};
/**
* @brief Fast Division
*/
#line 2 "Math/miller.hpp"
struct m64 {
using i64 = int64_t;
using u64 = uint64_t;
using u128 = __uint128_t;
static u64 mod;
static u64 r;
static u64 n2;
static u64 get_r() {
u64 ret = mod;
rep(_,0,5) ret *= 2 - mod * ret;
return ret;
}
static void set_mod(u64 m) {
assert(m < (1LL << 62));
assert((m & 1) == 1);
mod = m;
n2 = -u128(m) % m;
r = get_r();
assert(r * mod == 1);
}
static u64 get_mod() { return mod; }
u64 a;
m64() : a(0) {}
m64(const int64_t &b) : a(reduce((u128(b) + mod) * n2)){};
static u64 reduce(const u128 &b) {
return (b + u128(u64(b) * u64(-r)) * mod) >> 64;
}
u64 get() const {
u64 ret = reduce(a);
return ret >= mod ? ret - mod : ret;
}
m64 &operator*=(const m64 &b) {
a = reduce(u128(a) * b.a);
return *this;
}
m64 operator*(const m64 &b) const { return m64(*this) *= b; }
bool operator==(const m64 &b) const {
return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a);
}
bool operator!=(const m64 &b) const {
return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a);
}
m64 pow(u128 n) const {
m64 ret(1), mul(*this);
while (n > 0) {
if (n & 1) ret *= mul;
mul *= mul;
n >>= 1;
}
return ret;
}
};
typename m64::u64 m64::mod, m64::r, m64::n2;
bool Miller(ll n){
if(n<2 or (n&1)==0)return (n==2);
m64::set_mod(n);
ll d=n-1; while((d&1)==0)d>>=1;
vector<ll> seeds;
if(n<(1<<30))seeds={2, 7, 61};
else seeds={2, 325, 9375, 28178, 450775, 9780504};
for(auto& x:seeds){
if(n<=x)break;
ll t=d;
m64 y=m64(x).pow(t);
while(t!=n-1 and y!=1 and y!=n-1){
y*=y;
t<<=1;
}
if(y!=n-1 and (t&1)==0)return 0;
} return 1;
}
/**
* @brief Miller-Rabin
*/
#line 2 "Utility/random.hpp"
namespace Random {
mt19937_64 randgen(chrono::steady_clock::now().time_since_epoch().count());
using u64 = unsigned long long;
u64 get() {
return randgen();
}
template <typename T> T get(T L) { // [0,L]
return get() % (L + 1);
}
template <typename T> T get(T L, T R) { // [L,R]
return get(R - L) + L;
}
double uniform() {
return double(get(1000000000)) / 1000000000;
}
string str(int n) {
string ret;
rep(i, 0, n) ret += get('a', 'z');
return ret;
}
template <typename Iter> void shuffle(Iter first, Iter last) {
if (first == last)
return;
int len = 1;
for (auto it = first + 1; it != last; it++) {
len++;
int j = get(0, len - 1);
if (j != len - 1)
iter_swap(it, first + j);
}
}
template <typename T> vector<T> select(int n, T L, T R) { // [L,R]
if (n * 2 >= R - L + 1) {
vector<T> ret(R - L + 1);
iota(ALL(ret), L);
shuffle(ALL(ret));
ret.resize(n);
return ret;
} else {
unordered_set<T> used;
vector<T> ret;
while (SZ(used) < n) {
T x = get(L, R);
if (!used.count(x)) {
used.insert(x);
ret.push_back(x);
}
}
return ret;
}
}
void relabel(int n, vector<pair<int, int>> &es) {
shuffle(ALL(es));
vector<int> ord(n);
iota(ALL(ord), 0);
shuffle(ALL(ord));
for (auto &[u, v] : es)
u = ord[u], v = ord[v];
}
template <bool directed, bool simple> vector<pair<int, int>> genGraph(int n) {
vector<pair<int, int>> cand, es;
rep(u, 0, n) rep(v, 0, n) {
if (simple and u == v)
continue;
if (!directed and u > v)
continue;
cand.push_back({u, v});
}
int m = get(SZ(cand));
vector<int> ord;
if (simple)
ord = select(m, 0, SZ(cand) - 1);
else {
rep(_, 0, m) ord.push_back(get(SZ(cand) - 1));
}
for (auto &i : ord)
es.push_back(cand[i]);
relabel(n, es);
return es;
}
vector<pair<int, int>> genTree(int n) {
vector<pair<int, int>> es;
rep(i, 1, n) es.push_back({get(i - 1), i});
relabel(n, es);
return es;
}
}; // namespace Random
/**
* @brief Random
*/
#line 4 "Math/pollard.hpp"
vector<ll> Pollard(ll n) {
if (n <= 1)
return {};
if (Miller(n))
return {n};
if ((n & 1) == 0) {
vector<ll> v = Pollard(n >> 1);
v.push_back(2);
return v;
}
for (ll x = 2, y = 2, d;;) {
ll c = Random::get(2LL, n - 1);
do {
x = (__int128_t(x) * x + c) % n;
y = (__int128_t(y) * y + c) % n;
y = (__int128_t(y) * y + c) % n;
d = __gcd(x - y + n, n);
} while (d == 1);
if (d < n) {
vector<ll> lb = Pollard(d), rb = Pollard(n / d);
lb.insert(lb.end(), ALL(rb));
return lb;
}
}
}
/**
* @brief Pollard-Rho
*/
#line 4 "Math/primitive.hpp"
ll mpow(ll a, ll t, ll m) {
ll res = 1;
FastDiv im(m);
while (t) {
if (t & 1)
res = __int128_t(res) * a % im;
a = __int128_t(a) * a % im;
t >>= 1;
}
return res;
}
ll minv(ll a, ll m) {
ll b = m, u = 1, v = 0;
while (b) {
ll t = a / b;
a -= t * b;
swap(a, b);
u -= t * v;
swap(u, v);
}
u = (u % m + m) % m;
return u;
}
ll getPrimitiveRoot(ll p) {
vector<ll> ps = Pollard(p - 1);
sort(ALL(ps));
rep(x, 1, inf) {
for (auto &q : ps) {
if (mpow(x, (p - 1) / q, p) == 1)
goto fail;
}
return x;
fail:;
}
assert(0);
}
ll extgcd(ll a, ll b, ll &p, ll &q) {
if (b == 0) {
p = 1;
q = 0;
return a;
}
ll d = extgcd(b, a % b, q, p);
q -= a / b * p;
return d;
}
pair<ll, ll> crt(const vector<ll> &vs, const vector<ll> &ms) {
ll V = vs[0], M = ms[0];
rep(i, 1, vs.size()) {
ll p, q, v = vs[i], m = ms[i];
if (M < m)
swap(M, m), swap(V, v);
ll d = extgcd(M, m, p, q);
if ((v - V) % d != 0)
return {0, -1};
ll md = m / d, tmp = (v - V) / d % md * p % md;
V += M * tmp;
M *= md;
}
V = (V % M + M) % M;
return {V, M};
}
ll ModLog(ll a, ll b, ll p) {
ll g = 1;
for (ll t = p; t; t >>= 1)
g = g * a % p;
g = __gcd(g, p);
ll t = 1, c = 0;
for (; t % g; c++) {
if (t == b)
return c;
t = t * a % p;
}
if (b % g)
return -1;
t /= g, b /= g;
ll n = p / g, h = 0, gs = 1;
for (; h * h < n; h++)
gs = gs * a % n;
unordered_map<ll, ll> bs;
for (ll s = 0, e = b; s < h; bs[e] = ++s)
e = e * a % n;
for (ll s = 0, e = t; s < n;) {
e = e * gs % n, s += h;
if (bs.count(e)) {
return c + s - bs[e];
}
}
return -1;
}
ll mod_root(ll k, ll a, ll m) {
if (a == 0)
return k ? 0 : -1;
if (m == 2)
return a & 1;
k %= m - 1;
ll g = gcd(k, m - 1);
if (mpow(a, (m - 1) / g, m) != 1)
return -1;
a = mpow(a, minv(k / g, (m - 1) / g), m);
FastDiv im(m);
auto _subroot = [&](ll p, int e, ll a) -> ll { // x^(p^e)==a(mod m)
ll q = m - 1;
int s = 0;
while (q % p == 0) {
q /= p;
s++;
}
int d = s - e;
ll pe = mpow(p, e, m),
res = mpow(a, ((pe - 1) * minv(q, pe) % pe * q + 1) / pe, m), c = 1;
while (mpow(c, (m - 1) / p, m) == 1)
c++;
c = mpow(c, q, m);
map<ll, ll> mp;
ll v = 1, block = sqrt(d * p) + 1,
bs = mpow(c, mpow(p, s - 1, m - 1) * block % (m - 1), m);
rep(i, 0, block + 1) mp[v] = i, v = v * bs % im;
ll gs = minv(mpow(c, mpow(p, s - 1, m - 1), m), m);
rep(i, 0, d) {
ll err = a * minv(mpow(res, pe, m), m) % im;
ll pos = mpow(err, mpow(p, d - 1 - i, m - 1), m);
rep(j, 0, block + 1) {
if (mp.count(pos)) {
res = res *
mpow(c,
(block * mp[pos] + j) * mpow(p, i, m - 1) %
(m - 1),
m) %
im;
break;
}
pos = pos * gs % im;
}
}
return res;
};
for (ll d = 2; d * d <= g; d++)
if (g % d == 0) {
int sz = 0;
while (g % d == 0) {
g /= d;
sz++;
}
a = _subroot(d, sz, a);
}
if (g > 1)
a = _subroot(g, 1, a);
return a;
}
ull floor_root(ull a, ull k) {
if (a <= 1 or k == 1)
return a;
if (k >= 64)
return 1;
if (k == 2)
return sqrtl(a);
constexpr ull LIM = -1;
if (a == LIM)
a--;
auto mul = [&](ull &x, const ull &y) {
if (x <= LIM / y)
x *= y;
else
x = LIM;
};
auto pw = [&](ull x, ull t) -> ull {
ull y = 1;
while (t) {
if (t & 1)
mul(y, x);
mul(x, x);
t >>= 1;
}
return y;
};
ull ret = (k == 3 ? cbrt(a) - 1 : pow(a, nextafter(1 / double(k), 0)));
while (pw(ret + 1, k) <= a)
ret++;
return ret;
}
/**
* @brief Primitive Function
*/
#line 2 "FPS/multievalgeom.hpp"
template<typename T>vector<T> MultievalGeomSeq(vector<T>& f,T a,T w,int m){
int n=f.size();
vector<T> ret(m);
if(w==0){
T base=1;
rep(i,0,n)ret[0]+=base*f[i],base*=a;
rep(i,1,m)ret[i]=f[0];
return ret;
}
vector<T> tri(n+m-1),itri(n+m-1);
tri[0]=itri[0]=1;
T iw=w.inv(),pww=1,pwiw=1;
for(int i=1;i<n+m-1;i++,pww*=w,pwiw*=iw){
tri[i]=tri[i-1]*pww;
itri[i]=itri[i-1]*pwiw;
}
Poly<T> y(n),v(n+m-1);
T pwa=1;
for(int i=0;i<n;i++,pwa*=a){
y[i]=f[i]*itri[i]*pwa;
}
rep(i,0,n+m-1)v[i]=tri[i];
reverse(ALL(y));
y*=v;
rep(i,0,m)ret[i]=y[n-1+i]*itri[i];
return ret;
}
/**
* @brief Multipoint Evaluation on Geometric Sequence
*/
#line 4 "Convolution/multivariatecyclic.hpp"
template<typename T>vector<T> MultivariateCyclic
(vector<T> f,vector<T> g,vector<int>& a){
int MO=T::get_mod();
int k=a.size(),n=1;
for(auto& x:a)n*=x;
T pr=getPrimitiveRoot(MO),ipr=T(pr).inv();
int offset=1;
rep(x,0,k){
assert((MO-1)%a[x]==0);
T w=pr.pow((MO-1)/a[x]);
rep(i,0,n)if(i%(offset*a[x])<offset){
vector<T> na(a[x]),nb(a[x]);
rep(j,0,a[x]){
na[j]=f[i+offset*j];
nb[j]=g[i+offset*j];
}
na=MultievalGeomSeq(na,T(1),w,a[x]);
nb=MultievalGeomSeq(nb,T(1),w,a[x]);
rep(j,0,a[x]){
f[i+offset*j]=na[j];
g[i+offset*j]=nb[j];
}
}
offset*=a[x];
}
rep(i,0,n)f[i]*=g[i];
offset=1;
rep(x,0,k){
T iw=ipr.pow((MO-1)/a[x]);
rep(i,0,n)if(i%(offset*a[x])<offset){
vector<T> na(a[x]);
rep(j,0,a[x])na[j]=f[i+offset*j];
na=MultievalGeomSeq(na,T(1),iw,a[x]);
rep(j,0,a[x])f[i+offset*j]=na[j];
}
offset*=a[x];
}
T ninv=T(n).inv();
rep(i,0,n)f[i]*=ninv;
return f;
}
/**
* @brief Multivarate Convolution Cyclic
*/