JZOJ 4813 Running

注意到每个 \(a_i\) 的贡献是 \(\gcd(a_i,n)\) 的所有倍数(在模 \(n\) 意义下这是一个周期)。

考虑枚举某个 \(\gcd(a_i,n)\) 的倍数为 \(d\) 然后用特定条件来去重。
即,枚举 \(d | n\)(显然是 \(n\) 的约数),如果 \(\exists i,\gcd(a_i,n) | d\),它的贡献为 \(\sum\limits_{i = 1}^n [\gcd(i,n) = d] = \varphi(\dfrac n d)\)

为什么呢?如果直接枚举 \(\gcd(a_i,n)\) 的倍数的话,显然会算重。
但是如果我们改为枚举它的倍数 \(d\),而只算 \(\gcd(i,n) = d\) 的,就不会重复和遗漏,而且复杂度可以降下来。
首先,由于每个数与 \(n\)\(\gcd\) 的结果是唯一的,所以每个数只可能被一个 \(d\) 算到贡献。
第二,因为 \(d\)\(\gcd(a_i,n)\) 的倍数,所以就算某个数应该有贡献但是还没算到,它肯定也会被另一个 \(d\) 算到。
最后,显然 \(d\) 得是 \(n\) 的约数才有贡献,而众所周知,\(n\) 的约数上界是 \(O(\sqrt n)\) 个。

然后注意到 \(\dfrac n d\) 可能达到 \(10^9\),所以没法线性筛。但是有用的 \(\dfrac n d\) 只有 \(\sqrt n\) 个,可以考虑分解质因数然后算。
具体来说,设 \(n = \prod\limits_{i = 1}^m {p_i}^{a_i}\),其中 \(p_i\) 为互不相同的质数(即 \(n\) 的不同质因子)。
\(\varphi(n) = \prod\limits_{i = 1}^m ({p_i}^{a_i} - {p_i}^{a_i-1})\)
这个懒得证明了,会用就行。

然鹅我突然脑抽写了杜教筛还给同学这样讲了,大家 \(\infty\) 脸懵逼……
好尴尬……

对了如果要杜教筛的话要先线性筛到 \(10^7\) 不然过不掉。

跟来自 GDDG 同在 JZ 集训的 DJT 同学讨论了一下,发现了一种基于容斥的做法。
这位同学是根据欧拉函数的计算式启发得到的。
\(a_i' = \gcd(a_i,n)\),那么答案为 \(\sum\limits_{S \subseteq [1,m]}(-1)^{|S|}n\prod\limits_{i \in S} \dfrac1{a_i'}\)
但是这样子是 \(O(2^n)\) 的。
不过,我们会发现这个等价于 \(n\prod\limits_{i=1}^m(1-\dfrac1{a_i'})\)
展开后我们会发现这两个式子是等价的,而后面这个式子同时也很像欧拉函数的计算式(把前文的公式中所有 \({p_i}^{a_i}\) 项拖到 \(\prod\) 外面就是了)。

对不起,这个算法是错的。
由于 \(\dfrac na,\dfrac nb\) 重复的并非 \(\dfrac n{ab}\),而是 \(\dfrac n{\text{lcm}(a,b)}\),所以这个算法不对。

代码:

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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#include <cstdio>
#include <algorithm>
#include <tr1/unordered_map>
using namespace std;
using namespace tr1;
const int N = 1e9;
const int MX = 1e7;
const int M = 50;
int n,m,a[M + 5];
int prime[MX + 5],cnt,vis[MX + 5];
int ans;
long long phi[MX + 5];
unordered_map<int,long long> mem;
long long calc(int n)
{
if(n <= MX)
return phi[n];
if(mem.count(n))
return mem[n];
long long ret = (long long)n * (n + 1) / 2;
for(register int l = 2,r;l <= n;l = r + 1)
{
r = n / (n / l);
ret -= (r - l + 1) * calc(n / l);
}
return ret;
}
void solve(int d)
{
for(register int i = 1;i <= m;++i)
if(!(d % __gcd(a[i],n)))
{
ans += calc(n / d) - calc(n / d - 1);
break;
}
}
int main()
{
freopen("running.in","r",stdin),freopen("running.out","w",stdout);
phi[1] = 1;
for(register int i = 2;i <= MX;++i)
{
if(!vis[i])
phi[prime[++cnt] = i] = i - 1;
for(register int j = 1;j <= cnt && i * prime[j] <= MX;++j)
{
vis[i * prime[j]] = 1;
if(!(i % prime[j]))
{
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
else
phi[i * prime[j]] = phi[i] * (prime[j] - 1);
}
phi[i] += phi[i - 1];
}
scanf("%d%d",&n,&m);
for(register int i = 1;i <= m;++i)
scanf("%d",a + i);
for(register int i = 1;i * i <= n;++i)
if(!(n % i))
{
solve(i);
if(i ^ n / i)
solve(n / i);
}
printf("%d\n",n - ans);
}