P3768 简单的数学题(莫比乌斯反演+杜教筛)

最新推荐文章于 2022-04-02 21:01:31 发布

hongism

最新推荐文章于 2022-04-02 21:01:31 发布

阅读量183

点赞数

分类专栏：数论文章标签：莫比乌斯反演杜教筛

本文链接：https://blog.csdn.net/qq_40851016/article/details/102479697

版权

数论专栏收录该内容

3 篇文章 0 订阅

订阅专栏

P3768 简单的数学题

题意：

输入 $p(5\times 10^8\leq p\leq 1.1\times 10^9且p为质数),n(10^{10})$
求 $(\sum_{i=1}^{n}\sum_{j=1}^{n}ijgcd(i,j))\mod p$ 。

题解：

反演：
$\sum_{i=1}^{n}\sum_{j=1}^nijgcd(i,j)$
$=\sum_{d=1}^nd\sum_{i=1}^{n}\sum_{j=1}^{n}ij[gcd(i,j)==d]$
$=\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{n}{d}\rfloor}ij[gcd(i,j)==1]$
$=\sum_{d=1}^{n}d^3\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ij\sum_{k\mid gcd(i,j)}\mu(k)$
$=\sum_{d=1}^{n}d^3\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}k^2\sum_{i=1}^{\lfloor\frac{n}{dk}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{dk}\rfloor }ij\mu(k)$
$=\sum_{d=1}^{n}d^3\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor }k^2(sum(\lfloor\frac{n}{dk}\rfloor))^2\mu(k)$
令 $T = d k$ ；其中 $sum(i)=1+2+\dots+i$
$=\sum_{T=1}^{n}(sum(\lfloor\frac{n}{T}\rfloor))^2\sum_{d\mid T}d^3(\frac{T}{d})^2\mu(\frac{T}{d})$
$=\sum_{T=1}^{n}(sum(\lfloor\frac{n}{T}\rfloor))^2T^2\sum_{d\mid T}d\mu(\frac{T}{d})$
$=\sum_{T=1}^{n}(sum(\lfloor\frac{n}{T}\rfloor))^2T^2\phi(T)$
杜教筛：
$sum(\lfloor\frac{n}{T}\rfloor)$ 这一部分数论分块就好
$令f[i]=i^2\phi(i);g[i]=i$ ；
$f*g(n)=n^2$ 于是套用杜教筛：
$S[n]=\frac{\sum_{i=1}^nf*g(i)-\sum_{i=2}^{n}g[i]\times S[\lfloor\frac{n}{l}\rfloor]}{g[1]}$

注意： 除法尽量乘逆元，不然容易出错。

代码：

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define int ll
#define i128 __int128
ll p,n;
const int N=1e6+9;
int primes[N],tot;
bool prime[N];
ll f[N],he[N];
int inv2,inv6;
int poww(int a,int b,int c){
    int ans=1,base=a;
    while(b){
        if(b&1)ans=1ll*ans*base%c;
        base=1ll*base*base%c;
        b>>=1;
    }
    return ans;
}
void init(){
    f[1]=he[1]=1;
    for(int i=2;i<N;i++){
        if(!prime[i])primes[++tot]=i,f[i]=1ll*i*i%p*(i-1)%p;
        for(int j=1;j<=tot&&primes[j]*i<N;j++){
            prime[primes[j]*i]=1;
            if(i%primes[j]==0){
                f[i*primes[j]]=f[i]*primes[j]%p*primes[j]%p*primes[j]%p;
                break;
            }else f[i*primes[j]]=f[i]*f[primes[j]]%p;
        }
    }
    for(int i=2;i<N;i++)he[i]=(he[i-1]+f[i])%p;
}
unordered_map<ll,ll>mp;
ll gi(ll n){
    ll n1=n,n2=(n+1),n3=(2*n+1);
    if(n1%2==0)n1/=2;else n2/=2;
    if(n1%3==0)n1/=3;else if(n2%3==0)n2/=3;else n3/=3;
    return i128(n1)*n2%p*n3%p;
}
ll sum(ll x){
    if(x<N)return he[x];
    if(mp.count(x))return mp[x];
    ll t=i128(x)*(x+1)/2%p;t=t*t%p;
    for(ll l=2,r;l<=x;l=r+1){
        r=x/(x/l);
        ll u=(gi(r)+p-gi(l-1))%p;
        t=(t+p-sum(x/l)*u%p)%p;
    }
    return mp[x]=t;
}
signed main(){
  //  freopen("tt.in","r",stdin),freopen("tt.out","w",stdout);
    cin>>p>>n;
    init();
    inv2=poww(2,p-2,p),inv6=poww(6,p-2,p);
    ll ans=0;
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        ll t=n/l;
        t=i128(t)*(t+1)/2%p;
      //  t=(i128)t*(t+1)%p*inv2%p;
        ll y=(sum(r)+p-sum(l-1))%p;
        ans=(ans+t*t%p*y%p)%p;
    }
    cout<<ans<<endl;
    return 0;
}