Euclidean-like Algorithm

Introduction

The Euclidean-like algorithm was proposed by Hong Huadun during the 2016 Winter Camp camper exchange. It is commonly used to solve summation problems for sequences of the form

for sequences of this form (indexed by $i$). Its main idea is to use the recursive structure of fractions to reduce the problem to a smaller one and solve it recursively. Because the recursive structure of fractions has a direct connection with the Euclidean algorithm and a direct link to it in continued-fraction expansion, this summation method is also called the Euclidean-like algorithm.

Methods such as continued fractions and the Stern–Brocot tree also characterize the recursive structure of fractions, so problems solvable with the Euclidean-like algorithm can usually be solved with these methods as well. Compared with them, the Euclidean-like algorithm is generally easier to understand and its implementation is more concise.

Euclidean-like Algorithm

The simplest example is the following summation problem:

where $a,b,c,n$ are all positive integers.

Algebraic Solution

First, taking $a,b$ modulo $c$ simplifies the problem and reduces it to the case $0\le a,b<c$:

Now consider the transformed problem. Let

Then the original problem can be written as a double sum:

To reverse the order of summation, we need to compute, for each $j$, the range of $i$ satisfying the condition. To do so, transform the condition:

The properties of the floor and ceiling functions were used repeatedly during these transformations. Substituting the transformed condition, the original expression becomes:

Let $(a',b',c',n')=(c,c-b-1,a,m-1)$; this returns us to the case $a'>c'$ discussed above.

Combining these two transformations, we see that $(a,c)$ are repeatedly reduced modulo one another and swapped until $a=0$. This resembles applying the Euclidean algorithm to $(a,c)$, which is why the method is called the Euclidean-like algorithm. Its time complexity is $O(\log\min\{a,c\})$.

During the computation, the case $m=0$ may occur, in which case the inner recursion has $n=-1$. This does not affect the final result. However, if the algorithm is terminated immediately when $m=0$, its time complexity can be improved to $O(\log\min\{a,c,n\})$.

The similarity between this algorithm and the Euclidean algorithm makes it easy to show that its time complexity is $O(\log\min\{a,c\})$. Therefore, it only remains to show that if the algorithm terminates when $m=0$, its time complexity is also $O(\log n)$.

Let $m=\lfloor(an+b)/c\rfloor$, and write $S=mn$ and $k=m/n$. These correspond respectively to the area of the lattice-point diagram and the slope of the line in the geometric interpretation (see the next section). For sufficiently large $n$, we approximately have $k\doteq a/c$.

Consider how $S$ and $k$ change during the algorithm. In the first modulo operation, $n$ remains unchanged, and $k$ approximately changes from $a/c$ to $(a\bmod c)/c$, which corresponds to changing the slope from $k$ to $k-\lfloor k\rfloor$; $S$ also approximately becomes $(k-\lfloor k\rfloor)$ times its original value. In the second operation, which swaps the horizontal and vertical coordinates, $S$ remains approximately unchanged while $k$ becomes its reciprocal. Thus, if the pair $(k,S)$ becomes $(k',S')$ after the two operations, then $k'=(k-\lfloor k\rfloor)^{-1}$ and $S'=(k-\lfloor k\rfloor)S$.

Because $1\le\lfloor k'\rfloor\le k'<\lfloor k'\rfloor+1$, after two rounds of recursive computation, the minimum shrinkage factor of the product is

