I used heavy light decompositon combined with the segment tree.

Finding the solution online seems like the easiest approach :(

Initial idea: we can do this with matrix exponentiation. If $$$d(m)$$$ is the number of divisors of $$$m$$$ and $$$f(i)$$$ is the $$$i$$$-th divisor of $$$m$$$, then make a $$$d(m)\times d(m)$$$ matrix $$$A$$$ where

Then the answer for length $k$ is $$$e^T A^{k-1} e$$$ for $$$e$$$ as the vector of all 1s. Of course, this is way too slow (it would be $$$n[d(m)]^2$$$ or something like this), but there are some optimizations.

Optimization: if two numbers have the same prime factors, the transitions from that prime factor are the same. For example, if $$$m=12$$$ then $$$6 = 2\cdot 3$$$ and $$$12 = 2^2 \cdot 3$$$ can both transition to $$$2,3,4,6,12$$$. By multiplying $$$2\times 3 \times \cdots \times 47 = 0.6\times10^{18}$$$ we see that there can be at most 15 distinct prime factors. Each divisor of $$$m$$$ will be added into a state based on the bitmask of which prime factors it contains. Unfortunately, this is still to much (since this would be $$$2^{15}$$$ states)

Optimization: Consider the prime factorization of $$$m = \prod_{i=1}^q p_i^{a_i}$$$. If prime $$$p_i$$$ has power $$$a_i$$$ in this expansion of $$$m$$$, then call it a type-$$$a_i$$$ prime. Two divisors (of $$$m$$$) called $$$s$$$ and $$$t$$$ will have the same number of transitions if they have the same number of type-$$$k$$$ primes for each $$$k$$$.

For example, consider $$$m=180 = 2^2 \cdot 3^2 \cdot 5$$$. The number of transitions away from 10 or 20 (which both have $$$2$$$ and $$$5$$$ as prime factors) is equal to the number of transitions away from 15 or 45 (which both have $$$3$$$ and $$$5$$$ as prime factors), because 10, 15, 20, and 45 all have one type-1 prime and one type-2 prime.

Why is this true? Consider some number $$$x = p_1^{f_1}p_2^{f_2}\cdots$$$ (where some $$$f_j$$$ might be 0) and the set $$$S$$$ of all possible values that can appear immediately after $$$x$$$. If $$$p_i$$$ and $$$p_j$$$ are both type-$$$q$$$ primes, then we can swap $$$p_i$$$ and $$$p_j$$$ in $$$x$$$ and in each of the elements of $$$S$$$, and it will still be true that $$$S$$$ is the set of all possible values that can appear after $$$x$$$.

Based on this we see that we can further compress our states: each state can be represented as a list $$$[b_1, b_2, \cdots, b_k]$$$ indicating that numbers of this type have $$$b_1$$$ type-1 primes, $$$b_2$$$ type-2 primes, etc. For the example $$$m=180$$$ there are the following possible states:

• $$$[b_1 = 0, b_2 = 0]$$$ (i.e. the number 1)
• $$$[0, 1]$$$: $$$(2, 3, 4, 9)$$$
• $$$[0, 2]$$$: $$$(6, 12, 18, 36)$$$
• $$$[1, 0]$$$: $$$(5)$$$
• $$$[1, 1]$$$: $$$(10, 15, 20, 45)$$$
• $$$[1, 2]$$$: $$$(30, 60, 90, 180)$$$

We can build a matrix $$$A$$$ such that $$$a_{ij}$$$ is the number of possible state-$$$j$$$ values which can follow a state-$$$i$$$ value. Also build an initial vector $$$v$$$ where $$$v_i$$$ is the number of values in state $$$i$$$. Then, the answer for length $$$k$$$ is $$$e^T A^{k-1} v$$$.

In general if $$$m$$$ has $$$c_i$$$ type-$$$i$$$ primes, then the number of states will be $$$\prod (1+c_i)$$$. We found that this can be at most 192, which is achieved by $$$m= 2^6 \cdot 3^5 \cdot 5^4 \cdot 7^3 \cdot 11^2 \cdot 13^2 \cdot 17 \cdot 19 \cdot 23 = 506480603789160000$$$. Here there are 3 type-1 primes, 2 type-3 primes, and 1 each of type-3, type-4, type-5, and type-6 primes.

Unfortunately, this is still too slow to naively implement $$$e^T A^{k-1} v$$$, since this would require $$$nq^2$$$ for $$$q = 192$$$, which is $$$\approx 7 \times 10^{11}$$$. There is an additional optimization.

Faster matrix evaluation:

The first idea for how to evaluate $$$e^T A^{k-1} v$$$ for $$$k\in [1,n]$$$ might be to store $$$e^T A^{k-1}$$$ (which is a vector) at the $$$k$$$-th step, and then multiply by $$$A$$$ for each value of $$$k$$$. To obtain the answer for the $$$k$$$-th step we can just multiply $$$(e^T A^{k-1}) v$$$ in $$$O(n)$$$.

Instead, we can use an idea similar to baby-step giant-step. Let $$$m = \left\lceil \sqrt{n}\right\rceil$$$. First, compute $$$e^T A^{j}$$$ for $$$j = 0,1,\cdots,m$$$. Then, compute $$$A^{j} v$$$ for $$$j = 0, m, 2m, \cdots, m^2$$$. Now, we can compute the $$$k$$$-th step as $$$(e^T A^{(k-1)\%m})(A^{\left\lfloor(k-1)/m\right\rfloor m} v)$$$ as the inner product of two vectors, which is $$$O(q)$$$ when the size of the vectors are $$$q$$$.

There's one remaining step, which is: how do we compute $$$A^{j} v$$$ for $$$j=0,m,2m,\cdots$$$ ? This can be done with matrix exponentiation. Compute $$$B = A^{m}$$$ using matrix exponentiation (which takes $$$O(q^3\log \sqrt{n})$$$ if $$$q$$$ is the size of the matrix).

In total the matrix operations (after computing the initial matrix $$$A$$$) have complexity:

• $$$O(q^3 \log \sqrt{n})$$$ for matrix exponentiation to get $$$A^{\sqrt{n}}$$$
• $$$O(q^2 \sqrt{n})$$$ for computing $$$e^T A^j$$$ and $$$A^j v$$$
• $$$O(qn)$$$ for computing the final values