#include <vector>       // Required for std::vector
#include <numeric>      // Potentially useful, e.g. std::accumulate, but not used here
#include <cstdio>       // Required for fread, fwrite

// Fast I/O implementation using buffered input/output
namespace FastIO {
    // Input buffer (1MB)
    char buf[1 << 21], *p1 = buf, *p2 = buf; 
    // Output buffer (1MB)
    char obuf[1 << 21], *p = obuf; 

    // Function to get a character from the input buffer
    inline char gc() {
        // If buffer is empty, read more data from stdin
        if (p1 == p2) {
            p2 = (p1 = buf) + fread(buf, 1, 1 << 21, stdin);
            // If read failed or EOF reached, return EOF
            if (p1 == p2) return EOF;
        }
        // Return the next character from the buffer
        return *p1++;
    }

    // Function to read an integer value. Returns true on success, false on EOF.
    template<typename T>
    inline bool read(T &x) { 
        x = 0; // Initialize the value
        char c = gc(); // Get the first character
        // Check for EOF at the very beginning
        if (c == EOF) return false; 
        bool f = false; // Flag for negative numbers
        // Skip non-digit characters, check for negative sign and EOF
        while (c < '0' || c > '9') {
            if (c == '-') f = true;
            c = gc();
            // Check for EOF while skipping non-digits
            if (c == EOF) return true; // Return true if we potentially read a sign before EOF
        }
        // Read digits and build the number
        while (c >= '0' && c <= '9') {
            // Efficiently compute x = x * 10 + (c - '0')
            x = (x << 1) + (x << 3) + (c ^ 48); 
            c = gc();
        }
        // Apply negative sign if necessary
        if (f) x = -x; 
        return true; // Successfully read a number
    }

    // Variadic template function to read multiple arguments. Returns true if all reads succeed.
    template<typename T, typename ...Args>
    inline bool read(T &x, Args &...args) {
       // If the first read fails (e.g., EOF), return false
       if (!read(x)) return false; 
       // Recursively read the remaining arguments
       return read(args...); 
    }

    // Function to flush the output buffer to stdout
    inline void flush() {
        fwrite(obuf, p - obuf, 1, stdout); 
        p = obuf; // Reset the output buffer pointer
    }

    // Function to put a character into the output buffer
    inline void pc(char c) {
        // If the buffer is full, flush it
        if (p == obuf + (1 << 21)) flush();
        // Put the character into the buffer
        *p++ = c;
    }

    // Function to write an integer value to the output buffer
    template<typename T>
    void write(T x) {
        // Handle negative numbers
        if (x < 0) {
            pc('-');
            x = -x;
        }
        // Use a static buffer for converting number to string
        static char num_buf[20]; 
        int len = 0;
        // Handle the case x = 0 separately
        if (x == 0) {
           num_buf[len++] = '0';
        } else {
           // Convert positive number to string in reverse order
           while (x > 0) {
                num_buf[len++] = x % 10 + '0';
                x /= 10;
            }
        }
        // Output characters in the correct order
        while (len > 0) pc(num_buf[--len]);
    }
    
    // RAII class to ensure the output buffer is flushed at the end of the program scope
    struct flusher {
        ~flusher() { flush(); }
    } _flusher; // Global instance ensures flush() is called automatically

    // Base case for variadic write function (write single value)
    template<typename T>
    void write_sep(T x) {
        write(x);
    }

    // Variadic template function to write multiple values separated by spaces
    template<typename T, typename ...Args>
    void write_sep(T x, Args ...args) {
        write(x); // Write the first value
        pc(' ');  // Write a space separator
        write_sep(args...); // Recursively write the remaining values
    }

    // Function to write a value followed by a newline character
    template<typename T>
    void writeln(T x) {
        write(x);
        pc('\n');
    }
    
    // Variadic template function to write multiple values separated by spaces, ending with a newline
    template<typename T, typename ...Args>
    void writeln(T x, Args ...args) {
        write_sep(x, args...); // Write values with spaces
        pc('\n'); // Write the final newline
    }
    // Function to write just a newline character
     void writeln() { pc('\n'); }

} // namespace FastIO
using namespace FastIO; // Make FastIO functions available without namespace prefix

