Divisor Paths (Codeforces 1334E)
Tags: Greedy, combinatorics

You are given a positive integer $$$D$$$. Let's build the following graph from it:   each vertex is a divisor of $$$D$$$ (not necessarily prime, $$$1$$$ and $$$D$$$ itself are also included);  two vertices $$$x$$$ and $$$y$$$ ($$$x > y$$$) have an undirected edge between them if $$$x$$$ is divisible by $$$y$$$ and $$$\frac x y$$$ is a prime;  the weight of an edge is the number of divisors of $$$x$$$ that are not divisors of $$$y$$$. For example, here is the graph for $$$D=12$$$:   Edge $$$(4,12)$$$ has weight $$$3$$$ because $$$12$$$ has divisors $$$[1,2,3,4,6,12]$$$ and $$$4$$$ has divisors $$$[1,2,4]$$$. Thus, there are $$$3$$$ divisors of $$$12$$$ that are not divisors of $$$4$$$ — $$$[3,6,12]$$$.There is no edge between $$$3$$$ and $$$2$$$ because $$$3$$$ is not divisible by $$$2$$$. There is no edge between $$$12$$$ and $$$3$$$ because $$$\frac{12}{3}=4$$$ is not a prime.Let the length of the path between some vertices $$$v$$$ and $$$u$$$ in the graph be the total weight of edges on it. For example, path $$$[(1, 2), (2, 6), (6, 12), (12, 4), (4, 2), (2, 6)]$$$ has length $$$1+2+2+3+1+2=11$$$. The empty path has length $$$0$$$.So the shortest path between two vertices $$$v$$$ and $$$u$$$ is the path that has the minimal possible length.Two paths $$$a$$$ and $$$b$$$ are different if there is either a different number of edges in them or there is a position $$$i$$$ such that $$$a_i$$$ and $$$b_i$$$ are different edges.You are given $$$q$$$ queries of the following form:   $$$v$$$ $$$u$$$ — calculate the number of the shortest paths between vertices $$$v$$$ and $$$u$$$. The answer for each query might be large so print it modulo $$$998244353$$$.

Input: The first line contains a single integer $$$D$$$ ($$$1 \le D \le 10^{15}$$$) — the number the graph is built from.The second line contains a single integer $$$q$$$ ($$$1 \le q \le 3 \cdot 10^5$$$) — the number of queries.Each of the next $$$q$$$ lines contains two integers $$$v$$$ and $$$u$$$ ($$$1 \le v, u \le D$$$). It is guaranteed that $$$D$$$ is divisible by both $$$v$$$ and $$$u$$$ (both $$$v$$$ and $$$u$$$ are divisors of $$$D$$$).

Output: Print $$$q$$$ integers — for each query output the number of the shortest paths between the two given vertices modulo $$$998244353$$$.

Note: In the first example:   The first query is only the empty path — length $$$0$$$;  The second query are paths $$$[(12, 4), (4, 2), (2, 1)]$$$ (length $$$3+1+1=5$$$), $$$[(12, 6), (6, 2), (2, 1)]$$$ (length $$$2+2+1=5$$$) and $$$[(12, 6), (6, 3), (3, 1)]$$$ (length $$$2+2+1=5$$$).  The third query is only the path $$$[(3, 1), (1, 2), (2, 4)]$$$ (length $$$1+1+1=3$$$). 

Sample Inputs:
12
3
4 4
12 1
3 4

1
1
1 1

288807105787200
4
46 482955026400
12556830686400 897
414 12556830686400
4443186242880 325

Sample Outputs:
1
3
1

1

547558588
277147129
457421435
702277623

Solution (Python):
def gcd(a, b):
    while b:
        (a, b) = (b, a % b)
    return a

def isPrimeMR(n):
    d = n - 1
    d = d // (d & -d)
    L = [2, 3, 5, 7, 11, 13, 17]
    for a in L:
        t = d
        y = pow(a, t, n)
        if y == 1:
            continue
        while y != n - 1:
            y = y * y % n
            if y == 1 or t == n - 1:
                return 0
            t <<= 1
    return 1

