前言

萬年不寫公開博客了，這次填個坑

題目相關(guān)

鏈接

題目大意

求$$\sum _{i=1}^N\sum _{j=1}^{N}sgcd(i,j{)}^k(\mathrm{m}\mathrm{o}\mathrm{d}{2}^{32})$$
sgcd(i,j)為i,j的第二大公約數(shù)
特殊的，若gcd(i,j)=1,則sgcd(i,j)=0

數(shù)據(jù)范圍

$n\le 10^9,k\le 50$

題解

設(shè)m(x)為x的次大因數(shù)，特殊的，m(1)=0
$$\begin{array}{rlr}\sum _{i=1}^N\sum _{j=1}^{N}sgcd(i,j{)}^k& =\sum _{i=1}^N\sum _{j=1}^{N}m(gcd(i,j){)}^k& 根據(jù)定義\\ & =\sum _{g=1}^{N}m(g{)}^k\sum _{i=1}^{?\frac{N}{g}?}\sum _{j=1}^{?\frac{N}{g}?}[gcd(i,j)==1]& 提出gcd\\ & =\sum _{g=1}^{N}m(g{)}^k?2?\sum _{i=1}^{?\frac{N}{g}?}\phi (i)??1?& 大小互換也可以所以要乘以2，并且(1,1)只能算一次所以減一\\ & =\sum _{t=1}^n\sum _{?\frac{N}{g}?=t}m(g{)}^k\left(2\left(\sum _{i=1}^t\phi (i)\right)?1\right)& \end{array}$$
我們發(fā)現(xiàn)，滿足條件的t只有$\mathcal{O}(\sqrt{n})$個，而最后一項可以通過杜教篩來做
容易發(fā)現(xiàn)，對于同一個t，其滿足條件的g是連續(xù)的一串，那么只要對于所有滿足條件的t求${\sum }_{g=1}^{maxg}m(g{)}^k$，這個用min_25篩算
因為本質(zhì)上min25篩會枚舉最小因數(shù)，所以把最小因數(shù)排除掉即可

在預(yù)處理上述運算的時候要求${\sum }_{g=1}^n{g}^k$
我們引入斯特林數(shù)
$$x^n=\sum _{i=0}^x\left\{\begin{array}{c}n\\ i\end{array}\right\}{x}^{i\downarrow }$$
證明可以點進那篇博客看
這樣的話，求sigma就是
$$\begin{array}{rl}\sum _{x=1}^n{x}^k& =\sum _{x=1}^n\sum _{i=0}^x\left\{\begin{array}{c}k\\ i\end{array}\right\}{x}^{i\downarrow }\\ & =\sum _{i=0}^n\left\{\begin{array}{c}k\\ i\end{array}\right\}\sum _{x=i}^n{x}^{i\downarrow }\\ & =\sum _{i=0}^n\left\{\begin{array}{c}k\\ i\end{array}\right\}\sum _{x=i}^n\left(\genfrac{}{}{0ex}{}{x}{i}\right)i!\\ & =\sum _{i=0}^{n}i!\left\{\begin{array}{c}k\\ i\end{array}\right\}\sum _{x=i}^n\left(\genfrac{}{}{0ex}{}{x}{i}\right)\\ & =\sum _{i=0}^{n}i!\left\{\begin{array}{c}k\\ i\end{array}\right\}\left(\genfrac{}{}{0ex}{}{n+1}{i+1}\right)\\ & =\sum _{i=0}^{n}i!\left\{\begin{array}{c}k\\ i\end{array}\right\}\frac{(n+1)!}{(i+1)!(n?i)!}\\ & =\sum _{i=0}^n\left\{\begin{array}{c}k\\ i\end{array}\right\}\frac{(n+1)!}{(i+1)(n?i)!}\\ & =\sum _{i=0}^n\left\{\begin{array}{c}k\\ i\end{array}\right\}\frac{(n+1{)}^{i+1\downarrow }}{i+1}\end{array}$$
此處復(fù)雜度$\mathcal{O}(k^2)$
總復(fù)雜度$\mathcal{O}(\frac{n^{\frac{3}{4}}}{logn}+\sqrt{n}{k}^2+n^{\frac{2}{3}})$

