P3768 簡單的數學題

題目描述

由於出題人懶得寫背景了，題目還是簡單一點好。

輸入一個整數\(n\)和一個整數\(p，\)你需要求出\((\sum_{i=1}^n\sum_{j=1}^n ijgcd(i,j)) \bmod p\)，其中\(gcd(a,b)\)表示\(a\)與\(b\)的最大公約數。

剛才題面打錯了，已修改

輸入輸出格式

輸入格式：

一行兩個整數\(p\)、\(n\)。

輸出格式：

一行一個整數\((\sum_{i=1}^n\sum_{j=1}^n ijgcd(i,j))\bmod p\)。

說明

對於\(20\%\)的數據，\(n \leq 1000\)。

對於\(30\%\)的數據，\(n \leq 5000\)

對於\(60\%\)的數據，\(n \leq 10^6\)，時限\(1s\)。

對於另外\(20\%\)的數據，\(n \leq 10^9\)，時限\(3s\)。

對於最後\(20\%\)的數據，\(n \leq 10^{10}\)，時限\(6s\)。

對於\(100\%\)的數據，\(5 \times 10^8 \leq p \leq 1.1 \times 10^9\)且\(p\)為質數。

從各種方向推推式子，你會差不多發現有

\[\sum_{i=1}^n\sum_{j=1}^n ijgcd(i,j)\]

\[=\sum_{T=1}^nF(\lfloor\frac{n}{T}\rfloor)^2T^2\varphi(T)\]

其中\(F(n)=\sum\limits_{i=1}^ni\)

然後上杜教篩設\(\mathbf f(n)=n^2\varphi(n)\)，則有

\[\mathbf {Id}^3=\mathbf f *\mathbf {Id}^2\]

帶進去杜教篩得到

\[\mathbf s(n)=\sum_{i=1}^n i^3-\sum_{i=2}^n i^2 \mathbf s(\lfloor\frac{n}{i}\rfloor)\]

然後小學奧數一波算前綴和就行了

小心爆\(long \ long\)

Code:

#include <cstdio>
#include <unordered_map>
#define ll long long
const int N=5e6;
ll n,mod,phi[N+10],inv2,inv6;
int pri[N+10],ispri[N+10],cnt;
ll qp(ll d,ll k){ll re=1;while(k){if(k&1)re=re*d%mod;d=d*d%mod,k>>=1;}return re;}
ll f(ll x){x%=mod;return x*(x+1)%mod*inv2%mod;}
ll g(ll x){x%=mod;return x*(x+1)%mod*(2*x%mod+1)%mod*inv6%mod;}
void init()
{
    phi[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!ispri[i])
        {
            phi[i]=i-1;
            pri[++cnt]=i;
        }
        for(int j=1;j<=cnt&&i*pri[j]<=N;j++)
        {
            ispri[i*pri[j]]=1;
            if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j]%mod;break;}
            else phi[i*pri[j]]=phi[i]*(pri[j]-1)%mod;
        }
    }
    for(int i=1;i<=N;i++)
        phi[i]=(phi[i]*i%mod*i%mod+phi[i-1])%mod;
}
std::unordered_map <ll,ll> Phi;
ll calphi(ll n)
{
    if(n<=N) return phi[n];
    if(Phi.find(n)!=Phi.end()) return Phi[n];
    ll ret=f(n)*f(n)%mod;
    for(ll l=2,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        (ret-=(calphi(n/l)*(g(r)-g(l-1))%mod))%=mod;
    }
    ret=(ret%mod+mod)%mod;
    return Phi[n]=ret;
}
int main()
{
    scanf("%lld%lld",&mod,&n);
    init();
    ll ans=0;inv6=qp(6,mod-2);inv2=qp(2,mod-2);
    for(ll l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        (ans+=f(n/l)*f(n/l)%mod*(calphi(r)-calphi(l-1))%mod)%=mod;
    }
    ans=(ans%mod+mod)%mod;
    printf("%lld\n",ans);
    return 0;
}
 

2018.11.26

洛谷 P3768 簡單的數學題 解題報告