// Constants used in the problem
const int MOD = 998244353; // The specified modulus
// Maximum possible value for n + m (2e6 + 2e6 = 4e6). Add a small buffer.
const int MAX_N_PLUS_M = 4000005; 

// Arrays to store precomputed factorials and inverse factorials
long long fact[MAX_N_PLUS_M];
long long invFact[MAX_N_PLUS_M];

// Function for modular exponentiation (calculates base^exp % MOD)
long long power(long long base, long long exp) {
    long long res = 1; // Initialize result
    base %= MOD; // Reduce base modulo MOD
    while (exp > 0) {
        // If exponent is odd, multiply result by base
        if (exp % 2 == 1) res = (res * base) % MOD; 
        // Square the base
        base = (base * base) % MOD; 
        // Halve the exponent
        exp /= 2; 
    }
    return res;
}

// Function for calculating modular inverse using Fermat's Little Theorem
// Calculates n^(MOD-2) % MOD, which is the inverse of n modulo MOD (for prime MOD)
long long modInverse(long long n) {
    return power(n, MOD - 2);
}

// Function to precompute factorials and their modular inverses up to N
void precompute_combinations(int N) {
    // Basic check for non-negative N
    if (N < 0) return; 
    fact[0] = 1; // Base case: 0! = 1
    // Compute factorials iteratively: fact[i] = fact[i-1] * i % MOD
    for (int i = 1; i <= N; i++) {
        fact[i] = (fact[i - 1] * i) % MOD;
    }
    // Calculate inverse factorial of N using modular inverse
    // Ensure fact[N] is not 0 (should not happen for N < P)
     if (fact[N] == 0) {
         // This case is unlikely given P and N limits, but good practice to consider.
     }
    invFact[N] = modInverse(fact[N]);
    // Compute inverse factorials iteratively in reverse order:
    // invFact[i] = invFact[i+1] * (i+1) % MOD
    for (int i = N - 1; i >= 0; i--) {
        invFact[i] = (invFact[i + 1] * (i + 1)) % MOD;
    }
}

// Function to compute nCr mod P using precomputed factorials and inverses
// Calculates C(n, r) = n! / (r! * (n-r)!) % MOD
long long nCr_mod(int n, int r) {
    // Handle invalid arguments: r < 0 or r > n
    if (r < 0 || r > n) {
        return 0;
    }
    // Check if arguments are within bounds of precomputed arrays
    // This check is mostly for safety; based on MAX_N_PLUS_M calculation, it should be fine.
    if (n >= MAX_N_PLUS_M || r >= MAX_N_PLUS_M || (n-r) < 0 || (n-r) >= MAX_N_PLUS_M) {
       // Error: Arguments out of precomputed range.
       // Depending on requirements, could return an error code or throw exception.
       // For this problem setup, this path should not be hit.
       return -1; 
    }
    // Calculate nCr using the formula: fact[n] * invFact[r] * invFact[n-r] % MOD
    return (((fact[n] * invFact[r]) % MOD) * invFact[n - r]) % MOD;
}

// Global vector to store powers of k, reused across test cases for optimization
std::vector<long long> k_powers;
// Static variables to track the state of k_powers precomputation
static int current_k = -1; // Stores the k value for which powers are currently calculated
static int current_max_power = -1; // Stores the maximum power index (p in k^p) currently calculated

