[BZOJ4407]于神之怒加强版
BZOJ挂了。。。
先把程序放上来,如果A了在写题解吧。
1 #include<cstdio> 2 #include<algorithm> 3 #define N 5000010 4 #define ll long long 5 #define mod (int)(1e9+7) 6 using namespace std; 7 inline int read() 8 { 9 int x=0,f=1;char ch=getchar(); 10 while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();} 11 while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();} 12 return x*f; 13 } 14 inline int ksm(ll x,ll y) 15 { 16 ll p=1; 17 while(y) 18 { 19 if(y&1)p=p*x%mod; 20 x=x*x%mod;y>>=1; 21 } 22 return p; 23 } 24 int prime[N],tot,f[N],p[N],k; 25 bool vis[N]; 26 void pre(int t) 27 { 28 f[1]=1; 29 for(int i=2;i<=t;i++) 30 { 31 if(!vis[i]) 32 { 33 prime[++tot]=i; 34 p[tot]=ksm(i,k); 35 f[i]=p[tot]-1; 36 } 37 for(int j=1;j<=tot&&i*prime[j]<=t;j++) 38 { 39 vis[i*prime[j]]=1; 40 if(!(i%prime[j])) 41 { 42 f[i*prime[j]]=(ll)f[i]*p[j]%mod; 43 break; 44 } 45 f[i*prime[j]]=(ll)f[i]*f[prime[j]]%mod; 46 } 47 } 48 for(int i=1;i<=t;i++) 49 f[i]=(f[i]+f[i-1])%mod; 50 } 51 int n[2005],m[2005],t,maxn; 52 int main() 53 { 54 t=read();k=read(); 55 for(int i=1;i<=t;i++) 56 { 57 n[i]=read(); 58 m[i]=read(); 59 if(n[i]>m[i])swap(n[i],m[i]); 60 maxn=max(maxn,n[i]); 61 } 62 pre(maxn); 63 for(int z=1;z<=t;z++) 64 { 65 int ans=0; 66 for(int i=1,j;i<=n[z];i=j+1) 67 { 68 j=min(n[z]/(n[z]/i),m[z]/(m[z]/i)); 69 ans=(ans+((ll)(f[j]+mod-f[i-1])*(n[z]/i)%mod*(m[z]/i)%mod))%mod; 70 } 71 printf("%d\n",ans); 72 } 73 }
update:
$\sum\limits_{i=1}^n\sum\limits_{j=1}^n\gcd(i,j)^k$
$=\sum\limits_{d=1}^nd^k*\sum\limits_{i=1}^{[\frac nd]}\sum\limits_{j=1}^{[\frac md]}\gcd(i,j)==1$
$=\sum\limits_{d=1}^nd^k*\sum\limits_{i=1}^{[\frac nd]}\sum\limits_{j=1}^{[\frac md]}\sum\limits_{t|i,t|j}\mu(t)$
$=\sum\limits_{d=1}^nd^k*\sum\limits_{t=1}^n\mu(t)*[\frac n{td}]*[\frac m{td}]$
令$x=d*t$
原式=$\sum\limits_{x=1}^n[\frac nx][\frac mx]\sum\limits_{d|x}\mu(\frac xd)*d^k$
后面那个显然是积性函数,可以线性筛。线性筛预处理+分块 时间复杂度$O(n+T*\sqrt n)$
就让我永远不在这里写什么有意义的话--月下孤狼