1. 程式人生 > >洛谷3768:簡單的數學題——題解

洛谷3768:簡單的數學題——題解

show itl 作者 uri tab IT .org 什麽 sdn

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:簡單的數學題——題解