// Function to ensure that k_powers vector contains powers of k up to k^max_power_needed
void ensure_k_powers(int k, int max_power_needed) {
    // If requested power index is negative, nothing to do.
    if (max_power_needed < 0) return;

    // Check if an update (recomputation or extension) is necessary
    // Update if k has changed OR if we need a higher power than currently stored
    if (k != current_k || max_power_needed > current_max_power) {
        // Determine the required size of the vector (index 0 to max_power_needed)
        int required_size = max_power_needed + 1;
        int current_size = k_powers.size();
        
        // Resize the vector if k changed OR if the required size is larger than current capacity
        if (k != current_k || required_size > current_size) {
             k_powers.resize(required_size);
        }

        // Determine the starting index for computation
        int start_idx = 1; // Default start index is 1 (for k^1)
        if (k != current_k) {
             // If k changed, we must recompute from k^0
             k_powers[0] = 1; // k^0 = 1
             start_idx = 1;
        } else {
             // If k is the same, we extend from the last computed power
             start_idx = current_max_power + 1;
        }

        // Handle the special case where k = 0
        if (k == 0) {
            // k^0 = 1 is already set if k changed.
            // Set k^p = 0 for all p > 0.
             for(int i = start_idx; i <= max_power_needed; ++i) k_powers[i] = 0;
        } else {
            // General case for k > 0
            // Compute powers iteratively: k^i = k^(i-1) * k % MOD
            for (int i = start_idx; i <= max_power_needed; ++i) {
                 // Need k_powers[i-1] which should be available
                 k_powers[i] = (k_powers[i-1] * k) % MOD;
            }
        }

        // Update the state variables tracking the current k and max power stored
        current_k = k;
        // Update max power only if it increased or if k changed (resetting state)
        if(max_power_needed > current_max_power || k != current_k) {
           current_max_power = max_power_needed;
        } else {
            // If k is the same and required power is not higher, max_power remains same.
            // This line is technically redundant due to the outer 'if' condition,
            // but clarifies the state logic.
            current_max_power = current_max_power; 
        }
    }
     // If no update needed, the existing k_powers vector is sufficient.
}

int main() {
    // Precompute factorials and inverse factorials up to the maximum possible N = n + m
    // MAX_N_PLUS_M-1 because array indices go up to N.
    precompute_combinations(MAX_N_PLUS_M - 1); 

    int T; // Number of test cases
    read(T); // Read T using FastIO
    
    // Process each test case
    while (T--) {
        int n, m, k; // Input variables for the current test case
        read(n, m, k); // Read n, m, k using FastIO

        // Handle special case: k = 0
        // The shininess k^v is 0 for all v >= 1. Total shininess is 0.
        if (k == 0) {
            writeln(0); // Output 0
            continue; // Move to the next test case
        }
        
        // Handle special case: k = 1
        // The shininess k^v is 1 for all v >= 1. Total shininess is the total number of clusters.
        // The total number of clusters after m operations, starting with n, is F_m(1).
        // F_m(1) = C(n-1+m, m) = C(n+m-1, m).
        if (k == 1) {
             long long ans = nCr_mod(n + m - 1, m); 
             writeln(ans); // Output the result
             continue; // Move to the next test case
        }
        
        // General case: k > 1
        // Ensure powers of k up to k^(n-1) are precomputed or available in k_powers.
        // The maximum power needed is for the term k^(v-1) where v=n.
        ensure_k_powers(k, n - 1); 

        long long total_shininess = 0; // Initialize total shininess for this test case

        // Use the derived formula: Total Shininess = F_m(1)*k + (k-1) * sum_{v=2..n} F_m(v)*k^(v-1) (mod P)

        // Calculate the first component: F_m(1) * k
        // F_m(1) = C(n-1+m, m)
        long long Fm1 = nCr_mod(n - 1 + m, m);
        // Add F_m(1) * k % MOD to the total shininess
        total_shininess = (Fm1 * k) % MOD;

        // Calculate the second component: (k-1) * sum_{v=2..n} F_m(v)*k^(v-1)
        // Calculate (k-1) mod P safely (handles k=1 implicitly, but we already branched)
        long long k_minus_1 = (k - 1 + MOD) % MOD; 
        long long current_sum_term = 0; // Initialize the sum part
        // Loop through v from 2 to n
        for (int v = 2; v <= n; ++v) {
            // Calculate F_m(v) = C(n-v+m, m)
            long long Fmv = nCr_mod(n - v + m, m);
            // Get k^(v-1) from the precomputed/cached vector
            // Index needed is v-1. Max index is n-1.
            long long kv_minus_1 = k_powers[v-1]; 
            // Calculate the term F_m(v) * k^(v-1) % MOD
            long long term = (Fmv * kv_minus_1) % MOD;
            // Add the term to the running sum
            current_sum_term = (current_sum_term + term) % MOD;
        }
        // Add the second component to the total shininess: (sum * (k-1)) % MOD
        total_shininess = (total_shininess + (current_sum_term * k_minus_1) % MOD) % MOD;

        // Output the final calculated total shininess for the test case
        writeln(total_shininess);
    }

    return 0; // Indicate successful execution
}

