【BZOJ4407】於神之怒加強版
題面
題目分析
\[ \begin{split} \sum\limits_{i=1}^n\sum\limits_{j=1}^mgcd(i,j)^k&=\sum\limits_{d=1}^nd^k\sum\limits_{i=1}^n\sum\limits_{j=1}^m[gcd(i,j)==d]\\end{split} \]
設\(f(x)\)表示\(gcd(i,j)=x\),\(g(x)\)表示\(gcd(i,j)==kx,k\in Z\)。
\[
\begin{split}
g(x)&=\sum\limits_{x|d}^nf(d)\&=\sum\limits_{i=1}^n\sum\limits_{j=1}^m[x|gcd(i,j)]\&=\sum\limits_{i=1}^{\lfloor\frac n x\rfloor}\sum\limits_{j=1}^{\lfloor\frac m x\rfloor}\lfloor\frac n x\rfloor\lfloor\frac m x\rfloor\f(x)&=\sum\limits_{x|d}^n\mu(\frac dx)g(d)=\sum\limits_{x|d}^n\mu(\frac dx)\lfloor\frac n d\rfloor\lfloor\frac m d\rfloor
\end{split}
\]
\[ \begin{split} ans&=\sum\limits_{d=1}^nd^k\cdot f(d)\&=\sum\limits_{d=1}^nd^k\sum\limits_{d|T}^n\mu(\frac Td)\lfloor\frac n T\rfloor\lfloor\frac m T\rfloor\&=\sum\limits_{T=1}^n\lfloor\frac n T\rfloor\lfloor\frac m T\rfloor\sum\limits_{d|T}\mu(\frac Td)d^k \end{split} \]
由於\(\mu\)和\(d^k\)
前面部分用整除分塊加速。
代碼如下
#include<iostream> #include<cstring> #include<cmath> #include<algorithm> #include<cstdio> #include<iomanip> #include<cstdlib> #define MAXN 0x7fffffff typedef long long LL; const int N=5000005,mod=1e9+7; using namespace std; inline int Getint(){register int x=0,f=1;register char ch=getchar();while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}return x*f;} int g[N],mu[N],prime[N]; bool vis[N]; LL ksm(LL x,LL k){ LL ret=1; while(k){ if(k&1)ret=ret*x%mod; x=x*x%mod; k>>=1; } return ret; } int low[N]; int main(){ int T=Getint(),K=Getint(); mu[1]=g[1]=1; for(int i=2;i<=5e6;i++){ if(!vis[i]){ prime[++prime[0]]=i,mu[i]=-1; low[i]=i,g[i]=ksm(i,K)-1; } for(int j=1;j<=prime[0]&&1ll*prime[j]*i<=5e6;j++){ vis[i*prime[j]]=1; if(i%prime[j]==0){ low[i*prime[j]]=low[i]*prime[j]; if(low[i*prime[j]]==i*prime[j]) g[i*prime[j]]=g[i]*ksm(prime[j],K)%mod; else g[i*prime[j]]=(1ll*g[low[i*prime[j]]]*g[i*prime[j]/low[i*prime[j]]])%mod; break; } low[i*prime[j]]=prime[j]; g[i*prime[j]]=(1ll*g[i]*g[prime[j]])%mod; mu[i*prime[j]]=-mu[i]; } } for(int i=1;i<=5e6;i++)g[i]=(g[i]+g[i-1])%mod; while(T--){ int n=Getint(),m=Getint(); if(n>m)swap(n,m); int ans=0; for(int l=1,r;l<=n;l=r+1){ r=min(n/(n/l),m/(m/l)); ans=(ans+1ll*(n/l)*(m/l)%mod*(g[r]-g[l-1])%mod+mod)%mod; } cout<<ans<<'\n'; } return 0; }
【BZOJ4407】於神之怒加強版