【學術篇】分析礦洞 杜教篩
數論什麽的都去死吧！
看著題解我都能化式子用完4頁草紙。。。
另外吐槽一句出題人的拼音學的是真好, 不知道是不是故意的.
其實題解已經寫得挺詳細的了.
我就是提一些出題人覺得太easy沒必要提但是做題還是需要的一些東西....(因為這些東西我基本都是現學的)
然而之前剛剛學完mobius反演就暫時性脫坑的我啥也不會啊..
看到前排dp和曲神在水luogu的歡(bao/du)樂(ling/liu)賽, 就想去看看.
然後就點了報名但是發現自己什麽都不會.
去看了看T1. 就是這道題.
然後成功的化出了第一步的式子.
這樣就可以水30pts了.
一眼看出應該是反演類型的題目, 但是真的tmd不會啊,, 80pts都水不到.
30pts:
通過一眼看出法可以得到激光掃到的第一個點的坐標是
\[(\frac x{gcd(x,y)},\frac y{gcd(x,y)})\]
所以
\[(\frac{x+y}{\frac x{gcd(x,y)}+\frac y{gcd(x,y)}})^2=gcd^2(x,y)\]
於是很顯然地就是要求出
\[\sum_{i=1}^N\sum_{j=1}^N\varphi(gcd^2(i,j))\cdots①\]
這個東西, 於是就變成了一道純數論題. (本來就是一道純數論題不是?!)
然後就\(O(n^2logn)\)
80pts:
繼續化式子.
對於\(\varphi\)我們有這麽一種操作:
\[\varphi(n)=\prod_ip_i^{k_i-1}(p_i-1)\]
所以可以得到
\[\varphi(n^m)=\prod_ip_i^{mk_i-1}(p_i-1)=\prod_ip_i^{k_i-1}(p_i-1)*p_i^{(m-1)k_i}=\varphi(n)*n^{m-1}\]
我們就可以把①中的gcd拆出來, 把平方去掉, 變成
\[①=\sum_{i=1}^{N}\sum_{j=1}^{N}\varphi(gcd(i,j))*gcd(i,j)\]
然後反演的常見套路之枚舉公因數
\[=\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{d|i,d|j}\varphi(d)*d*b[gcd(i,j)=d]\]
其中\(b[x]\)表示\(x\)的真假性, \(x\)真則\(b[x]=1\), 否則\(b[x]=0\)
也就等價於在一個邊長為\(\left \lfloor\frac nd\right \rfloor\)的方陣中找互質的\((i,j)\), 然後對\(d*\varphi(d)\)求和.
\[=\sum_{d=1}^N\varphi(d)*d\sum_{i=1}^{\left \lfloor\frac Nd\right \rfloor}\sum_{j=1}^{\left \lfloor\frac Nd\right \rfloor} b[gcd(i,j)=1]\cdots②\]
後半邊式子有點眼熟?? 我們好像在哪裏見過這種形式.
儀仗隊!!!!!!
我們可以得出
\[\sum_{i=1}^N\sum_{j=1}^Nb[gcd(i,j)=1]\]
這個式子的結果是\(2*\sum_{i=1}^N\varphi(N)-1\).
所以
\[②=\sum_{d=1}^N\varphi(d)*d*(2*\sum_{i=1}^{\left \lfloor\frac Nd\right \rfloor}\varphi(i)-1)\cdots③\]
這樣這個式子就化到頭了.
而此時我們枚舉\(d\)就可以做到\(O(n\sqrt n)\)求解了..
這樣就能水到80pts. (個人感覺80pts部分分給得略高了)
100pts
滿分做法就要用到一種高端科技了..
杜教篩!
顧名思義是一種篩法. 但是要比線篩快一些.
舉個栗子, 我們來求一下
\[\sum_{i=1}^N\varphi(i)\] (其實跟上面的式子是有聯系的OvO)
那麽我們看數據範圍想算法:
- \(N<=1000\)?
- 枚舉, 每次從頭掃求一遍歐拉函數都能過.
- \(N<=10000000\)?
- \(\varphi(x)\)是個積性函數, 直接線篩一下就好咯.
- \(N<=10000000000\)?
- 這個...... \(O(n)\)過不了啊..
這就是說我們必須要想別的方法了.
比如杜教篩..
我們先來化一波式子, 盡量把N變小到能做的範圍.
對於\(\varphi\)函數有一條:
\[\sum_{d|n}\varphi(d)=n\]
那麽
\[\sum_{d|n,d<n}\varphi(d)+\varphi(n)=n\ \varphi(n)=n-\sum_{d|n,d<n}\varphi(d)\]
我們\[\phi(i)= \sum_{i=1}^N\varphi(i)=\sum_{i=1}^N(i-\sum_{d|i,d<i}\varphi(d))=\sum_{i=1}^Ni-\sum_{i=2}^N\sum_{d|i,d<i}\varphi(d) =\sum_{i=1}^Ni-\sum_{\frac id=2}^N\sum_{d=1}^{\left\lfloor\frac n{\frac id}\right\rfloor}\varphi(d)\ 令j=\frac id,則\phi(i)\sum_{i=1}^N\varphi(i)=\sum_{i=1}^Ni-\sum_{j=2}^N\sum_{d=1}^{\left\lfloor\frac nj\right\rfloor}\varphi(d)=\sum_{i=1}^N-\sum_{j=2}^N\phi(\left\lfloor\frac nj\right\rfloor)\]
其中減號前面的顯然是可以\(O(1)\)計算的(別說你不會), 後面的值是不會超過\(\sqrt n\)個的, 我們枚舉因數遞歸計算即可.
代碼太長而且基本是這題代碼的子集就不糊在這裏了..留個傳送門自己去看吧.
然後我們回到這道題. 我們化出了③式, 為了防止忘掉, 我們再貼一遍.
\[\sum_{d=1}^N\varphi(d)*d*(2*\sum_{i=1}^{\left \lfloor\frac Nd\right \rfloor}\varphi(i)-1)\]
這樣括號裏面的剛剛學習了怎麽篩(所以說是子集嘛), 所以問題就集中在了前面的
\[\sum_{d=1}^N\varphi(d)*d\]
怎麽快速的篩出來. 而這個題解已經說的挺清楚了的~~(說你是不是懶得繼續化了←_←)~~
我們令\(f(i)=\sum_{d=1}^N\varphi(d)*d\)
這個\(\varphi(d)*d=\varphi(d^2)\), 我們就猜測和\(\sum_{i=1}^Ni^2\)(這個式子可以用平方和公式\(O(1)\)求喲~)有什麽聯♂系.
(這個理由是蒙的, 比賽的時候怎麽湊知道該怎麽湊出來 不是很急但是會在線等= =)
那我們就化一下
\[\sum_{i=1}^Ni^2=\sum_{i=1}^N(i*\sum_{d|i}\varphi(d))=\sum_{i=1}^N(\frac id*\sum_{d|n}\varphi(d)*d)\]
這樣湊出了\(\varphi(d)*d\)的形式. 但是還沒有\(\sum_{i=1}^n\)的形式, 考慮枚舉.
我們令\(t=\frac id\), 然後枚舉\(t\). 因為要求1..N的和, 所以如果有\(t*d>N\)的\(d\)顯然不能對答案做出貢獻, 所以我們枚舉\(d\)的時候枚舉到\(\left\lfloor\frac Nt\right\rfloor\)即可. 也就是說
\[\sum_{i=1}^Ni^2=\sum_{t=1}^Nt*\sum_{d=1}^{\left\lfloor\frac Nt\right\rfloor}\varphi(d)*d=\sum_{d=1}^N\varphi(d)*d+\sum_{t=2}^Nt*\sum_{d=1}^{\left\lfloor\frac Nt\right\rfloor}\varphi(d)*d\]
那我們就可以得到
\[f(i)=\sum_{i=1}^Ni^2-\sum_{t=2}^Nt*\sum_{d=1}^{\left\lfloor\frac Nt\right\rfloor}\varphi(d)*d=\sum_{i=1}^Ni^2-\sum_{t=2}^Nf(\left\lfloor\frac Nt\right\rfloor)\]
就這麽得到了一個杜教篩的形式, 就可以仿照上面做咯~
md化式子化到吐系列……
聽說這題卡常數? 但也沒怎麽卡嘛 感覺隨便一跑就輕松第一頁了?
(好像出題人改了題, 所以只算改題之後的話應該就rank1了= =之前q1~q4數組開大了memset廢掉好多時間= =
用了一些小trick 比如全程能開int絕對不開long long... (所以因為少%了一次第一次交80pts嚶嚶嚶...
比如把結果記憶化一下. 看到std這一步是用map做的.
但是因為都是\(N\)的因數, 我們就可以把這些因數分比一個數大的和比這個數小的兩類分
這個數大約取\(\sqrt N\)即可. (大點也能存下, 但是小點好像可能會出現沖突)
開int的話一定要記得:步步取模, 強轉long long!
Emmmm還有一種非常無良的針對數據的壓常trick就是n<=1e5的時候線篩可以少篩一點...(這樣就穩穩rank1了)
const SQ=1e6/7; //142857
int p1[SQ],p2[SQ];
int& getaddr(LL x){
if(x<SQ) return p1[x];
return p2[n/x];
} //返回儲存地址這樣就能少個log, 就會非常快...
代碼(為什麽不去看std呢)
#include <cstdio>
#include <cstring>
const int N=5e6+5;
const int P=1e9+7;
const int SQ=1e6/7;
typedef long long LL;
LL n;
int p1[SQ],p2[SQ],p3[SQ],p4[SQ];
int prime[N],euler[N],mu[N],eusum[N],eumul[N],tot;
bool notp[N];
void shai(){
notp[1]=euler[1]=mu[1]=eusum[1]=eumul[1]=1;
for(int i=2;i<=5e6;++i){
if(!notp[i]){
prime[++tot]=i;
euler[i]=i-1;
mu[i]=-1;
}
for(int j=1;j<=tot&&i*prime[j]<5e6;++j){
int k=i*prime[j];
notp[k]=1;
if(i%prime[j]==0){
mu[k]=0;
euler[k]=euler[i]*prime[j];
break;
}
else{
mu[k]=-mu[i];
euler[k]=euler[i]*(prime[j]-1);
}
}
eusum[i]=(eusum[i-1]+euler[i])%P;
eumul[i]=(eumul[i-1]+1LL*euler[i]*i%P)%P;
}
}
inline int qpow(int a,int b,int s=1){
for(;b;b>>=1,a=1LL*a*a%P)
if(b&1) s=1LL*s*a%P;
return s;
}
int inv2=qpow(2,P-2),inv6=qpow(6,P-2);
inline int& getaddr(LL x,bool flag){
if(flag){
if(x<SQ) return p1[x];
return p2[n/x];
}
if(x<SQ) return p3[x];
return p4[n/x];
}
int eulersum(LL x){
if(x<=5e6) return eusum[x];
int& addr=getaddr(x,1);
if(addr!=-1) return addr;
int ans; LL last;
ans=1LL*(x%P)*(x%P+1)%P*inv2%P;
for(LL i=2;i<=x;i=last+1){
last=x/(x/i);
ans=(ans-1LL*(last-i+1)*eulersum(x/i)%P)%P;
}
return addr=(ans%P+P)%P;
}
int eumulsum(LL x){
if(x<=5e6) return eumul[x];
int& addr=getaddr(x,0);
if(addr!=-1) return addr;
int ans; LL last;
ans=1LL*(x%P)*(x%P+1)%P*((x+x+1)%P)%P*inv6%P;
for(LL i=2;i<=x;i=last+1){
last=x/(x/i);
ans=(ans-1LL*(i+last)%P*(last-i+1)%P*inv2%P*eumulsum(x/i)%P)%P;
}
return addr=(ans%P+P)%P;
}
int main(){
memset(p1,-1,sizeof(p1));
memset(p2,-1,sizeof(p2));
memset(p3,-1,sizeof(p3));
memset(p4,-1,sizeof(p4));
shai(); scanf("%lld",&n);
int s=0; LL last;
for(LL i=1;i<=n;i=last+1){
last=n/(n/i);
s=(s+1LL*(eumulsum(last)-eumulsum(i-1)+P)*((2*eulersum(n/i)%P-1+P)%P)%P)%P;
}
printf("%d",s);
}真是要吐了 完結撒花~
【學術篇】分析礦洞 杜教篩