def findFactorRho(n):
    m = 1 << n.bit_length() // 8
    for c in range(1, 99):
        f = lambda x: (x * x + c) % n
        (y, r, q, g) = (2, 1, 1, 1)
        while g == 1:
            x = y
            for i in range(r):
                y = f(y)
            k = 0
            while k < r and g == 1:
                ys = y
                for i in range(min(m, r - k)):
                    y = f(y)
                    q = q * abs(x - y) % n
                g = gcd(q, n)
                k += m
            r <<= 1
        if g == n:
            g = 1
            while g == 1:
                ys = f(ys)
                g = gcd(abs(x - ys), n)
        if g < n:
            if isPrimeMR(g):
                return g
            elif isPrimeMR(n // g):
                return n // g
            return findFactorRho(g)

def primeFactor(n):
    i = 2
    ret = {}
    rhoFlg = 0
    while i * i <= n:
        k = 0
        while n % i == 0:
            n //= i
            k += 1
        if k:
            ret[i] = k
        i += 1 + i % 2
        if i == 101 and n >= 2 ** 20:
            while n > 1:
                if isPrimeMR(n):
                    (ret[n], n) = (1, 1)
                else:
                    rhoFlg = 1
                    j = findFactorRho(n)
                    k = 0
                    while n % j == 0:
                        n //= j
                        k += 1
                    ret[j] = k
    if n > 1:
        ret[n] = 1
    if rhoFlg:
        ret = {x: ret[x] for x in sorted(ret)}
    return ret

def divisors(n):
    res = [1]
    prime = primeFactor(n)
    for p in prime:
        newres = []
        for d in res:
            for j in range(prime[p] + 1):
                newres.append(d * p ** j)
        res = newres
    res.sort()
    return res
mod = 998244353
import heapq

class Dijkstra:

    class Edge:

        def __init__(self, _to, _cost):
            self.to = _to
            self.cost = _cost

    def __init__(self, V):
        self.G = [[] for i in range(V)]
        self._E = 0
        self._V = V

    @property
    def E(self):
        return self._E

    @property
    def V(self):
        return self._V

    def add(self, _from, _to, _cost):
        self.G[_from].append(self.Edge(_to, _cost))

    def shortest_path(self, s):
        import heapq
        que = []
        d = [10 ** 15] * self.V
        d[s] = 0
        cnt = [0] * self.V
        cnt[s] = 1
        heapq.heappush(que, (0, s))
        while len(que) != 0:
            (cost, v) = heapq.heappop(que)
            if d[v] < cost:
                continue
            for i in range(len(self.G[v])):
                e = self.G[v][i]
                if d[e.to] > d[v] + e.cost:
                    d[e.to] = d[v] + e.cost
                    heapq.heappush(que, (d[e.to], e.to))
                    cnt[e.to] = cnt[v] % mod
                elif d[e.to] == d[v] + e.cost:
                    cnt[e.to] += cnt[v]
                    cnt[e.to] %= mod
        return (d, cnt)
import sys
from heapq import heappush, heappop, heapify
from math import gcd
input = sys.stdin.buffer.readline
D = int(input())
prime = primeFactor(D)
div = divisors(D)
div_dic = {e: i for (i, e) in enumerate(div)}
div_factors = [1 for i in range(len(div))]
n = len(div)
for i in range(n):
    d = div[i]
    for p in prime:
        tmp_d = d
        cnt = 0
        while tmp_d % p == 0:
            cnt += 1
            tmp_d //= p
        div_factors[i] *= cnt + 1
G = Dijkstra(n)
for i in range(n):
    for p in prime:
        nv = div[i] * p
        if nv in div_dic:
            j = div_dic[nv]
            G.add(i, j, div_factors[j] - div_factors[i])
            G.add(j, i, div_factors[j] - div_factors[i])
base = G.shortest_path(0)[1]

def dist(a, b):
    assert b % a == 0
    b //= a
    return base[div_dic[b]]
for _ in range(int(input())):
    (u, v) = map(int, input().split())
    g = gcd(u, v)
    print(dist(g, u) * dist(g, v) % mod)

Solution (Python):
from collections import Counter

def gcd(x, y):
    """"""
    while y:
        (x, y) = (y, x % y)
    return x

def memodict(f):
    """"""

    class memodict(dict):

        def __missing__(self, key):
            ret = self[key] = f(key)
            return ret
    return memodict().__getitem__

def pollard_rho(n):
    """"""
    if n & 1 == 0:
        return 2
    if n % 3 == 0:
        return 3
    s = (n - 1 & 1 - n).bit_length() - 1
    d = n >> s
    for a in [2, 325, 9375, 28178, 450775, 9780504, 1795265022]:
        p = pow(a, d, n)
        if p == 1 or p == n - 1 or a % n == 0:
            continue
        for _ in range(s):
            prev = p
            p = p * p % n
            if p == 1:
                return gcd(prev - 1, n)
            if p == n - 1:
                break
        else:
            for i in range(2, n):
                (x, y) = (i, (i * i + 1) % n)
                f = gcd(abs(x - y), n)
                while f == 1:
                    (x, y) = ((x * x + 1) % n, (y * y + 1) % n)
                    y = (y * y + 1) % n
                    f = gcd(abs(x - y), n)
                if f != n:
                    return f
    return n

@memodict
def prime_factors(n):
    """"""
    if n <= 1:
        return Counter()
    f = pollard_rho(n)
    return Counter([n]) if f == n else prime_factors(f) + prime_factors(n // f)

def distinct_factors(n, primes):
    """"""
    factors = [1]
    for (p, exp) in primes.items():
        factors += [p ** i * factor for factor in factors for i in range(1, exp + 1)]
    return factors

def all_factors(n):
    """"""
    (small, large) = ([], [])
    for i in range(1, int(n ** 0.5) + 1, 2 if n & 1 else 1):
        if not n % i:
            small.append(i)
            large.append(n // i)
    if small[-1] == large[-1]:
        large.pop()
    large.reverse()
    small.extend(large)
    return small
import os
import sys
from io import BytesIO, IOBase
BUFSIZE = 8192

class FastIO(IOBase):
    newlines = 0

    def __init__(self, file):
        self._fd = file.fileno()
        self.buffer = BytesIO()
        self.writable = 'x' in file.mode or 'r' not in file.mode
        self.write = self.buffer.write if self.writable else None

    def read(self):
        while True:
            b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
            if not b:
                break
            ptr = self.buffer.tell()
            (self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
        self.newlines = 0
        return self.buffer.read()

    def readline(self):
        while self.newlines == 0:
            b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
            self.newlines = b.count(b'\n') + (not b)
            ptr = self.buffer.tell()
            (self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
        self.newlines -= 1
        return self.buffer.readline()

    def flush(self):
        if self.writable:
            os.write(self._fd, self.buffer.getvalue())
            (self.buffer.truncate(0), self.buffer.seek(0))

class IOWrapper(IOBase):

    def __init__(self, file):
        self.buffer = FastIO(file)
        self.flush = self.buffer.flush
        self.writable = self.buffer.writable
        self.write = lambda s: self.buffer.write(s.encode('ascii'))
        self.read = lambda : self.buffer.read().decode('ascii')
        self.readline = lambda : self.buffer.readline().decode('ascii')
(sys.stdin, sys.stdout) = (IOWrapper(sys.stdin), IOWrapper(sys.stdout))
input = lambda : sys.stdin.readline().rstrip('\r\n')
from math import gcd, ceil

def prod(a, mod=10 ** 9 + 7):
    ans = 1
    for each in a:
        ans = ans * each % mod
    return ans

def lcm(a, b):
    return a * b // gcd(a, b)

def binary(x, length=16):
    y = bin(x)[2:]
    return y if len(y) >= length else '0' * (length - len(y)) + y
for _ in range(int(input()) if not True else 1):
    n = int(input())
    mod = 998244353
    primes = prime_factors(n)
    all = set(distinct_factors(n, primes))
    primes = [k for k in primes]
    factorial = [1]
    invfact = [1]
    for i in range(1, 10 ** 5 + 69):
        factorial += [factorial[-1] * i % mod]
        invfact += [pow(factorial[i], mod - 2, mod)]

    def f(x, y):
        count = 0
        x = y // x
        denominator = 1
        for prime in primes:
            uniques = 0
            while x % prime == 0:
                x = x // prime
                count += 1
                uniques += 1
            if uniques:
                denominator *= invfact[uniques]
                denominator %= mod
        return factorial[count] * denominator % mod
    for i in range(int(input())):
        (x, y) = map(int, input().split())
        (x, y) = (max(x, y), min(x, y))
        if x % y == 0:
            print(f(y, x))
        else:
            print(f(gcd(x, y), y) * f(gcd(x, y), x) % mod)

Solution (Python):
mod = 998244353
eps = 10 ** (-9)

def main():
    import sys
    input = sys.stdin.buffer.readline
    nmax = 100
    fac = [0] * nmax
    finv = [0] * nmax
    inv = [0] * nmax
    fac[0] = 1
    fac[1] = 1
    finv[0] = 1
    finv[1] = 1
    inv[1] = 1
    for i in range(2, nmax):
        fac[i] = fac[i - 1] * i % mod
        inv[i] = mod - inv[mod % i] * (mod // i) % mod
        finv[i] = finv[i - 1] * inv[i] % mod

    def comb(n, r):
        if n < r:
            return 0
        else:
            return fac[n] * (finv[r] * finv[n - r] % mod) % mod

    def PrimeDecomposition(N):
        ret = {}
        n = int(N ** 0.5)
        for d in range(2, n + 1):
            while N % d == 0:
                if d not in ret:
                    ret[d] = 1
                else:
                    ret[d] += 1
                N //= d
            if N == 1:
                break
        if N != 1:
            ret[N] = 1
        return ret
    D = int(input())
    PD = PrimeDecomposition(D)
    P = list(PD.keys())
    np = len(P)
    ans_all = []
    for _ in range(int(input())):
        (u, v) = map(int, input().split())
        pos = []
        neg = []
        for p in P:
            up = 0
            while u % p == 0:
                up += 1
                u //= p
            vp = 0
            while v % p == 0:
                vp += 1
                v //= p
            if up > vp:
                pos.append(up - vp)
            elif up < vp:
                neg.append(vp - up)
        ans = 1
        if pos:
            ans = ans * fac[sum(pos)] % mod
            for x in pos:
                ans = ans * finv[x] % mod
        if neg:
            ans = ans * fac[sum(neg)] % mod
            for x in neg:
                ans = ans * finv[x] % mod
        ans_all.append(ans)
    print(*ans_all, sep='\n')
if __name__ == '__main__':
    main()