[luogu3768] 简单的数学题 [杜教筛]

题面：

传送门

实际上就是求：

思路：

看到gcd就先反演一下，过程大概是这样：

明显的一步反演

这里设，S(x)等于1到x的和

然后把枚举d再枚举T变成先枚举T再枚举其约数d，变形：

后面其中两项展开，把T提出来

S那里可以数论分块，那么只要S后面那个东西可以筛出来，就可以O(sqrt(n))

发现后面的那部分可以狄利克雷卷积一波

这明显是一个积性函数，但是n有10^10，所以不能线筛

考虑使用杜教筛，令上述函数为f，函数S为f的前缀和

套用杜教筛模板式

现在问题就是选一个合适的g函数了

我们知道欧拉函数有一个卷积性质：

那么我们令g(x)=x^2

此时g与f的卷积变成了：

看起来真是赏心悦目

于是杜教筛的递推式变成了这样的：

一波记忆化搜索带走AC

Code：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<map>
#define ll long long
using namespace std;
inline ll read(){
    ll re=0,flag=1;char ch=getchar();
    while(ch>'9'||ch<'0'){
        if(ch=='-') flag=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9') re=(re<<1)+(re<<3)+ch-'0',ch=getchar();
    return re*flag;
}
ll n,MOD,inv2,inv6;
ll phi[8000010],pri[1000010],tot=0;bool vis[8000010];
void init(){
    ll i,j,k;phi[1]=1;phi[0]=0;
    for(i=2;i<=8000000;i++){
        if(!vis[i]){
            pri[++tot]=i;phi[i]=i-1;
        }
        for(j=1;j<=tot;j++){
            k=i*pri[j];if(k>8000000) break;
            vis[k]=1;
            if(i%pri[j]==0){
                phi[k]=1ll*phi[i]*pri[j]%MOD;
                break;
            }
            phi[k]=1ll*phi[i]*phi[pri[j]]%MOD;
        }
    }
    for(i=2;i<=8000000;i++) phi[i]=(phi[i-1]+1ll*i*i%MOD*phi[i]%MOD)%MOD;
}
ll sum1(ll x){x%=MOD;return x*(x+1)%MOD*inv2%MOD;}
ll sum2(ll x){x%=MOD;return x*(x+1)%MOD*((x<<1)+1)%MOD*inv6%MOD;}
map<ll,ll>m;
ll S(ll x){
    if(x<=8000000) return phi[x];
    if(m[x]) return m[x];
    ll re=sum1(x),tmp;re=re*re%MOD;ll i,j;
    for(i=2;i<=x;i=j+1){
        j=x/(x/i);
        tmp=sum2(j)-sum2(i-1);tmp=(tmp+MOD)%MOD;
        re-=S(x/i)*tmp%MOD;re%=MOD;
    }
    return m[x]=(re+MOD)%MOD;
}
ll fpow(ll a,ll b){
    ll re=1;a%=MOD;
    while(b){
        if(b&1) re=a*re%MOD;
        b>>=1;a=a*a%MOD;
    }
    return re;
}
int main(){
    MOD=read();n=read();ll i,j;ll ans=0,tmp,tt;
    inv2=fpow(2,MOD-2);inv6=fpow(6,MOD-2);init();
    for(i=1;i<=n;i=j+1){
        j=n/(n/i);
        tmp=sum1(n/i);
        tmp=(tmp+MOD)%MOD;tmp=(tmp*tmp)%MOD;
        tt=S(j)-S(i-1);
        tt=(tt+MOD)%MOD;
        ans=(ans+tmp*tt%MOD)%MOD;
    }
    printf("%lld
",(ans+MOD)%MOD);
}