您是对的，原始的 $O(N_{max} + \sum n)$ 复杂度分析中，虽然 $\sum n$ 可能在限制内（例如 $2 \times 10^7$），但如果 $T$ 很大（例如 $5 \times 10^5$）且某些测试用例的 $n$ 较大（例如 $10^5$），那么单次 $O(n)$ 的计算也可能导致超时。寻求一个与 $m$ 相关的复杂度是更优的策略。

让我们回顾之前的推导。我们得到了总闪耀度 $W_t = \sum_{v=1}^n f_t(v) k^v$ 的递推关系：
$W_t = C W_{t-1} - C \binom{n+t-2}{t-1}$
其中 $C = \frac{k}{k-1} \pmod{998244353}$，且此递推在 $k \neq 1$ 时成立。

我们解这个递推关系，得到：
$W_m = C^m W_0 - \sum_{i=1}^m C^{m-i+1} \binom{n+i-2}{i-1}$
其中 $W_0 = k^n$ (初始只有一个大小为 $n$ 的星团)。

所以，最终的闪耀度 $W_m$ 为：
$W_m = C^m k^n - C \sum_{i=1}^m C^{m-i} \binom{n+i-2}{i-1} \pmod{998244353}$

这个公式的计算过程如下：
1.  **处理特殊情况**:
    * 如果 $k=0$，答案是 $0$。
    * 如果 $k=1$，答案是总星团数 $F_m(1) = \binom{n+m-1}{m}$。

2.  **处理一般情况 $k > 1$**:
    * 计算 $C = k \cdot (k-1)^{-1} \pmod P$。这需要计算 $k-1$ 的模逆元。
    * 计算 $C^m \pmod P$ 和 $k^n \pmod P$。使用模快速幂，复杂度分别为 $O(\log m)$ 和 $O(\log n)$。
    * 计算求和项 $S = \sum_{i=1}^m C^{m-i} \binom{n+i-2}{i-1} \pmod P$。
        * 我们可以预先计算 $C^0, C^1, \dots, C^m$。这需要 $O(m)$ 时间。
        * 然后循环 $i$ 从 $1$ 到 $m$：
            * 取出预计算的 $C^{m-i}$。
            * 计算 $\binom{n+i-2}{i-1}$。使用预处理的阶乘和逆元，这需要 $O(1)$ 时间。注意所需的最大组合数上参数为 $n+m-2$，这在我们预处理的范围内 ($N_{max} \approx n+m$)。
            * 将 $C^{m-i} \times \binom{n+i-2}{i-1}$ 累加到 $S$ 中。
        * 计算这个和的总时间复杂度是 $O(m)$。
    * 计算最终结果 $W_m = (C^m k^n - C \cdot S) \pmod P$。注意处理减法可能产生的负数，要加上 $MOD$ 再取模。这需要 $O(1)$ 时间。

**时间复杂度分析 (新方法)**:
* **预处理组合数**: $O(N_{max}) = O(\max(n+m))$。
* **每个测试用例**:
    * 特殊情况 $k=0, k=1$：$O(1)$。
    * 一般情况 $k>1$：计算模逆元 $O(\log MOD)$，计算 $C$ 是 $O(1)$，计算 $C^m, k^n$ 是 $O(\log m + \log n)$，预计算 $C$ 的幂是 $O(m)$，计算求和项 $S$ 是 $O(m)$，最后组合结果是 $O(1)$。因此，每个测试用例的时间复杂度是 $O(m + \log n + \log m)$。
