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