标签： 数论

$$ \text{求}\sqrt{a}\left( mod\,\,p \right) ,\text{其中}p\text{是奇素数。} \\ \text{令}p-1=2^s\times t,\text{其中}t\text{是奇数。} \\ \text{如果}s=\text{1，且}a\text{是模}p\text{的二次剩余,那么}\sqrt{a}=\sqrt{a^{\frac{p-1}{2}}\times a}=\sqrt{a^t\times a}=a^{\frac{t+1}{2}},\text{又因为}t\text{为奇数，所以可以直接计算得到}\sqrt{a} \\ \text{对于}s>\text{1,设}x_{s-1}=a^{\frac{t+1}{2}} \\ a^{\frac{p-1}{2}}=a^{2^{s-1}\times t}=\left( a^{-1}\times \left( a^{\frac{t+1}{2}} \right) ^2 \right) ^{2^{s-1}}=\left( a^{-1}\times x_{s-1}^{2} \right) ^{2^{s-1}}=1\left( a\text{是模}p\text{意义下的二次剩余} \right) \\ \text{所以}a^{-1}\times x_{s-1}^{2}\text{是模}p\text{意义下的}2^{s-1}\text{次单位根。} \\ \text{我们的目标是计算出模}p\text{意义下的}2^0\text{单位根，即计算出}a^{-1}\times x_{0}^{2}=\text{1,在这个同时我们也可以计算出}x_0,\text{进而得到方程解。} \\ \text{设}e_k\text{为模}p\text{意义下的}2^k\text{次单位根，则有} \\ e_{s-1}=a^{-1}\times x_{s-1}^{2} \\ \text{假设现在我们知道了}e_{s-k}\text{和}x_{s-k}\text{，如果} \\ e_{s-k}^{\begin{array}{c} 2^{s-k-1}\\ \end{array}}\equiv 1\left( mod\,\,p \right) ,\text{那么}e_{s-k-1}=e_{s-k},x_{s-k-1}=x_{s-k} \\ \text{否则，如果}e_{s-k}^{\begin{array}{c} 2^{s-k-1}\\ \end{array}}\equiv -1\left( mod\,\,p \right) \\ \text{即}\left( a^{-1}\times x_{s-k}^{2} \right) ^{2^{s-k-1}}\equiv -\text{1,那么我们显然需要求一个数}q,\text{使得} \\ \left( a^{-1}\times \left( qx_{s-k} \right) ^2 \right) ^{2^{s-k-1}}\equiv 1 \\ \text{整理后得}q^{2^{s-k}}\equiv -1 \\ \text{寻找一个模}p\text{意义下的非二次剩余}b,\text{则有}b^{\frac{p-1}{2}}\equiv -1 \\ b^{2^{s-1}t}\equiv -1 \\ b^{2^{k-1}t2^{s-k}}\equiv -1 \\ \left( b^{2^{k-1}t} \right) ^{2^{s-k}}\equiv -\text{1,故}q=b^{2^{k-1}t} \\ \text{所以}x_{s-k-1}=x_{s-k}\times q \\ \text{一直递推到}x_0\text{即可。} \\ $$

模数固定为998244353:

int_t modSqrt(int_t a) {
    const int_t b = 3;
    //mod-1=2^s*t
    int_t s = 23, t = 119;
    int_t x = power(a, (t + 1) / 2);
    int_t e = power(a, t);
    int_t k = 1;
    while (k < s) {
        if (power(e, 1 << (s - k - 1)) != 1) {
            x = x * power(b, (1 << (k - 1)) * t) % mod;
        }
        e = power(a, -1) * x % mod * x % mod;
        k++;
    }
    return x;
}

2019-3-14更新：换成了拉格朗日插值

$$ \text{注意到}F_k\left( n \right) \text{是积性函数} \\ \text{对于质数幂}F_k\left( p^c \right) =\sum_{0\le i\le c}{F}_{k-1}\left( p^i \right) \\ \text{可以注意到}F_k\left( p^c \right) \text{的取值与}p\text{无关} \\ \text{设}F_k\left( p^n \right) \text{的生成函数为}G_k\left( x \right) =\sum_{n\ge 0}{F}_k\left( p^n \right) x^n \\ \text{由上述可得}G_k\left( x \right) =\frac{G_{k-1}\left( x \right)}{1-x} \\ \text{其中}G_0\left( x \right) =\frac{1}{1-x} \\ \text{故}G_k\left( x \right) =\frac{1}{\left( 1-x \right) ^{k+1}} \\ \text{由广义二项式定理可得} \\ F_k\left( p^n \right) =\left( \begin{array}{c} n+k\\ n\\ \end{array} \right) \\ \text{同时又注意到，}F_k\left( p_{1}^{c_1}p_{2}^{c_2}p_{3}^{c_3}….p_{t}^{c_t} \right) \text{的取值只与}\left\{ c_1,c_2….c_t \right\} \text{有关。} \\ \text{写个}DFS\text{看一下可行的}\left\{ c_1,c_2..c_t \right\} \text{的取值只有244种} \\ \text{又因为}\left( \begin{array}{c} c+x\\ c\\ \end{array} \right) =\frac{\left( c+x \right) !}{c!x!}\text{是一个关于}x\text{的}c\text{次多项式} \\ \text{所以,}\prod_i{\left( \begin{array}{c} c_i+x\\ c_i\\ \end{array} \right)}\text{是一个关于}x\text{的}\sum_i{c_i}\text{次多项式。} \\ \text{又因为}\sum_i{c_i}\text{的取值上界为}O\left( \log n \right) \\ \text{所以直接拉格朗日插值求出这个多项式，就可以}O\left( \log N \right) \text{的时间计算出}F_k\left( n \right) \\ $$

#pragma GCC optimize("unroll-loops")

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <map>
#include <numeric>
#include <set>
#include <vector>
using int_t = int;

using std::cin;
using std::cout;
using std::endl;

