题目地址：BZOJ：Problem 3561. — DZY Loves Math VI

题目描述

给定正整数n，m，求 \displaystyle \sum_{i=1}^n \sum_{i=1}^m \mathrm{lcm}(i, j)^{\mathrm{gcd}(i, j)} 。

输入输出格式

输入格式：
一行两个整数n,m。

输出格式：
一个整数，为答案模1000000007后的值。

输入输出样例

输入样例#1：

5 4

输出样例#1：

424

说明

数据规模：
1<=n,m<=500000，共有3组数据。

题解

参考资料：BZOJ3561 DZY Loves Math VI 数论 快速幂 莫比乌斯反演 – -zhouzhendong- – 博客园
我们知道有\mathrm{lcm}(i, j)\mathrm{gcd}(i, j) = ij，因此题目给的式子可以搞成
\sum_{i=1}^n \sum_{i=1}^m \left( \frac{ij}{\mathrm{gcd}(i, j)} \right)^{\mathrm{gcd}(i, j)}
类似反演经典题目那样，把gcd=k的问题转化成gcd=1的问题，这样就可以直接判定是否互质了
\sum_{d=1}^n \sum_{i=1}^{\left\lfloor \frac{n}{d} \right\rfloor} \sum_{j=1}^{\left\lfloor \frac{m}{d} \right\rfloor} \left( \frac{ijd^2}{d} \right)^d [\mathrm{gcd}(i, j) == 1]
然后你发现这个gcd判定可以用\mu(n)来换掉，有利于之后改变求和顺序
\sum_{d=1}^n \sum_{i=1}^{\left\lfloor \frac{n}{d} \right\rfloor} \sum_{j=1}^{\left\lfloor \frac{m}{d} \right\rfloor} \left( \frac{ijd^2}{d} \right)^d \sum_{p|i, p|j} \mu(p)
接着我们枚举gcd，改变求和顺序，在这个过程中还可以把不变的常数项提到前面来
\sum_{d=1}^n d^d \sum_{p=1}^{\left\lfloor \frac{n}{d} \right\rfloor} \mu(p) \sum_{i=1}^{\left\lfloor \frac{n}{pd} \right\rfloor} \sum_{j=1}^{\left\lfloor \frac{m}{pd} \right\rfloor} (ijp^2)^d
最后的整理
\sum_{d=1}^n d^d \sum_{p=1}^{\left\lfloor \frac{n}{d} \right\rfloor} \mu(p)p^{2d} \sum_{i=1}^{\left\lfloor \frac{n}{pd} \right\rfloor} i^d \sum_{j=1}^{\left\lfloor \frac{m}{pd} \right\rfloor} j^d
乘方求前缀和可以预处理，\mu(n)可以筛出来，d^d快速幂搞一下就好。总复杂度O(n \log n)。

代码

// Code by KSkun, 2018/6
#include <cstdio>
#include <cctype>

#include <algorithm>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2)
        ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1; register char c = fgc();
    for (; !isdigit(c); c = fgc()) if (c == '-') neg = -1;
    for (; isdigit(c); c = fgc()) res = (res << 1) + (res << 3) + c - '0';
    return res * neg;
}

const int MAXN = 500005, MO = 1e9 + 7;

bool notprime[MAXN];
int mu[MAXN], prime[MAXN], ptot;

inline void sieve() {
    notprime[1] = true; mu[1] = 1;
    for(int i = 2; i < MAXN; i++) {
        if(!notprime[i]) {
            prime[++ptot] = i; mu[i] = -1;
        }
        for(int j = 1; j <= ptot && prime[j] * i < MAXN; j++) {
            int k = prime[j] * i;
            notprime[k] = true;
            if(i % prime[j]) {
                mu[k] = -mu[i]; 
            } else {
                mu[k] = 0; break;
            }
        }
    }
}

inline LL fpow(LL x, LL k) {
    LL res = 1;
    for(; k; k >>= 1) {
        if(k & 1) res = res * x % MO;
        x = x * x % MO;
    }
    return res;
}

int n, m;
LL powd[MAXN], psum[MAXN];

int main() {
    sieve();
    n = readint(); m = readint();
    if(n > m) std::swap(n, m);
    for(int i = 1; i <= m; i++) {
        powd[i] = 1;
    }
    LL ans = 0;
    for(int d = 1; d <= n; d++) {
        LL res = 0;
        for(int i = 1; i <= m / d; i++) {
            powd[i] = powd[i] * i % MO; psum[i] = (psum[i - 1] + powd[i]) % MO;
        }
        for(int p = 1; p <= n / d; p++) {
            res = (res + mu[p] * powd[p] % MO * powd[p] % MO 
                * psum[n / p / d] % MO * psum[m / p / d] % MO) % MO;
        }
        res = (res % MO + MO) % MO;
        ans = (ans + res * fpow(d, d) % MO) % MO;
    }
    printf("%lld", ans);
    return 0;
}