仙气飘飘莫反题。
显现出了很多推式子习惯的问题。
思路
莫反 + 伯努利数 + Pr。
首先根据 \(\operatorname{lcm}(x, y) = \frac{xy}{\gcd(x, y)}\) 可以化简：\(\sum\limits_{i = 1}^n \gcd(i, n)^x (\frac{in}{\gcd(i, n)})^y\).
原式化简得：\(\sum\limits_{i = 1}^n \gcd(i, n)^{x - y} i^y n^y\).
- 与 \(\sum\) 无关的循环变量习惯性提到前面，作为常数处理。
于是原式等于 \(n^y \sum\limits_{i = 1}^n \gcd(i, n)^{x - y} i^y\).
- 习惯性枚举 \(\gcd\).
于是原式等于 \(n^y \sum\limits_{d = 1}^n d^{x - y} \sum\limits_{d \mid i} [\gcd(i, n) = d] i^y\).
等价于 \(n^y \sum\limits_{d = 1}^n d^{x - y} \sum\limits_{i = 1}^{\lfloor \frac{n}{d} \rfloor} [\gcd(i, \frac{n}{d}) = 1] (id)^y\).
化简得 \(n^y \sum\limits_{d = 1}^n d^x \sum\limits_{i = 1}^{\lfloor \frac{n}{d} \rfloor} [\gcd(i, \frac{n}{d}) = 1] i^y\).
-
形似 \(\gcd(x, y) = 1\) 的式子考虑 \(\mu\) 反演。
-
\(k \mid \gcd(x, y)\) 习惯性写成等价形式 \(k \mid x, k \mid y\).
套用莫反套路得 \(n^y \sum\limits_{d = 1}^n d^x \sum\limits_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \sum\limits_{k \mid i, k \mid \frac{n}{d}} \mu(k) i^y\).
- 合理将 \(\sum\) 以及与其有关的项整合在一起，必要时交换求和顺序。
整理得 \(n^y \sum\limits_{d = 1}^n d^x \sum\limits_{k \mid \frac{n}{d}} \mu(k) \sum\limits_{k \mid i} i^y\).
-
\(k \mid i\) 形式的式子考虑枚举 \(\frac{i}{k}\).
并且划分后要注意原本与 \(i\) 有关的项。
再次整理得 \(n^y \sum\limits_{d = 1}^n d^x \sum\limits_{k \mid \frac{n}{d}} \mu(k) k^y \sum\limits_{i = 1}^{\lfloor \frac{n}{kd} \rfloor} i^y\).
- 套路处理之外的部分考虑设出函数。
令 \(S(n, k) = \sum\limits_{i = 1}^n i^k\)，则原式等价于 \(n^y \sum\limits_{d = 1}^n d^x \sum\limits_{k \mid \frac{n}{d}} \mu(k) k^y S(\lfloor \frac{n}{kd} \rfloor, y)\).
有一个经典的拉插题结论：\(\sum\limits_{i = 1}^n i^k\) 是关于 \(n\) 的 \(k + 1\) 次多项式。
所以可以表示成 \(n^y \sum\limits_{i = 0}^{y + 1} f_i \sum\limits_{d = 1}^n d^x \sum\limits_{k \mid \frac{n}{d}} \mu(k) k^y (\frac{n}{kd})^i\).
-
考虑线筛 / PR 快速求满足积性的式子。
通常考虑其在质数的整次幂处的取值 以及 质因子对其的影响。
注意到后面的式子形似狄利克雷卷积，并且每一项都是积性函数，所以可以考虑用 PR 做。
因为含有质数的平方次幂必然有 \(\mu = 0\)，所以只需要考虑在 \(1\) 和质数的贡献。
设结果为 \(F(p^k) = \sum\limits_{j \mid p^k} \mu(j) j^y \sum\limits_{d \mid \frac{p^k}{j}} d^x (\frac{p^k}{jd})^i\)，分讨一下。
-
\(x = 1\)
\(F(p^k) = \sum\limits_{d \mid p^k} d^x (\frac{p^k}{d})^i\).
化简得 \(p^{ik} \sum\limits_{t = 0}^k p^{t(x - i)}\).
也就是 \(p^{ik} \sum\limits_{t = 0}^k p^{tx + (k - t)i}\)，躲掉逆元.
-
\(x = p\)
\(F(p^k) = -p^y \sum\limits_{d \mid p^{k - 1}} d^x (\frac{p^{k - 1}}{d})^i = -p^y F(p^{k - 1})\).
也就是说算 \(y\) 次 \(F\) 就行。
算幂次直接大力快速幂，信仰跑过。
考虑求 \(S(n, k)\).
但是这里拉插复杂度是假的，考虑使用伯努利数（伯努利数生成函数的第 \(k\) 项系数等于 \(\sum\limits_{i = 1}^n i^k\)）
考虑构造伯努利数的生成函数。
首先知道 \(\{1, c^1, c^2, \cdots\}\) 的 EGF 为 \(G_c(x) = \sum\limits_{i = 1}^{+ \infty} \frac{(cx)^i}{i!} = e^{cx}\).
考虑对 \(G_c(x)\) 求和：令 \(S_n(x) = \sum\limits_{c = 0}^{n - 1} G_c(x) = \sum\limits_{c = 0}^{n - 1} e^{cx} = \frac{e^{nx} - 1}{e^x - 1}\).
化简 \(S_n(x)\) 得：\(S_n(x) = \sum\limits_{k = 1}^{+ \infty} \frac{(\sum\limits_{i = 0}^{n - 1} c^k) x^k}{k!}\).
也就是 \(S_n(x) = \frac{1}{k!} \sum\limits_{i = 0}^{n - 1} i^k\).
分离出 \(B(x) = \frac{x}{e^x - 1}\).
所以 \(S_n(x) = B(x) \frac{e^{nx} - 1}{x}\).
设多项式 \(G(x) = \frac{x}{e^x - 1} = \sum\limits_{i = 0}^ \frac{n^{i + 1} x^i}{(i + 1)!}\).
则 \(\sum\limits_{i = 0}^{n - 1} i^k = k! S_n[k] = k! \sum\limits_{i + j = k} G[i] B[j] = k! \sum\limits_{i = 0}^k \frac{B[k - i] n^{i + 1}}{(i + 1)!}\).
于是得到了自然数幂和的多项式。
数据范围大的时候可以 poly 求，当然这里 \(x, y \leq 3 \times 10^3\) 可以直接大力算。
具体复杂度有 Pr 不太会算，结论是 \(O(T (\sqrt{4}{n} + y \log^2 n) + y^2)\).
代码
#include <cstdio>
#include <ctime>
#include <cstdlib>
#include <algorithm>
using namespace std;
typedef long long ll;
const int maxn = 3e3 + 10;
const int mod = 1e9 + 7;
int t, x, y, tot;
ll n;
ll fac[200], fpw[200], num[200], c[maxn];
const int lim = 200;
const int tpr[8] = {2, 3, 5, 7, 13, 19, 61, 233};
inline ll qpow(ll base, ll power = mod - 2)
{
ll res = 1;
base %= mod;
while (power)
{
if (power & 1) res = res * base % mod;
base = base * base % mod, power >>= 1;
}
return res;
}
namespace BNL //Bernoulli
{
ll B[maxn], fac[maxn], invf[maxn];
inline void init()
{
fac[0] = invf[0] = 1;
for (int i = 1; i <= 3e3 + 5; i++) fac[i] = fac[i - 1] * i % mod, invf[i] = invf[i - 1] * qpow(i) % mod;
B[0] = 1;
for (int i = 1; i <= 3e3 + 5; i++)
{
for (int j = 0; j < i; j++) B[i] = (B[i] + B[j] * invf[i - j + 1]) % mod;
B[i] = mod - B[i];
}
}
inline void getf(ll *f, int k)
{
for (int i = 1; i <= k + 1; i++) f[i] = B[k - i + 1] * invf[i] % mod * fac[k] % mod;
f[0] = 0;
if (k) f[k]++;
}
}
inline ll mul(ll a, ll b, ll m)
{
ll d = ((long double)a / m * b + 0.5), r = a * b - d * m;
return (r < 0 ? r + m : r);
}
inline ll qpow(ll base, ll power, ll mod)
{
ll res = 1;
while (power)
{
if (power & 1) res = mul(res, base, mod);
base = mul(base, base, mod), power >>= 1;
}
return res;
}
inline bool mr(ll n, ll a)
{
ll t = n - 1;
while (!(t & 1)) t >>= 1;
ll pw = qpow(a, t, n);
if ((pw == 1) || (pw == n - 1)) return false;
while ((t <<= 1) < n - 1)
{
pw = mul(pw, pw, n);
if (pw == n - 1) return false;
}
return true;
}
inline bool ptest(ll n)
{
if (n < 2) return false;
for (int i = 2; i * i <= min(n, (ll)1e4); i++)
if (n % i == 0) return false;
if (n <= 1e4) return true;
for (int i = 0; i < 8; i++)
if (mr(n, tpr[i])) return false;
return true;
}
inline ll gcd(ll a, ll b)
{
if (!a) return b;
if (a < 0) a = -a;
ll t;
while (b) t = a % b, a = b, b = t;
return a;
}
inline ll pr(ll n, ll c)
{
ll x1 = (c + 1) % n, x2 = (mul(x1, x1, n) + c) % n, val;
int tot = 0;
while (true)
{
val = mul(x1 - x2, val, n), num[++tot] = x1 - x2;
if (tot == lim)
{
if (gcd(val, n) > 1) break;
tot = 0;
}
x1 = (mul(x1, x1, n) + c) % n, x2 = (mul(x2, x2, n) + c) % n, x2 = (mul(x2, x2, n) + c) % n;
}
for (int i = 1; i <= tot; i++)
{
val = gcd(num[i], n);
if (val > 1) return val;
}
return n;
}
inline void get_fac(ll n)
{
if (ptest(n))
{
for (int i = 1; i <= tot; i++)
if (fac[i] == n) return fpw[i]++, void();
return fac[++tot] = n, fpw[tot] = 1, void();
}
ll val = pr(n, rand() % (n - 1) + 1);
while (val == n) val = pr(n, rand() % (n - 1) + 1);
get_fac(val), get_fac(n / val);
}
inline ll calc(ll p, int k, int i)
{
ll res = 0;
for (int t = 0; t <= k; t++) res = (res + qpow(p, t * x + (k - t) * i)) % mod;
return res;
}
inline void solve()
{
scanf("%lld%d%d", &n, &x, &y);
ll valn = n % mod;
BNL::getf(c, y);
for (int i = 2; (i <= 1e4) && (i * i <= n); i++)
if (n % i == 0)
{
fac[++tot] = i;
while (n % i == 0) n /= i, fpw[tot]++;
}
if (n > 1) get_fac(n);
ll ans = 0;
for (int i = 0; i <= y + 1; i++)
{
ll val = 1;
for (int j = 1; j <= tot; j++)
{
ll tmp = calc(fac[j], fpw[j], i);
tmp = (tmp - qpow(fac[j], y) * calc(fac[j], fpw[j] - 1, i)) % mod;
val = val * tmp % mod;
}
ans = (ans + c[i] * val) % mod;
}
printf("%lld\n", (ans * qpow(valn, y) % mod + mod) % mod);
for (int i = 1; i <= tot; i++) fpw[i] = 0;
tot = 0;
}
int main()
{
srand(time(0));
BNL::init();
scanf("%d", &t);
while (t--) solve();
return 0;
}