JZOJ 4813.Running

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

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

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

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

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

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

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

对不起,这个算法是错的。
由于 na,nb\dfrac na,\dfrac nb 重复的并非 nab\dfrac n{ab},而是 nlcm(a,b)\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);
}