51nod 1575 Gcd and Lcm
題目鏈接:http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1575
?
萬年巨坑終于填掉了……
首先是煞筆西瓜的做題歷程O_O。
原題要求$\sum_{i=1}^n\sum_{j=1}^i\sum_{k=1}^i?[(i,j),(i,k)]$
那么先推一波式子吧
balabala
我也忘記自己是怎么推的了(霧
總之最后推出來是這樣的
$ ans=\sum_{i=1}^{n} f(\left\lfloor\frac{n}{i}\right\rfloor)*g(i) $
其中?$ f(x)=\sum_{i=1}^{x} μ(i)*i^2*\frac{\frac{x}{i}(\frac{x}{i}+1)}{2} $ ?,$ g(x)=[\sum_{d|x} \frac{x}ze8trgl8bvbq*φ(d)]^2 $
然后接下來的問題就是怎么求f(x)的值和g(x)的前綴和了,求出來就能分段計(jì)算辣。
嗯……
想想怎么求……
$ μ(i)*i^2 $ 可以用杜教篩直接求嘛,然后f(x)就可以O(shè)(3/4)分段求,可是在總式里面再跑一次分段的話復(fù)雜度會……誒,不管了,先寫。
g(x)的前綴和可以洲閣篩求,嗯,碼碼碼……(這么復(fù)雜的題怎么才7級?
啊,T了,意料之中……
改改改
跑得越來越快了……
誒,極限數(shù)據(jù)要跑兩秒多,可是10組數(shù)據(jù)的話還是要10多秒
改不動,棄坑……
(51nod群上問了下,夾克老爺說他不會。
問問YJQ
他說$?\sum_{j=1}^i\sum_{k=1}^i?[(i,j),(i,k)] $ 這東西是個(gè)積性函數(shù)。
所以直接用洲閣篩對這個(gè)東西求前綴和就好了(掀桌……
也就是說,看到題目,你開始推式子,你就輸辣。
具體來說,對于一個(gè)質(zhì)數(shù)$ p $,當(dāng) $ i=p^k $ 時(shí),$?\sum_{j=1}^i\sum_{k=1}^i?[(i,j),(i,k)] =(2k+1)*(p^{2k}-p^{2k-1})+p^{k-1} $
優(yōu)越寫法才2.8k,第一種方法直接上5k……
代碼不貼辣。
?
upd at 2017.4.26好像這題爛大街了,我來發(fā)個(gè)洲閣篩模板吧
#include<cstdio> #include<cstring> #include<algorithm> #include<cmath> #define ui unsigned int #define ull unsigned long long #define MN 100001 #define SQ 64000 using namespace std; ui read_p,read_ca; inline ui read(){read_p=0;read_ca=getchar();while(read_ca<'0'||read_ca>'9') read_ca=getchar();while(read_ca>='0'&&read_ca<='9') read_p=read_p*10+read_ca-48,read_ca=getchar();return read_p; } const ui HA=1e6+7; ui T,n,MMH,p[MN/10],_p[MN/10],num=0,TTT,O,Num,f[SQ],_S[SQ],G_2[SQ],G_1[SQ],G_0[SQ],_T[SQ],P_P[MN/10],_G_2[MN/10],_G_1[MN/10],_O[SQ],_Num,_l[HA],NNN=0,LL;bool bo[MN]; ui gcd(ui x,ui y){return y?gcd(y,x%y):x;} ui _b_y[HA],_b_z[HA],_b_ne[HA]; inline void _in(ui x,ui y,ui z){_b_y[++_Num]=y;_b_z[_Num]=z;_b_ne[_Num]=_l[x];_l[x]=_Num;} inline ui Q_2(ui x){if (x%6==1)return (ull)(x+x+1)*(x+1)/6*x;else if (x%6==4) return (ull)(x+x+1)*x/6*(x+1);else return (ull)x*(x+1)/6*(x+x+1);} inline ui Q_1(ui x){return (ull)x*(x+1)>>1;} inline ui sqr(ui x){return x*x;} inline ui min(ui a,ui b){return a<b?a:b;} inline ui find(ui x){if (x<HA) return _b_z[_l[x]];register ui i;for (i=_l[x%HA];i;i=_b_ne[i]) if (_b_y[i]==x) return _b_z[i];return 0; } inline ui Mmh(ui n){register ui i,j,c;ui o=sqrt(n)+1e-9,MMH=0,k=0,R,Ls,Rs,P,SS=0;ull x,Q,O_O;while (p[LL+1]<=o&&LL<num) LL++;for (i=1;i<=n;i=j+1) j=n/(c=n/i),_S[++k]=c,_in(c%HA,c,k),_T[k]=f[k]=0,G_2[k]=Q_2(c),G_1[k]=Q_1(c),G_0[k]=c;for (i=k,j=0;i;_O[i--]=j) while (_S[i]>=p[j+1]) j++;for (f[i=1]=1,P=Ls=Rs=k;i<=LL;i++){if (i==1) while (_S[P]<_p[i]&&P) P--;else P=Rs;while (_S[Rs]<_p[i+1]&&Rs) SS+=f[Rs--];while (_S[Ls]<p[i+1]&&Ls) SS-=f[Ls--];if (i!=LL) MMH+=SS*P_P[i+1];for (j=P;j;j--)if (f[j]){for (x=p[i],Q=1,c=1;x<=_S[j];x*=p[i],Q*=p[i],c++){R=find(_S[j]/x);f[R]+=(O_O=f[j]*((2*c+1)*(p[i]-1)*sqr(Q)*p[i]+Q));if (Ls>=R&&R>Rs) SS+=O_O;if (i!=LL) if (_S[R]>=p[i+1]&&_S[R]<_p[i+1]) MMH+=O_O*P_P[i+1];}}for (j=1;j<=P;j++)if (_S[j]>=p[i]) c=min(i-1,_O[R=find(_S[j]/p[i])]),G_2[j]-=_p[i]*(G_2[R]-(_G_2[c]-_G_2[_T[R]]))+(_G_2[i-1]-_G_2[_T[j]]),G_1[j]-=p[i]*(G_1[R]-(_G_1[c]-_G_1[_T[R]]))+(_G_1[i-1]-_G_1[_T[j]]),G_0[j]-=(G_0[R]-(c-_T[R]))+(i-1-_T[j]),_T[j]=i;}for (j=1;j<=k;j++) c=min(LL,_O[j]),G_2[j]-=_G_2[c]-_G_2[_T[j]],G_1[j]-=_G_1[c]-_G_1[_T[j]],G_0[j]-=c-_T[j];for (i=1;i<=k;i++) MMH+=f[i]*((G_2[i]-G_1[i])*3+G_0[i]),_l[_S[k]%HA]=0;return _Num=0,MMH; } int main(){register ui i,j;for (i=2;i<MN;i++){if (!bo[i]) p[++num]=i,_p[num]=p[num]*p[num],_G_2[num]=_G_2[num-1]+_p[num],_G_1[num]=_G_1[num-1]+p[num],P_P[num]=3*p[num]*(p[num]-1)+1;for (j=1;j<=num&&(O=i*p[j])<MN;j++){bo[O]=1;if (!i%p[j]) break;}}T=read();while(T--){n=read();printf("%u\n",Mmh(n));} } View Code?
轉(zhuǎn)載于:https://www.cnblogs.com/Enceladus/p/5821862.html
總結(jié)
以上是生活随笔為你收集整理的51nod 1575 Gcd and Lcm的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: Spring-----多环境中加载资源配
- 下一篇: 封装javascript分页插件——可以