const int_t mod = 998244353;
const int_t LIMIT = 5e5;
const int_t LARGE = 5e5;
int fact[LIMIT + 1];
int factInv[LIMIT + 1];
bool isPrime[LIMIT + 1];
int_t factor[LIMIT + 1];
bool operator<(const std::vector<int_t>& A, const std::vector<int_t>& B) {
    if (A.size() != B.size())
        return A.size() < B.size();
    else {
        for (int_t i = 0; i < A.size(); i++) {
            if (A[i] != B[i]) return A[i] < B[i];
        }
        return false;
    }
}
bool operator==(const std::vector<int_t>& A, const std::vector<int_t>& B) {
    if (A.size() != B.size()) return false;
    return std::equal(A.begin(), A.end(), B.begin());
}
template <class T>
std::ostream& operator<<(std::ostream& os, const std::vector<T>& vec) {
    for (auto x : vec) os << x << " ";
    return os;
}
std::vector<std::vector<int>> states;
std::set<std::vector<int>> set;
int polys[244][25];
std::vector<int> state[LARGE + 1];
std::map<std::vector<int>, int> hash;
int prefix[LARGE + 1][244];
int C(int n, int m) {
    return 1ll * fact[n] * factInv[m] % mod * factInv[n - m] % mod;
}
int_t power(int_t base, int_t index) {
    int_t result = 1;
    base = (base % mod + mod) % mod;
    while (index) {
        if (index & 1) result = 1ll * result * base % mod;
        base = 1ll * base * base % mod;
        index >>= 1;
    }
    return result;
}
int maxn, maxk, q;
//将k代入状态
int subState(const std::vector<int>& state, int_t k) {
    int result = 1;
    for (int x : state) result = 1ll * result * C(x + k, x) % mod;
    return result;
}
int __attribute__((hot)) subPoly(int* poly, int x) {
    int result = 0;
    for (int i = poly[0]; i >= 1; i--) {
        result = (1ll * result * x + poly[i]) % mod;
    }
    return result;
};
void poly_mul(int* A, int n, const int* B, int m) {
    static int C[50];
    memset(C, 0, sizeof(C));
    for (int i = 0; i <= n; i++) {
        for (int j = 0; j <= m; j++) {
            C[i + j] = (C[i + j] + 1ll * A[i] * B[j] % mod) % mod;
        }
    }
    std::copy(C, C + n + m + 1, A);
}
//计算n次多项式A整除m次多项式B，忽略余数
void poly_div(const int* A, int n, const int* B, int m, int* C) {
    static int R[50];
    std::copy(A, A + n + 1, R);
    for (int i = n - m; i >= 0; i--) {
        int_t x = (int_t)R[i + m] * power(B[m], mod - 2) % mod;
        C[i] = x;
        for (int_t j = m; j >= 0; j--) {
            R[i + j] = (R[i + j] - 1ll * B[j] * x + 2ll * mod) % mod;
        }
    }
}
// state是一组{c_1,c_2...c_n}的取值
//求出对应的多项式
void interpolate(const std::vector<int>& state, int* output) {
    static int prod[50];
    static int result[50];
    memset(prod, 0, sizeof(prod));
    memset(result, 0, sizeof(result));
    prod[0] = 1;
    // k是要带进去的点的数量
    // k-1是结果多项式的次数
    int k = std::accumulate(state.begin(), state.end(), 0) + 1;
    for (int_t i = 0; i < k; i++) {
        static int A[20];
        A[0] = (mod - i) % mod;
        A[1] = 1;
        poly_mul(prod, k, A, 1);
    }
#ifdef DEBUG
    cout << "prod = ";
    for (int_t i = 0; i <= k; i++) cout << prod[i] << " ";
    cout << endl;
#endif
    for (int i = 0; i < k; i++) {
        int y = subState(state, i);
#ifdef DEBUG
        cout << "y of " << i << " = " << y << endl;
#endif
        static int A[50];
        static int B[50];
        B[0] = (mod - i) % mod;
        B[1] = 1;
        poly_div(prod, k, B, 1, A);
#ifdef DEBUG

        cout << "divide x-" << i << " = ";
        for (int_t i = 0; i <= k - 1; i++) cout << A[i] << " ";
        cout << endl;
#endif
        int x0 = 1;
        for (int_t j = 0; j < k; j++) {
            if (j != i) x0 = 1ll * x0 * (i - j + mod) % mod;
        }
#ifdef DEBUG
        cout << "x0 of " << i << " = " << x0 << endl;
#endif
        y = 1ll * y * power(x0, mod - 2) % mod;
        for (int_t j = 0; j < k; j++)
            result[j] = (1ll * result[j] + 1ll * y * A[j] % mod + mod) % mod;
    }
    output[0] = k;
    std::copy(result, result + k, output + 1);
}
int main() {
    // freopen("qwq.txt", "r", stdin);
    fact[0] = fact[1] = factInv[0] = factInv[1] = 1;
    for (int_t i = 2; i <= LIMIT; i++) {
        fact[i] = 1ll * fact[i - 1] * i % mod;
        factInv[i] = 1ll * (mod - mod / i) * factInv[mod % i] % mod;
    }
    for (int_t i = 2; i <= LIMIT; i++) {
        factInv[i] = 1ll * factInv[i - 1] * factInv[i] % mod;
    }
    memset(isPrime, true, sizeof(isPrime));
    isPrime[1] = 0;
    for (int_t i = 2; i <= LIMIT; i++) {
        if (isPrime[i]) {
            for (long long j = 1ll * i * i; j <= LIMIT; j += i) {
                isPrime[j] = false;
                factor[j] = i;
            }
            factor[i] = i;
        }
    }
    scanf("%d%d%d", &maxn, &maxk, &q);
    for (int_t i = 2; i <= maxn; i++) {
        int x = i;
        std::vector<int> index;
        while (x != 1) {
            int times = 0;
            int p = factor[x];
            while (x % p == 0) {
                times++;
                x /= p;
            }
            index.push_back(times);
        }
        std::sort(index.begin(), index.end());
        state[i] = index;
        set.insert(index);
    }
    states.resize(set.size());
    std::copy(set.begin(), set.end(), states.begin());

    for (const auto& vec : states) {
        hash[vec] = hash.size();
        interpolate(vec, polys[hash[vec]]);
#ifdef DEBUG
        cout << "poly of ";
        for (int_t x : vec) cout << x << " ";

        cout << "is ";
        for (int_t x : polys[hash[vec]]) cout << x << " ";
        cout << endl;
#endif
    }
    for (int i = 2; i <= maxn; i++) {
        memcpy(prefix[i], prefix[i - 1], sizeof(prefix[i]));
        auto hash0 = hash[state[i]];
        prefix[i][hash0] = (prefix[i][hash0] + 1) % mod;
#ifdef DEBUG
        cout << "count 1 at " << i << " for state " << state[i] << "("
             << hash[state[i]] << ")" << endl;
#endif
    }
    for (int_t i = 1; i <= q; i++) {
        int k, left, right;
        scanf("%d%d%d", &left, &right, &k);
        int result = 0;
        for (int j = 0; j < states.size(); j++) {
            int count = prefix[right][j] - prefix[left - 1][j];
            if (count > 0) {
                result =
                    (1ll * result + 1ll * count * subPoly(polys[j], k) % mod) %
                    mod;
            }
        }
        if (left == 1) result = (result + 1) % mod;
        printf("%d\n", result);
    }
    return 0;
}