* **总时间复杂度**: $O(N_{max} + \sum (m + \log n))$。
    * 因为 $m \ge 1$, $m$ 通常远大于 $\log n$ 和 $\log m$。
    * 所以总复杂度近似为 $O(N_{max} + \sum m)$。

**对比**:
* 旧方法复杂度: $O(N_{max} + \sum n)$
* 新方法复杂度: $O(N_{max} + \sum m)$

根据题目限制 $\sum m \le 2 \times 10^7$ 和 $N_{max} \approx 4 \times 10^6$，新方法的总复杂度约为 $O(4 \times 10^6 + 2 \times 10^7) = O(2.4 \times 10^7)$，这完全符合要求，并且避免了对 $\sum n$ 的依赖，从而解决了当单个 $n$ 很大时可能出现的问题。

下面是采用了这个 $O(m)$ 计算方法的 C++ 代码实现：

#include <vector>
#include <numeric> 
#include <cstdio>  

// --- Fast I/O Start ---
namespace FastIO {
    char buf[1 << 21], *p1 = buf, *p2 = buf; 
    char obuf[1 << 21], *p = obuf; 

    inline char gc() {
        if (p1 == p2) {
            p2 = (p1 = buf) + fread(buf, 1, 1 << 21, stdin);
            if (p1 == p2) return EOF;
        }
        return *p1++;
    }

    template<typename T>
    inline bool read(T &x) { 
        x = 0; 
        char c = gc(); 
        if (c == EOF) return false; 
        bool f = false; 
        while (c != EOF && (c < '0' || c > '9')) { 
            if (c == '-') f = true;
            c = gc();
        }
         if (c == EOF) return !f; 
        if (c < '0' || c > '9') return false; 

        while (c >= '0' && c <= '9') {
            x = (x << 1) + (x << 3) + (c ^ 48); 
            c = gc();
        }
        if (f) x = -x; 
        return true; 
    }

    template<typename T, typename ...Args>
    inline bool read(T &x, Args &...args) {
       if (!read(x)) return false; 
       return read(args...); 
    }

    inline void flush() {
        fwrite(obuf, p - obuf, 1, stdout); 
        p = obuf; 
    }

    inline void pc(char c) {
        if (p == obuf + (1 << 21)) flush();
        *p++ = c;
    }

    template<typename T>
    void write(T x) {
        if (x < 0) {
            pc('-');
            x = -x;
        }
        static char num_buf[20]; 
        int len = 0;
        if (x == 0) {
           num_buf[len++] = '0';
        } else {
           while (x > 0) {
                num_buf[len++] = x % 10 + '0';
                x /= 10;
            }
        }
        while (len > 0) pc(num_buf[--len]);
    }
    
    struct flusher {
        ~flusher() { flush(); }
    } _flusher; 

    template<typename T>
    void writeln(T x) {
        write(x);
        pc('\n');
    }
} // namespace FastIO
using namespace FastIO; 
// --- Fast I/O End ---

// Constants
const int MOD = 998244353; 
const int MAX_N_PLUS_M = 4000005; // Max N+M for combinations
const int MAX_M = 2000005; // Max M for C_powers array

// Precomputed arrays
long long fact[MAX_N_PLUS_M];
long long invFact[MAX_N_PLUS_M];
long long C_powers[MAX_M]; // Stores powers of C = k/(k-1)

// Modular exponentiation
long long power(long long base, long long exp) {
    long long res = 1;
    base %= MOD;
    if (base < 0) base += MOD; // Ensure base is non-negative
    while (exp > 0) {
        if (exp % 2 == 1) res = (res * base) % MOD; 
        base = (base * base) % MOD; 
        exp /= 2; 
    }
    return res;
}

// Modular inverse using Fermat's Little Theorem
long long modInverse(long long n) {
    n %= MOD;
    if (n < 0) n += MOD; 
    if (n == 0) return -1; // Inverse of 0 does not exist
    return power(n, MOD - 2);
}

