洛谷P3312 [SDOI2014]数表
标签
- 莫比乌斯反演
- 线段树
- 约数和筛
前言
- 我的csdn和博客园是同步的,欢迎来访danzh-博客园~
欢迎关注~ - 用数据结构维护的数论题…有点难,不过做完后收获很大~
简明题意
- 有一个 n ∗ m n *m n∗m的数表,每个点 ( i , j ) (i,j) (i,j)的值是能同时整除 i , j i,j i,j的自然数之和。再给定 a a a,需要你求表中所有值 > a >a >a的项之和。
思路
- 首先用数学公式表示出来:
∑ i = 1 n ∑ j = 1 m ( ∑ d ∣ i 且 d ∣ j d ) ∗ [ ∑ d ∣ i 且 d ∣ j d < = a ] \sum_{i=1}^n\sum_{j=1}^m(\sum_{d|i且d|j}d)*[\sum_{d|i且d|j}d<=a] i=1∑nj=1∑m(d∣i且d∣j∑d)∗[d∣i且d∣j∑d<=a] - 我们知道关于
g
c
d
gcd
gcd有一个很常用的性质:
d
∣
g
c
d
(
i
,
j
)
  
⟺
  
d
∣
i
且
d
∣
j
d|gcd(i,j) \iff d|i且d|j
d∣gcd(i,j)⟺d∣i且d∣j(这个太常用了,不会的同学一定要自行搞清楚)。然后我们替换掉,就成了:
∑ i = 1 n ∑ j = 1 m ( ∑ d ∣ g c d ( i , j ) d ) ∗ [ ∑ d ∣ g c d ( i , j ) d < = a ] \sum_{i=1}^n\sum_{j=1}^m(\sum_{d|gcd(i,j)}d)*[\sum_{d|gcd(i,j)}d<=a] i=1∑nj=1∑m(d∣gcd(i,j)∑d)∗[d∣gcd(i,j)∑d<=a] - 然后就是改为枚举
g
c
d
(
i
,
j
)
gcd(i,j)
gcd(i,j),显然上限是
n
n
n
∑ x = 1 n ∑ i = 1 n ∑ j = 1 m ∑ d ∣ x d ∗ [ ∑ d ∣ x d < = a ] [ g c d ( i , j ) = = x ] \sum_{x=1}^n\sum_{i=1}^n\sum_{j=1}^m\sum_{d|x}d*[\sum_{d|x}d<=a][gcd(i,j)==x] x=1∑ni=1∑nj=1∑md∣x∑d∗[d∣x∑d<=a][gcd(i,j)==x]
式子反而变复杂了?不要担心,通常改为枚举会使得式子最前面多一个 ∑ \sum ∑,被加式中多一个条件。但后面的步骤就会把他们都简化
- 接下来移项,成为:
∑ x = 1 n ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = = x ] ∑ d ∣ x d ∗ [ ∑ d ∣ x d < = a ] \sum_{x=1}^n\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)==x]\sum_{d|x}d*[\sum_{d|x}d<=a] x=1∑ni=1∑nj=1∑m[gcd(i,j)==x]d∣x∑d∗[d∣x∑d<=a] - 然后由莫比乌斯反演我们又有这样一个式子:
∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = = x ] = ∑ d = 1 [ n x ] ( μ ( d ) ∗ [ n d x ] ∗ [ m d x ] ) \sum_{i=1}^{n}\sum_{j=1}^m[gcd(i,j)==x]=\sum_{d=1}^{[\frac nx]}\left( \mu(d)*[\frac {n}{dx}]*[\frac {m}{dx}]\right) i=1∑nj=1∑m[gcd(i,j)==x]=d=1∑[xn](μ(d)∗[dxn]∗[dxm])这个式子很常用。不清楚的同学参见我的另一篇博客戳这里
- 把这个式子代进去,原式就是:
∑ x = 1 n ∑ d = 1 [ n x ] ( μ ( d ) ∗ [ n d x ] ∗ [ m d x ] ) ∑ d ∣ x d ∗ [ ∑ d ∣ x d < = a ] \sum_{x=1}^n\sum_{d=1}^{[\frac nx]}\left( \mu(d)*[\frac {n}{dx}]*[\frac {m}{dx}]\right)\sum_{d|x}d*[\sum_{d|x}d<=a] x=1∑nd=1∑[xn](μ(d)∗[dxn]∗[dxm])d∣x∑d∗[d∣x∑d<=a] - 然后观察这个式子:
∑ d ∣ x d \sum_{d|x}d d∣x∑d - 数论比较好的同学可以立马发现这就是
σ
\sigma
σ函数,方便起见,我们就用
σ
\sigma
σ替换
∑ x = 1 n ∑ d = 1 [ n x ] ( μ ( d ) ∗ [ n d x ] ∗ [ m d x ] ) σ ( x ) ∗ [ σ ( x ) < = a ] \sum_{x=1}^n\sum_{d=1}^{[\frac nx]}\left( \mu(d)*[\frac {n}{dx}]*[\frac {m}{dx}]\right)\sigma(x)*[\sigma(x)<=a] x=1∑nd=1∑[xn](μ(d)∗[dxn]∗[dxm])σ(x)∗[σ(x)<=a] - 然后式子比较长,我们令
f
(
x
)
=
σ
(
x
)
∗
[
σ
(
x
)
<
=
a
]
f(x)=\sigma(x)*[\sigma(x)<=a]
f(x)=σ(x)∗[σ(x)<=a]:
∑ x = 1 n f ( x ) ∑ d = 1 [ n x ] ( μ ( d ) ∗ [ n d x ] ∗ [ m d x ] ) \sum_{x=1}^nf(x)\sum_{d=1}^{[\frac nx]}\left( \mu(d)*[\frac {n}{dx}]*[\frac {m}{dx}]\right) x=1∑nf(x)d=1∑[xn](μ(d)∗[dxn]∗[dxm]) - 到了这里,我们可以用线性筛预处理出 f , μ f,\mu f,μ函数(不会约数个数和筛的同学请看看别的博客弄明白),再前缀和一下,后面的分以下块,整体复杂度就会是 O ( n + T n n ) O(n+Tn\sqrt{n}) O(n+Tnn),可以拿到60分。我们进一步优化。
- 这里我们令
k
=
d
x
k=dx
k=dx,然后改为枚举
k
k
k:
∑ k = 1 n [ n k ] [ m k ] ∑ d ∣ k f ( d ) μ ( [ k d ] ) \sum_{k=1}^n[\frac {n}{k}][\frac {m}{k}]\sum_{d|k}f(d)\mu([\frac kd]) k=1∑n[kn][km]d∣k∑f(d)μ([dk])
这一步一定很多同学会产生疑惑。我在这里解释一下,我们有这样一个结论:
∑ i = 1 n ∑ j = 1 [ n i ] f ( i j ) ∗ g ( i ) ∗ h ( j ) = ∑ k = 1 n f ( k ) ∗ ∑ d ∣ k g ( [ k d ] ) h ( d ) \sum_{i=1}^n\sum_{j=1}^{[\frac ni]}f(ij)*g(i)*h(j)=\sum_{k=1}^nf(k)*\sum_{d|k}g([\frac kd])h(d) i=1∑nj=1∑[in]f(ij)∗g(i)∗h(j)=k=1∑nf(k)∗d∣k∑g([dk])h(d)
只要是二重枚举式,第一重枚举随意,且第二重枚举枚举的是第一重的上限除以第一重的枚举项,就可以将被加式写成关于 i , j , i ∗ j i,j,i*j i,j,i∗j的三个函数,所以可以令 k = i j k=ij k=ij。这个式子要牢记。
- 现在令
g
(
k
)
=
∑
d
∣
k
f
(
d
)
μ
(
[
k
d
]
)
g(k) =\sum\limits_{d|k}f(d)\mu([\frac kd])
g(k)=d∣k∑f(d)μ([dk]),原式就是:
∑ k = 1 n [ n k ] [ m k ] g ( k ) \sum_{k=1}^n[\frac {n}{k}][\frac {m}{k}]g(k) k=1∑n[kn][km]g(k) - 这个式子,我们十分想对它进行整除分块,如果整除分块,显然应该预处理出
g
(
k
)
g(k)
g(k)的前缀和。然而,每组样例
a
a
a不同,会影响
g
(
k
)
g(k)
g(k)的取值,因此无法对
g
(
k
)
g(k)
g(k)前缀和。怎么办呢?所以现在的关键就是如何处理对不同的
a
a
a计算。现在先把
f
(
d
)
f(d)
f(d)代回去:
g ( k ) = ∑ d ∣ k σ ( d ) ∗ [ σ ( d ) < = a ] μ ( [ k d ] g(k) =\sum\limits_{d|k}\sigma(d)*[\sigma(d)<=a]\mu([\frac kd] g(k)=d∣k∑σ(d)∗[σ(d)<=a]μ([dk] - 我们先不考虑 g ( k ) g(k) g(k)的前缀和,而是考虑单个的 g ( k ) g(k) g(k)。实际上, g ( k ) g(k) g(k)的取值并不是只和 k k k有关,还与 a a a有关,严谨一点, g ( k ) g(k) g(k)应该是一个二元函数 g ( k , a ) g(k,a) g(k,a)。但是实际上 k k k并不重要,因为每一次分块都要求 g ( k ) g(k) g(k)的前缀和,所以实际上对于给定的 a a a,必须求出所有的 g ( k ) g(k) g(k),而不是只用求出单独一项 g ( k ) g(k) g(k)。所以,关键的是 a a a的取值。这是我们换一个角度思考这个问题
- 对于给定的 k k k,显然 a a a越大,就会有更多的 σ \sigma σ贡献到 g ( k ) g(k) g(k)。如果 a a a从 3 3 3变化到 10 10 10,对于同一个 k k k,会有 3 < σ ( d ) < = 10 的 σ ( d ) ∗ μ ( k d ) ( d ∣ k ) 3<\sigma(d)<=10的\sigma(d)*\mu(\frac kd)(d|k) 3<σ(d)<=10的σ(d)∗μ(dk)(d∣k)被加入到了原先的 g ( k ) g(k) g(k)。(多读几遍黑色的部分),是不是注意到了 a a a越大, g ( k ) g(k) g(k)越大,而且更大的 a a a得到的 g ( k ) g(k) g(k)可以由更小的 a a a推过来。
- 所以这里,我们先离线一下,把
a
a
a从小到大排序。然后对
g
(
k
)
g(k)
g(k)建一棵线段树,每一次
a
a
a增大,更新一下线段树。具体更新的方法是:把位于
(
上
一
次
的
a
,
当
前
a
]
(上一次的a,当前a]
(上一次的a,当前a]范围内的
σ
(
d
)
\sigma(d)
σ(d)找出来,符合要求的每一个
σ
(
d
)
\sigma(d)
σ(d)都会对
d
∣
k
d|k
d∣k的
g
(
k
)
g(k)
g(k)产生贡献。就对
d
d
d枚举它的所有倍数
k
k
k,让其值加上
σ
(
d
)
∗
μ
(
k
d
)
\sigma(d)*\mu(\frac kd)
σ(d)∗μ(dk)。每一次修改操作进行完,
g
(
k
)
g(k)
g(k)就是正确的,所以回想刚刚要求的式子:
∑ k = 1 n [ n k ] [ m k ] g ( k ) \sum_{k=1}^n[\frac {n}{k}][\frac {m}{k}]g(k) k=1∑n[kn][km]g(k) - 我们直接数论分块,每一块求出区间 [ l , r ] [l,r] [l,r]的 g ( k ) g(k) g(k)之和乘以 [ n l ] [ m k ] [\frac nl][\frac mk] [ln][km],全部累加起来就是这一次的答案。
注意事项
- 无
总结
- 对 2 31 2^{31} 231取模等价于最终的答案& 2 31 − 1 2^{31}-1 231−1
- 前缀和如果动态变化,可以尝试用线段树维护。
-
∑
i
=
1
n
∑
j
=
1
[
n
i
]
f
(
i
j
)
∗
g
(
i
)
∗
h
(
j
)
=
∑
k
=
1
n
f
(
k
)
∗
∑
d
∣
k
g
(
[
k
d
]
)
h
(
d
)
\sum_{i=1}^n\sum_{j=1}^{[\frac ni]}f(ij)*g(i)*h(j)=\sum_{k=1}^nf(k)*\sum_{d|k}g([\frac kd])h(d)
i=1∑nj=1∑[in]f(ij)∗g(i)∗h(j)=k=1∑nf(k)∗d∣k∑g([dk])h(d)
这个式子很常用! - 注意约数和筛的写法!
AC代码
// luogu-judger-enable-o2
#include<cstdio>
#include<algorithm>
using namespace std;
const int maxn = 1e5 + 10;
struct Sig
{
int sigma, id;
bool operator < (const Sig& a)const
{
return sigma < a.sigma;
}
};
Sig sigma[maxn];
bool no_prime[maxn];
int prime[maxn], s[maxn], ssum[maxn], mu[maxn], pre_mu[maxn];
int shai(int n)
{
int cnt = 0;
mu[1] = 1, s[1] = 1;
for (int i = 2; i <= n; i++)
{
if (!no_prime[i])
prime[++cnt] = i, mu[i] = -1, s[i] = i + 1, ssum[i] = i + 1;
for (int j = 1; j <= cnt && prime[j] * i <= n; j++)
{
no_prime[prime[j] * i] = 1;
mu[prime[j] * i] = i % prime[j] == 0 ? 0 : -mu[i];
s[prime[j] * i] = i % prime[j] == 0 ? s[i] / ssum[i] + s[i] * prime[j] : s[i] * (1 + prime[j]);
ssum[prime[j] * i] = i % prime[j] == 0 ? ssum[i] * prime[j] + 1 : prime[j] + 1;
if (i % prime[j] == 0) break;
}
}
for (int i = 1; i <= n; i++)
pre_mu[i] = pre_mu[i - 1] + mu[i], sigma[i].sigma = s[i], sigma[i].id = i;
sort(sigma + 1, sigma + 1 + n);
return cnt;
}
int n, m, a;
int f(int x)
{
if (s[x] <= a) return s[x];
return 0;
}
struct Query
{
int n, m, a, id;
bool operator < (const Query& b)const
{
return a < b.a;
}
};
Query query[maxn];
struct Node
{
int l, r, sum;
};
Node tree[maxn * 4];
void build(int o, int l, int r)
{
tree[o].l = l, tree[o].r = r;
if (l == r)
return;
int mid = (l + r) / 2;
build(o * 2, l, mid);
build(o * 2 + 1, mid + 1, r);
}
void change(int o, int x, int k)
{
if (tree[o].l == tree[o].r)
{
tree[o].sum += k;
return;
}
int mid = (tree[o].l + tree[o].r) / 2;
if (x <= mid)
change(o * 2, x, k);
else
change(o * 2 + 1, x, k);
tree[o].sum = (tree[o * 2].sum + tree[o * 2 + 1].sum);
}
int ask(int o, int l, int r)
{
if (tree[o].l == l && tree[o].r == r)
return tree[o].sum;
int mid = (tree[o].l + tree[o].r) / 2;
if (r <= mid)
return ask(o * 2, l, r);
else if (l > mid)
return ask(o * 2 + 1, l, r);
return ask(o * 2, l, mid) + ask(o * 2 + 1, mid + 1, r);
}
int cal(int n, int m)
{
int l = 1, r, ans = 0;
while (l <= n)
{
r = min(n / (n / l), m / (m / l));
ans += (n / l) * (m / l) * ask(1, l, r);
l = r + 1;
}
return ans;
}
int ans[maxn];
void solve()
{
shai(maxn - 10);
build(1, 1, maxn - 10);
int t;
scanf("%d", &t);
for (int i = 1; i <= t; i++)
{
scanf("%d%d%d", &query[i].n, &query[i].m, &query[i].a);
query[i].id = i;
if (query[i].n > query[i].m) swap(query[i].n, query[i].m);
}
sort(query + 1, query + 1 + t);
int pt = 1;
for (int i = 1; i <= t; i++)//i表示第i组样例
{
for (; sigma[pt].sigma <= query[i].a; pt++)//处理每一个sigma
{
for (int k = sigma[pt].id; k <= maxn - 10; k += sigma[pt].id)
change(1, k, sigma[pt].sigma * mu[k / sigma[pt].id]);
}
ans[query[i].id] = cal(query[i].n, query[i].m);
}
for (int i = 1; i <= t; i++)
printf("%d\n", ans[i] & 2147483647);
}
int main()
{
solve();
return 0;
}