$$ \text{求}F\left( n \right) \,\,mod\,\,p,n\le 10^{30000000},p\le 2^{31}-1 \\ \text{斐波那契数列在模}p\text{意义下的循环节长度为}O\left( p \right) \text{级别。} \\ \text{首先将}p\text{进行质因数分解，得到}p_{1}^{k_1}p_{2}^{k_2}p_{3}^{k_3}… \\ \text{设}T\left( k \right) \text{为模}k\text{意义下的循环节长度，则}T\left( p \right) =\prod{T\left( p_i^{k_i} \right)} \\ \text{同时对于}p^k,\text{有}T\left( p^k \right) =T\left( p \right) \times p^{k-1} \\ \text{考虑如何求出模质数的循环节长度。} \\ \text{对于2,3,5直接特判。} \\ \text{对于其他质数}p\text{，如果5是模}p\text{意义下的二次剩余}\left( \text{即}5^{\frac{p-1}{2}}\equiv 1\left( mod\,\,p \right) \right) ,\text{则模}p\text{意义下的循环节长度是}p-\text{1的约数。} \\ \text{部分证明：} \\ \text{斐波那契数列的通项公式为}F_n=\frac{\sqrt{5}}{5}\left( \left( \frac{\left( 1+\sqrt{5} \right)}{2} \right) ^n-\left( \frac{\left( 1-\sqrt{5} \right)}{2} \right) ^n \right) \\ \text{设5在模}p\text{意义下的二次剩余为}x\text{，则}F_n\text{一定可以表示为}k\left( x_0^n-x_1^n \right) \text{的形式。} \\ \text{其中}x_0^n\text{和}x_1^n\text{的周期一定是}\varphi \left( p \right) \text{的约数。} \\ \text{故}F_n\text{的周期是}\varphi \left( p \right) \text{的约数。} \\ \\ \text{反之，模}p\text{意义下的循环节长度是}2\left( p+1 \right) \text{的约数。} \\ \text{枚举约数验证即可。} \\ \text{这个我不会证} \\ \\ $$

#include <climits>
#include <cstring>
#include <iostream>
#include <map>
using int_t = long long int;
using std::cin;
using std::cout;
using std::endl;
struct Matrix mat_mul(const struct Matrix& a, const struct Matrix& b,
                      int_t mod);
int_t power(int_t base, int_t index, int_t mod);
struct Matrix {
    int_t data[4][4];
    Matrix() { memset(data, 0, sizeof(data)); }
    int_t& operator()(int_t r, int_t c) { return data[r][c]; }
    const int_t N = 2;
    Matrix& operator=(const Matrix& mat2) {
        memcpy(data, mat2.data, sizeof(data));
        return *this;
    }
    Matrix power(int_t index, int_t mod) {
        Matrix result;
        Matrix base = *this;
        result(1, 1) = result(2, 2) = 1;
        while (index) {
            if (index & 1) result = mat_mul(result, base, mod);
            base = mat_mul(base, base, mod);
            index >>= 1;
        }
        return result;
    }
};

Matrix mat_mul(const Matrix& a, const Matrix& b, int_t mod) {
    const auto N = a.N;
    Matrix result;
    for (int_t i = 1; i <= N; i++) {
        for (int_t j = 1; j <= N; j++) {
            for (int_t k = 1; k <= N; k++) {
                result.data[i][j] =
                    (result.data[i][j] + a.data[i][k] * b.data[k][j] % mod) %
                    mod;
            }
        }
    }
    return result;
}
Matrix F(int_t k, int_t mod) {
    Matrix mat;
    mat(1, 1) = 0;
    mat(1, 2) = 1;
    mat(2, 1) = 1;
    mat(2, 2) = 1;
    return mat.power(k, mod);
}
//计算模p意义下的循环周期，p是质数
int_t getPeriod(int_t p) {
    if (p == 2)
        return 3;
    else if (p == 3)
        return 8;
    else if (p == 5)
        return 20;
    int_t fac;
    if (power(5, (p - 1) / 2, p) == 1) {
        fac = p - 1;
    } else {
        fac = 2 * (p + 1);
    }
    const auto check = [&](int_t x) -> bool {
        Matrix res = ::F(x + 1, p);
        return res(1, 1) == 0 && res(1, 2) == 1;
    };
    int_t result = 1e12;
    for (int_t i = 1; i * i <= fac; i++) {
        if (fac % i == 0) {
            if (i != 1 && check(i)) {
                result = std::min(result, i);
            }
            if (fac / i != i && check(fac / i)) {
                result = std::min(result, fac / i);
            }
        }
    }
    return result;
}
std::map<int_t, int_t> divide(int_t x) {
    std::map<int_t, int_t> result;
    for (int_t i = 2; i * i <= x; i++) {
        if (x % i == 0) {
            result[i] = 0;
            while (x % i == 0) {
                x /= i;
                result[i]++;
            }
        }
    }
    if (x > 1) result[x] += 1;
    return result;
}
int main() {
    static char buf[30000000 + 10];
    scanf("%s", buf);
    int_t p;
    cin >> p;
    int_t length = 1;
    for (const auto& pair : divide(p)) {
        length *= ::getPeriod(pair.first) *
                  power(pair.first, pair.second - 1, int_t(1) << 40);
    }
    int_t sum = 0;
    for (auto ptr = buf; *ptr != '\0'; ptr++) {
        sum = (sum * 10 + *ptr - '0') % length;
    }
    cout << F(sum , p)(1, 2) << endl;
    return 0;
}

int_t power(int_t base, int_t index, int_t mod) {
    int_t result = 1;
    while (index) {
        if (index & 1) result = result * base % mod;
        index >>= 1;
        base = base * base % mod;
    }
    return result;
}

Min_25筛板子.

注意到对于质数$p$,$f(p)=p-1(p\geq 3)$,而对于$2$,$f(2)=3$.

所以先默认对于所有质数,$f(p)=p-1$,拆成两个函数$h(x)=x$和$g(x)=1$,分别按照Min25筛里求出来$g_0(n,k)$和$g_1(n,k)$即可。

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <inttypes.h>
using int_t = long long int;

