Euclidean-like Algorithm

Introduction

The Euclidean-like algorithm was introduced by Hong Huadun during the 2016 Winter Camp participant exchange. It is commonly used to solve summation problems involving sequences (indexed by \(i\)) of the form

\[ \left\lfloor\dfrac{ai+b}{c}\right\rfloor \]

The main idea is to utilize the recursive structure of fractions to reduce the problem to smaller subproblems and solve them recursively. Because the recursive structure of fractions is directly related to the Euclidean algorithm, this summation method is also called the Euclidean-like algorithm.

Since methods like continued fractions and the Stern–Brocot tree also describe the recursive structure of fractions, problems solvable by the Euclidean-like algorithm can usually also be solved by these methods. Compared to them, the Euclidean-like algorithm is often easier to understand and its implementation is more concise.

Euclidean-like Algorithm

The simplest example is the summation problem:

\[ f(a,b,c,n)=\sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor, \]

where \(a,b,c,n\) are positive integers.

Algebraic Solution

First, reduce \(a\) and \(b\) modulo \(c\) to simplify the problem, transforming it into the case where \(0\le a,b<c\):

\[ \begin{aligned} f(a,b,c,n)&=\sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor\\ &=\sum_{i=0}^n\left\lfloor \frac{\left(\left\lfloor\frac{a}{c}\right\rfloor c+(a\bmod c)\right)i+\left(\left\lfloor\frac{b}{c}\right\rfloor c+(b\bmod c)\right)}{c}\right\rfloor\\ &=\sum_{i=0}^n\left(\left\lfloor\frac{a}{c}\right\rfloor i+\left\lfloor\frac{b}{c}\right\rfloor+\left\lfloor\frac{\left(a\bmod c\right)i+\left(b\bmod c\right)}{c} \right\rfloor\right)\\ &=\frac{n(n+1)}{2}\left\lfloor\frac{a}{c}\right\rfloor +(n+1)\left\lfloor\frac{b}{c}\right\rfloor+f(a\bmod c,b\bmod c,c,n). \end{aligned} \]

Now consider the transformed problem. Let

\[ m = \left\lfloor \frac{an+b}{c} \right\rfloor. \]

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

\[ \sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor =\sum_{i=0}^n\sum_{j=0}^{m-1}\left[j<\left\lfloor \frac{ai+b}{c} \right\rfloor\right]. \]

Swapping the summation order requires determining the range of \(i\) satisfying the condition for each \(j\). To do this, transform the condition:

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

The transformation repeatedly uses properties of the floor function. Substituting the transformed condition, the expression becomes:

\[ \begin{aligned} f(a,b,c,n)&=\sum_{j=0}^{m-1} \sum_{i=0}^n\left[i>\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor \right]\\ &=\sum_{j=0}^{m-1}\left(n-\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor\right)\\ &=nm-f\left(c,c-b-1,a,m-1\right). \end{aligned} \]

Let \((a',b',c',n')=(c,c-b-1,a,m-1)\). This returns to the previously discussed case where \(a'>c'\).

Combining these two transformations, we observe that during the process, \((a,c)\) repeatedly swap positions after taking modulo, until \(a=0\). This resembles performing the Euclidean algorithm on \((a,c)\), hence the name "Euclidean-like algorithm". Its time complexity is \(O(\log\min\{a,c\})\).

During computation, the case \(m=0\) may occur, causing the inner recursion to have \(n=-1\). This does not affect the final result. However, if we require the algorithm to terminate immediately when \(m=0\), the time complexity can be improved to \(O(\log\min\{a,c,n\})\).

Explanation of Complexity

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

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

