码迷,mamicode.com
首页 > 其他好文 > 详细

洛谷3768:简单的数学题——题解

时间:2018-03-20 20:47:31      阅读:193      评论:0      收藏:0      [点我收藏+]

标签:com   phi   简单的   map   pac   博客   分块   blank   name   

https://www.luogu.org/problemnew/show/P3768

题面来自洛谷,因为没用markdown所以直接截的图。

剩余的图是我用markdown写完然后截的图。

技术分享图片

参考洛谷第一篇题解。

这个式子直观感受就需要莫比乌斯反演,大致的过程参考:BZOJ2693:jzptab

那么跳过暴力推式子,我们能够得到:

技术分享图片

(如果你疑问为什么是miu(k/d)而不是miu(d),其实二者皆可,但是为什么这么干请往下看)

显然技术分享图片可以分块O(sqrt(n))做,那么后面的那一串东西怎么做呢?

参考jzptab的想法的话我们知道后面的式子是积性函数可以O(n)筛出来。

emmm……虽然时间可能能过但是你空间也存不下啊,所以我们要舍弃这个想法。

我们有这样的式子:技术分享图片

于是我们可以把后面的式子转换成这个,那么就变成了:

技术分享图片

phi的前缀和可以杜教筛做,k^2也是分块求和。

具体的杜教筛式子可以看:http://blog.csdn.net/samjia2000/article/details/70147436

#include<cstdio>
#include<queue>
#include<map>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=4700050;
ll su[N],phi[N],sum[N],n,p,inv2,inv6;
bool he[N];
ll maxn=4700000;
map<ll,ll>mp;
inline ll s2(ll x){x%=p;return x*(x+1)%p*inv2%p;}
inline ll s3(ll x){x%=p;return x*(x+1)%p*(2*x+1)%p*inv6%p;}
inline ll sqr(ll x){x%=p;return x*x%p;}
ll pow(ll x,ll y){
    ll res=1;
    while(y){
    if(y&1)res=res*x%p;
    x=x*x%p;
    y>>=1;
    }
    return res;
}
void Euler(int n){
    int tot=0;
    phi[1]=1;he[1]=1;
    for(int i=2;i<=n;i++){
    if(!he[i]){
        su[++tot]=i;
        phi[i]=i-1;
    }
    for(int j=1;j<=tot;j++){
        if(i*su[j]>n)break;
        he[i*su[j]]=1;
        if(i%su[j]==0){
        phi[i*su[j]]=phi[i]*su[j]%p;
        break;
        }
        else phi[i*su[j]]=phi[i]*phi[su[j]]%p;
    }
    }
    for(int i=1;i<=n;i++){
    phi[i]=(phi[i-1]+phi[i]*sqr(i)%p)%p;
    }
    return;
}
ll f(ll x){
    if(x<=maxn)return phi[x];
    if(mp[x])return mp[x];
    ll res=sqr(s2(x));
    for(ll i=2,j;i<=x;i=j+1){
    j=x/(x/i);
    res-=f(x/i)*(s3(j)-s3(i-1))%p;
    res%=p;
    }
    return mp[x]=(res+p)%p;
    
}
int main(){
    scanf("%lld%lld",&p,&n);
    inv2=pow(2,p-2);inv6=pow(6,p-2);
    maxn=min(maxn,n);
    Euler(maxn);
    ll ans=0;
    for(ll i=1,j;i<=n;i=j+1){
    j=n/(n/i);
    ans+=(f(j)-f(i-1))%p*sqr(s2(n/i))%p;
    ans%=p;
    }
    printf("%lld\n",(ans+p)%p);
    return 0;
}

+++++++++++++++++++++++++++++++++++++++++++

+本文作者:luyouqi233。               +

+欢迎访问我的博客:http://www.cnblogs.com/luyouqi233/+

+++++++++++++++++++++++++++++++++++++++++++

洛谷3768:简单的数学题——题解

标签:com   phi   简单的   map   pac   博客   分块   blank   name   

原文地址:https://www.cnblogs.com/luyouqi233/p/8611961.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!