using std::cin;
using std::cout;
using std::endl;

const int_t mod = 1000000007;
const int_t LARGE = 1e6;

int_t hash1[LARGE + 1], hash2[LARGE + 1];
int_t sum0[LARGE + 1], sum1[LARGE + 1];
int_t g0[LARGE + 1], g1[LARGE + 1];
int_t primes[LARGE + 1];
int_t numbers[LARGE + 1];
int_t sqr;
int_t n;
void process(int_t n)
{
    static bool isPrime[LARGE + 1];
    memset(isPrime, 1, sizeof(isPrime));
    isPrime[1] = false;
    for (int_t i = 2; i <= n; i++)
    {
        if (isPrime[i])
        {
            primes[++primes[0]] = i;
            sum0[primes[0]] = (sum0[primes[0] - 1] + 1) % mod;
            sum1[primes[0]] = (sum1[primes[0] - 1] + i) % mod;
#ifdef DEBUG
            cout << "sum0 " << primes[0] << " = " << sum0[primes[0]] << endl;
            cout << "sum1 " << primes[0] << " = " << sum1[primes[0]] << endl;

#endif
            for (int_t j = i * i; j <= LARGE; j += i)
            {
                isPrime[j] = false;
            }
        }
    }
}

int_t S(int_t n, int_t k)
{
    if (n <= 1 || primes[k] > n)
        return 0;
    int_t result = 0;
    if (n <= sqr)
        result += (g1[hash1[n]] - g0[hash1[n]] + 2 * mod) % mod;
    else
        result += (g1[hash2[::n / n]] - g0[hash2[::n / n]] + 2 * mod) % mod;
    result = (result - (sum1[k - 1] - sum0[k - 1]) + 2 * mod) % mod;
#ifdef DEBUG
    cout << "S " << n << ",k(" << primes[k] << ") prime answer = " << result << endl;
#endif
    for (int_t i = k; i <= primes[0] && primes[i] * primes[i] <= n; i++)
    {
        for (int_t p = primes[i], j = 1; p * primes[i] <= n; p *= primes[i], j++)
        {
            //应该是去查质因子大于等于P_{i+1},而不是k+1(因为已经把i+1的除掉了)
            result = (result + S(n / p, i + 1) * (primes[i] ^ j) + ((primes[i]) ^ (j + 1))) % mod;
        }
    }
#ifdef DEBUG
    cout << "S " << n << "," << k << "(" << primes[k] << ") = " << result << endl;
#endif
    return result;
}

int main()
{
    const auto S2 = [](int_t x) -> int_t { return (__int128_t)x * (x + 1) / 2 % mod; };
    cin >> n;
    if(n==1){
         cout<<1<<endl;
         return 0;
    }
    
    sqr = sqrt(n);
    process(sqr);
    for (int_t i = 1, next; i <= n; i = next + 1)
    {
        next = n / (n / i);
        numbers[++numbers[0]] = n / i;
        if (n / i <= sqr)
            hash1[n / i] = numbers[0];
        else
            hash2[n / (n / i)] = numbers[0];
        //numbers[0]写成numbers[i]
        g0[numbers[0]] = n / i - 1;
        g1[numbers[0]] = S2(n / i) - 1;
    }
    for (int_t j = 1; j <= primes[0]; j++)
    {
        for (int_t i = 1; i <= numbers[0] && primes[j] * primes[j] <= numbers[i]; i++)
        {
            int_t trans;
            //要判断numbers[i]/primes[j]而不是别的
            if (numbers[i] / primes[j] <= sqr)
                trans = hash1[numbers[i] / primes[j]];
            else
                trans = hash2[n / (numbers[i] / primes[j])];
            g0[i] = (g0[i] - (g0[trans] - sum0[j - 1]) + mod * 2) % mod;
            //primes[j]没写
            g1[i] = (g1[i] - primes[j] * (g1[trans] - sum1[j - 1]) + mod * 2) % mod;
#ifdef DEBUG
            cout << "G0 " << numbers[i] << "," << j << "(" << primes[j] << ") = " << g0[i] << endl;
            cout << "G1 " << numbers[i] << "," << j << "(" << primes[j] << ") = " << g1[i] << endl;

#endif
        }
    }
    //忘了加上1的情况
    cout << (S(n, 1) + 3) % mod << endl;
    return 0;
}

$$ \\ \text{设素数集合}P,P_i\text{表示第}i\text{个素数,}P_1=2 \\ \text{设积性函数}F\left( x \right) ,\text{在}x\text{为质数时满足}F\left( x \right) =f\left( x \right) \\ \text{设}g\left( n,k \right) \text{为}n\text{范围内，质数和最小质因子大于}P_k\text{的数的}f\left( x \right) \text{之和} \\ \text{即把}f\left( 1 \right) ,f\left( 2 \right) …f\left( n \right) \text{这些数写出来，运行埃氏筛}k\text{次后，没被筛掉的数}+\text{素数的}f\text{函数的和} \\ g\left( n,k \right) \text{可以以类似于递推的形式进行转移。} \\ \text{如果}P_k^2>n \\ \text{此时}g\left( n,k \right) =g\left( n,k-1 \right) \\ \text{这时候运行筛法筛不掉任何东西。} \\ \text{如果}P_k^2\le n \\ \text{那么}g\left( n,k \right) =g\left( n,k-1 \right) -\left( \text{最小质因子恰好为}P_k\text{的数的影响} \right) \\ \text{我们筛掉所有最小质因子为}P_k\text{的数} \\ g\left( n,k \right) =\begin{cases} g\left( n,k-1 \right) \left( P_{k}^{2}>n \right)\\ g\left( n,k-1 \right) -f\left( p_k \right) \left( g\left( \frac{n}{p_k},k-1 \right) -\sum_{1\le i\le k-1}{f\left( p_i \right)} \right) \left( P_{k}^{2}\le k \right)\\ \end{cases} \\ -\sum_{1\le i\le k-1}{f\left( P_i \right)}\text{表示把最小质因子不是}P_i\text{的质数减掉，这些数不应该参与这次筛。} \\ \text{最终,}g\left( n,|P| \right) \text{表示}n\text{范围内，所有质数的}f\left( x \right) \text{之和。} \\ \text{即运行完埃氏筛后筛出来的数之和。} \\ \text{这部分的复杂度是}O\left( \frac{n^{\frac{3}{4}}}{\log\text{ }n} \right) \\ \text{然后如何求积性函数前缀和？} \\ \text{设}S\left( n,k \right) \text{表示}n\text{之内，最小质因子大于等于}P_k\text{的}F\left( x \right) \text{之和} \\ \text{首先统计}n\text{范围内,大于等于}P_k\text{质数的贡献为}g\left( n,|P| \right) -\sum_{1\le i\le k-1}{F\left( P_i \right)} \\ \text{然后暴力枚举大于等于}P_k\text{的质因子及其出现次数,根据积性函数的性质统计答案}. \\ \text{注意要额外加上}f\left( P_i^{e+1} \right) ,\text{因为}S\left( \frac{n}{P_i^e},i+1 \right) \text{不包括}f\left( P_i^e \right) \\ S\left( n,k \right) =g\left( n,|P| \right) -\sum_{1\le i\le k-1}{F\left( P_i \right)}+\sum_{i\ge k\,\,and\,\,P_i^2\le n}{\sum_{e\ge \text{1 }and\,\,P_i^{e+1}\le n}{\left( S\left( \frac{n}{P_i^e},i+1 \right) f\left( P_i^e \right) +f\left( P_i^{e+1} \right) \right)}} \\ \text{则}\sum_{1\le i\le n}{F\left( i \right)}=F\left( 1 \right) +S\left( n,1 \right) \\ $$

