-
Notifications
You must be signed in to change notification settings - Fork 0
/
DivisorZetaMoebiusTransform.cpp
45 lines (45 loc) · 2.16 KB
/
DivisorZetaMoebiusTransform.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
/*
* @title DivisorZetaMoebiusTransform - 約数のゼータ・メビウス変換 (gcd/lcm畳み込み)
* @docs md/math/DivisorZetaMoebiusTransform.md
*/
class DivisorZetaMoebiusTransform {
template<class T> inline static void zeta_multiple(vector<T>& v, const Eratosthenes& eratosthenes) {
assert(v.size()); int N = v.size();
for(const int& p:eratosthenes.prime) for(int i=(N-1)/p; i; --i) v[i] += v[i*p];
}
template<class T> inline static void mobius_multiple(vector<T>& v, const Eratosthenes& eratosthenes) {
assert(v.size()); int N = v.size();
for(const int& p:eratosthenes.prime) for(int i=1; i*p<N; ++i) v[i] -= v[i*p];
}
template<class T> inline static void zeta_divisor(vector<T>& v, const Eratosthenes& eratosthenes) {
assert(v.size()); int N = v.size();
for(const int& p:eratosthenes.prime) for(int i=1; i*p<N; ++i) {v[i*p] += v[i];}
}
template<class T> inline static void mobius_divisor(vector<T>& v, const Eratosthenes& eratosthenes) {
assert(v.size()); int N = v.size();
for(const int& p:eratosthenes.prime) for(int i=(N-1)/p; i; --i) {v[i*p] -= v[i];}
}
public:
template<class T> inline static vector<T> gcd_convolution(const vector<T>& a,const vector<T>& b, const Eratosthenes& eratosthenes) {
int N = max(a.size(),b.size());
assert(N <= eratosthenes.size());
vector<T> f(N,0),g(N,0);
for(int i=0;i<N;++i) f[i] = a[i];
for(int i=0;i<N;++i) g[i] = b[i];
zeta_multiple(f,eratosthenes);zeta_multiple(g,eratosthenes);
for(int i=0;i<N;++i) f[i] = f[i]*g[i];
mobius_multiple(f,eratosthenes);
return f;
}
template<class T> inline static vector<T> lcm_convolution(const vector<T>& a,const vector<T>& b, const Eratosthenes& eratosthenes) {
int N = max(a.size(),b.size());
assert(N <= eratosthenes.size());
vector<T> f(N,0),g(N,0);
for(int i=0;i<N;++i) f[i] = a[i];
for(int i=0;i<N;++i) g[i] = b[i];
zeta_divisor(f,eratosthenes);zeta_divisor(g,eratosthenes);
for(int i=0;i<N;++i) f[i] = f[i]*g[i];
mobius_divisor(f,eratosthenes);
return f;
}
};