#include<cstdio> #include<cctype> #include<cstring> #include<algorithm> #include<cmath> #include<map> #define rg register typedef long long LL; typedef unsigned int ULL; template <typename T> inline T max(const T a,const T b){return a>b?a:b;} template <typename T> inline T min(const T a,const T b){return a<b?a:b;} template <typename T> inline void mind(T&a,const T b){a=a<b?a:b;} template <typename T> inline void maxd(T&a,const T b){a=a>b?a:b;} template <typename T> inline T abs(const T a){return a>0?a:-a;} template <typename T> inline void swap(T&a,T&b){T c=a;a=b;b=c;} template <typename T> inline T gcd(const T a,const T b){if(!b)return a;return gcd(b,a%b);} template <typename T> inline T lcm(const T a,const T b){return a/gcd(a,b)*b;} template <typename T> inline T square(const T x){return x*x;}; template <typename T> inline void read(T&x) {char cu=getchar();x=0;bool fla=0;while(!isdigit(cu)){if(cu=='-')fla=1;cu=getchar();}while(isdigit(cu))x=x*10+cu-'0',cu=getchar();if(fla)x=-x; } template <typename T> inline void printe(const T x) {if(x>=10)printe(x/10);putchar(x%10+'0'); } template <typename T> inline void print(const T x) {if(x<0)putchar('-'),printe(-x);else printe(x); } namespace du {const int MAX=2000000;bool isprime[MAX+1];int n,prime[MAX+1],primesize;LL phi[MAX+1],phi1[MAX+1];inline void getlist(const LL listsize){memset(isprime,1,sizeof(isprime));isprime[1]=0,phi[1]=1;;for(rg int i=2;i<=listsize;i++){if(isprime[i])prime[++primesize]=i,phi[i]=i-1;for(rg int j=1;j<=primesize&&(LL)i*prime[j]<=listsize;j++){isprime[i*prime[j]]=0;if(i%prime[j]==0){phi[i*prime[j]]=phi[i]*prime[j];break;}phi[i*prime[j]]=phi[i]*(prime[j]-1);}}}inline LL Val(const int x){if(x<=MAX)return phi[x];return phi1[n/x];}void pre(int ttt){n=ttt;getlist(MAX);for(rg int i=1;i<=MAX;i++)phi[i]+=phi[i-1];for(rg int i=min((int)sqrt(n),n/(MAX+1));i>=1;i--){const int u=n/i;LL res=(LL)u*(u+1)/2;int limit=sqrt(u);for(rg int j=1;j<=limit;j++)res-=Val(u/j);limit++;for(rg int j=1;limit*j<=u;j++)res-=phi[j]*((u/j)-(u/(j+1)));phi1[i]=res;}} } int n,k,part,tot,qz[700001],prime[700001],sq[700001],primesize,sum0[700001],id[700001]; inline int ck(const int x){return x<=part?x:tot-n/x+1;} ULL S[51][51],sum2[700001],sum3[700001]; LL sum1[700001]; void pre() {S[0][0]=1;for(rg int i=1;i<=k;i++)for(rg int j=1;j<=i;j++)S[i][j]=S[i-1][j-1]+j*S[i-1][j]; } inline ULL qzktime(int x) {ULL ans=0,dq=0;for(rg int i=0;i<=k;i++)for(rg int j=x-i+1;j<=x+1;j++)if(j%(i+1)==0){dq=S[k][i];for(rg int l=x-i+1;l<=x+1;l++)if(j==l)dq*=l/(i+1);else dq*=l;ans+=dq;break;}return ans; } int val[700001]; ULL R[700001],ans,QZ[700001]; ULL pow(ULL x,int y) {ULL res=1;for(;y;y>>=1,x=x*x)if(y&1)res*=x;return res; } int main() {read(n),read(k);du::pre(n),pre();part=sqrt(n);tot=primesize=0;for(rg int i=1;i<=part;i++)id[++tot]=i,sum0[tot]=i-1,sum1[tot]=((((LL)i+1)*i)>>1)-1,sum2[tot]=qzktime(i)-1,sum3[tot]=i-1;id[++tot]=n/part;if(id[tot]!=part)sum0[tot]=id[tot]-1,sum1[tot]=((LL)(((LL)id[tot]+1)*id[tot])>>1)-1,sum2[tot]=qzktime(id[tot])-1,sum3[tot]=id[tot]-1;else tot--;for(rg int i=part-1;i>=1;i--)id[++tot]=n/i,sum0[tot]=id[tot]-1,sum1[tot]=((((LL)id[tot]+1)*id[tot])>>1)-1,sum2[tot]=qzktime(id[tot])-1,sum3[tot]=id[tot]-1;for(rg int i=2;i<=part;i++)if(sum0[i]!=sum0[i-1]){const int limit=i*i;ULL I=pow(i,k);for(rg int j=tot;id[j]>=limit;j--){const int t=ck(id[j]/i);sum3[j]+=sum2[t]-QZ[primesize]-sum0[t]+primesize;sum2[j]-=(sum2[t]-QZ[primesize])*I;sum1[j]-=(sum1[t]-qz[primesize])*i;sum0[j]-=sum0[t]-primesize;}prime[++primesize]=i,sq[primesize]=i*i,qz[primesize]=qz[primesize-1]+i,QZ[primesize]=QZ[primesize-1]+I;}for(rg int i=1;i<=tot;i++)sum1[i]-=sum0[i];//printf("%lld\n",sum1[i]);//sum1為phi的前綴和for(rg int i=1;i<=part;i++)val[i]=n/i,R[i]=sum3[ck(val[i])];//printf("%d %u %lld\n",val[i],R[i],sum1[i]+1);for(rg int i=part+1;i<=tot;i++)val[i]=tot-i+1,R[i]=sum3[ck(val[i])];//printf("%d %u %lld\n",val[i],R[i],sum1[i]+1);for(rg int i=1;i<=part;i++)ans+=(R[i]-R[i+1])*(2*(ULL)du::Val(i)-1);//printf("%lld\n",sum1[i]);for(rg int i=part+1;i<=tot;i++)ans+=(R[i]-R[i+1])*(2*(ULL)du::Val(n/(tot-i+1))-1);print(ans),putchar('\n');return 0; }

后記

本題可以用min_25而不用杜教篩
但我太菜了沒有學(xué)過非遞歸min_25
QAQ

總結(jié)

以上是生活随笔為你收集整理的[51nod1847][算法马拉松23（飞越愚人节）F]奇怪的数学题的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。