https://www.luogu.org/problemnew/show/P4213

#include <iostream>
#include <cstring>
#include <cstdio>
#include <vector>
#include <unordered_map>
#include <cmath>
#include <map>
#include <assert.h>
using int_t = long long int;

using std::cin;
using std::cout;
using std::endl;

const int_t LARGE = 1e6;
int_t primes[LARGE + 1];

int_t primeCount = 0;
int_t hash1[LARGE + 1], hash2[LARGE + 1];
int_t numbers[LARGE + 1];
//g0- f(x)=1
//g1 -f(x)=x
int_t g0[LARGE + 1], g1[LARGE + 1];
int_t sum0[LARGE + 1], sum1[LARGE + 1];
int_t sqr;
int_t n;
int_t S2(int_t x)
{
    return x * (x + 1) / 2;
}

void process(int_t n)
{
    static bool isPrime[LARGE + 1];
    memset(isPrime, 1, sizeof(isPrime));
    isPrime[1] = false;
    for (int_t i = 2; i <= n; i++)
    {
        for (int_t j = i * i; j <= n; j += i)
        {
            isPrime[j] = false;
        }
    }
    for (int_t i = 1; i <= n; i++)
    {
        if (isPrime[i])
        {
            primes[++primeCount] = i;
            sum0[primeCount] = sum0[primeCount - 1] + 1;
            sum1[primeCount] = sum1[primeCount - 1] + i;
        }
    }
}

//返回n以内，最小质因子大于等于p_k的数的欧拉函数之和
int_t Sphi(int_t n, int_t k)
{
    if (n <= 1 || primes[k] > n)
        return 0;

    int_t result = 0;
    if (n <= sqr)
        result += g1[hash1[n]] - g0[hash1[n]];
    else
        result += g1[hash2[::n / n]] - g0[hash2[::n / n]];
    //统计质数贡献的答案
    result -= sum1[k - 1] - sum0[k - 1];
    //枚举质因子
    for (int_t i = k; primes[i] * primes[i] <= n && i <= primeCount; i++)
    {
        for (int_t p = primes[i], e = 1; p * primes[i] <= n; p = p * primes[i], e++)
        {
            result += Sphi(n / p, i + 1) * (p - p / primes[i]) + (p * primes[i] - p);
        }
    }
    return result;
}

int_t Smu(int_t n, int_t k)
{
    if (n <= 1 || primes[k] > n)
        return 0;

    int_t result = 0;
    if (n <= sqr)
        result += -g0[hash1[n]];
    else
        result += -g0[hash2[::n / n]];
    //统计质数贡献的答案
    result -= -sum0[k - 1];
    //枚举质因子
    for (int_t i = k; primes[i] * primes[i] <= n && i <= primeCount; i++)
    {
        result += Smu(n / primes[i], i + 1) * -1;
    }
    return result;
}

void process()
{
    cin >> n;
    sqr = ::sqrt(n);
    numbers[0] = 0;
    for (int_t i = 1, next = 1; i <= n; i = next + 1)
    {
        next = n / (n / i);
        numbers[++numbers[0]] = n / i;
        if (n / i <= sqr)
        {
            hash1[n / i] = numbers[0];
        }
        else
        {
            hash2[n / (n / i)] = numbers[0];
        }
        g0[numbers[0]] = (n / i) - 1;
        g1[numbers[0]] = S2(n / i) - 1;
    }
    for (int_t j = 1; j <= primeCount; j++)
    {
        for (int_t i = 1; i <= numbers[0] && primes[j] * primes[j] <= numbers[i]; i++)
        {
            int_t trans = 0;
            if (numbers[i] / primes[j] <= sqr)
            {
                trans = hash1[numbers[i] / primes[j]];
            }
            else
            {
                trans = hash2[n / (numbers[i] / primes[j])];
            }
            g0[i] = g0[i] - (g0[trans] - sum0[j - 1]);
            g1[i] = g1[i] - primes[j] * (g1[trans] - sum1[j - 1]);
        }
    }
    cout << Sphi(n, 1) + 1 << " " << Smu(n, 1) + 1 << endl;
}
int main()
{
    process(1e6);
    int_t T;
    cin >> T;
    for (int_t i = 1; i <= T; i++)
    {
        process();
    }
    return 0;
}

