bzoj 4407: 於神之怒加強版【莫比烏斯反演+線性篩】
阿新 • • 發佈:2019-04-17
isp space names bre esp clas ios getch [1]
\[ =\sum_{d=1}^{n}d^k\sum_{g=1}^{\lfloor \frac{n}{d} \rfloor }\mu(g)\lfloor\frac{n}{dg}\rfloor\lfloor\frac{m}{dg}\rfloor \]
一般這樣就行了,但是這裏T很大,所以看看有沒有能預處理的東西,枚舉p=dg
\[ =\sum_{d=1}^{n}d^k\sum_{d|p}\mu(\frac{p}{d})\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor \]
\[ =\sum_{d=1}^{n}\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor\sum_{d|p}\mu(\frac{p}{d})d^k \]
前面那段和nm有關,分塊來做;考慮怎麽預處理後面的
顯然是個積性的,所以考慮線性篩出來前綴和即可
看著就像反演,所以先推式子(默認n<m):
\[
\sum_{d=1}^{n}d^k\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)==d]
\]
\[
=\sum_{d=1}^{n}d^k\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor }\sum_{j=1}^{\lfloor \frac{m}{d} \rfloor }[gcd(i,j)==d]
\]
\[
=\sum_{d=1}^{n}d^k\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor }\sum_{j=1}^{\lfloor \frac{m}{d} \rfloor }\sum_{g|i,g|j}\mu(g)
\]
\[ =\sum_{d=1}^{n}d^k\sum_{g=1}^{\lfloor \frac{n}{d} \rfloor }\mu(g)\lfloor\frac{n}{dg}\rfloor\lfloor\frac{m}{dg}\rfloor \]
一般這樣就行了,但是這裏T很大,所以看看有沒有能預處理的東西,枚舉p=dg
\[ =\sum_{d=1}^{n}d^k\sum_{d|p}\mu(\frac{p}{d})\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor \]
\[ =\sum_{d=1}^{n}\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor\sum_{d|p}\mu(\frac{p}{d})d^k \]
前面那段和nm有關,分塊來做;考慮怎麽預處理後面的
顯然是個積性的,所以考慮線性篩出來前綴和即可
#include<iostream> #include<cstdio> using namespace std; const int N=5000005,mod=1e9+7; int T,k,n,m,p[N],tot,s[N],f[N],sm[N],ans; bool v[N]; int read() { int r=0,f=1; char p=getchar(); while(p>'9'||p<'0') { if(p=='-') f=-1; p=getchar(); } while(p>='0'&&p<='9') { r=r*10+p-48; p=getchar(); } return r*f; } int ksm(int a,int b) { int r=1; while(b) { if(b&1) r=1ll*r*a%mod; a=1ll*a*a%mod; b>>=1; } return r; } int main() { T=read(),k=read(); f[1]=1; for(int i=2;i<=5000000;i++) { if(!v[i]) { p[++tot]=i; s[i]=ksm(i,k); f[i]=s[i]-1; } for(int j=1;j<=tot&&i*p[j]<=5000000;j++) { v[i*p[j]]=1; if(i%p[j]==0) { f[i*p[j]]=1ll*f[i]*s[p[j]]%mod; break; } f[i*p[j]]=1ll*f[i]*f[p[j]]%mod; } } for(int i=1;i<=5000000;i++) sm[i]=(sm[i-1]+f[i])%mod; while(T--) { n=read(),m=read(),ans=0; if(n>m) swap(n,m); for(int i=1,la;i<=n;i=la+1) { la=min(n/(n/i),m/(m/i)); ans=(ans+1ll*(n/i)*(m/i)%mod*(sm[la]-sm[i-1])%mod)%mod; } printf("%d\n",(ans+mod)%mod); } return 0; }
bzoj 4407: 於神之怒加強版【莫比烏斯反演+線性篩】