Examine the changes in \(S\) and \(k\) during the algorithm. In the first modulo step, \(n\) remains unchanged, and \(k\) approximately changes from \(a/c\) to \((a\bmod c)/c\), equivalent to the slope changing from \(k\) to \(k-\lfloor k\rfloor\), while \(S\) also approximately becomes \((k-\lfloor k\rfloor)\) times its original value. In the second step swapping coordinates, \(S\) approximately remains unchanged, while \(k\) becomes its reciprocal. Therefore, if after two steps the pair \((k,S)\) becomes \((k',S')\), we have \(k'=(k-\lfloor k\rfloor)^{-1}\) and \(S'=(k-\lfloor k\rfloor)S\).

Since \(1\le\lfloor k'\rfloor\le k'<\lfloor k'\rfloor+1\), after two rounds of recursion, the product reduction factor is at least

\[ (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}. \]

Thus, after at most \(O(\log S)\) rounds, the algorithm must terminate. Because from the second round onward, the \(S\) at the start of each round is at most the \(S\) after the modulo step of the previous round, and the latter is roughly \(kn^2\) with \(k<1\), so \(O(\log S)\subseteq O(\log n)\). This yields the above conclusion.

A reference implementation for the template problem is as follows:

Template Problem Implementation (Library Checker - Sum of Floor of Linear)

#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 Interpretation

This algorithm can also be understood geometrically. The Euclidean-like algorithm mainly solves problems of counting lattice points under a line.

As shown in the leftmost part of the figure below, the summation corresponds to counting lattice points below the line

\[ y = \dfrac{ax+b}{c} \]

and above the x-axis (excluding the x-axis), with x-coordinates in \([0,n]\).

First, remove the integer parts of the slope and intercept. This step corresponds to separately counting the number of blue points in the middle part of the figure. When both slope and intercept are integers, the blue points necessarily form a trapezoidal array; that is, lattice points in different vertical columns form an arithmetic progression, so their count is easy to compute. After removing these points, the remaining lattice points correspond to the red points in the rightmost part of the figure. The problem reduces to the case where both slope and intercept are less than one. Since the trapezoid has height \(n+1\) and base lengths \(\lfloor b/c\rfloor\) and \((\lfloor a/c\rfloor n+\lfloor b/c\rfloor)\), using the trapezoid area formula, this step can be summarized as

\[ f(a,b,c,n) = f(a\bmod c,b\bmod c,c,n) + \dfrac{1}{2}(n+1)\left(\left\lfloor\dfrac{b}{c}\right\rfloor+\left(\left\lfloor\dfrac{a}{c}\right\rfloor n+\left\lfloor\dfrac{b}{c}\right\rfloor\right)\right). \]

Then, swap the coordinate axes. As shown in the leftmost part of the figure below, the red and blue points form a rectangular lattice with horizontal length \(n\) and vertical length \(m=\lfloor(an+b)/c\rfloor\). To count the red points, we can count the blue points and subtract from the total lattice points. After swapping, the blue lattice in the left part becomes a red lattice under some line. Moreover, after swapping, the slope becomes greater than one, returning to the case already handled earlier.

The key is to compute the equation of the line above the new red lattice. Swapping the axes of the leftmost part yields the middle part. The line above the red lattice after swapping (solid line in the middle) does not correspond to the line before swapping (solid line on the left), but to the line before swapping shifted slightly upward and leftward (dashed line on the left). This is because directly swapping the line (solid line on the left) yields the dashed line in the middle, but by definition, lattice points exactly on the line would be counted twice. To avoid this, we need to shift the swapped line \(y=(ax+b)/c\) to \(y=(cx-b)/a\) downward slightly to \(y=(cx-b-1)/a\), so that the lattice below it exactly corresponds to the blue lattice before swapping.

Another detail needs handling. The intercept of the line in the middle part is negative, meaning we haven't returned to the initial case. To make the intercept non-negative, we simply shift the line (solid line in the middle) left by one unit. This does not miss any lattice points because the blue lattice before swapping had no points with y-coordinate zero, so after swapping there are no points with x-coordinate zero. Finally, the line equation becomes \(y=(cx+c-b-1)/a\); simultaneously, the upper bound of the x-coordinate changes from \(m\) to \(m-1\). This step can be summarized as

\[ f(a,b,c,n) = mn - f(c,c-b-1,a,m-1). \]

This recursive algorithm works for two main reasons:

• First, the slope of the line repeatedly takes the fractional part and then the reciprocal, which is equivalent to computing the continued fraction expansion of the slope \(k=a/c\). Since the length of the continued fraction expansion of a rational number is \(O(\log\min\{a,c\})\), this process must terminate after \(O(\log\min\{a,c\})\) steps.
• Second, because each time we swap axes, the slope is less than one, intuitively we should have \(m<n\), meaning after each iteration the x-coordinate range shrinks. The earlier complexity analysis rigorously shows that after every two iterations, \(n\) is at most halved, so this process must terminate after \(O(\log n)\) steps.

This is also why the complexity of the Euclidean-like algorithm for rational slopes is \(O(\log\min\{a,c,n\})\).

Using similar geometric intuition, the Euclidean-like algorithm can be extended to irrational slopes; see the examples later for details.

Examples

【Template】Euclidean-like Algorithm

Multiple queries. Given positive integers \(a,b,c,n\), compute

\[ \begin{aligned} f(a,b,c,n) &= \sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor,\\ g(a,b,c,n) &= \sum_{i=0}^ni\left\lfloor \frac{ai+b}{c} \right\rfloor,\\ h(a,b,c,n) &= \sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor^2. \end{aligned} \]

Solution 1

Similar to the derivation for \(f\), we can obtain recursive expressions for \(g\) and \(h\).

First, use modulo to reduce to the case \(0\le a,b<c\):

\[ \begin{aligned} g(a,b,c,n) &=g(a\bmod c,b\bmod c,c,n)+\left\lfloor\frac{a}{c}\right\rfloor\frac{n(n+1)(2n+1)}{6}+\left\lfloor\frac{b}{c}\right\rfloor\frac{n(n+1)}{2}, \\ h(a,b,c,n)&=h(a\bmod c,b\bmod c,c,n)\\ &\quad+2\left\lfloor\frac{b}{c}\right\rfloor f(a\bmod c,b\bmod c,c,n) +2\left\lfloor\frac{a}{c}\right\rfloor g(a\bmod c,b\bmod c,c,n)\\ &\quad+\left\lfloor\frac{a}{c}\right\rfloor^2\frac{n(n+1)(2n+1)}{6}+\left\lfloor\frac{b}{c}\right\rfloor^2(n+1) +\left\lfloor\frac{a}{c}\right\rfloor\left\lfloor\frac{b}{c}\right\rfloor n(n+1). \end{aligned} \]

Then, swap summation order to further transform. Similarly, let

\[ m = \left\lfloor \frac{an+b}{c} \right\rfloor. \]

For the sum \(g\), we have

\[ \begin{aligned} g(a,b,c,n)&=\sum_{i=0}^ni\left\lfloor \frac{ai+b}{c} \right\rfloor\\ &=\sum_{i=0}^n \sum_{j=0}^{m-1}i \left[j<\left\lfloor\frac{ai+b}{c}\right\rfloor\right] \\ &=\sum_{j=0}^{m-1}\sum_{i=0}^n i\left[i>\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor \right]\\ &=\sum_{j=0}^{m-1}\dfrac{1}{2}\left(\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor+n+1\right)\left(n-\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor\right)\\ &=\dfrac{1}{2}mn(n+1) - \dfrac{1}{2}\sum_{j=0}^{m-1}\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor - \dfrac{1}{2}\sum_{j=0}^{m-1}\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor^2\\ &=\dfrac{1}{2}mn(n+1) - \dfrac{1}{2}f(c,c-b-1,a,m-1) - \dfrac{1}{2}h(c,c-b-1,a,m-1). \end{aligned} \]

For the sum \(h\), we have

\[ \begin{aligned} h(a,b,c,n)&=\sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor^2\\ &=\sum_{i=0}^n\sum_{j=0}^{m-1}(2j+1)\left[j<\left\lfloor\frac{ai+b}{c}\right\rfloor\right]\\ &=\sum_{j=0}^{m-1}\sum_{i=0}^n(2j+1)\left[i>\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor \right]\\ &=\sum_{j=0}^{m-1}(2j+1)\left(n-\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor\right)\\ &=nm^2 - \sum_{j=0}^{m-1}\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor - 2\sum_{j=0}^{m-1}j\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor\\ &=nm^2 - f(c,c-b-1,a,m-1) - 2g(c,c-b-1,a,m-1). \end{aligned} \]

From a geometric perspective, these nonlinear sums correspond to assigning weights \(w(i,j)\) to each point \((i,j)\) in the region. Apart from these weights, the rest of the computation process is identical. For weight choices, generally

\[ \sum_{i=0}^ni^r\left\lfloor \frac{ai+b}{c} \right\rfloor^s = \sum_{i=0}^n\sum_{j=0}^{m-1} i^r\left((j+1)^s-j^s\right)\left[j<\left\lfloor\frac{ai+b}{c}\right\rfloor\right]. \]

Another feature of this problem is that \(g\) and \(h\) intertwine during recursive computation. Therefore, we need to recursively compute the triple \((f,g,h)\) simultaneously.

#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;
}

