看了看唐老師的blog，照貓畫虎的做了幾道題目，感覺對杜教篩有些感覺了
但是稍微有一點難度的題目還是做不出來，放假的時候爭取都A掉(挖坑ing)
這篇文章以后等我A掉那些題目之后再UPD上去就好啦
由於懶得去寫怎么用編輯器寫公式，所以公式就准備直接copy唐老師的啦
首先積性函數和完全積性函數什么的就不再多說了
列舉常見的積性函數：
1、約數個數函數和約數個數和函數
2、歐拉函數phi
3、莫比烏斯函數mu
4、元函數e 其中e(n)=[n==1]
5、恆等函數I 其中I(n)=1
6、單位函數id 其中id(n)=n
有關狄利克雷卷積
定義：
我們簡單記為f*g
注意，狄利克雷卷積滿足交換律，結合律，對於加法滿足分配率
兩個很顯然的事情是：1、f*e=e*f=f
2、設h=f*g，若f和g是積性函數，則h是積性函數
積性函數的一些性質：
1、
2、e=mu*I
3、id=phi*I
之后我們由第二條性質可以得到一些很有趣的結論：
設存在函數g=f*I (通過f求g）
那么我們可以得到g*mu=f*I*mu
進而化簡得到f=g*mu (通過g求f)
唐老師的博客里還舉了一些更為有趣的例子
好了，前置技能都說完了，我們來說正題
杜教篩基本可表示成如下形式：
求F(n)=sigma(f(i))
存在g=f*I，定義G(n)=simga(g(i))
就可以得到F(n)=G(n)-sigma(F(n/i)) (這里i從2開始)
如果G(n)是可以在一定時間內求解(不一定要求O(1))
那么我們就可以做到O(n^(3/4))求解F(n)
如果加一些預處理我們可以做到O(n^(2/3))求解F(n)
下面通過幾個例子來說一說為什么是這樣的：
1、求sigma(phi(i))
由於我們知道id=phi*I
而且id的前綴和是可以O(1)的計算的
所以我們可以做以下推導
對，沒錯，推完了之后就是我們上面提到過的形式，那么至於證明上面的形式的正確性就可以用這種推導方式啦
時間復雜度的分析:
這個結論非常的神奇，因為看上去復雜度像是sqrt(n)*sqrt(n)=n的，有興趣的話可以自己化簡一下，時間復雜度的確是上述所說
2、求sigma(mu(i))
有了上面的經驗，我們又知道e=mu*I
那么我們就可以直接寫了
這樣我們就機智的做出來了
唐老師的blog上還有一道題目，不過公式太長了，就不說了
有了上面學習的一些東西，我們就可以刷一些水題了
(感覺還是很好理解和學習的呢）
51Nod 莫比烏斯函數之和
直接套用上面的公式即可，但是真正寫起來還是要注意一些事情的
由於杜教篩的時間復雜度的原因，通常情況下給定的n是超過int范圍內，這就要求我們在計算過程中要時刻小心算術溢出
第二點是預處理多少比較合適，自己動手實驗了一下，唐老師說的不錯，大概預處理n^(2/3)較優
第三點就是在寫的過程中要把每次遞歸計算好的函數用哈希表存進來實現記憶化，不然時間復雜度難以保證
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<algorithm>
#include<cstring>
using namespace std;
const int mod=1333331;
typedef long long LL;
LL a,b;
LL mu[5000010];
int p[500010],cnt=0;
bool vis[5000010];
struct HASHMAP{
int h[mod+10],cnt;
int next[100010];
LL st[100010],S[100010];
void push(LL k,LL v){
int key=k%mod;
for(int i=h[key];i;i=next[i]){
if(S[i]==k)return;
}
++cnt;next[cnt]=h[key];h[key]=cnt;
S[cnt]=k;st[cnt]=v;
}
LL ask(LL k){
int key=k%mod;
for(int i=h[key];i;i=next[i]){
if(S[i]==k)return st[i];
}return -1;
}
}H;
void Get_Prime(){
mu[1]=1;
for(int i=2;i<=5000000;++i){
if(!vis[i]){
p[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt;++j){
if(1LL*p[j]*i>5000000)break;
int ip=i*p[j];
vis[ip]=true;
if(i%p[j]==0)break;
mu[ip]=-mu[i];
}
}
for(int i=2;i<=5000000;++i)mu[i]+=mu[i-1];
}
LL mu_sum(LL n){
if(n<=5000000)return mu[n];
LL tmp=H.ask(n);
if(tmp!=-1)return tmp;
LL la,A=1;
for(LL i=2;i<=n;i=la+1){
LL now=n/i;
la=n/now;
A=A-(la-i+1)*mu_sum(n/i);
}H.push(n,A);return A;
}
int main(){
scanf("%lld%lld",&a,&b);
Get_Prime();
printf("%lld\n",mu_sum(b)-mu_sum(a-1));
return 0;
}
51Nod 歐拉函數之和
直接寫就可以了
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<algorithm>
#include<cstring>
using namespace std;
const int mod=1333331;
typedef long long LL;
LL a,b;
LL mu[5000010];
int p[500010],cnt=0;
bool vis[5000010];
struct HASHMAP{
int h[mod+10],cnt;
int next[100010];
LL st[100010],S[100010];
void push(LL k,LL v){
int key=k%mod;
for(int i=h[key];i;i=next[i]){
if(S[i]==k)return;
}
++cnt;next[cnt]=h[key];h[key]=cnt;
S[cnt]=k;st[cnt]=v;
}
LL ask(LL k){
int key=k%mod;
for(int i=h[key];i;i=next[i]){
if(S[i]==k)return st[i];
}return -1;
}
}H;
void Get_Prime(){
mu[1]=1;
for(int i=2;i<=5000000;++i){
if(!vis[i]){
p[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt;++j){
if(1LL*p[j]*i>5000000)break;
int ip=i*p[j];
vis[ip]=true;
if(i%p[j]==0)break;
mu[ip]=-mu[i];
}
}
for(int i=2;i<=5000000;++i)mu[i]+=mu[i-1];
}
LL mu_sum(LL n){
if(n<=5000000)return mu[n];
LL tmp=H.ask(n);
if(tmp!=-1)return tmp;
LL la,A=1;
for(LL i=2;i<=n;i=la+1){
LL now=n/i;
la=n/now;
A=A-(la-i+1)*mu_sum(n/i);
}H.push(n,A);return A;
}
int main(){
scanf("%lld%lld",&a,&b);
Get_Prime();
printf("%lld\n",mu_sum(b)-mu_sum(a-1));
return 0;
}
BZOJ 3944 Sum
上面兩道題目的混合版，唯一的坑點就是如果你想用int讀入，在i=la+1的時候會有爆掉的風險
#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cstdlib>
using namespace std;
typedef long long LL;
const int mod=1333331;
int T;
int p[500010],cnt=0;
bool vis[2000010];
LL mu[2000010],phi[2000010];
LL n;
struct HASHMAP{
int h[mod+10],cnt;
int next[3000010];
LL S[3000010],M[3000010],P[3000010];
void push(LL k,LL v1,LL v2){
int key=k%mod;
if(key<0)key+=mod;
for(int i=h[key];i;i=next[i]){
if(S[i]==k)return;
}
++cnt;next[cnt]=h[key];h[key]=cnt;
S[cnt]=k;M[cnt]=v1;P[cnt]=v2;
}
pair<LL,LL>ask(LL k){
int key=k%mod;
if(key<0)key+=mod;
for(int i=h[key];i;i=next[i]){
if(S[i]==k)return make_pair(M[i],P[i]);
}return make_pair(-1,-1);
}
}H;
void Get_Prime(){
mu[1]=1;phi[1]=1;
for(int i=2;i<=2000000;++i){
if(!vis[i]){
p[++cnt]=i;
phi[i]=i-1;mu[i]=-1;
}
for(int j=1;j<=cnt;++j){
if(1LL*p[j]*i>2000000)break;
int ip=i*p[j];
vis[ip]=true;
if(i%p[j]==0){
phi[ip]=phi[i]*p[j];
break;
}mu[ip]=-mu[i];phi[ip]=phi[i]*(p[j]-1);
}
}
for(int i=2;i<=2000000;++i)phi[i]+=phi[i-1],mu[i]+=mu[i-1];
}
pair<LL,LL>sum(LL n){
if(n<=2000000)return make_pair(mu[n],phi[n]);
pair<LL,LL>tmp=H.ask(n);
if(tmp.first!=-1)return tmp;
LL A=1,B=1LL*n*(n+1)/2;LL la;
for(LL i=2;i<=n;i=la+1){
LL now=n/i;
la=n/now;
tmp=sum(now);
A=A-1LL*tmp.first*(la-i+1);
B=B-1LL*tmp.second*(la-i+1);
}H.push(n,A,B);
return make_pair(A,B);
}
int main(){
scanf("%d",&T);
Get_Prime();
while(T--){
scanf("%lld",&n);
pair<LL,LL>ans=sum(n);
printf("%lld %lld\n",ans.second,ans.first);
}return 0;
}
hdu 5608
QAQ 標准的不能在標准的杜教篩形式
首先n^2-3n+2的前綴和我們可以O(1)的求出來(表告訴我你不會平方前綴和公式）
之后我們考慮如果預處理f，我們設g(n)=n^2-3n+2
可以得到g=f*I，進而得到g*mu=f*I*mu,則f=g*mu
我自己比較懶，並不知道f是不是積性函數(貌似不是)
然后直接O(nlogn)的貢獻過去的，這樣就可以預處理前100w的f了
之后運用杜教篩就可以O(n^(2/3))求出答案了
#include<cstdio>
#include<cstring>
#include<iostream>
#include<cstdlib>
#include<algorithm>
using namespace std;
const int mod=1e9+7;
const int MOD=1333331;
typedef long long LL;
int T,n,inv1,inv2;
bool vis[1000010];
int p[1000010],cnt=0;
int mu[1000010];
int f[1000010];
struct HASHMAP{
int h[MOD+10],cnt;
int next[3000010];
int S[3000010],st[3000010];
void push(int k,int v){
int key=k%MOD;
for(int i=h[key];i;i=next[i]){
if(S[i]==k)return;
}
++cnt;next[cnt]=h[key];h[key]=cnt;
S[cnt]=k;st[cnt]=v;
}
int ask(int k){
int key=k%MOD;
for(int i=h[key];i;i=next[i]){
if(S[i]==k)return st[i];
}return -1;
}
}H;
int pow_mod(int v,int p){
int tmp=1;
while(p){
if(p&1)tmp=1LL*tmp*v%mod;
v=1LL*v*v%mod;p>>=1;
}return tmp;
}
void Get_Prime(){
mu[1]=1;
for(int i=2;i<=1000000;++i){
if(!vis[i]){
p[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt;++j){
if(1LL*p[j]*i>1000000)break;
int ip=i*p[j];
vis[ip]=true;
if(i%p[j]==0)break;
mu[ip]=-mu[i];
}
}return;
}
void Get_f(){
for(int i=1;i<=1000000;++i){
int now=(1LL*i*i-3*i+2)%mod;
int cnt=0;
for(int j=i;j<=1000000;j+=i){
cnt++;
f[j]=f[j]+now*mu[cnt];
if(f[j]<0)f[j]+=mod;
if(f[j]>=mod)f[j]-=mod;
}
}
for(int i=2;i<=1000000;++i){
f[i]+=f[i-1];
if(f[i]>=mod)f[i]-=mod;
}
return;
}
int sum(int n){
if(n<=1000000)return f[n];
int tmp=H.ask(n);
if(tmp!=-1)return tmp;
int A=1LL*n*(n+1)%mod*(2*n+1)%mod*inv1%mod;
A=A-3LL*n*(n+1)%mod*inv2%mod;
if(A<0)A+=mod;
A=A+2*n%mod;
if(A>=mod)A-=mod;
int la;
for(int i=2;i<=n;i=la+1){
int now=n/i;
la=n/now;
A=A-1LL*sum(n/i)*(la-i+1)%mod;
if(A<0)A+=mod;
}H.push(n,A);
return A;
}
int main(){
scanf("%d",&T);
Get_Prime();Get_f();
inv1=pow_mod(6,mod-2);
inv2=pow_mod(2,mod-2);
while(T--){
scanf("%d",&n);
printf("%d\n",sum(n));
}return 0;
}
51Nod 最大公約數之和
貌似網上並沒有題解，那我就來發一份吧，沒准能騙騙訪問量
注意到這道題是n和n，不是n和m
我們不妨枚舉最大公約數k
又因為我們知道(i,j)=k成立當且僅當(i/k,j/k)=1成立
然后我們就可以把原來的式子轉化一下變成：
設P(n)=sigma(phi(i))
原式=2*sigma(k*P(n/k))-n*(n+1)/2
式子轉化的正確性是顯然的
這樣我們利用杜教篩就可以求出P(n/k)的值
由於我們知道杜教篩求解所有n的約數k的所有的P(k)的值只需要O(n^(2/3))
所以我們在外面套上一個分塊也不會影響時間復雜度
總時間復雜度O(sqrt(n)+n^(2/3))
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cstdlib>
using namespace std;
const int mod=1e9+7;
const int MOD=1333331;
typedef long long LL;
int inv;
int p[500010],cnt=0;
int phi[5000010];
bool vis[5000010];
LL n;
struct HASHMAP{
int h[MOD+10],cnt;
int next[100010],st[100010];
LL S[100010];
void push(LL k,int v){
int key=k%MOD;
for(int i=h[key];i;i=next[i]){
if(S[i]==k)return;
}
++cnt;next[cnt]=h[key];h[key]=cnt;
S[cnt]=k;st[cnt]=v;
}
int ask(LL k){
int key=k%MOD;
for(int i=h[key];i;i=next[i]){
if(S[i]==k)return st[i];
}return -1;
}
}H;
int pow_mod(int v,int p){
int tmp=1;
while(p){
if(p&1)tmp=1LL*tmp*v%mod;
v=1LL*v*v%mod;p>>=1;
}return tmp;
}
void Get_Prime(){
phi[1]=1;
for(int i=2;i<=5000000;++i){
if(!vis[i]){
p[++cnt]=i;
phi[i]=i-1;
}
for(int j=1;j<=cnt;++j){
if(1LL*p[j]*i>5000000)break;
int ip=i*p[j];
vis[ip]=true;
if(i%p[j]==0){
phi[ip]=phi[i]*p[j];
break;
}phi[ip]=phi[i]*(p[j]-1);
}
}
for(int i=2;i<=5000000;++i){
phi[i]+=phi[i-1];
if(phi[i]>=mod)phi[i]-=mod;
}return;
}
int sum(LL n){
if(n<=5000000)return phi[n];
int tmp=H.ask(n);
if(tmp!=-1)return tmp;
LL la,k=n%mod;
int A=k*(k+1)%mod*inv%mod;
for(LL i=2;i<=n;i=la+1){
LL now=n/i;
la=n/now;
A=A-(la-i+1)%mod*sum(now)%mod;
if(A<0)A+=mod;
}H.push(n,A);
return A;
}
int main(){
scanf("%lld",&n);
inv=pow_mod(2,mod-2);
Get_Prime();
LL la,ans=0;
for(LL i=1;i<=n;i=la+1){
LL now=n/i;
la=n/now;
ans=ans+(i+la)%mod*(la-i+1)%mod*inv%mod*sum(now)%mod;
if(ans>=mod)ans-=mod;
}ans<<=1;if(ans>=mod)ans-=mod;
LL k=n%mod;
ans=ans-k*(k+1)%mod*inv%mod;
if(ans<0)ans+=mod;
printf("%lld\n",ans);
return 0;
}
BZOJ 4176
貌似是SDOI2015某道題的加強版？
接下來是一個悲傷的故事，故事的名稱叫做不會用博客園寫公式的代價
首先我們有一個結論：
好吧，湊合着看吧，結論的證明我知道的有三種，其中最簡單的一種是利用數學歸納法
然后式子就變成了(把e函數等價代換成mu函數)
然后我們枚舉k，式子等價轉換成
我們很容易發現中間那兩個sigma連一起的奇怪的東西其實就是這個玩意
好了，終於不用在寫這歪歪扭扭的公式了，結束了QAQ
g函數是約數個數函數，我們求約數個數函數的前綴和可以做到O(sqrt(n))
實際做法是枚舉每個數i，(n/i)就是這個數對答案的貢獻，分塊求和即可
然后我們整體枚舉k，做分塊求和，某一段mu的和我們可以用杜教篩搞定
由於段的起點和終點都是n的約數，所以杜教篩的時間復雜度是O(n^(2/3))的
至於兩次分塊求和疊加看似是sqrt(n)*sqrt(n)=n的，但是根據上面的時間復雜度分析
時間復雜度是O(n^(3/4))的，到此這道題我們就做完了
(我一定要學會如何用blog寫公式(太丑了，捂臉。。挖坑ing))
寫這道題的時候我發現一個很逗的事情
就是我之前的題目的哈希表的push其實不用看之前有沒有push過
直接push進去就好了，因為之前已經ask過發現並沒有放進去了
代碼在BZOJ上跑的還是很快的
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cstdlib>
using namespace std;
const int mod=1e9+7;
const int MOD=1333331;
typedef long long LL;
int n,inv;
int mu[1000010];
bool vis[1000010];
int p[1000010],cnt=0;
struct HASHMAP{
int h[MOD+10],cnt;
int next[100010];
int S[100010],st[100010];
void push(int k,int v){
int key=k%MOD;
++cnt;next[cnt]=h[key];h[key]=cnt;
S[cnt]=k;st[cnt]=v;
}
int ask(int k){
int key=k%MOD;
for(int i=h[key];i;i=next[i]){
if(S[i]==k)return st[i];
}return -1;
}
}H;
int pow_mod(int v,int p){
int tmp=1;
while(p){
if(p&1)tmp=1LL*tmp*v%mod;
v=1LL*v*v%mod;p>>=1;
}return tmp;
}
void Get_Prime(){
mu[1]=1;
for(int i=2;i<=1000000;++i){
if(!vis[i]){
p[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt;++j){
if(1LL*p[j]*i>1000000)break;
int ip=i*p[j];
vis[ip]=true;
if(i%p[j]==0)break;
mu[ip]=-mu[i];
}
}
for(int i=2;i<=1000000;++i)mu[i]+=mu[i-1];
}
int sum(int n){
if(n<=1000000)return mu[n];
int tmp=H.ask(n);
if(tmp!=-1)return tmp;
int la,A=1;
for(int i=2;i<=n;i=la+1){
int now=n/i;
la=n/now;
A=A-1LL*(la-i+1)*sum(now)%mod;
if(A<0)A+=mod;
}H.push(n,A);
return A;
}
int g(int n){
int la,A=0;
for(int i=1;i<=n;i=la+1){
int now=n/i;
la=n/now;
A=A+1LL*(la-i+1)*now%mod;
if(A>=mod)A-=mod;
}return 1LL*A*A%mod;
}
int main(){
scanf("%d",&n);
Get_Prime();inv=pow_mod(2,mod-2);
int la,ans=0;
for(int i=1;i<=n;i=la+1){
int now=n/i;
la=n/now;
ans=ans+1LL*(sum(la)-sum(i-1))*g(now)%mod;
if(ans>=mod)ans-=mod;
if(ans<0)ans+=mod;
}printf("%d\n",(ans%mod+mod)%mod);
return 0;
UPD:最近又做了幾道51Nod的題目，公式推導不知道比上面的長到哪里去了
鑒於本人並不會用編輯器寫公式，所以就只放代碼了
51Nod 最小公倍數之和
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long LL;
const int mod=1e9+7;
const int MOD=1333331;
LL n,inv2,inv6;
bool vis[5000010];
int p[500010],cnt=0;
LL phi[5000010];
struct HASHMAP{
int h[MOD+10],cnt;
int next[100010];
LL S[100010],st[100010];
LL ask(LL k){
int key=k%MOD;
for(int i=h[key];i;i=next[i]){
if(S[i]==k)return st[i];
}return -1;
}
void push(LL k,LL v){
int key=k%MOD;
++cnt;next[cnt]=h[key];h[key]=cnt;
S[cnt]=k;st[cnt]=v;
}
}H;
void Get_Prime(){
phi[1]=1;
for(int i=2;i<=5000000;++i){
if(!vis[i]){
p[++cnt]=i;
phi[i]=i-1;
}
for(int j=1;j<=cnt;++j){
if(1LL*i*p[j]>5000000)break;
int ip=i*p[j];
vis[ip]=true;
if(i%p[j]==0){
phi[ip]=phi[i]*p[j];
break;
}phi[ip]=phi[i]*(p[j]-1);
}
}
for(int i=2;i<=5000000;++i){
phi[i]=phi[i]*i%mod*i+phi[i-1];
phi[i]%=mod;
}return;
}
LL pow_mod(LL v,LL p){
LL tmp=1;
while(p){
if(p&1)tmp=tmp*v%mod;
v=v*v%mod;p>>=1;
}return tmp;
}
LL Get_mod(LL v){return v%mod;}
LL S(LL n){
if(n<=5000000)return phi[n];
LL tmp=H.ask(n);
if(tmp!=-1)return tmp;
LL la,A=Get_mod(n);
A=A*(A+1)%mod*inv2%mod;
A=A*A%mod;
for(LL i=2;i<=n;i=la+1){
LL now=n/i;
la=n/now;
LL t1=Get_mod(la),t2=Get_mod(i-1);
t1=t1*(t1+1)%mod*((t1<<1)|1)%mod*inv6;
t2=t2*(t2+1)%mod*((t2<<1)|1)%mod*inv6;
t1=t1-t2;t1%=mod;if(t1<0)t1+=mod;
A=A-t1*S(now)%mod;
if(A<0)A+=mod;
}H.push(n,A);
return A;
}
int main(){
scanf("%lld",&n);
Get_Prime();
LL la,ans=0;
inv2=pow_mod(2LL,mod-2);
inv6=pow_mod(6LL,mod-2);
for(LL i=1;i<=n;i=la+1){
LL now=n/i;
la=n/now;
ans=ans+Get_mod(i+la)*Get_mod(la-i+1)%mod*inv2%mod*S(now);
ans%=mod;
}printf("%lld\n",ans);
return 0;
}
51Nod 平均最小公倍數
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long LL;
const int mod=1e9+7;
const int MOD=1333331;
int p[100010],cnt=0;
bool vis[1000010];
LL phi[1000010];
LL L,R,inv2,inv6;
struct HASHMAP{
int h[MOD+10],cnt;
int next[100010];
LL S[100010],st[100010];
LL ask(LL k){
int key=k%MOD;
for(int i=h[key];i;i=next[i]){
if(S[i]==k)return st[i];
}return -1;
}
void push(LL k,LL v){
int key=k%MOD;
++cnt;next[cnt]=h[key];h[key]=cnt;
S[cnt]=k;st[cnt]=v;
}
}H;
void Get_Prime(){
phi[1]=1;
for(int i=2;i<=1000000;++i){
if(!vis[i]){
p[++cnt]=i;
phi[i]=i-1;
}
for(int j=1;j<=cnt;++j){
if(1LL*i*p[j]>1000000)break;
int ip=i*p[j];
vis[ip]=true;
if(i%p[j]==0){
phi[ip]=phi[i]*p[j];
break;
}
phi[ip]=phi[i]*(p[j]-1);
}
}
for(int i=2;i<=1000000;++i){
phi[i]=phi[i]*i%mod;
phi[i]+=phi[i-1];
if(phi[i]>=mod)phi[i]-=mod;
}return;
}
LL pow_mod(LL v,LL p){
LL tmp=1;
while(p){
if(p&1)tmp=tmp*v%mod;
v=v*v%mod;p>>=1;
}return tmp;
}
LL sum(LL n){
if(n<=1000000)return phi[n];
LL tmp=H.ask(n);
if(tmp!=-1)return tmp;
LL la,A=n*(n+1)%mod*((n<<1)|1)%mod*inv6%mod;
for(LL i=2;i<=n;i=la+1){
LL now=n/i;
la=n/now;
A=A-(i+la)*(la-i+1)%mod*inv2%mod*sum(now)%mod;
if(A<0)A+=mod;
}H.push(n,A);
return A;
}
LL Get_ans(LL n){
LL la,ans=0;
for(LL i=1;i<=n;i=la+1){
LL now=n/i;
la=n/now;
ans=ans+(la-i+1)*sum(now)%mod;
if(ans>=mod)ans-=mod;
}ans+=n;if(ans>=mod)ans-=mod;
return ans*inv2%mod;
}
int main(){
scanf("%lld%lld",&L,&R);
inv2=pow_mod(2LL,mod-2);
inv6=pow_mod(6LL,mod-2);
Get_Prime();
printf("%lld\n",(Get_ans(R)-Get_ans(L-1)+mod)%mod);
return 0;
}
51Nod 最小公倍數計數
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long LL;
LL L,R;
int mu[1000010];
int p[1000010],cnt=0;
bool vis[1000010];
void Get_Prime(){
mu[1]=1;
for(int i=2;i<=1000000;++i){
if(!vis[i]){
p[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt;++j){
if(1LL*p[j]*i>1000000)break;
int ip=i*p[j];
vis[ip]=true;
if(i%p[j]==0)break;
mu[ip]=-mu[i];
}
}return;
}
LL S(LL n){
LL t1=0,t2=0,t3=0;
LL p1=(int)(pow(n,1.0/3.0)+0.00000001);
LL p2=(int)(sqrt(n));
for(LL i=1;i<=p1;++i){
t3++;
LL lim=n/i;
LL p3=(int)(sqrt(lim));
for(int j=i;j<=p3;++j){
t1=t1+(lim/j)-j+1;
}
}
for(LL i=1;i<=p2;++i)t2=t2+(n/(i*i));
t1=6LL*t1-3LL*t2-2LL*t3;
return t1;
}
LL Solve(LL n){
LL d=1,ans=0,sum,la;
for(LL i=1;i<=n;i=la+1){
LL now=n/i;
la=n/now;sum=0;
while(d*d<=la){sum+=mu[d];d++;}
if(sum)ans=ans+sum*S(now);
}ans+=n;ans>>=1;
return ans;
}
int main(){
scanf("%lld%lld",&L,&R);
Get_Prime();
printf("%lld\n",Solve(R)-Solve(L-1));
return 0;
}
51Nod 約數之和
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cstdlib>
using namespace std;
typedef long long LL;
const int maxn=1000010;
const int mod=1e9+7;
const int MOD=1333331;
LL n,inv;
int mu[maxn];
int p[maxn],cnt=0;
bool vis[maxn];
struct HASHMAP{
int h[MOD+10],cnt;
int next[100010];
LL S[100010],st[100010];
LL ask(LL k){
int key=k%MOD;
for(int i=h[key];i;i=next[i]){
if(S[i]==k)return st[i];
}return -1;
}
void push(LL k,LL v){
int key=k%MOD;
++cnt;next[cnt]=h[key];h[key]=cnt;
st[cnt]=v;S[cnt]=k;
}
}H;
void Get_Prime(){
mu[1]=1;
for(int i=2;i<=1000000;++i){
if(!vis[i]){
p[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt;++j){
if(1LL*p[j]*i>1000000)break;
int ip=i*p[j];
vis[ip]=true;
if(i%p[j]==0)break;
mu[ip]=-mu[i];
}
}
for(int i=2;i<=1000000;++i){
mu[i]=1LL*mu[i]*i%mod;
mu[i]+=mu[i-1];
if(mu[i]>=mod)mu[i]-=mod;
}return;
}
LL pow_mod(LL v,LL p){
LL tmp=1;
while(p){
if(p&1)tmp=tmp*v%mod;
v=v*v%mod;p>>=1;
}return tmp;
}
LL sum(LL n){
if(n<=1000000)return mu[n];
LL tmp=H.ask(n);
if(tmp!=-1)return tmp;
LL la,A=1;
for(LL i=2;i<=n;i=la+1){
LL now=n/i;
la=n/now;
LL t1=(i+la)*(la-i+1)%mod*inv%mod;
A=A-t1*sum(now)%mod;
if(A<0)A+=mod;
}H.push(n,A);
return A;
}
LL f(LL n){
LL la,ans=0;
for(LL i=1;i<=n;i=la+1){
LL now=n/i;
la=n/now;
ans=ans+(la-i+1)*now%mod*(now+1)%mod*inv%mod;
if(ans>=mod)ans-=mod;
}return ans;
}
LL g(LL n){
LL la,ans=0;
for(LL i=1;i<=n;i=la+1){
LL now=n/i;
la=n/now;
ans=ans+(i+la)*(la-i+1)%mod*inv%mod*now%mod;
if(ans>=mod)ans-=mod;
}return ans;
}
int main(){
scanf("%lld",&n);
Get_Prime();
inv=pow_mod(2LL,mod-2);
LL la,ans=0;
for(LL i=1;i<=n;i=la+1){
LL now=n/i;
la=n/now;
LL t1=sum(la)-sum(i-1);
t1%=mod;if(t1<0)t1+=mod;
ans=ans+t1*f(now)%mod*g(now)%mod;
if(ans>=mod)ans-=mod;
}printf("%lld\n",ans);
return 0;
}
Em，總結一下吧
杜教篩具有非常明顯的特征：
1、數據范圍特征:10^9,10^10,10^11
2、能解決的問題是求解某個函數的前綴和或者一段區間的和
但是前提是其和I函數的狄利克雷卷積的前綴和是比較好求的
3、可以在預處理的情況下O(n^(2/3))的時間復雜度內求出O(sqrt(n))個值，也就是n的約數的所有值
因為有記憶化的存在，所以在計算時間復雜度的時候杜教篩不能簡單的跟其他東西相乘來計算
4、一定不要忘記記憶化，一定不要忘記預處理
5、由於數據范圍的特殊性，所以在計算過程中非常容易整數溢出
在寫的時候一定要仔細分析之后在碼，碼完之后要通過對拍或者測極限數據等方法來檢驗
留下的一些坑：
51Nod 加權約數之和
BZOJ 3512