JZOJ 4813 Running

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

考虑枚举某个 gcd(ai,n) 的倍数为 d 然后用特定条件来去重。
即，枚举 d|n（显然是 n 的约数），如果 ∃i,gcd(ai,n)|d，它的贡献为 n∑i=1[gcd(i,n)=d]=φ(nd)。

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

然后注意到 nd 可能达到 109，所以没法线性筛。但是有用的 nd 只有 √n 个，可以考虑分解质因数然后算。
具体来说，设 n=m∏i=1piai，其中 pi 为互不相同的质数（即 n 的不同质因子）。
有 φ(n)=m∏i=1(piai−piai−1)。
这个懒得证明了，会用就行。

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

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

跟来自 GDDG 同在 JZ 集训的 DJT 同学讨论了一下，发现了一种基于容斥的做法。
这位同学是根据欧拉函数的计算式启发得到的。
设 a′i=gcd(ai,n)，那么答案为 ∑S⊆[1,m](−1)|S|n∏i∈S1a′i。
但是这样子是 O(2n) 的。
不过，我们会发现这个等价于 nm∏i=1(1−1a′i)。
展开后我们会发现这两个式子是等价的，而后面这个式子同时也很像欧拉函数的计算式（把前文的公式中所有 piai 项拖到 ∏ 外面就是了）。

对不起，这个算法是错的。
由于 na,nb 重复的并非 nab，而是 nlcm(a,b)，所以这个算法不对。

代码：

#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);
}