【LuoguP3768】简单的数学题-杜教筛+数论分块

测试地址：简单的数学题
做法：本题需要用到杜教筛+数论分块。
首先推式子，因为∑d|nφ(d)=n$\sum _{d|n}\phi (d)=n$，所以有：
∑ni=1∑nj=1ijgcd(i,j)$\sum _{i=1}^n\sum _{j=1}^{n}ijgcd(i,j)$
=∑ni=1i∑nj=1j∑d|gcd(i,j)φ(d)$=\sum _{i=1}^{n}i\sum _{j=1}^{n}j\sum _{d|gcd(i,j)}\phi (d)$
=∑nd=1φ(d)∑⌊nd⌋i=1id∑⌊nd⌋j=1jd$=\sum _{d=1}^n\phi (d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }id\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }jd$
=∑nd=1φ(d)⋅d2⋅(⌊nd⌋(⌊nd⌋+1)2)2$=\sum _{d=1}^n\phi (d)\cdot d^2\cdot (\frac{\lfloor \frac{n}{d}\rfloor (\lfloor \frac{n}{d}\rfloor +1)}{2}{)}^2$
显然如果我们能求f(n)=φ(n)⋅n2$f(n)=\phi (n)\cdot n^2$的前缀和的话，就可以数论分块做了。由于f$f$是积性函数，考虑杜教筛，令g(n)=n2$g(n)=n^2$，那么有(f∗g)(n)=n3$(f\ast g)(n)=n^3$，这两个函数都能O(1)$O(1)$算出前缀和，那么我们就解决了这一题。
要注意的一点是，因为数据范围很大，达到了1010${10}^{10}$，有的时候直接计算乘法会溢出，一定要在乘上一个1010${10}^{10}$大小的数前先对这个数取模。
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll hashsiz=8000009;
ll n,mod,limit,hashlist[hashsiz+10]={0},hashval[hashsiz+10];
ll phi[5000010],sum[5000010],prime[5000010],inv[7];
bool vis[5000010]={0};

void calc()
{
    phi[1]=1;
    prime[0]=0;
    for(ll i=2;i<=limit;i++)
    {
        if (!vis[i])
        {
            prime[++prime[0]]=i;
            phi[i]=i-1;
        }
        for(ll j=1;j<=prime[0]&&i*prime[j]<=limit;j++)
        {
            vis[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
    }
    sum[0]=0;
    for(ll i=1;i<=limit;i++)
        sum[i]=(sum[i-1]+i*i%mod*phi[i]%mod)%mod;
}

ll sumg(ll n)
{
    return (n%mod)*((n+1)%mod)%mod*((2*n+1)%mod)%mod*inv[6]%mod;
}

ll sumfg(ll n)
{
    return (n%mod)*(n%mod)%mod*((n+1)%mod)%mod*((n+1)%mod)%mod*inv[4]%mod;
}

ll hashfind(ll x)
{
    ll pos=x%hashsiz;
    while(hashlist[pos]&&hashlist[pos]!=x) pos++;
    if (hashlist[pos]) return pos;
    else return -1;
}

void hashinsert(ll x,ll v)
{
    ll pos=x%hashsiz;
    while(hashlist[pos]&&hashlist[pos]!=x) pos++;
    hashlist[pos]=x;
    hashval[pos]=v;
}

ll solve(ll n)
{
    ll pos=hashfind(n);
    if (n<=limit) return sum[n];
    if (pos!=-1) return hashval[pos];
    ll ans=sumfg(n);
    for(ll i=n;i>=2;i=n/(n/i+1))
    {
        ll l=max(n/(n/i+1)+1,2ll),r=i;
        ans-=solve(n/i)*(sumg(r)-sumg(l-1));
        ans=(ans+mod)%mod;
    }
    hashinsert(n,ans);
    return ans;
}

ll F(ll x)
{
    ll ans=(x%mod)*((x+1)%mod)%mod*inv[2]%mod;
    ans=ans*ans%mod;
    return ans;
}

int main()
{
    scanf("%lld%lld",&mod,&n);
    inv[1]=1;
    for(ll i=2;i<=6;i++)
        inv[i]=(mod-mod/i)*inv[mod%i]%mod;
    for(ll i=1;i*i*i<=n;i++)
        if ((i+1)*(i+1)*(i+1)>n)
        {
            limit=i*i;
            break;
        }
    calc();

    ll ans=0;
    for(ll i=n;i>=1;i=n/(n/i+1))
    {
        ll l=n/(n/i+1)+1,r=i;
        ans=((ans+F(n/i)*(solve(r)-solve(l-1)))%mod+mod)%mod;
    }
    printf("%lld",ans);

    return 0;
}