The content below was generated entirely by machine translation. Please verify its accuracy. If anything is unclear, consult the Chinese source version.

Euclidean Algorithm

Find the greatest common divisor $(a,b)$ , where $a>b$ .

We have:

$a \div b = q \ldots r\\ a = bq + r \\ r = a - bq$

The Euclidean algorithm states that $\gcd(a,b)=\gcd(b,r)$ .

Let $d$ be any common divisor of $a$ and $b$ , and let:

$m = a \div d,\ n = b \div d$

Then:

$a = dm,\ b = dn\\ \begin{align*} r &= a - bq\\ &= dm - (dn)q\\ &= d(m - nq) \end{align*}$

Since $r=d(m-nq)$ , if $a$ and $b$ have any common divisor $d$ , then $d$ is also a common divisor of $r$ and $b$ (except when $r=0$ ; in that case the greatest common divisor is $a$ ). If $AB$ is the set of common divisors of $a$ and $b$ , and $RB$ is the set of common divisors of $r$ and $b$ , then $AB \in RB$ when $r\ne0$ .

This alone does not prove $\gcd(a,b)=\gcd(b,r)$ , because $RB$ might contain a number larger than $\gcd(a,b)$ . If we prove $RB\in AB$ , then $AB=RB$ , so $RB$ cannot contain a number larger than $\gcd(a,b)$ .

Let $e$ be any number in $RB$ . Then:

$b=me,\qquad r=ne$

Substitute $b$ and $r$ back into $a=bq+r$ :

$\begin{align*} a &= bq + r\\ a &= (me)q + (ne)r\\ a &= (mq + nr)e \end{align*}$

Thus $e\mid a$ , so every $e$ in $RB$ is also in $AB$ ; that is, $RB\in AB$ . Therefore $\gcd(a,b)=\gcd(b,r)$ .

The Euclidean algorithm is very concise in code:

int gcd(int a, int b){
    if(b)
        return gcd(b, a % b);
    else
        return a;
}

References:

Extended Euclidean Algorithm (exgcd)

The extended Euclidean algorithm finds one solution to:

$ax+by=\gcd(a,b).$

For example, the following Euclidean-algorithm calculation computes $\gcd(1180,482)=2$ :

$\begin{align*} \gcd(a,b)&=\gcd(b,r)\\ a&=bq+r\\ 1180&=482(2)+216\\ 482&=216(2)+50\\ 216&=50(4)+16\\ 50&=16(3)+2\\ 16&=2(8)+0 \end{align*}$

We can derive a solution to $2=1180x+482y$ . Start with the penultimate step, $50=16(3)+2$ , and rewrite it as:

$2=50+16(-3).$

Applying $a=bq+r\to r=a+b(-q)$ to the preceding steps gives:

$\begin{align*} 216=50(4)+16&\to16=216+50(-4)\\ 482=216(2)+50&\to50=482+216(-2)\\ 1180=482(2)+216&\to216=1180+482(-2) \end{align*}$

Substitute these expressions into $2=50+16(-3)$ . We can first replace the $16$ in the expression with the sum of $216$ and $50(-4)$ .

The expression now has the form $2=216x+50y$ , where $x=-3$ and $y=13$ .

Next, replace $50$ with the sum of $482$ and $216(-2)$ , turning the expression into $2=482x+216y$ .

Finally, replace $216$ with the sum of $1180$ and $482(-2)$ . The resulting expression is:

$2=1180x+482y$

where $x=-29$ and $y=71$ .

This is exactly the answer we wanted.

It is apparent that exgcd is somewhat like running the Euclidean algorithm in reverse. It uses the calculation process of the Euclidean algorithm to derive one solution to $ax+by=\gcd(a,b)$ .

Now let us generalize the pattern we just observed. What we want to solve first is:

$ax+by=\gcd(a,b)$

Because $\gcd(a,b)=\gcd(b,a\bmod b)$ , $\gcd(b,a\bmod b)$ can also be written in the form $ax+by$ :

$\gcd(b,a\bmod b)=bx_2+(a\bmod b)y_2.$

Although $\gcd(b,a\bmod b)$ is equal to $\gcd(a,b)$ here, giving

$ax+by=bx_2+(a\bmod b)y_2,$

the two sides use the same general form but have different values in the positions of $a$ and $b$ . Their solutions $x,y$ and $x_2,y_2$ are therefore different. Suppose we already know $x_2$ and $y_2$ . If we can determine how to calculate $x$ and $y$ from them, we can solve $x$ and $y$ recursively.

We can simplify $bx_2+(a\bmod b)y_2$ :
Then:

$\begin{align*} ax+by&=bx_2+(a\bmod b)y_2\\ &=bx_2+(a-\lfloor a/b\rfloor b)y_2\\ &=ay_2+b(x_2-\lfloor a/b\rfloor y_2). \end{align*}$

Thus, if we have already found the solution $(x_2,y_2)$ to $bx_2+(a\bmod b)y_2=\gcd(b,a\bmod b)$ , then in the original expression $ax+by=\gcd(a,b)$ we have $x=y_2$ and $y=x_2-\lfloor a/b\rfloor y_2$ . This gives a recursive solution.

