分析矿洞题解

暴力枚举激光的预期得分为$0$。

观察发现每个点$(x,y)$的价值是$gcd(x,y)^2$

其实要求的就是这个:

$\sum ^{N}_{i=1} \sum^{N}_{j=1} \varphi(gcd(i,j)^2)$

直接暴力求期望得分$30$分。

考虑反演：

$\sum ^{N}_{i=1} \sum^{N}_{j=1} \varphi(gcd(i,j)^2)$

$=\sum ^{N}_{i=1} \sum^{N}_{j=1} \varphi(gcd(i,j))*gcd(i,j)$

$=\sum ^{N}_{i=1} \sum^{N}_{j=1} \sum_{d|i,d|j} \varphi(d)*d*\varepsilon(gcd(i,j)==d)$

$=\sum ^{N}_{d=1} \varphi(d)*d*\sum ^{\lfloor N/d \rfloor}_{i=1} \sum^{\lfloor N/d \rfloor}_{j=1} \varepsilon(gcd(i,j))$

$=\sum ^{N}_{d=1} \varphi(d)*d*(2*\sum ^{\lfloor N/d \rfloor}_{i=1} \varphi(i)-1)$

直接分块做可以到$80$分。

$N=10^{10}$需要用杜教筛。那么$\sum^{N}_{d=1}  \varphi(d)*d$怎么快速筛？

$\sum^{N}_{i=1} i^2$

$=\sum^{N}_{i=1} i*\sum_{d|i} \varphi(d)$

$=\sum^{N}_{i=1} \sum_{d|i} \varphi(d)*d*\frac{i}{d}$

令$t=\frac{i}{d}$并先枚举$t$，考虑到$t*d=i$只有当小于等于$N$时才会对$\varphi(d)$产生贡献，所以被一个$t$贡献的$\varphi(d)$中的$d$必须要小于等于$\lfloor \frac{N}{t} \rfloor$（这里不是很好理解，可以理解为计算$t$对$\varphi(d)$的贡献。）：

$=\sum^{N}_{t=1} t *\sum^{\lfloor N/t \rfloor}_{d=1} \varphi(d)*d$

$=\sum^{N}_{d=1} \varphi(d)*d+\sum^{N}_{t=2} t \sum^{\lfloor N/t \rfloor}_{d=1} \varphi(d)*d$

这样我们得到

$\sum^{N}_{d=1} \varphi(d)*d=\sum^{N}_{i=1}i^2-\sum^{N}_{t=2} t \sum^{\lfloor N/t \rfloor}_{d=1} \varphi(d)*d$

这样就可以递归往下筛了。

线性筛预处理的范围的话，在$10^{10*\frac{2}{3}}=4.7*10^6$左右都没有问题。

好像有点卡常，调高了点时限，注意$(10^{10})^{2}$是爆$long\ long$的，所以求$\sum^{N}_{i=1} i^2$一类的值要先取模再代公式。

更坑的是$10^{10}*(10^9+7)$也正好爆$long \ long$（不爆$unsigned$），所以一路乘下去一边取模也会错，虽然几率有点小。（没有特殊构造卡这个的数据）

代码：

``` cpp
#include<bits/stdc++.h>
#define maxn 5000010
#define ll long long
using namespace std;

int pri[maxn],phi[maxn],pricnt=0;
ll sphi[maxn],s2[maxn],n,inv2,inv6;
const ll p=1e9+7;
bool notpri[maxn];
map<ll,ll> mp,mp2;
ll pw(ll a,ll b){
    ll ret=1;
    while(b){
        if(b&1)
            ret=ret*a%p;
        a=a*a%p;
        b>>=1;
    }
    return ret;
}
void init(){
    phi[1]=1;
    for(int i=2;i<maxn;i++){
        if(!notpri[i]){
            pri[++pricnt]=i;
            phi[i]=i-1;
        }
        for(int j=1;j<=pricnt;j++){
            if(i*pri[j]>=maxn)
                break;
            notpri[i*pri[j]]=true;
            if(i%pri[j]==0){
                phi[i*pri[j]]=phi[i]*pri[j];
                break;
            }
            phi[i*pri[j]]=phi[i]*(pri[j]-1);
        }
    }
    for(int i=1;i<maxn;i++)
        sphi[i]=(sphi[i-1]+phi[i])%p,s2[i]=(s2[i-1]+1LL*i*phi[i])%p;
    inv2=pw(2,p-2),inv6=pw(6,p-2);
}
ll c1(ll x){
    x%=p;
    return 1LL*x*(x+1)/2%p;
}
ll c2(ll x){
    x%=p;
    return 1LL*x*(x+1)%p*(2*x+1)%p*inv6%p;
}
ll _phi(ll x){
    if(x<maxn)
        return sphi[x];
    if(mp[x])
        return mp[x];
    ll pos,ret=c1(x);
    for(ll i=2;i<=x;i=pos+1){
        pos=x/(x/i);
        (ret-=1LL*(pos-i+1)%p*(_phi(x/i)%p))%=p;
    }
    return mp[x]=(ret+p)%p;
}
ll calc(ll x){
    if(x<maxn)
        return s2[x];
    if(mp2[x])
        return mp2[x];
    ll pos,ret=c2(x);
    for(ll i=2;i<=x;i=pos+1){
        pos=x/(x/i);
        (ret-=1LL*(pos+i)%p*((pos-i+1)%p*inv2%p)%p*calc(x/i))%=p;
    }
    return mp2[x]=(ret+p)%p;
}
ll solve(){
    ll pos,ret=0;
    for(ll i=1;i<=n;i=pos+1){
        pos=n/(n/i);
        (ret+=1LL*(calc(pos)-calc(i-1)+p)*(2*_phi(n/i)-1+p)%p)%=p;
    }
    return (ret+p)%p;
}
int main(){
    scanf("%lld",&n);
    init();
    printf("%lld\n",solve());
    return 0;
} 
```

正在渲染内容...

分析矿洞 题解