BZOJ 4407: 于神之怒加强版

Posted 2020-09-09 NeighThorn

4407: 于神之怒加强版

Time Limit: 80 Sec  Memory Limit: 512 MB
Submit: 560  Solved: 271
[Submit][Status][Discuss]

Description

给下N,M,K.求

 

Input

输入有多组数据，输入数据的第一行两个正整数T,K，代表有T组数据，K的意义如上所示，下面第二行到第T+1行，每行为两个正整数N,M，其意义如上式所示。

Output

如题

Sample Input

1 2
3 3

Sample Output

20

HINT

1<=N,M,K<=5000000,1<=T<=2000

题解：JudgeOnline/upload/201603/4407.rar

Source

命题人：成都七中张耀楠，鸣谢excited上传。

分析：

定义$n<=m$...

$\sum _{i=1}^{n}\sum _{i=1}^{m} gcd(i,j)^k$

$=\sum _{d=1}^{n} d^k \sum _{x=1}^{n} \mu (x) \left \lfloor \frac{n}{xd} \right \rfloor \left \lfloor \frac{m}{xd} \right \rfloor$

令$y=xd$

$\sum _{d=1}^{n} d^k \sum _{d\mid y} \mu (\frac{y}{d}) \left \lfloor \frac{n}{y} \right \rfloor \left \lfloor \frac{m}{y} \right \rfloor$

$=\sum _{y=1}^{n} \left \lfloor \frac{n}{y} \right \rfloor \left \lfloor \frac{m}{y} \right \rfloor \sum _{d\mid y} \mu( \frac{y}{d} ) d^k$

令$f(n)=\sum _{d\mid n} \mu (\frac{n}{d}) d^k$，可以通过线性筛线性求出$f(x)$...

因为$f(x)$是一个积性函数，所以有以下转化：

$f(n)=\Pi _{i=1}^{t} f(p_i^{x_i})$

$=\Pi _{i=1}^{t} \sum _{d\mid p_i^{x_i}} \mu (\frac{p_i^{x_i}}{d}) d^k$

因为当一个数含有平方因子的时候$\mu (x)=0$，所以可以发现之后下面的两项是有用的...

$=\Pi _{i=1}^{t} \mu (1) p_i^{x_i} +\mu (p_i) p_i^{k(x_i-1)}$

$=\Pi _{i=1}^{t} p_i^{kx_i}-p_i^{k(x_i-1)}$

$=\Pi _{i=1}^{t} p_i^{k(x_i-1)}(p_i^k-1)$

我们发现，当$p$是一个质数的时候，$f(p)=p^k-1$，所以我们可以得到一下等式：

$f[x*p]=f[x]f[p]$   ---   $x$%$p\neq 0$

$f[x*p]=f[x]*p^k$ ---   $x$%$p = 0$

就可以愉快地分块计算了...

代码：

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
//by NeighThorn 
using namespace std;

const int maxn=5000000+5,mod=1e9+7;

int n,m,k,ans,cas,cnt,f[maxn],g[maxn],pri[maxn],vis[maxn];

inline int power(int x,int y){
	int res=1;
	while(y){
		if(y&1)
			res=1LL*res*x%mod;
		x=1LL*x*x%mod,y>>=1;
	}
	return res;
}

signed main(void){
	scanf("%d%d",&cas,&k);g[1]=1;
	for(int i=2;i<=5000000;i++){
		if(!vis[i])
			pri[++cnt]=i,vis[i]=1,f[i]=power(i,k),g[i]=f[i]-1;
		for(int j=1;j<=cnt&&1LL*i*pri[j]<=5000000;j++){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0){
				g[i*pri[j]]=1LL*g[i]*f[pri[j]]%mod;
				break;
			}
			g[i*pri[j]]=1LL*g[i]*g[pri[j]]%mod;
		}
	}
	for(int i=2;i<=5000000;i++)
		g[i]=(g[i-1]+g[i])%mod;
	while(cas--){
		scanf("%d%d",&n,&m);ans=0;
		if(n>m) swap(n,m);
		for(int i=1,r;i<=n;i=r+1){
			r=min(n/(n/i),m/(m/i));
			ans=(ans+1LL*(n/i)*(m/i)%mod*((g[r]-g[i-1]+mod)%mod)%mod)%mod;
		}
		printf("%d\n",ans);
	}
	return 0;
}

　　

---

By NeighThorn