Solution for https://codeforces.com/contest/1334/problem/E

sol.cc

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using ull = unsigned long long;
#define FOR(i, j, k, in) for (int i = j; i < k; i += in)
#define RFOR(i, j, k, in) for (int i = j; i >= k; i -= in)
#define REP(i, j) FOR(i, 0, j, 1)
#define RREP(i, j) RFOR(i, j, 0, 1)

namespace prime {

using u64 = uint64_t;
using u32 = uint32_t;
// using u128 = __uint128_t;

u64 binpower(u64 base, u64 e, u64 mod) {
    u64 result = 1;
    base %= mod;
    while (e) {
        if (e & 1)
            result = (u64)result * base % mod;
        base = (u64)base * base % mod;
        e >>= 1;
    }
    return result;
}

bool check_composite(u64 n, u64 a, u64 d, int s) {
    u64 x = binpower(a, d, n);
    if (x == 1 || x == n - 1)
        return false;
    for (int r = 1; r < s; r++) {
        x = (u64)x * x % n;
        if (x == n - 1)
            return false;
    }
    return true;
};

bool miller_rabin(u64 n, int iter = 5) {  // returns true if n is probably prime, else returns false.
    if (n < 4)
        return n == 2 || n == 3;

    int s = 0;
    u64 d = n - 1;
    while ((d & 1) == 0) {
        d >>= 1;
        s++;
    }

    for (int i = 0; i < iter; i++) {
        int a = 2 + rand() % (n - 3);
        if (check_composite(n, a, d, s))
            return false;
    }
    return true;
}

vector<u32> primes, lp;

void init_primes(u64 N) {
    lp.resize(N + 1, 0), primes.clear();
    for (u64 i = 2; i <= N; ++i) {
        if (lp[i] == 0) {
            lp[i] = i;
            primes.push_back(i);
        }
        for (int j = 0; j < (int)primes.size() && primes[j] <= lp[i] && i * primes[j] <= N; ++j)
            lp[i * primes[j]] = primes[j];
    }
}

unordered_map<u64, u64> factorize(u64 x) {
    unordered_map<u64, u64> ret;
    if (x > lp.size()) {
        u64 pc = primes.size();
        for (u64 j = 1; j < pc + 1; j++) {
            while (x % primes[pc - j] == 0) {
                ret[primes[pc - j]]++;
                x /= primes[pc - j];
            }
            if (x <= lp.size()) break;
        }
    }
    if (x > lp.size()) {
        ret[x]++;
        return ret;
    }
    while (x != 1) {
        ret[lp[x]]++;
        x /= lp[x];
    }

    return ret;
}

void _rec_divisors(const vector<pair<u64, u64>> &factors, int index, u64 curr, vector<u64> &ret) {
    int n = factors.size();
    if (index == n) {
        ret.push_back(curr);
        return;
    }

    u64 base = factors[index].first, exp = factors[index].second;
    for (u64 i = 0; i <= exp; i++) {
        _rec_divisors(factors, index + 1, curr, ret);
        curr *= base;
    }
}

vector<u64> divisors(u64 x) {
    if (x == 1) return {1};
    auto factor_map = factorize(x);
    vector<pair<u64, u64>> factors;
    for (const auto &el : factor_map) {
        factors.emplace_back(el.first, el.second);
    }
    vector<u64> ret;
    _rec_divisors(factors, 0, 1, ret);
    return ret;
}

}  // namespace prime

namespace modulo {
using ull = unsigned long long;

vector<ull> fact, inv_fact;
ull mod;

ull pow(ull base, ull exp) {
    ull ret = 1;
    while (exp) {
        if (exp % 2 == 1) ret = (ret * base) % mod;
        base = (base * base) % mod;
        exp >>= 1;
    }

    return ret;
}

ull invert(ull x) {
    return pow(x, mod - 2);
}

void init_factorials(ull n, ull p) {
    fact.resize(n + 1);
    inv_fact.resize(n + 1);
    mod = p;
    fact[0] = 1;
    for (ull i = 1; i <= n; i++) fact[i] = (i * fact[i - 1]) % mod;
    inv_fact[n] = invert(fact[n]);
    assert((inv_fact[n] * fact[n]) % mod == 1);
    for (ull i = n - 1; i >= 1; i--) inv_fact[i] = ((i + 1) * inv_fact[i + 1]) % mod;
}

ull multi_comb(const ull nom, const vector<ull> &denoms) {
    ull ret = fact[nom];
    for (const auto denom : denoms) {
        ret = (ret * inv_fact[denom]) % mod;
    }

    return ret;
}
}  // namespace modulo

constexpr prime::u64 p = 998244353, limit = 31622777;

void quick_factor(const vector<ull> &prime_factors, vector<ull> &exps, ull x) {
    ull pfs = prime_factors.size();
    REP(i, pfs) {
        auto prime_factor = prime_factors[i];
        while (x % prime_factor == 0) {
            exps[i]++;
            x /= prime_factor;
        }
        if (x == 1) return;
    }
}

int main() {
    prime::init_primes(limit);
    modulo::init_factorials(200000, p);
    prime::u64 d;
    scanf("%llu\n", &d);
    auto d_factors = prime::factorize(d);
    vector<ull> prime_factors;
    for (const auto &el : d_factors) prime_factors.push_back(el.first);
    sort(begin(prime_factors), end(prime_factors));
    ull pfs = prime_factors.size();
    // cin >> d;
    ull q;
    scanf("%llu\n", &q);
    // cin >> q;
    REP(i, q) {
        ull one, other;
        scanf("%llu %llu\n", &one, &other);
        // cin >> one >> other;
        if (one == other) {
            printf("1\n");
            // cout << 1 << endl;
            continue;
        }
        auto gcd = __gcd(one, other);
        vector<ull> one_exps(pfs), other_exps(pfs), gcd_exps(pfs);
        quick_factor(prime_factors, one_exps, one);
        quick_factor(prime_factors, other_exps, other);
        quick_factor(prime_factors, gcd_exps, gcd);
        ull res;
        if (min(one, other) == gcd) {
            ull nom = 0;
            vector<ull> denoms;
            if (other > one) {
                REP(i, pfs) {
                    auto diff = other_exps[i] - one_exps[i];
                    if (diff == 0) continue;
                    nom += diff;
                    denoms.push_back(diff);
                }
            } else {
                REP(i, pfs) {
                    auto diff = one_exps[i] - other_exps[i];
                    if (diff == 0) continue;
                    nom += diff;
                    denoms.push_back(diff);
                }
            }
            res = modulo::multi_comb(nom, denoms);
        } else {
            ull nom = 0;
            vector<ull> denoms;
            REP(i, pfs) {
                auto diff = other_exps[i] - gcd_exps[i];
                if (diff == 0) continue;
                nom += diff;
                denoms.push_back(diff);
            }
            res = modulo::multi_comb(nom, denoms);
            nom = 0;
            denoms.clear();
            REP(i, pfs) {
                auto diff = one_exps[i] - gcd_exps[i];
                if (diff == 0) continue;
                nom += diff;
                denoms.push_back(diff);
            }
            res = (res * modulo::multi_comb(nom, denoms)) % p;
        }

        printf("%llu\n", res);
        // cout << res << endl;
    }

    return 0;
}