$$(k'-\lfloor k'\rfloor)(k-\lfloor k\rfloor) = 1-\dfrac{\lfloor k'\rfloor}{k'} < 1-\dfrac{\lfloor k'\rfloor}{\lfloor k'\rfloor+1} = \dfrac{1}{\lfloor k'\rfloor+1}\le \dfrac{1}{2}.$$

Therefore, the algorithm must terminate after at most $O(\log S)$ rounds. Starting from the second round, the $S$ at the beginning of each round is never greater than the $S$ after the modulo operation in the preceding round; the latter is approximately $kn^2$ with $k<1$, so $O(\log S)\subseteq O(\log n)$. This proves the claim above.

The reference implementation for the template problem is as follows:

#include <iostream>

long long solve(long long a, long long b, long long c, long long n) {
  long long n2 = n * (n + 1) / 2;
  if (a >= c || b >= c)
    return solve(a % c, b % c, c, n) + (a / c) * n2 + (b / c) * (n + 1);
  long long m = (a * n + b) / c;
  if (!m) return 0;
  return m * n - solve(c, c - b - 1, a, m - 1);
}

int main() {
  int t;
  std::cin >> t;
  for (; t; --t) {
    int a, b, c, n;
    std::cin >> n >> c >> a >> b;
    std::cout << solve(a, b, c, n - 1) << '\n';
  }
  return 0;
}

Geometric Intuition

This algorithm can also be understood geometrically. The main problems solvable by the Euclidean-like algorithm are lattice-point counting problems below a line.

As shown in the leftmost part of the figure below, the sum is equivalent to counting the lattice points below the line

above the $x$-axis (excluding the $x$-axis), whose horizontal coordinates lie in $[0,n]$.

First remove the integer parts of the slope and intercept. This is equivalent to separately counting the blue points in the middle part of the figure above. When the slope and intercept are both integers, the blue points form a trapezoidal array: the lattice points in different columns form arithmetic progressions, so their number is easy to compute. After removing these points, the remaining lattice points match the number of red points in the rightmost part of the figure. The problem is reduced to the case where both the slope and intercept are less than one. Since the trapezoid has height $n+1$ and its two bases have lengths $\lfloor b/c\rfloor$ and $(\lfloor a/c\rfloor n+\lfloor b/c\rfloor)$, the trapezoid area formula gives

Next, flip the horizontal and vertical axes. As shown in the leftmost part of the figure below, the red and blue points form a rectangular lattice of horizontal length $n$ and vertical length $m=\lfloor(an+b)/c\rfloor$. To count the red points, count the blue points and subtract that number from the number of points in the rectangular lattice. After the flip, the blue lattice in the left half of the figure becomes the red lattice below a line. Moreover, after the flip, the slope is greater than one, returning to the case handled above.

The key is how to calculate the equation of the line above the new red lattice. Flipping the coordinate axes in the leftmost part of the figure gives the middle part. The line above the flipped red lattice (the solid line in the middle) does not correspond to the original line (the solid line on the left); instead, it is obtained by shifting the original line slightly to the upper left (the dashed line on the left). This is because directly flipping the line (the solid line on the left) produces the dashed line in the middle, but by definition the lattice points below it include those exactly on the line, causing points on the line to be counted twice. To avoid this, after flipping $y=(ax+b)/c$ to obtain $y=(cx-b)/a$, shift the resulting line slightly downward to obtain $y=(cx-b-1)/a$; the lattice below it is then exactly the blue lattice before the flip.

There is one more detail to handle. The intercept of the line in the middle part of the figure is negative, so we have not yet returned to the original situation. To make the intercept nonnegative, shift the line (the solid line in the middle) one unit to the left. This does not omit any lattice point because the blue lattice before the flip contains no point with vertical coordinate zero, so after the flip there is no point with horizontal coordinate zero either. The line equation finally becomes $y=(cx+c-b-1)/a$; at the same time, the upper bound of the lattice's horizontal coordinate changes from $m$ to $m-1$. This step can be summarized as

This recursive algorithm works mainly for two reasons:

• First, the slope of the line repeatedly takes its fractional part and then its reciprocal, which is equivalent to computing the continued-fraction expansion of the line's slope $k=a/c$. Since the continued-fraction expansion of a rational fraction has length $O(\log\min\{a,c\})$, this process must terminate after $O(\log\min\{a,c\})$ steps;
• Second, because the slope of the line is less than one whenever the axes are flipped, intuition suggests that $m<n$; that is, the range of horizontal coordinates keeps shrinking after each such iteration. The complexity analysis above rigorously shows that after every two iterations, $n$ is at most half its previous value, so this process must terminate after $O(\log n)$ steps.

This is also why the Euclidean-like algorithm has complexity $O(\log\min\{a,c,n\})$ when the slope is rational.

The Euclidean-like algorithm can be extended to the case of an irrational slope using similar geometric intuition; see the examples later in the article for the detailed analysis.

Examples

#include <iostream>

struct Data {
  int f, g, h;
};

Data solve(long long a, long long b, long long c, long long n) {
  constexpr long long M = 998244353;
  constexpr long long i2 = (M + 1) / 2;
  constexpr long long i6 = (M + 1) / 6;
  long long n2 = (n + 1) * n % M * i2 % M;
  long long n3 = (2 * n + 1) * (n + 1) % M * n % M * i6 % M;
  Data res = {0, 0, 0};
  if (a >= c || b >= c) {
    auto tmp = solve(a % c, b % c, c, n);
    long long aa = a / c, bb = b / c;
    res.f = (tmp.f + aa * n2 + bb * (n + 1)) % M;
    res.g = (tmp.g + aa * n3 + bb * n2) % M;
    res.h = (tmp.h + 2 * bb * tmp.f % M + 2 * aa * tmp.g % M +
             aa * aa % M * n3 % M + bb * bb % M * (n + 1) % M +
             2 * aa * bb % M * n2 % M) %
            M;
    return res;
  }
  long long m = (a * n + b) / c;
  if (!m) return res;
  auto tmp = solve(c, c - b - 1, a, m - 1);
  res.f = (m * n - tmp.f + M) % M;
  res.g = (m * n2 + (M - tmp.f) * i2 + (M - tmp.h) * i2) % M;
  res.h = (n * m % M * m - tmp.f - tmp.g * 2 + 3 * M) % M;
  return res;
}

int main() {
  int t;
  std::cin >> t;
  for (; t; --t) {
    int n, a, b, c;
    std::cin >> n >> a >> b >> c;
    auto res = solve(a, b, c, n);
    std::cout << res.f << ' ' << res.h << ' ' << res.g << '\n';
  }
  return 0;
}

To apply the Euclidean-like algorithm, first transform the sum into a familiar form:

$$\begin{aligned} \sum_{d=1}^n(-1)^{\lfloor d\sqrt{r}\rfloor} &= \sum_{d=1}^n\left(1 - 2(\lfloor d\sqrt{r}\rfloor\bmod 2)\right)\\ &= n - 2\sum_{d=1}^n\left(\lfloor d\sqrt{r}\rfloor - 2\left\lfloor\dfrac{\lfloor d\sqrt{r}\rfloor}{2}\right\rfloor\right) \\ &= n - 2\sum_{d=1}^n\lfloor d\sqrt{r}\rfloor + 4 \sum_{d=1}^n\left\lfloor\dfrac{d\sqrt{r}}{2}\right\rfloor\\ &= n - 2f(n,1,0,1) + 4f(n,1,0,2) \end{aligned}$$

The function $f$ has the form

$$f(a,b,c,n) = \sum_{i=1}^n\left\lfloor\dfrac{a\sqrt{r}+b}{c}i\right\rfloor.$$

Unlike the algorithm in the main text, the slope here is no longer rational. Let the slope be

$$k = \dfrac{a\sqrt{r}+b}{c}.$$

We again discuss two cases. If $k\ge 1$, then

$$\begin{aligned} f(a,b,c,n) &= \sum_{i=1}^n \lfloor ki\rfloor = \sum_{i=1}^n \lfloor(k-\lfloor k\rfloor)i\rfloor + \lfloor k\rfloor \sum_{i=1}^ni\\ &= \lfloor k\rfloor\dfrac{n(n+1)}{2} + f(a,b-c\lfloor k\rfloor,c,n). \end{aligned}$$

The problem is reduced to the case where the slope is less than one. If $k<1$, let $m=\lfloor nk\rfloor$. Then

$$\begin{aligned} f(a,b,c,n) &= \sum_{i=1}^n \lfloor ki\rfloor = \sum_{i=1}^n\sum_{j=1}^m[j\le\lfloor ki\rfloor]\\ &= \sum_{j=1}^m\sum_{i=1}^n[i>\lfloor k^{-1}j\rfloor] = nm - \sum_{j=1}^m\sum_{i=1}^n[i\le\lfloor k^{-1}j\rfloor]. \end{aligned}$$

In this derivation, the condition for exchanging $i$ and $j$ is simpler than in the main text because the line $y=kx$ has no lattice points other than the origin. The key is to write the exchanged sum in the form of $f(a,b,c,n)$, which requires $a',b',c'$ to satisfy

$$k^{-1} = \dfrac{a'\sqrt{r}+b'}{c'}.$$

This is not difficult: rationalizing the denominator gives

$$k^{-1} = \dfrac{c}{a\sqrt{r}+b} = \dfrac{ca\sqrt{r}-cb}{a^2r-b^2}.$$

Therefore,

$$a'=ca,~b'=-cb,~c'=a^2r-b^2.$$

This gives

$$f(a,b,c,n) = nm - f(ca,-cb,a^2r-b^2,m).$$

To avoid integer overflow, divide $a,b,c$ by their greatest common divisor at every step. Because this computation is exactly the same as computing the continued fraction of $k$, the theory of continued fractions shows that, as long as $\gcd(a,b,c)=1$, they necessarily remain within the integer range during the computation. In addition, although $(a,b,c,n)$ does not overflow, under this problem's data limits $f(a,b,c,n)$ may exceed the range of $64$-bit integers. Natural overflow is sufficient and no extra handling is needed; the final result is always in $[-n,n]$.

Although the slope never becomes zero, the algorithm still has complexity $O(\log n)$, as can be seen from the preceding complexity argument.

#include <cmath>
#include <iostream>

long long r;
long double sqrt_r;

long long gcd(long long a, long long b) { return b ? gcd(b, a % b) : a; }

unsigned long long f(long long a, long long b, long long c, long long n) {
  if (!n) return 0;
  auto d = gcd(a, gcd(b, c));
  a /= d;
  b /= d;
  c /= d;
  unsigned long long k = (a * sqrt_r + b) / c;
  if (k) {
    return n * (n + 1) / 2 * k + f(a, b - c * k, c, n);
  } else {
    unsigned long long m = n * (a * sqrt_r + b) / c;
    return n * m - f(c * a, -c * b, a * a * r - b * b, m);
  }
}

unsigned long long solve(long long n, long long r) {
  long long sqr = sqrt_r = sqrtl(r);
  if (r == sqr * sqr) return r % 2 ? (n % 2 ? -1 : 0) : n;
  return n - 2 * f(1, 0, 1, n) + 4 * f(1, 0, 2, n);
}

int main() {
  int t;
  std::cin >> t;
  for (; t; --t) {
    int n;
    std::cin >> n >> r;
    long long res = solve(n, r);
    std::cout << res << '\n';
  }
  return 0;
}

#include <iostream>

void solve(int a, int b, int& p, int& q, int c, int d) {
  if ((a / b + 1) * d < c) {
    p = a / b + 1;
    q = 1;
  } else {
    solve(d, c - d * (a / b), q, p, b, a % b);
    p += q * (a / b);
  }
}

int main() {
  int a, b, c, d, p, q;
  while (std::cin >> a >> b >> c >> d) {
    solve(a, b, p, q, c, d);
    std::cout << p << '/' << q << '\n';
  }
  return 0;
}

Universal Euclidean Algorithm

The derivation of the Euclidean-like algorithm discussed in the previous section is usually rather cumbersome, and the sums it can solve are mainly those reducible to (weighted) lattice-point counting below a line. This section presents a more general method that abstracts the process above and can solve more problems. It is therefore also called the universal Euclidean algorithm. It likewise uses the recursive structure of fractions, but its approach to reducing the problem differs slightly from that of the Euclidean-like algorithm.

We still consider the most classical summation problem:

where $a,b,c,n$ are all positive integers.

Problem Transformation

Let the line with parameters $(a,b,c,n)$ be

For this line segment, define a string $S$ consisting of $U$ and $R$, also called the operation sequence, as follows:

• The string consists of exactly $n$ $R$'s and $m=\lfloor(an+b)/c\rfloor$ $U$'s;
• The number of $U$'s before the $i$-th $R$ is exactly $\lfloor(ai+b)/c\rfloor$, where $i=1,\cdots,n$.

Geometrically, this is roughly equivalent to starting at the origin: write an $R$ each time a vertical grid line is crossed to the right, and write a $U$ each time a horizontal grid line is crossed upward. As shown below:

Naturally, this definition must account for several special cases:

• When an integer lattice point is passed (that is, an upward crossing and a rightward crossing occur simultaneously), write $U$ before $R$;
• At the start of the string, in addition to the number of grid lines crossed upward in the interval $(0,1]$, add $\lfloor b/c\rfloor$ extra $U$'s;
• At the end of the string, there must be no extra $U$'s.

If any part of the geometric description is unclear, refer to the algebraic definition above for help. The geometric description helps explain the algorithm below.

The basic idea of the universal Euclidean algorithm is to regard the $U$ and $R$ in the operation sequence as elements of a certain monoid, regard the entire operation sequence as the product of elements in that monoid, and relate the final answer to this product.

For example, in this problem, define a state vector $v = (1,y,\sum y)$ to represent the current state after crossing some horizontal and vertical grid lines from the origin. The first component is a constant, the second is the vertical coordinate $y$, and the third is the sum being sought. Initially, $v=(1,0,0)$. Each upward crossing increases the vertical coordinate by one, which is equivalent to right-multiplying the state vector by the matrix

Each rightward crossing adds the current vertical coordinate to the sum, which is equivalent to right-multiplying the state vector by the matrix

Therefore, the final state is the product $(1,0,0)S$, where $S$ is understood as the product of the matrices above. The answer is the third component of the final state.

Besides defining elements of the monoid as matrices, they can also be defined as the contribution of an operation sequence to the final result, with the product of operations defined as the combination of the contributions of two operation sequences.

For this problem, define the contribution of each operation sequence as $(x,y,\sum y)$. To explain this notation rigorously, regard each component as a function of the operation sequence: for an operation sequence $S$, its contribution is $(x(S),y(S),(\sum y)(S))$. Here, $x(S)$ and $y(S)$ are respectively the numbers of $R$ and $U$ in $S$, that is, the numbers of rightward and upward grid-line crossings. In general, the summation symbol in the last term is defined as follows: for a function $f(S)$ on operation sequences, define $(\sum f)(S)$, or write it as $\sum_S f$, by the expression below:

Here, $S_r$ is the $r$-th character of $S$, and $S_{[1,r]}$ is the prefix consisting of the first $r$ characters of $S$. In other words, this summation notation can be viewed as summing over all prefixes of $S$ that end in $R$. For example,

Also, $\sum y$ is the accumulated number of previous upward grid-line crossings each time a rightward grid line is crossed. For the entire operation sequence, the values of $y$ at all prefixes ending in $R$ are the values of $\lfloor(ai+b)/c\rfloor$ for $i=1,\cdots,n$. Therefore, $\sum y$ computed for the entire operation sequence is exactly the quantity required by the problem.

Initially, $U=(0,1,0)$ and $R=(1,0,0)$. We can further define the product of two elements $(x_1,y_1,s_1)$ and $(x_2,y_2,s_2)$ as

The result of combining the contribution in the last term can be obtained as follows:

It is easy to verify that this multiplication is associative and that the identity is $(0,0,0)$, so these elements form a monoid under this multiplication. The answer is the third component of the product.

Both methods produce the correct result. However, because the matrix method retains more redundant information and therefore has larger constant factors, the second method is more practical for actual problems.

Algorithm

Unlike the global reduction of the Euclidean-like algorithm, the universal Euclidean algorithm reduces the problem by combining these operations in batches. Denote the product of the operations corresponding to the string by

The reduction proceeds as follows:

• When $b\ge c$, the operation sequence begins with $\lfloor b/c\rfloor$ $U$'s. Compute their product directly and remove these $U$'s from the sequence. The number of $U$'s before the $i$-th $R$ is then

  $$\left\lfloor\dfrac{ai+b}{c}\right\rfloor - \left\lfloor\dfrac{b}{c}\right\rfloor = \left\lfloor\dfrac{ai+(b\bmod c)}{c}\right\rfloor.$$

  Therefore, this is equivalent to changing the line parameters from $(a,b,c,n)$ to $(a,b\bmod c,c,n)$. Thus, in this case,

  $$F(a,b,c,n,U,R) = U^{\lfloor b/c\rfloor}F(a,b\bmod c,c,n,U,R).$$

• When $a\ge c$, there are at least $\lfloor a/c\rfloor$ $U$'s before every $R$ in the operation sequence, so they can be merged into $R$. In other words, replace $R$ by $U^{\lfloor a/c\rfloor}R$. The number of $U$'s before the $i$-th $R$ in the merged string is

  $$\left\lfloor\dfrac{ai+b}{c}\right\rfloor - \left\lfloor\dfrac{a}{c}\right\rfloor i = \left\lfloor\dfrac{(a\bmod c)i+b}{c}\right\rfloor.$$

  Therefore, this is equivalent to changing the line parameters from $(a,b,c,n)$ to $(a\bmod c,b,c,n)$. Thus, in this case,

  $$F(a,b,c,n,U,R) = F(a\bmod c,b,c,n,U,U^{\lfloor a/c\rfloor}R).$$

• For the remaining case, flip the horizontal and vertical coordinates. This essentially swaps $U$ and $R$, but the parameters of the flipped line must be calculated carefully. From the definition of the operation sequence, determine coefficients $(a',b',c',n')$ such that the number of $R$'s before the $j$-th $U$ in the original operation sequence is $\lfloor(a'j+b')/c'\rfloor$, and there are $n'$ $U$'s in total. By definition,

  $$n'=\left\lfloor\dfrac{an+b}{c}\right\rfloor = m,$$

  The number of $R$'s before the $j$-th $U$ is the largest $i$ such that

  $$\begin{aligned} \left\lfloor\dfrac{ai+b}{c}\right\rfloor < j &\iff \dfrac{ai+b}{c} < j \iff i < \dfrac{cj-b}{a} \\ &\iff i < \left\lceil\dfrac{cj-b}{a}\right\rceil = \left\lfloor\dfrac{cj-b - 1}{a}\right\rfloor + 1. \end{aligned}$$

  Therefore, $i = \lfloor(cj-b-1)/a\rfloor$. This derivation is similar to that of the Euclidean-like algorithm above and likewise uses properties of floor and ceiling functions.

  There are two details to handle:

  • The intercept term $-(b+1)/a$ is negative. If the line is shifted one unit to the left, the intercept becomes nonnegative because $(c-b-1)/a\ge 0$. Therefore, extract the initial segment $R^{\lfloor(c-b-1)/a\rfloor}U$ before the swap, and swap only the remaining $U$'s and $R$'s;
  • After swapping $U$ and $R$, there are extra $U$'s at the end. Therefore, before swapping $U$ and $R$, first extract the final segment of $R$'s, and swap only the remaining $U$'s and $R$'s. The number of $R$'s in this segment is $n-\lfloor(cm-b-1)/a\rfloor$.

  After removing some characters from the beginning and end, the number of $R$'s before the $j$-th $U$ becomes:

  $$\left\lfloor\dfrac{c(j+1)-b-1}{a}\right\rfloor - \left\lfloor\dfrac{c-b-1}{a}\right\rfloor = \left\lfloor\dfrac{cj+(c-b-1)\bmod a}{a}\right\rfloor.$$

  Recall that the number of $U$'s before the swap is $m = \lfloor(an+b)/c\rfloor$. The one-unit left shift above requires at least one $U$ before the swap, that is, $m>0$. Under this condition, there are two cases:

  • When $m>0$, after handling the two points above, the operation sequence obtained by swapping $U$ and $R$ is the valid sequence corresponding to parameters $(c,(c-b-1)\bmod a,a,m-1)$. Thus,

    $$F(a,b,c,n,U,R) = R^{\lfloor(c-b-1)/a\rfloor}UF(c,(c-b-1)\bmod a,a,m-1,R,U)R^{n-\lfloor(cm-b-1)/a\rfloor}.$$

  • In the special case $m=0$, the operation sequence before the swap contains only $n$ $R$'s, so no swap is needed and we can return directly:

    $$F(a,b,c,n,U,R) = R^n.$$

    Unlike the Euclidean-like algorithm, this special case of the universal Euclidean algorithm must be handled separately; otherwise negative exponents would make the computation incorrect.

These observations allow the problem to be solved recursively.

Assume that multiplying two monoid elements takes $O(1)$ time. If powers of these elements are computed by binary exponentiation, the resulting algorithm has complexity $O(\log\max\{a,c\}+\log(b/c))$.¹

The universal Euclidean algorithm can be written as a unified template; when applying it to a specific problem, only the implementation of the template type T needs to be changed.

// Class T implements the monoid.
// Assume that it provides a multiplication operator
//     and a default constructor returning the unity in the monoid.

// Binary exponentiation.
template <typename T>
T pow(T a, int b) {
  T res;
  for (; b; b >>= 1) {
    if (b & 1) res = res * a;
    a = a * a;
  }
  return res;
}

// Universal Euclidean algorithm.
template <typename T>
T euclid(int a, int b, int c, int n, T U, T R) {
  if (b >= c) return pow(U, b / c) * euclid(a, b % c, c, n, U, R);
  if (a >= c) return euclid(a % c, b, c, n, U, pow(U, a / c) * R);
  auto m = ((long long)a * n + b) / c;
  if (!m) return pow(R, n);
  return pow(R, (c - b - 1) / a) * U *
         euclid(c, (c - b - 1) % a, a, m - 1, R, U) *
         pow(R, n - (c * m - b - 1) / a);
}

The implementation of the template problem obtained using the universal Euclidean algorithm is as follows:

#include <array>
#include <iostream>

// Switch between matrix and info merging approaches.
#define MATRIX 1

// --8<-- [start:euclidean]
// Class T implements the monoid.
// Assume that it provides a multiplication operator
//     and a default constructor returning the unity in the monoid.

// Binary exponentiation.
template <typename T>
T pow(T a, int b) {
  T res;
  for (; b; b >>= 1) {
    if (b & 1) res = res * a;
    a = a * a;
  }
  return res;
}

// Universal Euclidean algorithm.
template <typename T>
T euclid(int a, int b, int c, int n, T U, T R) {
  if (b >= c) return pow(U, b / c) * euclid(a, b % c, c, n, U, R);
  if (a >= c) return euclid(a % c, b, c, n, U, pow(U, a / c) * R);
  auto m = ((long long)a * n + b) / c;
  if (!m) return pow(R, n);
  return pow(R, (c - b - 1) / a) * U *
         euclid(c, (c - b - 1) % a, a, m - 1, R, U) *
         pow(R, n - (c * m - b - 1) / a);
}

// --8<-- [end:euclidean]
#if MATRIX

template <size_t N>
struct Matrix {
  std::array<long long, N * N> mat;

  auto loc(size_t i, size_t j) const { return mat[i * N + j]; }

  auto& loc(size_t i, size_t j) { return mat[i * N + j]; }

  Matrix() : mat{} {
    for (size_t i = 0; i != N; ++i) {
      loc(i, i) = 1;
    }
  }

  Matrix operator*(const Matrix& rhs) const {
    Matrix res;
    res.mat.fill(0);
    for (size_t i = 0; i != N; ++i) {
      for (size_t k = 0; k != N; ++k) {
        for (size_t j = 0; j != N; ++j) {
          res.loc(i, j) += loc(i, k) * rhs.loc(k, j);
        }
      }
    }
    return res;
  }
};

long long solve(int a, int b, int c, int n) {
  if (!n) return 0;
  Matrix<3> U, R;
  U.loc(0, 1) = R.loc(1, 2) = 1;
  auto res = euclid(a, b, c, n, U, R);
  return res.loc(0, 2);
}

#else

struct Info {
  long long x, y, s;

  Info() : x(0), y(0), s(0) {}

  Info operator*(const Info& rhs) const {
    Info res;
    res.x = x + rhs.x;
    res.y = y + rhs.y;
    res.s = s + rhs.s + rhs.x * y;
    return res;
  }
};

long long solve(int a, int b, int c, int n) {
  if (!n) return 0;
  Info U, R;
  U.y = 1;
  R.x = 1;
  auto res = euclid(a, b, c, n, U, R);
  return res.s;
}

#endif

int main() {
  int t;
  std::cin >> t;
  for (; t; --t) {
    int a, b, c, n;
    std::cin >> n >> c >> a >> b;
    std::cout << solve(a, b, c, n - 1) << '\n';
  }
  return 0;
}

Examples