【题解】P4464 [国家集训队] JZPKIL

仙气飘飘莫反题。

显现出了很多推式子习惯的问题。

思路

莫反 + 伯努利数 + 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;
}