First, the number of operations replacing \(A\) with \(f(A)\), \(N\), is not very important. Let us assume that sufficiently many operations will be done. From now on, we do not refer to the number of operations.

---

◆ The coefficient sequence

The inserting operations change the sequence \(A\) as follows.

• \(a, b\)
• \(a, a+b, b\)
• \(a, 2a+b, a+b, a+2b, b\)
• \(a, 3a+b, 2a+b, 3a+2b, a+b, 2a+3b, a+2b, a+3b, b\)
• \(\vdots\)

As seen here, every term is in the form \(xa + yb\), so let us consider the sequences of coefficients \((x, y)\) for these elements.

• \((1, 0), (0, 1)\)
• \((1, 0), (1,1), (0, 1)\)
• \((1, 0), (2,1), (1,1), (1,2), (0, 1)\)
• \((1, 0), (3,1), (2,1), (3,2), (1,1), (2,3), (1,2), (1,3), (0, 1)\)
• \(\vdots\)

Below, we call these sequences of pairs coefficient sequences.

The properties of the coefficient sequence

We can show the following.

[Property (1)] If \((x_1, y_1)\) and \((x_2,y_2)\) are adjacent in a coefficient sequence after some number of inserting operations, \(x_1y_2 - x_2y_1 = 1\).

We can easily prove it by induction. Conversely, we can also show the following.

[Property (2)] If \(x_1y_2 - x_2y_1 = 1\) holds for \(0\leqq x_1,y_1,x_2,y_2\), \((x_1, y_1)\) and \((x_2,y_2)\) will be adjacent in a coefficient sequence after some number of inserting operations.

We can easily prove it by, for example, induction on \(x_2 + y_1\). If \(x_2 + y_1 > 0\) and \(x_1y_2 - x_2y_1 = 1\), we can prove that one of the following holds.

• \(x_1\leqq x_2\) and \(y_1\leqq y_2\);
• \(x_1\geqq x_2\) and \(y_1\geqq y_2\).

In the former case, if we let \((x_3, y_3) = (x_2 - x_1, y_2 - y_1)\), we have \(x_3, y_3\geqq 0\) and \(x_1y_3 - x_3y_1 = 1\). From the inductive hypothesis, these mean that \((x_1,y_1)\) and \((x_3,y_3)\) will be adjacent in a coefficient sequence after some number of inserting operations. Then, the subsequent operation will insert \((x_2, y_2)\) between them. Similarly goes the proof for the latter case.

From these properties, we can also see the following.

[Property (3)] A pair \((x,y)\) of non-negative integers such that \(\gcd(x,y) = 1\) will appear in a coefficient sequence after some number of inserting operations.

[Property (4)] In a coefficient sequence at any point, pairs \((x,y)\) will be in ascending order of \(y/x\). (For \((x,y) = (0,1)\), \(y/x\) is treated as \infty$.)

---

◆ The sequence \(B\) in other words

In the end, the sequence \(B\) turns out to be the following.

Among the pairs of non-negative integers \(x, y\) such that \(ax + by\leqq N\), consider the ones that satisfy \(\gcd(x,y) = 1\). Sort these pairs \((x,y)\) in ascending order of \(y/x\). Then, \(B_i = ax_i + by_i\) holds, where \((x_i, y_i)\) is the \(i\)-th pair.

---

◆ Computing \(B_n\)

Let us first solve the partial problem of computing one \(B_n\). In combination with binary search, it will be reduced to the following counting problem.

Given are positive integers \(a, b, N\) and a rational number \(c\). Count the pairs of non-negative integers \((x,y)\neq (0,0)\) that satisfy all of the following conditions:

• \(ax + by\leqq N\);
• \(\gcd(x,y) = 1\);
• \(y/x \leqq c\).

Let us fix \(a, b, c\) and let \(f(N)\) denote the answer to the above counting problem. Additionally, let \(F(N)\) denote the answer when ignoring the condition \(\gcd(x, y) = 1\). By classifying pairs \((x,y)\) according to \(\gcd(x,y) = d\), we have \(F(N) = \sum_{d=1}^N f(\lfloor N/d\rfloor)\). (We have omitted the pair \((x,y) = (0,0)\) to make this classification simpler.)