$$ \sum_{1\le i\le A}{\sum_{1\le j\le B}{\sum_{1\le k\le C}{\sum_{x|i}{\sum_{y|j}{\sum_{z|k}{\left[ gcd\left( x,y \right) ==1 \right] \left[ gcd\left( x,z \right) ==1 \right] \left[ gcd\left( y,z \right) ==1 \right]}}}}}} \\ \sum_{1\le i\le A}{\sum_{1\le j\le B}{\sum_{1\le k\le C}{\sum_{x|i}{\sum_{y|j}{\sum_{z|k}{\left[ gcd\left( x,y \right) ==1 \right] \left[ gcd\left( x,z \right) ==1 \right] \left[ gcd\left( y,z \right) ==1 \right]}}}}}} \\ \sum_{1\le x\le A}{\sum_{1\le i\le \lfloor \frac{A}{x} \rfloor}{\sum_{1\le y\le B}{\sum_{1\le j\le \lfloor \frac{B}{y} \rfloor}{\sum_{1\le z\le C}{\sum_{1\le k\le \lfloor \frac{C}{z} \rfloor}{\left[ gcd\left( x,y \right) ==1 \right] \left[ gcd\left( x,z \right) ==1 \right] \left[ gcd\left( y,z \right) ==1 \right]}}}}}} \\ \sum_{1\le x\le A}{\sum_{1\le y\le B}{\sum_{1\le z\le C}{\lfloor \frac{A}{x} \rfloor \lfloor \frac{B}{y} \rfloor \lfloor \frac{C}{z} \rfloor \sum_{u|x\ and\ u|y}{\mu \left( u \right)}\sum_{v|x\ and\ v|z}{\mu \left( v \right)}\sum_{w|y\ and\ w|z}{\mu \left( w \right)}}}} \\ \sum_{1\le u\le \min \left( A,B \right)}{\sum_{1\le v\le \min \left( A,C \right)}{\sum_{1\le w\le \min \left( B,C \right)}{\mu \left( u \right) \mu \left( v \right) \mu \left( w \right)}}}\sum_{u|x\ and\ u|y}{\sum_{v|x\ and\ v|z}{\sum_{w|y\ and\ w|z}{\lfloor \frac{A}{x} \rfloor \lfloor \frac{B}{y} \rfloor \lfloor \frac{C}{z} \rfloor}}} \\ \sum_{1\le u\le \min \left( A,B \right)}{\sum_{1\le v\le \min \left( A,C \right)}{\sum_{1\le w\le \min \left( B,C \right)}{\mu \left( u \right) \mu \left( v \right) \mu \left( w \right)}}}\sum_{lcm\left( u,v \right) |x}{\lfloor \frac{A}{x} \rfloor \sum_{lcm\left( u,w \right) |y}{\lfloor \frac{B}{y} \rfloor \sum_{lcm\left( w,v \right) |z}{\lfloor \frac{C}{z} \rfloor}}} \\ \text{设}f\left( p,A \right) =\sum_{p|x\ and\ x\le A}{\lfloor \frac{A}{x} \rfloor} \\ \sum_{1\le u\le \min \left( A,B \right)}{\sum_{1\le v\le \min \left( A,C \right)}{\sum_{1\le w\le \min \left( B,C \right)}{\mu \left( u \right) \mu \left( v \right) \mu \left( w \right) f\left( lcm\left( u,v \right) ,A \right) f\left( lcm\left( u,w \right) ,B \right) f\left( lcm\left( w,v \right) ,C \right)}}} \\ \text{可以注意到不为0的}\mu \left( x \right) \text{的个数不多。} \\ \text{可以构建一张图，图中的点是满足}\mu \left( x \right) \text{不为0的}x\text{，两个点}u,v\text{之间有边当且仅当}\mu \left( u \right) \ne \text{0且}\mu \left( v \right) \ne \text{0且}lcm\left( u,v \right) <\max \left( A,B,C \right) \\ \text{可以断定满足}lcm\left( u,v \right) <10^5\text{的}u,v\text{并不多}. \\ \text{所以直接枚举}lcm\left( u,v \right) \text{建图。} \\ \text{如果要满足}\mu \left( u \right) \ne \text{0且}\mu \left( v \right) \ne \text{0,那么一定满足}\mu \left( uv \right) \ne 0. \\ \text{所以枚举所有}\mu \left( x \right) \ne \text{0且}x\le \max \left( A,B,C \right) \text{的}x,\text{枚举}x\text{的质因子集合}S\text{的一个子集}u\text{，然后再枚举}u\text{的一个子集}v \\ \text{这样}lcm\left( v\times \frac{x}{u},u \right) =x\text{显然成立}. \\ \text{设}S=2\times 3\times 5\times 7 \\ u=2\times 5 \\ \text{则如果一个}v\text{满足}lcm\left( u,v \right) =S,\text{那么}v\text{中一定包含质因子5,7,同时可以包含}u\text{的任何一个子集。} \\ \text{建图之后，去掉所有重边和自环。} \\ \text{然后考虑如何枚举三元环。} \\ \text{直接暴力枚举显然是有问题的。} \\ \text{首先把所有无向边重定向为有向边，对于边}\left( u,v \right) ,\text{如果}deg\left( u \right) <deg\left( v \right) ,\text{则让}u\text{指向}v,\text{如果}deg\left( u \right) =deg\left( v \right) , \\ \text{则让编号小的指向编号大的} \\ \text{可以证明，这样每个点的出度不会超过}\sqrt{n} \\ \text{然后直接对于每个点}u\text{，暴力枚举两条出边}\left( u,v \right) \ \left( u,w \right) ,\text{然后看是否满足}lcm\left( v,w \right) \le \max\text{,如果满足，则说明存在无序三元环} \\ \left( u,v,w \right) ,\text{统计进答案。} \\ \text{这只考虑了三元环三个数均不相同的情况，除此之外还需要考虑两个数相同和三个数都相同的情况。} \\ \\ \\ \\ \\ \ $$

#pragma GCC optimize("O3")
#include <iostream>
#include <algorithm>
#include <vector>
#include <cstring>
#include <assert.h>
using int_t = long long int;
using std::cin;
using std::cout;
using std::endl;

const int_t mod = 1e9 + 7;

const int_t LARGE = 1e5;

char mu[LARGE + 1];
int_t factor[LARGE + 1];
std::vector<int_t> pfactor[LARGE + 1];
int f[4][LARGE + 1];