[Tsinghua Training 2014] Sum

Multiple queries. Given positive integers \(n\) and \(r\), compute

\[ \sum_{d=1}^n(-1)^{\lfloor d\sqrt{r}\rfloor}. \]

Solution 1

If \(r\) is a perfect square, then when \(\sqrt{r}\) is even, the sum is \(n\); otherwise, the sum alternates between \(0\) and \(-1\) depending on the parity of \(n\). Now consider the case where \(r\) is not a perfect square.

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} \]

where 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}. \]

Similarly, 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 reduces 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, swapping \(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 swapped sum in the form \(f(a,b,c,n)\), which requires finding \(a',b',c'\) such that

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

This is not difficult; rationalizing the denominator yields

\[ 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 implies

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

To avoid integer overflow, we need to divide \(a,b,c\) by their greatest common divisor each time. Since this computation process is identical to computing the continued fraction of \(k\), according to continued fraction theory, as long as \(\gcd(a,b,c)=1\), they will remain within integer range during computation. Additionally, although \((a,b,c,n)\) will not overflow, under the problem's constraints, \(f(a,b,c,n)\) may exceed 64-bit integer range. Natural overflow suffices; no extra handling is needed, and the final result will be within \([-n,n]\).

Although the slope never becomes zero, the algorithm's complexity remains \(O(\log n)\), as can be seen from the earlier 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;
}

Fraction

Given positive integers \(a,b,c,d\), find the fraction \(p/q\) with minimal \((q,p)\) in lexicographic order among all reduced fractions satisfying \(a/b<p/q<c/d\).

Solution

This problem is also a classic application of the Stern–Brocot tree; a related solution can be found here. Because it only depends on the recursive structure of fractions, it can also be solved using a method similar to the Euclidean algorithm, hence it can be considered an application of the Euclidean-like algorithm.

If there exists at least one integer between \(a/b\) and \(c/d\) (excluding endpoints), we can directly take \((q,p)=(1,\lfloor a/b\rfloor+1)\). Otherwise, we must have

\[ \left\lfloor\dfrac{a}{b}\right\rfloor \le \dfrac{a}{b} <\dfrac{p}{q} <\dfrac{c}{d}\le\left\lfloor\dfrac{a}{b}\right\rfloor+1. \]

From this inequality, we can determine that the integer part of \(p/q\) is \(\lfloor a/b\rfloor\). Remove this integer part, then take the reciprocal to determine the fractional part. This is exactly the basic method for determining the continued fraction of \(p/q\). If the final answer is \(p/q\), the algorithm's time complexity is \(O(\log\min\{p,q\})\).

Here, there is a subtlety: whether the fraction with minimal lexicographic order after taking the reciprocal is also the fraction with minimal lexicographic order before taking the reciprocal. In other words, among fractions \(p/q\) satisfying \(a/b<p/q<c/d\), is the one with minimal \((q,p)\) also the one with minimal \((p,q)\)? Suppose not; let \(p/q\) be minimal in \((q,p)\), but \(r/s\neq p/q\) be minimal in \((r,s)\). Then necessarily \(r<p\) and \(q<s\). But this implies

\[ \dfrac{a}{b} < \dfrac{r}{s} < \dfrac{r}{q} < \dfrac{p}{q} < \dfrac{c}{d}. \]

Therefore, \(r/q\) is strictly smaller than the current solution under any lexicographic order. This contradicts the assumption. Hence, the algorithm is correct.

#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 derivations of the Euclidean-like algorithm discussed in the previous section are often cumbersome, and the sums it can solve are mainly those that can be transformed into (weighted) lattice point counting problems under a line. This section discusses a more general method that further abstracts the above process, enabling it to solve more problems. Therefore, this method is also called the universal Euclidean algorithm. It also utilizes the recursive structure of fractions to solve problems, but the approach to reducing problems differs slightly from the Euclidean-like algorithm.

Consider again the classic summation problem:

\[ f(a,b,c,n)=\sum_{i=1}^n\left\lfloor \frac{ai+b}{c} \right\rfloor, \]

where \(a,b,c,n\) are positive integers.

Problem Transformation

Define a segment with parameters \((a,b,c,n)\) as

\[ y = \frac{ax+b}{c},~0< x\le n. \]

For this segment, define a string \(S\) consisting of \(U\) and \(R\), called the operation sequence, as follows:

• The string has 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 roughly corresponds to starting from the origin, writing an \(R\) each time we cross a vertical grid line to the right, and writing a \(U\) each time we cross a horizontal grid line upward. As shown below:

Of course, this definition needs to handle several special cases:

• When passing through a lattice point (i.e., crossing both upward and right simultaneously), write \(U\) first, then \(R\);
• At the start of the string, besides the number of upward crossings in the interval \((0,1]\), we need to additionally prepend \(\lfloor b/c\rfloor\) \(U\)'s;
• At the end of the string, there should be no extra \(U\)'s.

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

The basic idea of the universal Euclidean algorithm is to treat each \(U\) and \(R\) in the operation sequence as elements of some monoid, treat the entire operation sequence as the product of monoid elements, and relate the final answer to this product.

For example, in this problem, we can define a state vector \(v = (1,y,\sum y)\), representing the current state after starting from the origin and undergoing several upward and right crossings. The first component is constant, the second is the y-coordinate \(y\), and the third is the required sum. Initially, \(v=(1,0,0)\). Each upward crossing increments the y-coordinate by one, equivalent to right-multiplying the state vector by the matrix

\[ U = \begin{pmatrix}1 & 1 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1\end{pmatrix}. \]

Each right crossing adds the current y-coordinate to the sum, equivalent to right-multiplying by

\[ R = \begin{pmatrix}1 & 0 & 0 \\ 0 & 1 & 1 \\ 0 & 0 & 1\end{pmatrix}. \]

Thus, the final state is the product \((1,0,0)S\), where \(S\) is interpreted as the product of these matrices. The answer is the third component of the final state.

Alternatively, instead of defining monoid elements as matrices, we can define them as the contribution of a segment of the operation sequence to the final result, and define the product of operations as merging the contributions of two segments.

In this problem, we can define the contribution of an operation sequence segment as \((x,y,\sum y)\). To explain these notations rigorously, treat each component as a function of the operation sequence. That is, for an operation sequence \(S\), its contribution can be written as \((x(S),y(S),(\sum y)(S))\). Here, \(x(S)\) and \(y(S)\) are the counts of \(R\) and \(U\) in \(S\), i.e., the number of right and upward crossings. The summation symbol in the last term is generally defined as follows: for a function \(f(S)\) on the operation sequence, define \((\sum f)(S)\), or \(\sum_S f\), as

\[ \sum_S f := \sum\{f(S_{[1,r]}):S_r=R\}. \]

Here, \(S_r\) is the \(r\)-th character of \(S\), and \(S_{[1,r]}\) is the prefix of \(S\) consisting of the first \(r\) characters. That is, this summation is over all prefixes of \(S\) ending with \(R\). For example,

\[ \sum_S 1 = x,~ \sum_S x = \dfrac{1}{2}x(x+1). \]

Also, \(\sum y\) is the sum of the number of upward crossings before each right crossing. For the entire sequence, the values of \(y\) at prefixes ending with \(R\) are exactly all \(\lfloor(ai+b)/c\rfloor\) for \(i=1,\cdots,n\). Therefore, \(\sum y\) computed for the entire sequence is the desired quantity.

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

\[ (x_1,y_1,s_1)\cdot (x_2,y_2,s_2) = (x_1+x_2,y_1+y_2,s_1+s_2+x_2y_1). \]

The last term's merging result can be derived as:

\[ \sum_{S_1+S_2}y = \sum_{S_1}y + \sum_{S_2}(y+y_1) = \sum_{S_1}y + \sum_{S_2}y + y_1\sum_{S_2}1 = s_1+s_2+x_2y_1. \]

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

Both methods yield correct results. However, because matrix operations retain more redundant information, their constant factor is larger, so the second method is more practical for actual problems.

Algorithm Process

Unlike the Euclidean-like algorithm which reduces the problem as a whole, the universal Euclidean algorithm reduces the problem by merging operations in batches. Denote the product of operations corresponding to the string as

\[ F(a,b,c,n,U,R). \]

The reduction process is as follows:

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

  \[ \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. \]

  Thus, this is equivalent to changing segment parameters from \((a,b,c,n)\) to \((a,b\bmod c,c,n)\). So for 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\), each \(R\) in the sequence is preceded by at least \(\lfloor a/c\rfloor\) \(U\)'s. We can merge them into the \(R\). That is, replace \(R\) with \(U^{\lfloor a/c\rfloor}R\). After merging, the number of \(U\)'s before the \(i\)-th \(R\) becomes

  \[ \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. \]

  Thus, this is equivalent to changing parameters from \((a,b,c,n)\) to \((a\bmod c,b,c,n)\). So for 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 cases, we need to swap the coordinates, essentially exchanging \(U\) and \(R\), but the parameters of the segment after swapping need careful calculation. From the definition of the operation sequence, we need to determine coefficients \((a',b',c',n')\) such that before swapping, the number of \(R\)'s before the \(j\)-th \(U\) is exactly \(\lfloor(a'j+b')/c'\rfloor\), and there are \(n'\) \(U\)'s total. From the definition,

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

  and 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} \]

  Hence, \(i = \lfloor(cj-b-1)/a\rfloor\). This derivation is similar to that in the Euclidean-like algorithm, again using properties of floor and ceiling functions.

  Two details need handling:

  • The intercept term \(-(b+1)/a\) is negative. Note that shifting the segment left by one unit makes the intercept non-negative, since \((c-b-1)/a\ge 0\). Therefore, we can extract the first segment \(R^{\lfloor(c-b-1)/a\rfloor}U\) before swapping, and only swap the remaining \(U\) and \(R\) in the sequence;
  • After swapping \(U\) and \(R\), there may be extra \(U\)'s at the end. Thus, before swapping, we first extract the last segment of \(R\)'s, and only swap the remaining \(U\) and \(R\). The number of these \(R\)'s is \(n-\lfloor(cm-b-1)/a\rfloor\).

  After removing the head and tail characters, 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 swapping is $m = \

---

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用