From the Möbius inversion formula, we have \(f(N) = \sum_{d=1}^N \mu(d)F(\lfloor N/d \rfloor)\), using which we can compute \(f(N)\).

Since we can compute \(F(1), F(2), \ldots, F(N)\) in a total of \(O(N)\) time using accumulated sums and so on, we can compute \(f(N)\) in \(O(N)\) time assuming that the Möbius function is precomputed.

By combining it with binary search, we have found that we can compute one term \(B_n\) in \(O(N\log N)\) time.

---

As a side note, assuming that the Möbius function is precomputed, we can also compute \(B_n\) in \(O(N^{1/2}\log^2 N)\) time, by computing each \(F(n)\) in \(O(\log n)\) time using [AtCoder Library] floor_sum.

---

◆ Computing \(B_L, \ldots, B_R\)

Considering how the sequence \(B\) is constructed, when \(B_{n}\) and \(B_{n+1}\) are known, we can find \(B_{n+2}\) in \(O(1)\) time. Thus, if we just compute \(B_L\) and \(B_{L+1}\), we can compute the remaining \(B_i\) in a total of \(O(R-L)\), solving the problem.

We can also solve the problem by jointly listing the \(L\)-th through \(R\)-th coefficients \((x, y)\) by combining binary searches on \(L\) and \(R\), and then sorting these coefficients.

---

◆ Reference: The connection to the Farey sequences

A Farey sequence is the sequence of completely reduced fractions between \(0\) and \(1\) whose denominators are at most \(N\), in ascending order. Below are some of these sequences.

• For \(n = 1\): \(\frac01, \frac11\)
• For \(n = 2\): \(\frac01, \frac12, \frac11\)
• For \(n = 3\): \(\frac01, \frac13, \frac12, \frac23, \frac11\)
• For \(n = 4\): \(\frac01, \frac14, \frac13, \frac12, \frac23, \frac34, \frac11\)
• For \(n = 5\): \(\frac01, \frac15, \frac14, \frac13, \frac25, \frac12, \frac35, \frac23, \frac34, \frac45, \frac11\)

It is known that the Farey sequences can be calculated as follows. (We can apply our arguments in “The properties of the coefficient sequence” almost as is.)

• Begin with \(\frac{0}{1}, \frac{1}{1}\).
• Repeatedly insert \(\frac{y_1+y_2}{x_1+x_2}\) between \(\frac{y_1}{x_1}\) and \(\frac{y_2}{x_2}\).

This method is extremely similar to the computation of coefficient sequences in this article. Actually, we can associate the coefficient sequences with the Farey sequences by making a pair of non-negative integers \((x_1,y_1)\) correspond to a fraction \(\frac{x_1}{x_1+y_1}\).

---

We will show another way to see the problem, with more awareness of the connection to the Farey sequences. First, since we can easily reduce it to the case \(\gcd(a, b) = 1\), below we assume that it holds.

From \(\gcd(a,b) = 1\), we can take non-negative integers \(a', b'\) such that \(a'b - ab' = 1\).

Instead of \((a, b)\), let us consider two fractions \((a'/a, b'/b)\). Then, if by inserting \((y_1+y_2)/(x_1+x_2)\) between \(y_1/x_1\) and \(y_2/x_2\), the inserting operation in the problem can be identified with the insertion in the method to calculate the Farey sequences. Thus, the sequence \(B\) in the problem can be rephrased into the following.

• Arrange in ascending order the completely reduced fractions \(y/x\) between \(a'/a\) and \(b'/b\) whose denominators are at most \(N\). Extracting just the denominators \(x\) from this will produce the sequence \(B\).

For example, if \(a = 2, b = 1, N = 6\), we have the following.

• The sequence \(B\): \((2, 5, 3, 4, 5, 1) \)
• The list of completely reduced fractions between \(1/2\) and \(1/1\) whose denominators are at most \(6\), sorted in ascending order: \(\left(\frac12, \frac25, \frac23, \frac34, \frac45, \frac11\right)\)

In the end, it has turned out that the problem is equivalent to the following.

Of the list of completely reduced fractions between \(\frac{a'}{a}\) and \(\frac{b'}{b}\) whose denominators are at most \(N\), sorted in ascending order, what are the denominators of the \(L\)-th through \(R\)-th terms?