[SMOJ2004]^=^笑脸

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

如果对一些小范围的数据，将所有情况输出，可以提出这样一个猜想：任意的 ab=cd 都可以化成 (xy)b=(xz)d 的形式。
下面就来证明这个猜想。
考虑我们对 ab 进行质因数分解，得到形如 (d1)k1×(d2)k2×…×(dm)km 的结果，则显而易见 cd 质因数分解后得到的结果和 ab 相同。
根据幂的乘方的逆运算法则，可以将上述结果化为 [(d1)k1÷b]b×[(d2)k2÷b]b×…×[(dm)km÷b]b 的形式，由此不难看出任意的 ki 都为 b 的倍数。
把公共指数 b 提出来，就得到

a=(d1)k1÷b×(d2)k2÷b×…×(dm)km÷b

类似地，有

c=(d1)k1÷d×(d2)k2÷d×…×(dm)km÷d

于是我们设 x=∏di，显然可以将 a 和 c 分别化成 xy 和 xz 的形式。

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

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

b=g×lcm(y,z)y

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

⎧⎩⎨g×lcm(y,z)y≤ng×lcm(y,z)z≤n

变形得到

⎧⎩⎨g≤nylcm(y,z)g≤nzlcm(y,z)

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

g≤min⎧⎩⎨nylcm(y,z)nzlcm(y,z)

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

g≤n×min{y,z}lcm(y,z)

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

g=⌊n×min{y,z}lcm(y,z)⌋

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

参考代码：

//2004.cpp#include <algorithm>#include <cmath>#include <cstdio>#include <cstdlib>#include <cstring>#include <iostream>using namespace std;#define aspe MODconst 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;}