// Precompute factorials and their modular inverses up to N
void precompute_combinations(int N) {
    if (N < 0) return; 
    if (N >= MAX_N_PLUS_M) N = MAX_N_PLUS_M - 1; 
    fact[0] = 1; 
    for (int i = 1; i <= N; i++) {
        fact[i] = (fact[i - 1] * i) % MOD;
    }
    long long invN = modInverse(fact[N]);
    if (invN == -1) { return; } // Error handling if needed
    invFact[N] = invN;
    for (int i = N - 1; i >= 0; i--) {
        invFact[i] = (invFact[i + 1] * (i + 1)) % MOD;
    }
}

// Compute nCr mod P using precomputed values
long long nCr_mod(int n, int r) {
    if (r < 0 || r > n) {
        return 0;
    }
    if (n >= MAX_N_PLUS_M) { return 0; } 
    
    // Ensure n, r, n-r indices are valid for precomputed arrays
    if (r >= MAX_N_PLUS_M || (n-r) < 0 || (n-r) >= MAX_N_PLUS_M) {
        // This should not happen if n < MAX_N_PLUS_M and r <= n
        return 0; 
    }

    return (((fact[n] * invFact[r]) % MOD) * invFact[n - r]) % MOD;
}


int main() {
    // Precompute combinations up to n+m-2. Needs index up to MAX_N_PLUS_M-1.
    precompute_combinations(MAX_N_PLUS_M - 1); 

    int T; 
    read(T); 
    
    while (T--) {
        int n;
        long long m_ll; 
        int k; 
        read(n); read(m_ll); read(k); 
        
        int m = (int)m_ll; // m fits in int

        // Handle k=0
        if (k == 0) {
            writeln(0);
            continue; 
        }
        
        // Handle k=1
        if (k == 1) {
             if (n <= 0) { writeln(0); continue; } // Assuming n>=1
             long long ans = nCr_mod(n + m - 1, m); 
             writeln(ans);
             continue; 
        }
        
        // General case k > 1
        
        long long k_ll = k; 
        long long k_minus_1 = (k_ll - 1 + MOD) % MOD;
        long long inv_k_minus_1 = modInverse(k_minus_1);
        if (inv_k_minus_1 == -1) { continue; } // Should not happen

        // C = k / (k-1) mod P
        long long C = (k_ll * inv_k_minus_1) % MOD;

        // Check bounds for m
        if (m < 0) { /* handle appropriately */ continue; } // m >= 1 in problem
        if (m >= MAX_M) { /* Error: m too large for C_powers array */ continue;} 

        // Precompute powers C^0 to C^m
        C_powers[0] = 1;
        for(int i=1; i <= m; ++i) {
             C_powers[i] = (C_powers[i-1] * C) % MOD;
        }

        // Calculate k^n
        long long k_pow_n = power(k_ll, n);
        // Get C^m
        long long Cm = C_powers[m]; 

        // Calculate Sum = sum_{i=1..m} C^{m-i} * C(n+i-2, i-1)
        long long sum_term = 0;
        for (int i = 1; i <= m; ++i) {
             long long C_pow_m_minus_i = C_powers[m-i];
             // Calculate binomial term C(n+i-2, i-1)
             long long binom_term = nCr_mod(n + i - 2, i - 1);
             // Calculate term = C^{m-i} * BinomialTerm % MOD
             long long term = (C_pow_m_minus_i * binom_term) % MOD;
             // Add to the sum
             sum_term = (sum_term + term) % MOD;
        }

        // Calculate W_m = C^m * k^n - C * Sum (mod P)
        long long final_sum_part = (C * sum_term) % MOD;
        long long first_part = (Cm * k_pow_n) % MOD;

        // Calculate Result = (first_part - final_sum_part) mod P, ensuring non-negative result
        long long result = (first_part - final_sum_part + MOD) % MOD;

        // Output the result for the test case
        writeln(result);
    }

    return 0;
}

