483 字
2 分钟
BZOJ 4407 于神之怒加强版

【题目描述】#

给定n,m,k,计算 i=1nj=1mgcd(i,j)k\sum_{i=1}^{n}{\sum_{j=1}^{m}{gcd(i,j)^k}} 对1000000007取模的结果

【输入格式】#

多组数据。 第一行是两个数T,K; 之后的T行,每行两个整数n,m;

【输出格式】#

K行,每行一个结果

【样例输入】#

1 2
3 3

【样例输出】#

20

【提示】#

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

题解#

第一步推式子

i=1nj=1mgcd(i,j)k\sum_{i=1}^{n}{\sum_{j=1}^{m}{gcd(i,j)^k}}
=d=1ndki=1nj=1m[gcd(i,j)==d]=\sum_{d=1}^{n}{d^k\sum_{i=1}^{n}{\sum_{j=1}^{m}{[gcd(i,j)==d]}}}
=d=1ndki=1ndj=1md[gcd(i,j)==1]=\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)==1]}}}
=d=1ndki=1ndj=1mdci,cjμc=\sum_{d=1}^{n}{d^k\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}{\sum_{j=1}^{\lfloor \frac{m}{d} \rfloor}{\sum_{c|i,c|j}{\mu{c}}}}}
T=dcT = dc
=T=1ndTdkμTdnTmT=\sum_{T=1}^{n}{\sum_{d|T}{d^k \mu{\frac{T}{d}} \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor}}
=T=1nnTmTdTdkμTd=\sum_{T=1}^{n}{\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor \sum_{d|T}{d^k \mu{\frac{T}{d}}}}

然后我们需要线性筛出dTdkμTd\sum_{d|T}{d^k \mu{\frac{T}{d}}}就可以了
g(T)=dTdkμTdg(T) = \sum_{d|T}{d^k \mu{\frac{T}{d}}}
显然是积性函数 别问我怎么知道的
当互质时 g(TP)=g(T)g(p)g(T*P) = g(T) * g(p)
当不互质时 又要推狮子
根据μ\mu函数的定义
当且仅当不含平方因子时不为零
所以质因子只有选和不选两种状态,有用的状态数是不变的
每一个有用的μ\mu 都乘了一个p^k; 所以g(Tp)=g(T)pkg(T*p) = g(T) * p^k

#define LL long long
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
int prime[500005], tot;
const int MOD = 1000000007;
bool isnprime[5000005];
LL g[5000005];
LL primeK[500005], n, k, m;
LL N = 5000000;
LL pow_mod(LL a, LL b)
{
LL ans = 1;
while (b)
{
if (b & 1)
ans = ans * a % MOD;
b >>= 1;
a = a * a % MOD;
}
return ans;
}
void Get_g()
{
g[1] = 1;
for (int i = 2; i <= N; i++)
{
if (!isnprime[i])
{
prime[++tot] = i;
g[i] = (primeK[tot] = pow_mod(i, k)) - 1;
}
for (int j = 1; j <= tot; j++)
{
if (i * prime[j] > N)
break;
isnprime[i * prime[j]] = 1;
if (i % prime[j] == 0)
{
g[i * prime[j]] = g[i] * primeK[j] % MOD;
break;
}
g[i * prime[j]] = g[i] * g[prime[j]] % MOD;
}
}
for (int i = 2; i <= N; i++)
g[i] += g[i - 1], g[i] %= MOD;
}
int main()
{
freopen("bzoj_4407.in","r",stdin);
freopen("bzoj_4407.out","w",stdout);
int t;
scanf("%d%lld", &t, &k);
Get_g();
while (t--)
{
scanf("%lld%lld", &n, &m);
int last;
LL ret = 0;
if (n > m)
swap(n, m);
for (int i = 1; i <= n; i = last + 1)
{
last = min(n / (n / i), m / (m / i));
(ret += (n / i) * (m / i) % MOD * (g[last] - g[i - 1] + MOD) % MOD)%=MOD;
}
printf("%lld\n", ret);
}
}
BZOJ 4407 于神之怒加强版
https://www.nekomio.com/posts/93/
作者
NekoMio
发布于
2017-08-13
许可协议
CC BY-NC-SA 4.0