BZOJ4407: 于神之怒加强版

Time Limit: 80 Sec Memory Limit: 512 MB
Submit: 1220 Solved: 550
[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上传。

题解

[sum_{i=1}^{n}sum_{j=1}^{m}gcd(i,j)^k ]

[=sum_{d=1}^{min(n,m)}sum_{i=1}^{n}sum_{j=1}^{m}d^k imes[gcd(i,j)=d] ]

[=sum_{d=1}^{min(n,m)}d^ksum_{i=1}^{lfloor frac{n}{d} floor}sum_{j=1}^{lfloor frac{m}{d} floor}[gcd(i,j)=1] ]

[=sum_{d=1}^{min(n,m)}d^ksum_{i=1}^{lfloor frac{n}{d} floor}sum_{j=1}^{lfloor frac{m}{d} floor}sum_{t|gcd(i,j)}mu(t) ]

[=sum_{d=1}^{min(n,m)}d^ksum_{t=1}^{min(lfloor frac{n}{d} floor,lfloor frac{m}{d} floor)}mu(t)sum_{i=1}^{lfloor frac{n}{d} floor}sum_{j=1}^{lfloor frac{m}{d} floor}sum_{t|i}sum_{t|j}1 ]

[=sum_{d=1}^{min(n,m)}d^ksum_{t=1}^{min(lfloor frac{n}{d} floor,lfloor frac{m}{d} floor)}mu(t)sum_{i=1}^{lfloor frac{n}{d} floor}sum_{t|i}1sum_{j=1}^{lfloor frac{m}{d} floor}sum_{t|j}1 ]

[=sum_{d=1}^{min(n,m)}d^ksum_{t=1}^{min(lfloor frac{n}{d} floor,lfloor frac{m}{d} floor)}mu(t)lfloorfrac{n}{dt} floorlfloorfrac{m}{dt} floor ]

然后呢？

然后复杂度就爆炸了。

其实化式子熟练了一边写一遍就出来了，不用推这么多步

然后有一个有点套路的东西（这个至少我是第三次见了）

令(T=dt)，有

[sum_{i=1}^{n}sum_{j=1}^{m}gcd(i,j)^k ]

[=sum_{T=1}^{min(n,m)}lfloor frac{n}{T} floor lfloorfrac{m}{T} floorsum_{d|T}mu(frac{T}{d})d^k ]

不难发现(g(x)=d^k)是一个积性函数（还是一个完全积性函数）

(sum_{d|T}mu(frac{T}{d})d^k)这东西其实就是(mu)和(g)的卷积

设(f(T)=(mu*g)(T)=sum_{d|T}mu(frac{T}{d})d^k)

这东西显然可以线筛

怎么筛呢？

线筛三步走：

第一步，T为质数时，(f(T)=T^k-1)

第二步，T由他的最小质因子(a)（次数唯一）与另一个数(b)相乘得到时（此时两个数互质），有(f(T)=f(a)f(b))

第三步，T是由他的最小质因子(a)与另一个数(b)相乘得到，且(a|b)时（此时两个数不互质），我们设(ab=a^cq,q=frac{b}{a^{c-1}}),有(f(T)=f(ab)=f(a^cq)=f(a^c)f(q))

考虑从前面递推过来

[f(frac{T}{a})=f(a^{c-1})f(q)=f(q)(a^{k(c-1)}-a^{k(c-2)}) ]

[f(T)=f(a^c)f(q)=f(q)(a^{kc}-a^{k(c-1)}) ]

不难发现

[f(T)=f(frac{T}{a})a^k ]

完工，分块算即可。（上面式子有个地方我写的(-),zyf神犇（%%%%%）的博客里写的是(+)，如有错误还请指正）

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <map>
#include <cmath>
inline long long max(long long a, long long b){return a > b ? a : b;}
inline long long min(long long a, long long b){return a < b ? a : b;}
inline void swap(long long &x, long long &y){long long  tmp = x;x = y;y = tmp;}
inline void read(long long &x)
{
    x = 0;char ch = getchar(), c = ch;
    while(ch < '0' || ch > '9') c = ch, ch = getchar();
    while(ch <= '9' && ch >= '0') x = x * 10 + ch - '0', ch = getchar();
    if(c == '-') x = -x;
}
const long long INF = 0x3f3f3f3f3f3f3f3f;
const long long MAXN = 5000000;
const long long MOD = 1e9 + 7; 

long long pow(long long a, long long b)
{
	long long r = 1, base = a;
	for(;b;b >>= 1)
	{
		if(b & 1) r *= base, r %= MOD;
		base *= base, base %= MOD;
	}
	return r;
}

long long t, k, p[MAXN + 10], bp[MAXN + 10], tot, f[MAXN], n, m;

void make_f()
{
	f[1] = 1;
	for(long long i = 2;i <= MAXN;++ i)
	{
		if(!bp[i]) p[++ tot] = i, f[i] = (pow(i, k) - 1 + MOD) % MOD;
		for(long long j = 1;j <= tot && p[j] * i <= MAXN;++ j)
		{
			bp[i * p[j]] = 1;
			if(i % p[j] == 0)
			{
				f[i * p[j]] = f[i] * pow(p[j], k) % MOD;
				break; 
			}
			f[i * p[j]] = f[i] * f[p[j]] % MOD;
		}
	}
	for(long long i = 1;i <= MAXN;++ i) f[i] += f[i - 1];
} 

int main()
{
	read(t), read(k);
	make_f();
	for(;t;-- t)
	{
		read(n), read(m);
		long long ans = 0, r, mi = min(n, m);
		for(long long T = 1;T <= mi;++ T)
		{
			r = min(n/(n/T), min(m/(m/T), mi));
			ans += ((n/T) * (m/T)) % MOD * (f[r] - f[T - 1]) % MOD;
			if(ans >= MOD) ans -= MOD;
			T = r;
		}
		printf("%lld
", ans);
	}
 	return 0;
}