为了解决这个问题，我们需要计算小K施展魔法后的星团闪耀度之和。通过分析，我们发现每次操作会将每个星团分解为较小的星团，并保留原有星团。我们需要高效地计算这些分解操作后的闪耀度之和。

### 方法思路
1. **预处理阶乘和逆阶乘**：为了快速计算组合数，我们需要预处理阶乘和逆阶乘。这样可以快速计算出组合数，避免重复计算。
2. **特殊情况处理**：处理k=0的情况，此时所有闪耀度之和为0。
3. **等比数列求和**：使用等比数列求和公式快速计算星团数量的贡献。
4. **组合数计算**：对于每个可能的星团大小，使用预处理的阶乘和逆阶乘快速计算组合数，从而确定其贡献。

### 解决代码
#include <bits/stdc++.h>
using namespace std;

const int MOD = 998244353;
const int MAX_FACT = 4e6 + 10;

vector<int> fact(MAX_FACT), inv_fact(MAX_FACT);

int pow_mod(int a, int b) {
    int res = 1;
    while (b) {
        if (b & 1) res = 1LL * res * a % MOD;
        a = 1LL * a * a % MOD;
        b >>= 1;
    }
    return res;
}

void precompute() {
    fact[0] = 1;
    for (int i = 1; i < MAX_FACT; ++i) {
        fact[i] = 1LL * fact[i-1] * i % MOD;
    }
    inv_fact[MAX_FACT-1] = pow_mod(fact[MAX_FACT-1], MOD-2);
    for (int i = MAX_FACT-2; i >= 0; --i) {
        inv_fact[i] = 1LL * inv_fact[i+1] * (i+1) % MOD;
    }
}

int comb(int a, int b) {
    if (a < 0 || b < 0 || a < b) return 0;
    return 1LL * fact[a] * inv_fact[b] % MOD * inv_fact[a - b] % MOD;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    precompute();

    int T;
    cin >> T;
    while (T--) {
        int n, m, k;
        cin >> n >> m >> k;

        if (k == 0) {
            cout << "0\n";
            continue;
        }

        int term1 = pow_mod(k, n);

        if (n == 1) {
            cout << term1 << '\n';
            continue;
        }

        int S1;
        if (k == 1) {
            S1 = (n - 1) % MOD;
        } else {
            int kn = pow_mod(k, n-1);
            S1 = 1LL * k * (kn - 1) % MOD;
            S1 = 1LL * S1 * pow_mod(k-1, MOD-2) % MOD;
        }
        int term2 = 1LL * m * S1 % MOD;

        int term3 = 0;
        int max_d = n - 2;
        for (int d = 1; d <= max_d; ++d) {
            int a = m + d - 1;
            int b = d + 1;
            int c = comb(a, b);

            int rem = n - d - 1;
            if (rem <= 0) {
                continue;
            }

            int s;
            if (k == 1) {
                s = rem % MOD;
            } else {
                int kr = pow_mod(k, rem);
                s = 1LL * k * (kr - 1) % MOD;
                s = 1LL * s * pow_mod(k - 1, MOD-2) % MOD;
            }

            term3 = (term3 + 1LL * c * s) % MOD;
        }

        int total = (term1 + term2) % MOD;
        total = (total + term3) % MOD;
        cout << total << '\n';
    }

    return 0;
}

### 代码解释
1. **预处理阶乘和逆阶乘**：使用快速幂计算阶乘和逆阶乘，以便快速计算组合数。
2. **处理输入**：读取输入数据，处理每个测试用例。
3. **特殊情况处理**：当k=0时直接输出0，因为任何星团的闪耀度都为0。
4. **计算初始闪耀度**：初始星团大小为n的闪耀度。
5. **等比数列求和**：计算星团分解后的闪耀度之和。
6. **组合数计算**：遍历每个可能的分解大小，使用预处理数据计算组合数，累加其贡献。

通过这种方法，我们能够高效地计算出多次魔法操作后的闪耀度之和，满足题目中的时间和空间复杂度要求。