int gcd(int a, int b)
{
    if (b == 0)
        return a;
    else
        return gcd(b, a % b);
}
int lcm(int a, int b)
{
    return a / gcd(a, b) * b;
}
void calc_f(int_t id, int_t A)
{
    for (int_t i = 1; i <= A; i++)
    {
        for (int_t j = 1; j * i <= A; j++)
        {
            f[id][i] = (f[id][i] + A / (j * i)) % mod;
        }
    }
}
int prod(int set, int n)
{
    int prod = 1;
    for (size_t i = 0; i < pfactor[n].size(); i++)
    {
        if (set & (1 << i))
            prod *= pfactor[n][i];
    }
    return prod;
}
int_t process(int_t A, int_t B, int_t C)
{
    struct Pair
    {
        int v1, v2;
    };
    static std::vector<int> graph[LARGE + 1];
    static int degree[LARGE + 1];
    std::vector<Pair> edges;
    memset(f, 0, sizeof(f));
    calc_f(1, A);
    calc_f(2, B);
    calc_f(3, C);
    int max = std::max({A, B, C});
    for (int i = 1; i <= max; i++)
    {
        graph[i].clear();
        degree[i] = 0;
    }
    for (int i = 1; i <= max; i++)
    {
        for (int j = 0; j < (1 << pfactor[i].size()); j++)
        {
            int set = j;
            int c = ((1 << pfactor[i].size()) - 1) & (~set);
            for (int k = 0; k <= j; k++)
            {
                if ((k & j) == k)
                {
                    int v1 = prod(j, i), v2 = prod(k | c, i);
                    edges.push_back({v1, v2});
                    degree[v1] += 1;
                    degree[v2] += 1;
                }
            }
        }
    }
    for (const auto &edge : edges)
    {
        if (degree[edge.v1] < degree[edge.v2] || (degree[edge.v1] == degree[edge.v2] && edge.v1 < edge.v2))
        {
            graph[edge.v1].push_back(edge.v2);
        }
    }

    int_t result = 0;
    ///三个点互异的三元环
    const auto f1 = [&](int a, int b, int c) -> int_t {
        int_t res = 0;
        int p = mu[a] * mu[b] * mu[c];
        int ab = lcm(a, b);
        int bc = lcm(b, c);
        int ca = lcm(c, a);
        //(ab)(bc)(ca)
        res = (res + 1ll * ::f[1][ab] * ::f[2][bc] * ::f[3][ca]) % mod;
        //(ab)(ca)(bc)
        res = (res + 1ll * ::f[1][ab] * ::f[2][ca] * ::f[3][bc]) % mod;
        //(bc)(ab)(ca)
        res = (res + 1ll * ::f[1][bc] * ::f[2][ab] * ::f[3][ca]) % mod;
        //(bc)(ca)(ab)
        res = (res + 1ll * ::f[1][bc] * ::f[2][ca] * ::f[3][ab]) % mod;
        //(ca)(ab)(bc)
        res = (res + 1ll * ::f[1][ca] * ::f[2][ab] * ::f[3][bc]) % mod;
        //(ca)(bc)(ab)
        res = (res + 1ll * ::f[1][ca] * ::f[2][bc] * ::f[3][ab]) % mod;
        return (res * p % mod + mod) % mod;
    };
    for (int i = 1; i <= max; i++)
    {
        auto &curr = graph[i];
        for (int j = 0; j < curr.size(); j++)
        {
            for (int k = j + 1; k < curr.size(); k++)
            {
                if (lcm(curr[j], curr[k]) <= max)
                    result = (result + f1(i, curr[j], curr[k])) % mod;
            }
        }
    }
    //有两个点相同的三元环
    const auto f2 = [&](int a, int b) -> int_t {
        // return 0;
        int_t res = 0;
        int p1 = mu[a] * mu[a] * mu[b];
        int p2 = mu[a] * mu[b] * mu[b];
        int ab = lcm(a, b);
        //(ba) (aa) (ab)
        res = (res + 1ll * p1 * ::f[1][ab] * ::f[2][a] * ::f[3][ab]) % mod;
        //(ab) (ba) (aa)
        res = (res + 1ll * p1 * ::f[1][ab] * ::f[2][ab] * ::f[3][a]) % mod;
        //(aa) (ab) (ba)
        res = (res + 1ll * p1 * ::f[1][a] * ::f[2][ab] * ::f[3][ab]) % mod;
        //(bb) (ba) (ab)
        res = (res + 1ll * p2 * ::f[1][b] * ::f[2][ab] * ::f[3][ab]) % mod;
        //(ba) (ab) (bb)
        res = (res + 1ll * p2 * ::f[1][ab] * ::f[2][ab] * ::f[3][b]) % mod;
        //(ab) (bb) (ba)
        res = (res + 1ll * p2 * ::f[1][ab] * ::f[2][b] * ::f[3][ab]) % mod;
        return (res % mod + mod) % mod;
    };
    for (int i = 1; i <= max; i++)
    {
        auto &curr = graph[i];
        for (int j = 0; j < curr.size(); j++)
        {
            result = (result + f2(i, curr[j])) % mod;
        }
    }
    //三个点都相同的三元环
    const auto f3 = [&](int x) -> int_t {
        return (1ll * mu[x] * mu[x] * mu[x] * f[1][x] * f[2][x] % mod * f[3][x] % mod + mod) % mod;
    };
    for (int_t i = 1; i <= max; i++)
    {
        if (mu[i] != 0)
            result = (result + f3(i)) % mod;
    }
    return result;
}
void init()
{

    for (int_t i = 1; i <= LARGE; i++)
    {
        factor[i] = i;
    }
    for (int_t i = 2; i <= LARGE; i++)
    {
        if (factor[i] == i)
        {
            for (int_t j = i * i; j <= LARGE; j += i)
            {
                factor[j] = i;
            }
        }
    }
    mu[1] = 1;
    for (int_t i = 2; i <= LARGE; i++)
    {
        int_t fact = factor[i];
        if (i / fact % fact == 0)
            mu[i] = 0;
        else
            mu[i] = -mu[i / fact];
    }
    for (int_t i = 2; i <= LARGE; i++)
    {
        if (mu[i] == 0)
            continue;
        int_t curr = i;
        while (curr != 1)
        {
            pfactor[i].push_back(factor[curr]);
            curr /= factor[curr];
        }
    }
}
int main()
{
    init();
    int_t T;
    cin >> T;
    for (int_t i = 1; i <= T; i++)
    {
        int_t A, B, C;
        cin >> A >> B >> C;
        cout << process(A, B, C) << endl;
    }
    return 0;
}