The boundary condition is similar to the ordinary Euclidean algorithm: $b=0$ . Then:

$\begin{align*} ax+by &= \gcd(a,b)\\ ax+(0)y &= a\\ x &= 1 \end{align*}$

Although $y$ can take any value in this case, we normally return 0.

The following is the code, using the C++20 standard:

template<typename T>
concept Integral = std::is_integral<T>::value;
// gcd, x, y
template<Integral T>
tuple<T, T, T> ex_gcd(T a, T b){
    if (b == 0) {
        return {a, 1, 0};
    }
    auto [gcd, x2, y2] = ex_gcd(b, a % b);
    // Derive x and y from x2 and y2.
    T x = y2;
    T y = x2 - (a / b) * y2;
    return {gcd, x, y};
}

References:

Modular Multiplicative Inverse

The modular multiplicative inverse of $a\bmod p$ is the solution $x$ of $ax\equiv1\pmod b$ .

A modular inverse is somewhat like an additive inverse under a modulus.

exgcd

When $a$ and $b$ are coprime, exgcd can solve this problem. Since $\gcd(a,b)=1$ , extended Euclid solves $ax+by=1$ .

Rewrite $ax\equiv1\pmod b$ as:

$ax\equiv1\pmod b\\ ax\bmod b=1\\ ax-bk=1.$

Setting $y=-k$ gives $ax+by=1$ .

One of $x$ and $y$ in $ax+by=1$ may be negative. A negative $y$ causes no problem, but if $x$ is negative, the answer we obtain is not the smallest positive integer among all feasible values of $x$ .

Looking at $ax+by=1$ , we can add a multiple of $b$ to $x$ , transforming the expression into $a(x+bn)+b(y+an)=1$ (note that $b$ is negative, so the $abn$ terms cancel). This lets us make $x$ positive without changing the equality $ax+by=1$ .

We can therefore write:

1	x = (x % b + b) % b;

Assume first that $x$ is negative.

The first x % b adds some multiples of $b$ to $x$ , turning it into the largest negative number that still satisfies the condition. For example, suppose $b$ is 13 and $x$ is -25. After x = x % b, $x$ becomes -12, which is equivalent to adding 13 to $x$ .

The following +b turns this largest valid negative number into the smallest valid positive number. For example, $x+b=-12+13=1$ . What, then, is the purpose of the final % b?

It handles the case where $x$ was positive to begin with. By subtracting suitable multiples of $b$ from $x$ , it reduces $x$ to the smallest positive number satisfying the condition.

For the modular-inverse template problem:

int n, p;
template<typename T>
concept Integral = std::is_integral<T>::value;
template<Integral T>
tuple<T, T, T> ex_gcd(T a, T b){
    if (b == 0) return {a, 1, 0};
    auto[gcd, x2, y2] = ex_gcd(b, a % b);
    T x = y2;
    T y = x2 - (a / b) * y2;
    return {gcd, x, y};
}
int main() { 
    cin>>n>>p;
    for(int i = 1; i <= n; i++){
        auto[gcd, x, y] = ex_gcd(i, p);
        x = (x % p + p) % p;
        cout<<x<<endl;
    }
}

Because the data size is $3e6$ and the time limit is 500 ms, an $n\log p$ algorithm is too slow; use the linear algorithm below.

References:

Linear Recurrence

The linear recurrence computes the modular inverses of all integers from $1$ to $n$ modulo a prime $p$ in $O(n)$ time. $p$ must be prime to ensure all values in the range are coprime to it.

Because this is a recurrence algorithm, it needs an initial condition. It is easy to see that the inverse of 1 modulo any integer is 1 itself, because $1\times1=1$ . This gives us the initial condition.

Suppose the recurrence has now reached the number $i$ . Write $p\div i=k\ldots r$ , so $p=ki+r$ . Converting this to a congruence gives:

$ki+r\equiv0\pmod p.$

Let $i^{-1}$ and $r^{-1}$ denote the modular inverses of $i$ and $r$ modulo $p$ , respectively. Multiplying both sides of the congruence by $i^{-1}r^{-1}$ and expanding gives:

$kr^{-1}+i^{-1}\equiv0\pmod p,$

so:

$i^{-1}\equiv-k r^{-1}\pmod p.$

The simplification uses $i^{-1}i\equiv1\pmod p$ , with the same relationship holding for $r$ and $r^{-1}$ . Since $k=\lfloor p/i\rfloor$ and $r=p\bmod i$ :

$i^{-1}\equiv-\lfloor p/i\rfloor\times(p\bmod i)^{-1}\pmod p.$

Normalize the possibly negative value in the same way:

1	x = (x % b + b) % b;

The template code is:

ll inv[MAXN];
template <typename T>
concept Int_t = is_integral<T>::value;
template <Int_t T>
inline T mod_norm(T val, T m) {
    return (val % m + m) % m;
}
int main() {
    ios::sync_with_stdio(0);
    cin.tie(0);
    ll n, p;
    cin >> n >> p;
    inv[1] = 1;
    cout << inv[1] <<'\n';
    for (int i = 2; i <= n; i++) {
        inv[i] = mod_norm(-p / i * inv[p % i] % p, p);
        cout << inv[i] <<'\n';
    }
}