[SMOJ2004]^=^笑脸

首先一看题目就是个数学题，而且数据范围比较大，需要 $O(n)$ 或 $O(n\log_2n)$ 级别的算法才能通过。直接从枚举相同的幂的角度来入手，显然是不可行的。
我们需要换个角度，同样是枚举，怎样枚举才能快速考虑所有可能的情况呢？这就需要我们通过观察或分析，得到 $a^b=c^d$ 的一些性质。

如果对一些小范围的数据，将所有情况输出，可以提出这样一个猜想：任意的 $a^b=c^d$ 都可以化成 $(x^y)^b=(x^z)^d$ 的形式。
下面就来证明这个猜想。
考虑我们对 $a^b$ 进行质因数分解，得到形如 $(d_1)^{k_1}\times(d_2)^{k_2}\times\ldots\times(d_m)^{k_m}$ 的结果，则显而易见 $c^d$ 质因数分解后得到的结果和 $a^b$ 相同。
根据幂的乘方的逆运算法则，可以将上述结果化为 $[(d_1)^{k_1\div b}]^b\times[(d_2)^{k_2\div b}]^b\times\ldots\times[(d_m)^{k_m\div b}]^b$ 的形式，由此不难看出任意的 $k_i$ 都为 $b$ 的倍数。
把公共指数 $b$ 提出来，就得到

$$a=(d_1)^{k_1\div b}\times(d_2)^{k_2\div b}\times\ldots\times(d_m)^{k_m\div b}$$

类似地，有

$$c=(d_1)^{k_1\div d}\times(d_2)^{k_2\div d}\times\ldots\times(d_m)^{k_m\div d}$$

于是我们设 $x=\prod d_i$，显然可以将 $a$ 和 $c$ 分别化成 $x^y$ 和 $x^z$ 的形式。

明确了这条基本事实之后，就可以枚举 $x$、$y$ 和 $z$。由于 的增长速度是指数级别的，且按照题目要求，其大小要控制在 $n$ 以内，因此实际上枚举得很少。
现在问题是，对于已知的 $x^y$ 和 $x^z$，存在多少对 $(b, d)$ 使得 $(x^y)^b=(x^z)^d$。
可以变形为 $x^{yb}=x^{zd}$，亦即 $yb=zd=f$。显然最小的 $f$ 就是 $\text{lcm}(y,z)$。而 $2f$，$3f$…… 也是能够满足的。最多能有多少对呢？

这时应想想题目的限制条件。在 $a^b=c^d$ 中，要求 $a, b, c, d$ 均不超过 $n$。对 $a$ 和 $c$ 的限制可以在枚举 $x, y$, 和 $z$ 的过程中实现。
而对 $b, d$ 的限制，就是这里求出合法对数的关键。不妨设有 $yb=dz=g\times\text{lcm}(y,z)$，则可以移项得到

$$b=\frac{g\times\text{lcm}(y,z)}{y}$$

$d$ 也是类似的。题目要求满足 $b$ 和 $d$ 都小于 $n$，即

$$\large \begin{cases}\frac{g\times\text{lcm}(y,z)}{y}\leq n\\ \frac{g\times\text{lcm}(y,z)}{z}\leq n\end{cases}$$

变形得到

$$\large \begin{cases}g\leq \frac{ny}{\text{lcm}(y,z)}\\ g\leq \frac{nz}{\text{lcm}(y,z)}\end{cases}$$

根据不等式组“同小取小”的原则，可得

$$\large g\leq\min\begin{cases}\frac{ny}{\text{lcm}(y,z)}\\\frac{nz}{\text{lcm}(y,z)}\end{cases}$$

不难发现两个式子的分母其实是一样的，则只需考虑分子当中的 $y$ 和 $z$。也就是

$$g\leq\frac{n\times\min\{y,z\}}{\text{lcm}(y,z)}$$

答案必须是整数，而且应该向下取整。这样我们就得到了，对于确定的 $x$、$y$ 和 $z$，满足 $yb=dz=g\times\text{lcm}(y,z)<n$ 的 $g$ 的数量：

$$g=\left \lfloor \frac{n\times\min\{y,z\}}{\text{lcm}(y,z)} \right \rfloor$$

这样就完美解决了问题。但是要注意 $x=1$ 是比较特殊的情况，应该一开始就单独考虑，之后令 从 2 开始枚举。

参考代码：

//2004.cpp
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>

using namespace std;

#define aspe MOD

const long long MOD = 1e9 + 7;
const long long MAXN = 1e6 + 1;

inline long long gcd(long long x, long long y) { return y ? gcd(y, x % y) : x; }

bool flag[MAXN];

int main(void) {
    freopen("2004.in", "r", stdin);
    freopen("2004.out", "w", stdout);
    for (long long i = 2; i * i <= MAXN; i++)
        for (long long j = i * i; j <= MAXN; j *= i)
            flag[j] = true;
    long long T; cin >> T; //printf("%lld\n", T);
    while (T--) {
        long long N; cin >> N; //cout << N << endl; //printf("%lld\n", N);
        long long ans = (long long)N * N % aspe;
        for (long long k = 2; k <= N; k++) {
            if (flag[k]) continue; //cout << k << endl;
            for (long long x = 1, kx = k; kx <= N; x++, kx *= k)
                for (long long y = 1, ky = k; ky <= N; y++, ky *= k) {
                    (ans += min(x, y) * N / (x * y % aspe / gcd(x, y))) %= aspe;
//                  printf("k = %d, x = %d, y = %d, ans = %d\n", k, x, y, ans);
                }
        }
        cout << ans << endl;
    }
    return 0;
}