$$ \text{计算}\left( \begin{array}{c} n\\ m\\ \end{array} \right) \ mod\ p,\text{满足}p\le 10^6,n,m\le 10^{18},\text{不保证}p\text{是质数} \\ \text{首先可以考虑将}p\text{进行分解，得到}p_1^{k_1}p_2^{k_2}p_{3}^{k_3}…. \\ \text{然后考虑分别计算出}\left( \begin{array}{c} n\\ m\\ \end{array} \right) \ mod\ p_i^{k_i},\text{然后用}CRT\text{进行合并} \\ \text{现在的问题在于如何计算}\frac{n!}{m!\left( n-m \right) !}\ mod\ p_i^{k_i} \\ \text{因为}p_i^{k_i}\text{不是质数，所以不能用卢卡斯定理，但是我们考虑到}p_i^{k_i}\text{中只有}p_i\text{一个质因子，} \\ \text{所以我们可以把}n!\text{中所有的}p_i\text{提出来} \\ \text{例如现在计算11!对}3^2\text{取模} \\ \text{11!}=1\times 2\times 3\times 4\times 5\times 6\times 7\times 8\times 9\times 10\times \text{11\ } \\ =\left( 1\times 2 \right) \times \left( 3 \right) \times \left( 4\times 5 \right) \times \left( 2\times 3 \right) \times \left( 7\times 8 \right) \times \left( 3\times 3 \right) \times \left( 10\times 11 \right) \\ =3^3\times \left( 1\times 2\times 3 \right) \times \left( 1\times 2\times 4\times 5\times 7\times 8\times 10\times 11 \right) \\ =3^3\times \left( 1\times 2\times 3 \right) \times \left( 1\times 2\times 4\times 5\times 7\times 8\times 1\times 2 \right) \\ \text{对于}3^3,\text{额外记录下来，对于}\left( 1\times 2\times 3 \right) ,\text{递归计算,对于}\left( 1\times 2\times 4\times 5\times 7\times 8 \right) \times \left( 1\times 2 \right) ,\text{求出单个循环节和不足一个循环节的部分} \\ \text{每个循环节的长度是}p_i^{k_i}-p_i^{k_i-1},\text{一共循环了}\lfloor \frac{n-\lfloor \frac{n}{p_i} \rfloor}{p_i^{k_i}-p_i^{k_i-1}} \rfloor \text{次,不足一个循环节部分的长度是}\left( n-\lfloor \frac{n}{p_i} \rfloor \right) \ mod\left( p_i^{k_i}-p_i^{k_i-1} \right) \ \\ \text{对于不足一个循环节的部分暴力计算} \\ \text{最后可以得到11!中的3的幂次和除掉3之后的数在模}3^2\text{下的值} \\ \text{复杂度}O\left( p_i^{k_i}\log n \right) \\ \text{这样对于}\left( \begin{array}{c} n\\ m\\ \end{array} \right) =\frac{n!}{m!\left( n-m \right) !},\text{我们可以分别计算出}n!,m!,\left( n-m \right) !\text{中}p_i\text{的幂次和除掉}p_i\text{的幂次的部分，} \\ \text{然后就可以计算了} \\ \\ $$

#include <iostream>
#include <algorithm>
#include <utility>
#include <inttypes.h>
using int_t = int64_t;
using pair_t = std::pair<int_t, int_t>;

using std::cin;
using std::cout;
using std::endl;

const int_t LARGE = 1e6;
int_t power(int_t base, int_t index, int_t mod);
int_t gcd(int_t a, int_t b)
{
    if (b == 0)
        return a;
    return gcd(b, a % b);
}

pair_t exgcd(int_t a, int_t b)
{
    if (b == 0)
        return pair_t(1, 0);
    auto ret = exgcd(b, a % b);
    return pair_t(ret.second, ret.first - a / b * ret.second);
}

int_t excrt(int64_t *coef, int64_t *mod, int64_t n)
{
    int_t prex = coef[1], M = mod[1];
    for (int_t i = 2; i <= n; i++)
    {
        int_t A = M, B = mod[i];
        int_t C = coef[i] - prex;
        C = (C % B + B) % B;
        int_t gcd = ::gcd(A, B);
        auto x = exgcd(A, B).first;
        x = (x % (B / gcd) + (B / gcd)) % (B / gcd) * (C / gcd);
        int_t preM = M;
        M = M / gcd * mod[i];
        prex = ((prex + x * preM) % M + M) % M;
    }
    return prex;
}
//计算n!中,p的次数和除掉p的幂之后的数对p^k取模的值
void fact(int_t n, int_t p, int_t k, int_t *idx, int_t *rem, int_t *facts)
{
    if (n < p)
    {
        *idx = 0;
        *rem = facts[n];
        return;
    }
    fact(n / p, p, k, idx, rem, facts);
    *idx += n / p;
    int_t prod = 1, remaining = 1;
    int_t count = 0;
    const int_t mod = power(p, k, 998244353);
    int_t length = mod - mod / p;
    for (int_t i = 1; i < mod; i++)
    {
        if (i % p != 0)
        {
            prod = prod * i % mod;
            if (count < (n - n / p) % length)
            {
                remaining = remaining * i % mod;
                count++;
            }
        }
    }
    (*rem) = (*rem) * power(prod, (n - n / p) / (length), mod) % mod * remaining % mod;
}

//计算C(n,m)对p^k取模的结果
int_t C(int_t n, int_t m, int_t p, int_t k)
{
    static int_t fact[LARGE + 1];
    int_t mod = power(p, k, 998244353);
    const int_t phi = mod - (mod / p);
    fact[0] = 1;
    for (int_t i = 1; i < p; i++)
        fact[i] = fact[i - 1] * i % mod;
    int_t idx1, idx2, idx3, rem1, rem2, rem3;
    ::fact(n, p, k, &idx1, &rem1, fact);
    ::fact(m, p, k, &idx2, &rem2, fact);
    ::fact(n - m, p, k, &idx3, &rem3, fact);
    int_t index = idx1 - (idx2 + idx3);
    int_t remaining = rem1 * power(rem2 * rem3 % mod, phi - 1, mod) % mod;
    return power(p, index, mod) * remaining % mod;
}

int main()
{
    int_t n, m, p;
    cin >> n >> m >> p;
    static int_t index[LARGE + 1];
    static int_t prime[LARGE + 1];
    static int_t mod[LARGE + 1];
    static int_t coef[LARGE + 1];
    int_t used = 0;
    for (int_t i = 2; i <= p && p != 1; i++)
    {
        if (p % i == 0)
        {
            used++;
            prime[used] = i;
            while (p % i == 0)
            {
                index[used]++;
                p /= i;
            }
        }
    }
    for (int_t i = 1; i <= used; i++)
    {
        mod[i] = power(prime[i], index[i], 998244353);
        coef[i] = C(n, m, prime[i], index[i]);
    }
    cout << excrt(coef, mod, used) << endl;
    return 0;
}
int_t power(int_t base, int_t index, int_t mod)
{
    int_t result = 1;
    while (index)
    {
        if (index & 1)
            result = result * base % mod;
        base = base * base % mod;
        index >>= 1;
    }
    return result;
}