P3768 简单的数学题

洛谷题目链接

毒瘤题

弄模数弄了我一天，结果少模了一个地方

做法

我们要求：$$(sumlimits_{i=1}^n sumlimits_{j=1}^n i imes j imes gcd(i,j))mod p$$

我们先不管后面的模(p)，考虑前面枚举(gcd(i,j))：$$sumlimits_{d=1}^n d sumlimits_{i=1}^n sumlimits_{j=1}^n i imes j [gcd(i,j)=d]$$

套路的除(d)：$$sumlimits_{d=1}^n d^3 sumlimits_{i=1}^{leftlfloor frac{n}{d} ight floor} sumlimits_{j=1}^{leftlfloor frac{n}{d} ight floor} i imes j [gcd(i,j)=1]$$

后面显然用莫比乌斯反演，我们用性质来变换：$$sumlimits_{d=1}^n d^3 sumlimits_{i=1}^{leftlfloor frac{n}{d} ight floor} sumlimits_{j=1}^{leftlfloor frac{n}{d} ight floor} i imes j sumlimits_{k|gcd(i,j)} mu (k)$$

套路的枚举(k)：

[sumlimits_{d=1}^n d sumlimits_{k=1}^{leftlfloor frac{n}{d} ight floor} mu (k) sumlimits_{i=1}^{leftlfloor frac{n}{d} ight floor} sumlimits_{j=1}^{leftlfloor frac{n}{d} ight floor} i imes j [k|gcd(i,j)] ]

除(k)把后面去掉：$$sumlimits_{d=1}^n d^3 sumlimits_{k=1}^{leftlfloor frac{n}{d} ight floor} mu (k) k^{2sumlimits_{i=1}}{leftlfloor frac{n}{dk} ight floor} sumlimits_{j=1}^{leftlfloor frac{n}{dk} ight floor} i imes j$$

我们用(T=dk)带入：$$sumlimits_{T=1}^n sumlimits_{d|T} d^3 mu (frac{T}{d}) {(frac{T}{d})}^{2sumlimits_{i=1}}{leftlfloor frac{n}{T} ight floor} sumlimits_{j=1}^{leftlfloor frac{n}{T} ight floor} i imes j$$

前面的(d^3)可以和({(frac{T}{d})}^2)消去分母：$$sumlimits_{T=1}^n sumlimits_{d|T} d mu (frac{T}{d}) T^{2sumlimits_{i=1}}{leftlfloor frac{n}{T} ight floor} sumlimits_{j=1}^{leftlfloor frac{n}{T} ight floor} i imes j$$

换个位置：$$sumlimits_{T=1}^n T^2 sumlimits_{d|T} d mu (frac{T}{d}) sumlimits_{i=1}^{leftlfloor frac{n}{T} ight floor} sumlimits_{j=1}^{leftlfloor frac{n}{T} ight floor} i imes j$$

我们设函数(S)：$$S(n)=(sumlimits_{i=1}^n i)^2$$

那么原式可化为：$$sumlimits_{T=1}^n T^2 sumlimits_{d|T} d mu (frac{T}{d}) S(leftlfloor frac{n}{T} ight floor)$$

因为(φ=id imes mu)，所以：$$sumlimits_{T=1}^n T^2 φ(T) S(leftlfloor frac{n}{T} ight floor)$$

那么后面我们可以用整除分块来搞，可是前面呢？

杜教筛大法吼啊~~~

我们要求：$$L(n)=sumlimits_{i=1}^n i^2φ(i)$$

那么就要考虑怎么配(g,h)函数：(h(n)=sumlimits_{d|n} d^2 φ(d) g(frac{n}{d}))$

我们发现前面的(d^2)很讨厌，自然要想着消去，于是我们设(g(i)=i^2)：$$h(n)=sumlimits_{d|n} d^2 φ(d) (frac{n}{d})^2$$

[h(n)=sumlimits_{d|n} n^2φ(d) ]

[h(n)=n^2sumlimits_{d|n} φ(d) ]

[h(n)=n^3 ]

那么我们发现这样配对的话可以很容易求出(g,h)，（别说你没学过小学奥数）

把杜教筛的套路式弄出来：$$g(1)L(n)=sumlimits_{i=1}^n h(i) -sumlimits_{d=2}^n g(d)L(frac{n}{d})$$

带入(g,h)函数：$$L(n)=sumlimits_{i=1}^n i^3 -sumlimits_{d=2}^n d^2L(frac{n}{d})$$

而我们知道：$$sumlimits_{i=1}^n i^2 =frac{n(1+n)(2n+1)}{6}$$

[sumlimits_{i=1}^n i^3=(sumlimits_{i=1}^n i)^2=frac{n^2(1+n)^2}{4} ]

所以很容易用杜教筛求出来啦

要注意要模(p)，所以要算出(4,6)的逆元

接下来是美滋滋的代码时间~~~

#include<iostream>
#include<cstdio>
#include<cmath>
#include<tr1/unordered_map>
#define int long long
#define N 8000007
using namespace std;
int n,mod,cnt,inv4,inv6;
int phi[N],prime[N];
bool isp[N];
tr1::unordered_map<int,int> map_;
void Get_phi(int x)
{
    phi[1]=1;
    for(int i=2;i<=x;++i)
    {
        if(!isp[i])
        {
            prime[++cnt]=i;
            phi[i]=i-1;
        }
        for(int j=1;j<=cnt&&(i*prime[j])<=x;++j)
        {
            isp[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                phi[i*prime[j]]=(phi[i]*prime[j])%mod;
                break;
            }
            else
                phi[i*prime[j]]=(phi[i]*phi[prime[j]])%mod;
        }
    }
    for(int i=1;i<=x;++i)
        phi[i]=(phi[i-1]+phi[i]*i%mod*i%mod)%mod;
}
int qpow(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)
            ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}
int L(int x)
{
    x%=mod;
    int num=x%mod*(x+1)%mod;
    num=num*num%mod*inv4%mod;
    return num;
}
int L1(int x)
{
    x%=mod;
    int num=x*(x+1)%mod*(x+x+1)%mod*inv6%mod;
    return num;
}
int Ans(int x)
{
    if(x<=N-7)
        return phi[x];
    if(map_[x])
        return map_[x];
    int ans=L(x);
    for(int l=2,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        ans-=(Ans(x/l)*(L1(r)-L1(l-1)+mod)%mod)%mod;
        ans%=mod;
    }
    //printf("%lld
",ans);
    return map_[x]=(ans+mod)%mod;
}
signed main()
{
    scanf("%lld%lld",&mod,&n);
    inv4=qpow(4,mod-2),inv6=qpow(6,mod-2);
    Get_phi(N-7);
    int ans=0;
    for(int l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans+=(((Ans(r)-Ans(l-1)+mod)%mod)*(L(n/l)%mod))%mod;
        ans%=mod;
        //printf("%lld
",ans);
    }
    printf("%lld",(ans+mod)%mod);
    return 0;
}
/*
1000000007 9786510294 

*/