Skip to content

Higher Residues & Roots of Unity

This article discusses higher-degree residues and roots of unity modulo a modulus, and introduces algorithms for extracting roots modulo a modulus.

Higher-Degree Residues

Higher-degree residues modulo a modulus essentially discuss the feasibility of extracting higher-degree roots modulo a modulus. It is a generalization of quadratic residues.

\(k\)-th Residue

Let integer \(k\geq 2\), integer \(a\) and positive integer \(m\) be coprime. If there exists an integer \(x\) such that

\[ x^k\equiv a\pmod m, \]

then \(a\) is called a \(k\)-th residue modulo \(m\), and \(x\) is called a \(k\)-th root of \(a\) modulo \(m\); otherwise, \(a\) is called a \(k\)-th nonresidue modulo \(m\).

That is, a \(k\)-th root of \(a\) modulo \(m\) exists if and only if \(a\) is a \(k\)-th residue modulo \(m\).

Properties

Similar to quadratic residues, one can discuss the criteria, count, and number of residue classes for \(k\)-th residues. Like other congruence equation problems, they can be reduced to prime power moduli via the Chinese Remainder Theorem. Depending on the existence of primitive roots, this further divides into odd prime power moduli and powers of \(2\) as moduli.

The case of odd prime powers is relatively simple. In fact, for all cases where primitive roots exist, we have the following result:

Theorem

Let integer \(k\geq 2\), integer \(a\) and positive integer \(m\) be coprime. Assume a primitive root modulo \(m\) exists, and let \(g\) be a primitive root modulo \(m\). Denote \(d=\gcd(k,\varphi(m))\) and \(d'=\dfrac{\varphi(m)}{d}\), where \(\varphi(m)\) is the Euler totient function. Then:

  1. \(a\) is a \(k\)-th residue modulo \(m\) if and only if

    \[ a^{d'} \equiv 1 \pmod m. \]
  2. When \(a\) is a \(k\)-th residue modulo \(m\), modulo equivalence, \(a\) has exactly \(d\) distinct \(k\)-th roots modulo \(m\), and they are of the form

    \[ x \equiv g^{y_0+id'}\pmod{\varphi(m)},~0\le y_0 < d',~i=0,1,\cdots,d-1. \]
  3. The number of \(k\)-th residue classes modulo \(m\) is \(d'\), and their complete set is

    \[ \{g^{di}\bmod m : 0 \le i < d'\}. \]
Proof

Since \(a\perp m\), we have \(x\perp m\). Since \(g\) is a primitive root modulo \(m\), both \(x\) and \(a\) are congruent to some power of \(g\). Let \(x\equiv g^y\pmod m\), then the equation \(x^k\equiv a\pmod m\) is equivalent to

\[ g^{ky} \equiv g^{\operatorname{ind}_g a}\pmod m. \]

Here, \(\operatorname{ind}_g a\) is the discrete logarithm. By the properties of order and \(\delta_m(g)=\varphi(m)\), this is equivalent to the congruence equation

\[ ky \equiv \operatorname{ind}_g a \pmod{\varphi(m)}. \]

This is a linear congruence equation in \(y\). Applying the analysis of its solution structure from that page, we know the equation has a solution if and only if \(d\mid\operatorname{ind}_g a\), and the general solution is of the form

\[ y = y_0 + id' \pmod{\varphi(m)},~0\le y_0 < d',~i=0,1,\cdots,d-1. \]

From this, almost all content of the theorem follows; the only extra part to justify is the criterion \(a^{d'} \equiv 1 \pmod m\). By property 3 of order,

\[ \delta_m(a) = \delta_m(g^{\operatorname{ind}_g a}) = \dfrac{\varphi(m)}{\gcd(\varphi(m),\operatorname{ind}_g a)} = \dfrac{\varphi(m)}{\operatorname{ind}_g a}. \]

And we know the equation has a solution if and only if \(d\mid \operatorname{ind}_g a\), i.e., \(\delta_m(a)\mid d'\). By property 2 of order, this is equivalent to the criterion.

The case where the modulus is a power of \(2\) is special. To handle this case, we need a result about the structure of the reduced residue system modulo \(2^e\): every odd integer \(a\) is uniquely congruent to an integer of the form \((-1)^s5^r\bmod 2^e\), where \(s\in\{0,1\}\) and \(0\le r < 2^{e-2}\). Using this result, we have:

Theorem

Let integer \(k\ge 2\), odd integer \(a\), and positive integer \(m=2^e\) with \(e \ge 2\). When \(k\) is odd:

  1. \(a\) is always a \(k\)-th residue modulo \(m\).
  2. \(a\) has exactly one \(k\)-th root modulo \(m\).
  3. The number of \(k\)-th residue classes modulo \(m\) is \(2^{e-1}\), and they are exactly all reduced residue classes.

When \(k\) is even, denote \(d=\gcd(k,2^{e-2})\) and \(d'=\dfrac{2^{e-2}}{d}\):

  1. \(a\) is a \(k\)-th residue modulo \(m\) if and only if \(a\equiv 1\pmod 4\) and \(a^{d'}\equiv 1\pmod m\).
  2. When \(a\) is a \(k\)-th residue modulo \(m\), modulo equivalence, \(a\) has exactly \(2d\) distinct \(k\)-th roots modulo \(m\), and they are of the form

    \[ x \equiv \pm 5^{y_0 + id'} \pmod{2^{e-1}},~ 0 \le y_0 < d',~i = 0, 1,\cdots,d-1. \]
  3. The number of \(k\)-th residue classes modulo \(m\) is \(d'\), and their complete set is

    \[ \{5^{di}\bmod m : 0 \le i < d'\}. \]
Proof

Since \(a\perp m\), we have \(x\perp m\). Since \(x\) and \(a\) are odd, by the aforementioned result, we can set \(a\equiv (-1)^s5^r\pmod{2^e}\) and \(x=(-1)^z5^{y}\pmod{2^e}\). Because the representation is unique, the congruence equation \(x^k\equiv a\pmod{2^e}\) is equivalent to the system of linear congruence equations:

\[ \begin{aligned} kz &\equiv s \pmod{2},\\ ky &\equiv r \pmod{2^{e-2}}. \end{aligned} \]

Combining with the analysis of solutions to linear congruence equations from that page, we obtain the solution structure of \(x^k\equiv a\pmod{2^e}\). Depending on the parity of \(k\), there are two cases:

  • When \(k\) is odd, since \(\gcd(k,2)=\gcd(k,2^{e-2})=1\), both linear congruence equations have solutions for all \(s,r\), hence the original congruence equation always has a solution for all odd \(a\).
  • When \(k\) is even, the first equation has a solution if and only if \(2\mid s\), and the second equation has a solution if and only if \(d=\gcd(k,2^{e-2})\mid r\). Combining these gives the complete form of \(k\)-th residue classes. Direct calculation shows the first condition is equivalent to \(a\equiv 1\pmod 4\); repeating the analysis for odd prime power cases shows the second condition is equivalent to \(a^{d'}=1\). Combining these yields the criterion in the theorem. The general solutions of the two linear congruence equations are known:

    \[ \begin{aligned} z &\equiv0,1\pmod 2, \\ y &\equiv y_0 + id' \pmod{2^{e-2}},~ 0\le y_0 < 2^{e-2}. \end{aligned} \]

    Combining them gives the general solution of the original equation.

This completely solves the problem of determining \(k\)-th residues for different moduli. Concepts like the Legendre symbol and quadratic reciprocity in quadratic residues can also be generalized to higher-degree residues, but this is non-trivial and requires concepts like cyclotomic fields. In algebraic number theory, quadratic reciprocity ultimately generalizes to Artin reciprocity.

Roots of Unity

As a special case of \(k\)-th roots, this section discusses the concept of \(k\)-th (primitive) roots of unity. It can be seen as the counterpart in the reduced residue system \(\mathbf Z_m^*\) modulo \(m\) of the concept of \(k\)-th roots of unity in the complex field \(\mathbf C\). When the modulus \(m\) is suitable, using a \(k\)-th primitive root of unity modulo \(m\) in place of the complex root \(\omega_k\) can accelerate computations.

Similar to the complex case, we have the following definition:

\(k\)-th Root of Unity Modulo \(m\)

For modulus \(m\), a \(k\)-th root of the element \(1\) is called a \(k\)-th root of unity modulo \(m\). In particular, if \(x\) is a \(k\)-th root of unity modulo \(m\) and it is not any \(k' < k\) root of unity modulo \(m\), then \(x\) is also called a \(k\)-th primitive root of unity modulo \(m\).

Comparing with the definition of primitive root, a primitive root \(g\) is exactly a \(\varphi(m)\)-th primitive root of unity modulo \(m\), where \(\varphi(m)\) is the Euler totient function.

When a \(k\)-th primitive root of unity modulo \(m\) exists, its algebraic properties are consistent with the \(k\)-th primitive complex root of unity \(\omega_k\), and it can replace \(\omega_k\) in various computations. For example, applying it to the Fast Fourier Transform yields the Number Theoretic Transform over finite fields[^fnnt].

Properties

In the complex field, roots of unity (primitive) of any order exist. However, this is not the case in number theory.

Properties

For modulus \(m\), let \(\lambda(m)\) be its Carmichael function. Then:

  1. All integers \(a\) coprime to \(m\) are \(\delta_m(a)\)-th primitive roots of unity modulo \(m\), where \(\delta_m(a)\) is the order of \(a\) modulo \(m\).
  2. If element \(a\) is a \(k\)-th root of unity modulo \(m\), and \(k'\) is any multiple of \(k\), then \(a\) is also a \(k'\)-th root of unity modulo \(m\).
  3. If element \(a\) is a \(k\)-th (primitive) root of unity modulo \(m\), then element \(a^{\ell}\) is a \(\dfrac{k}{\gcd(k,\ell)}\)-th (primitive, respectively) root of unity modulo \(m\).
  4. As \(k'\) runs through the divisors of \(k\), all \(k'\)-th primitive roots of unity modulo \(m\) exactly constitute a partition of the \(k\)-th roots of unity modulo \(m\). Moreover, for \(\ell\perp k\), the mapping \(x\mapsto x^\ell\) gives a bijection among \(k\)-th roots of unity, preserving the partition: it maps \(k'\)-th primitive roots of unity (for \(k'\mid k\)) to \(k'\)-th primitive roots of unity.
  5. A \(k\)-th primitive root of unity modulo \(m\) exists if and only if \(k\mid\lambda(m)\). In particular, a \(\lambda(m)\)-th primitive root of unity modulo \(m\) exists, called a \(\lambda\)-primitive root modulo \(m\).
  6. Element \(a\) is a \(k\)-th root of unity modulo \(m\) if and only if \(a^k\equiv 1\pmod{m}\) and for every prime factor \(p\mid k\), \(a^{k/p}\not\equiv 1\pmod{m}\).
Proof

By the definition of order, all integers \(a\) coprime to \(m\) are \(\delta_m(a)\)-th primitive roots of unity modulo \(m\), where \(\delta_m(a)\) is the order of \(a\) modulo \(m\). Conversely, if \(a\) is a \(k\)-th root of unity modulo \(m\), then \(\gcd(a^k,m)=1\), so \(\gcd(a,m)=1\). Therefore, \(a\) is a (primitive) root of unity modulo \(m\) if and only if \(a\) is coprime to \(m\). This is property 1.

Direct verification of the definition shows that if \(k\mid k'\), then \(a^k\equiv 1\pmod m\) implies \(a^{k'}\equiv 1\pmod m\). This is property 2. By properties of order,

\[ \delta(a^\ell) = \dfrac{\delta_m(a)}{\gcd(\delta_m(a),\ell)}. \]

If \(a\) is a \(k\)-th primitive root of unity modulo \(m\), then \(\delta_m(a)=k\). Substituting into the above gives that \(a^\ell\) is a \(\dfrac{k}{\gcd(k,\ell)}\)-th primitive root of unity modulo \(m\). If \(a\) is only a \(k\)-th root of unity modulo \(m\), suppose it is a \(k'\mid k\)-th primitive root of unity, then \(a^\ell\) is a \(\dfrac{k'}{\gcd(k',\ell)}\)-th primitive root of unity modulo \(m\). Since \(k'\mid k\), we have

\[ \dfrac{k'}{\gcd(k',\ell)} \mid \dfrac{k}{\gcd(k,\ell)}, \]

and by property 2, \(a^\ell\) is a \(\dfrac{k}{\gcd(k,\ell)}\)-th root of unity modulo \(m\). This is property 3.

For \(k'\mid k\), by property 2, \(k'\)-th primitive roots of unity modulo \(m\) are necessarily \(k\)-th roots of unity modulo \(m\). They are pairwise disjoint, thus forming a partition. For \(\ell\perp k\), we always have \(\ell\perp k'\), so for a \(k'\)-th primitive root of unity \(a\) modulo \(m\), \(a^\ell\) is always a \(k'\)-th primitive root of unity modulo \(m\). Taking \(\ell'=\ell^{-1}\bmod k\), we can verify that \(x\mapsto x^\ell\) and \(x\mapsto x^{\ell'}\) are inverse mappings, hence \(x\mapsto x^\ell\) is a bijection. This is property 4.

By properties of the Carmichael function, a \(\lambda(m)\)-th primitive root of unity modulo \(m\) always exists; let it be \(a\), with \(\delta_m(a)=\lambda(m)\). For \(k\mid\lambda(m)\), let \(k'=\dfrac{\lambda(m)}{k}\). Then

\[ \delta_m(a^{k'}) = \dfrac{\lambda(m)}{(\lambda(m),k')} = \dfrac{\lambda(m)}{k'} = k. \]

Therefore, \(a^{k'}\) is a \(k\)-th primitive root of unity. And by the definition of the Carmichael function, the order of any \(x\perp m\) is a divisor of \(\lambda(m)\). This yields property 5.

Almost repeating the proof of the primitive root criterion gives property 6. This criterion essentially verifies \(\delta_m(a)=k\).

From these properties, we see that relative to cases where primitive roots exist, the \(\lambda\)-primitive root modulo \(m\) plays a similar foundational role. Unlike primitive roots, powers of a \(\lambda\)-primitive root cannot generate all roots of unity modulo \(m\). Nevertheless, since the density of \(\lambda\)-primitive roots is not low[^lambda-density], if one indeed needs to find a \(k\)-th primitive root of unity, one can first find a \(\lambda\)-primitive root via a randomized method, then obtain a \(k\)-th primitive root of unity by exponentiation.

If a \(k\)-th root of \(a\) modulo \(m\) is known, all \(k\)-th roots of \(a\) modulo \(m\) can be generated using all \(k\)-th roots of unity modulo \(m\).

Theorem

Let \(x\) be a \(k\)-th root of \(a\) modulo \(m\). As \(r\) runs through all \(k\)-th roots of unity modulo \(m\), \(xr\) runs through all \(k\)-th roots of \(a\) modulo \(m\).

Proof

For two \(k\)-th roots \(x,y\) of \(a\) modulo \(m\), let \(r=x^{-1}y\bmod m\). Then \(r\) satisfies \(r^k\equiv 1\pmod m\), so it is a \(k\)-th root of unity modulo \(m\). Conversely, if \(r\) is a \(k\)-th root of unity modulo \(m\), then \((xr)^{k}= x^kr^k\equiv a\pmod m\), meaning \(xr\) is a \(k\)-th root of \(a\) modulo \(m\).

Using \(k\)-th roots of unity to generate all \(k\)-th roots is analogous to using solutions of a homogeneous linear system to generate the general solution of a non-homogeneous linear system.

The above discusses the general case. For the case where primitive roots exist, the structure of roots of unity is simpler:

Theorem

For modulus \(m\), assume a primitive root modulo \(m\) exists, and let \(a\) be a \(k\)-th primitive root of unity modulo \(m\). Then \(b\) is a \(k\)-th root of unity modulo \(m\) if and only if it can be expressed as a power of \(a\).

Proof

Let \(g\) be a primitive root modulo \(m\). Then all elements coprime to \(m\) can be expressed as powers of \(g\). So \(a\) is a \(k\)-th primitive root of unity modulo \(m\) if and only if

\[ \delta_m(a) = \delta_m(g^{\operatorname{ind}_ga}) = \dfrac{\varphi(m)}{\gcd(\varphi(m),\operatorname{ind}_ga)} = k. \]

Similarly, \(b\) is a \(k\)-th root of unity modulo \(m\) if and only if

\[ \delta_m(b) = \delta_m(g^{\operatorname{ind}_gb}) = \dfrac{\varphi(m)}{\gcd(\varphi(m),\operatorname{ind}_gb)} = k' \mid k. \]

Thus,

\[ \gcd(\varphi(m),\operatorname{ind}_ga) \mid \gcd(\varphi(m),\operatorname{ind}_gb)\mid \operatorname{ind}_gb. \]

By the analysis of linear congruence equations, this condition is equivalent to the equation

\[ (\operatorname{ind}_ga) x \equiv \operatorname{ind}_gb \pmod{\varphi(m)} \]

having a solution. Taking powers with base \(g\) yields \(a^x\equiv b\pmod{m}\), i.e., \(b\) can be expressed as a power of \(a\).

This theorem shows that when primitive roots exist, all \(k\)-th roots of unity form a cyclic group structure, and the \(k\)-th primitive roots of unity are generators of this cyclic group. Later we will see that the Tonelli–Shanks algorithm exploits this to accelerate the discrete logarithm part in root extraction.

Root Extraction Modulo a Modulus

Finally, this article discusses methods for finding \(k\)-th roots. For \(k=2\), there are many efficient algorithms for extracting square roots modulo a modulus. However, for general \(k\), no known polynomial-time algorithm exists. This section introduces two common algorithms that can find a \(k\)-th root in \(O(m^{1/2})\) and \(O(m^{1/4+\varepsilon})\) time, respectively. Using the Chinese Remainder Theorem, the problem can always be reduced to prime power moduli, so this section focuses on solutions for prime power moduli.

Naive Algorithm

The analysis of properties of \(k\)-th residues earlier essentially points out a method for solving \(k\)-th roots modulo prime powers. Strictly speaking, the case solved earlier is when the radicand \(a\) is coprime to the modulus \(m\). The algorithm process is summarized as follows:

  • When \(m=p^e\) is an odd prime power, let \(g\) be a primitive root modulo \(m\). Then the equation \(x^k\equiv a\pmod m\) can be transformed into the linear congruence equation

    \[ ky \equiv \operatorname{ind}_g a \pmod{\varphi(m)}. \]

    Here, \(\operatorname{ind}_g a\) can be found via the BSGS algorithm, and the complete solution set of the linear congruence equation is easy to find. Thus, all \(k\)-th roots \(x\equiv g^y\pmod m\) of \(a\) can be obtained.

    Alternatively, there is a similar approach. Still set \(x\equiv g^y\pmod m\), but transform via

    \[ x^k \equiv (g^k)^y \equiv a \pmod m \]

    into finding the discrete logarithm of \(a\) with base \(g^k\). This can also be solved by the BSGS algorithm to find a particular solution. Its general solution can be obtained from the solution expression mentioned earlier, i.e., by multiplying the particular solution with all \(k\)-th roots of unity.

    Regardless of the approach, when a primitive root is known, the complexity of finding a single solution is \(O(m^{1/2})\). Since a primitive root can be found in \(o(m^{1/2})\) time, the overall time complexity remains \(O(m^{1/2})\).

  • When \(m=2^e\) and \(e\in\mathbf N_+\), first find \(s,r\) in \(a\equiv (-1)^s5^r\pmod m\). Among these exponents, \(s\) can be determined in \(O(1)\) time:

    \[ s = \begin{cases}0, & a\equiv 1\pmod 4, \\ 1, & a\equiv 3\pmod 4.\end{cases} \]

    And \(r=\operatorname{ind}_5((-1)^sa)\) can be found via the BSGS algorithm in \(O(m^{1/2})\) time. Next, solve the system of linear congruence equations:

    \[ \begin{aligned} kz &\equiv s \pmod{2},\\ ky &\equiv r \pmod{2^{e-2}}. \end{aligned} \]

    The general solution \((z,y)\) of this system is easy to find, and \(x=(-1)^z5^y\) is the desired root. The complexity of finding a single solution remains \(O(m^{1/2})\).

Of course, for cases with no solution, it can be quickly determined in \(O(\log m)\) time using the criteria discussed earlier, without needing to check during the solving process.

A reference implementation for finding \(k\)-th roots modulo a prime is as follows: (The code is for illustration only; due to high time complexity, it cannot pass the problem)

Reference Implementation for Problem Library Checker - Kth Root (Mod)
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
// Submission (TLE): https://judge.yosupo.jp/submission/320582
#include <algorithm>
#include <cmath>
#include <iostream>
#include <tuple>
#include <unordered_map>
#include <vector>

struct PrimePower {
  int p, e, pe;

  PrimePower(int p, int e, int pe) : p(p), e(e), pe(pe) {}
};

// Factorization.
auto factorize(int n) {
  std::vector<PrimePower> ans;
  for (int x = 2; x * x <= n; ++x) {
    int e = 0, pe = 1;
    for (; n % x == 0; n /= x, ++e, pe *= x);
    if (e) ans.emplace_back(x, e, pe);
  }
  if (n > 1) ans.emplace_back(n, 1, n);
  return ans;
}

// Binary exponentiation.
int pow(int a, int b, int m = 0) {
  int res = 1;
  for (; b; b >>= 1) {
    if (b & 1) res = m ? (long long)res * a % m : res * a;
    a = m ? (long long)a * a % m : a * a;
  }
  return res;
}

// Find a primitive root modulo prime.
int primitive_root(int p) {
  std::vector<int> exp;
  for (auto factor : factorize(p - 1)) {
    exp.push_back((p - 1) / factor.p);
  }
  int ans = 0;
  bool succ = false;
  while (!succ) {
    ++ans;
    succ = true;
    for (int b : exp) {
      if (pow(ans, b, p) == 1) {
        succ = false;
        break;
      }
    }
  }
  return ans;
}

// Discrete logarithm. (BSGS Algorithm)
int log(int g, int a, int m) {
  int b = std::sqrt(m + 0.25l) + 1;
  std::unordered_map<int, int> mp;
  int po0 = a % m, po1 = 1;
  for (int i = 0; i < b; ++i) {
    mp[po0] = i;
    po0 = (long long)po0 * g % m;
  }
  po0 = pow(g, b, m);
  for (int j = 1; j <= b; ++j) {
    po1 = (long long)po1 * po0 % m;
    if (mp.count(po1)) return j * b - mp[po1];
  }
  return -1;
}

// Extended Euclidean Algorithm.
int ex_gcd(int a, int b, int& x, int& y) {
  if (!b) {
    x = 1;
    y = 0;
    return a;
  } else {
    int d = ex_gcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
  }
}

// Solves the linear congruence equation: Ax = B mod N.
// Return the least nonnegative solution and the common difference.
std::pair<int, int> solve_linear(int a, int b, int n) {
  int x, y;
  int d = ex_gcd(a, n, x, y);
  if (b % d) return {-1, -1};
  n /= d;
  x = ((long long)x * (b / d) % n + n) % n;
  return {x, n};
}

// Subroutine: Find a K-th root with a primitive root G known.
int calc(int g, int k, int a, int p) {
  int ind = log(g, a, p);
  if (ind == -1) return -1;
  int y0, d;
  std::tie(y0, d) = solve_linear(k, ind, p - 1);
  if (y0 == -1) return -1;
  return pow(g, y0, p);
}

// Find a K-th root of A modulo prime P.
int kth_roots_mod_p(int k, int a, int p) {
  a %= p;
  if (k == 0) return a == 1 ? 0 : -1;
  if (a == 0) return 0;
  return calc(primitive_root(p), k, a, p);
}

void solve() {
  int k, y, p;
  std::cin >> k >> y >> p;
  std::cout << kth_roots_mod_p(k, y, p) << '\n';
}

int main() {
  int t;
  std::cin >> t;
  for (; t; --t) {
    solve();
  }
}

Improved Tonelli–Shanks Algorithm

By appropriately generalizing the Tonelli–Shanks algorithm for extracting square roots modulo a modulus, one can solve root extraction modulo prime powers. A straightforward generalization is the Adleman–Manders–Miller algorithm[^amm], but its complexity is still not optimal[^amm-comp]. This section introduces the improved Tonelli–Shanks algorithm proposed by sugarknri, Min_25, 37zigen, et al. It can find a \(k\)-th root in \(O(m^{1/4+\varepsilon})\) time.

The core idea of the Tonelli–Shanks algorithm is to solve the discrete logarithm in a group of order \(2^e\), thereby reducing time complexity. Similarly, for discrete logarithms in groups of prime power order \(p^e\), one can also solve them relatively efficiently, but the algorithm complexity is \(\Omega\left(\sqrt{p}\right)\). The Adleman–Manders–Miller algorithm splits the \(k\)-th root extraction into computing discrete logarithms in multiple prime power order groups, but limited by the size of the largest prime factor \(p_\text{max}(k)\) of \(k\), the algorithm complexity remains \(\Omega\left(\sqrt{p_\text{max}(k)}\right)\). This section's algorithm further improves this process, avoiding discrete logarithm computations for large prime factors, thereby controlling the overall complexity to \(O(m^{1/4+\varepsilon})\).

Process

Consider computing a \(k\)-th root of \(a\) modulo a prime power \(m\), i.e., solving the congruence equation:

\[ x^k \equiv a \pmod m. \]

Specifically, for the case \(m=2^e\), we also require \(a\equiv 1\pmod{4}\), so \(a\) can be written as a power of \(g=5\). As discussed earlier, computing \(k\)-th roots modulo \(2^e\) can always be reduced to this case. When handling modulo \(2^e\), \(\varphi(m)\) mentioned in this section should be replaced by \(\delta_m(5)=2^{e-2}\).

First, the problem can be reduced to the case where the root exponent divides \(\varphi(m)\). Let \(d=\gcd(k,\varphi(m))\). Then, by properties of \(k\)-th residues, when \(a\) is a \(k\)-th residue modulo \(m\), \(a\) is always a \(\dfrac{\varphi(m)}{d}\)-th root of unity modulo \(m\). By properties of roots of unity, for any \(\ell\perp\dfrac{\varphi(m)}{d}\), the mapping \(x\mapsto x^{\ell}\) is a bijection among \(\dfrac{\varphi(m)}{d}\)-th roots of unity. Therefore, we can take

\[ \ell = \left(\dfrac{k}{d}\right)^{-1}\bmod\dfrac{\varphi(m)}{d}. \]

Raising both sides of the original congruence equation to the power \(\ell\) yields

\[ x^d\equiv x^{k\ell} \equiv a^{\ell} =: b \pmod{m}. \]

The leftmost congruence uses Euler's theorem and the congruence relation: (for some \(c\in\mathbf Z\))

\[ k\ell = d\left(\frac{k}{d}\ell\right) = d\left(c\dfrac{\varphi(m)}{d}+1\right) \equiv d \pmod{\varphi(m)}. \]

For the transformed problem, consider the prime factorization of \(d\):

\[ d = \prod_{p\in\mathbf P}p^e. \]

Starting from \(b=a^\ell\), for each \(p^e\neq 1\), extract the \(p^e\)-th root successively; finally, we obtain a \(d\)-th root of \(b\), which is a \(k\)-th root of \(a\).

Finally, the problem reduces to solving the equation:

\[ x^{p^e} \equiv b \pmod m. \]

Assume \(\varphi(m)=p^sr\) with \(p\perp r\). Let \(q\in\mathbf N_+\) be a solution to \(qr\equiv -1\pmod{p^e}\). Then, since \(b\) is an \(rp^{s-e}\)-th root of unity, \(b^{qr}\) must be a \(p^{s-e}\)-th root of unity. Also, let \(\zeta\) be a \(p^s\)-th primitive root of unity modulo \(m\). Then \(\zeta^{p^e}\) is a \(p^{s-e}\)-th primitive root of unity, so there exists \(h\in\mathbf N\) such that \(b^{qr}\equiv \zeta^{hp^{e}}\pmod{m}\). Direct verification shows that

\[ x\equiv b^{(qr+1)/p^e}\zeta^{-h} \pmod{m} \]

is a \(p^e\)-th root of \(b\) modulo \(m\).

To compute \(x\), we need to find a \(p\)-th nonresidue \(\eta\) modulo \(m\). For this, by earlier properties, it suffices to randomly pick \(\eta\perp m\) and verify \(\eta^{\varphi(m)/p}\bmod{m}\neq 1\). The density of such numbers is

\[ \dfrac{\varphi(m)}{m}\left(1-\dfrac{1}{p}\right) \ge \dfrac{1}{4}. \]

Thus, we expect to find one within \(4\) random trials on average. Note that \(\eta^{rp^{s-1}}\not\equiv 1\pmod m\) and \(\eta^{rp^s}\equiv 1\pmod m\), so if we set \(\zeta=\eta^r\bmod m\) and \(\xi=\eta^{rp^{s-1}}\bmod m\), then they are \(p^s\)-th and \(p\)-th primitive roots of unity, respectively.

Finally, we need to compute \(h\in\mathbf N\). Clearly, we can take \(h < p^{s-e}\). Consider the base-\(p\) expansion of \(h\):

\[ h = \sum_{j=0}^{s-e-1}h_jp^j = h_0 + h_1p + h_2p^2 +\cdots. \]

Compute these digits sequentially. When the first \(j\) digits are computed, we must have

\[ \left(b^{qr}\zeta^{-p^e(h_0+h_1p+\cdots + h_{j-1}p^{j-1})}\right)^{p^{s-e-j-1}} \equiv \zeta^{h_jp^{s-1}} \equiv \xi^{h_j} \pmod{m}. \]

Thus, \(h_j\) can be found by computing the discrete logarithm with base \(\xi\). To achieve better time complexity, the BSGS algorithm is used. In total, \((s-e)\) discrete logarithms need to be computed. Suppose we precompute \(B\) powers of \(\xi\); then the time complexity for a single discrete logarithm query is \(O(p/B)\), and the total time complexity is

\[ O\left(B+(s-e)\dfrac{p}{B}\right). \]

When \(B=\sqrt{(s-e)p}\), the total time complexity is minimized to \(O\left(\sqrt{(s-e)p}\right)\). After obtaining \(h\), substitute into the expression for \(x\) to get a particular solution.

Time Complexity

The time complexity of this algorithm is \(O(m^{1/4+\varepsilon})\). In this complexity discussion, we assume a single multiplication takes \(O(1)\) time, and when computing powers, we always reduce exponents using Euler's theorem, so computing a single power involved always takes \(O(\log m)\) time.

First consider computing a single \(p^e\)-th root. Finding a \(p\)-th nonresidue requires verifying \(O(1)\) numbers on average, with total time \(O(\log m)\). Computing \(s,r,\zeta,\eta,b^{qr}\) each takes \(O(\log m)\) time. When computing \(h\), each digit requires \(O(\log m)\) time to compute a power, and there are \((s-e)\) digits, so total time is \(O((s-e)\log m)\). As mentioned earlier, the total time for preprocessing and \((s-e)\) discrete logarithm queries is \(O\left(\sqrt{(s-e)p}\right)\). Since \(s-e\in O(\log m)\), the time complexity for computing a single \(p^e\)-th root is \(O(p^{1/2+\varepsilon})\). In particular, when \(s=e\), the time complexity can be further reduced to \(O(\log m)\).

Now consider the overall algorithm time complexity. Computing \(\varphi(m),d,\ell\) each takes \(O(\log m)\) time. Next, we need to factorize \(d=\prod_p{p^e}\); using the Pollard Rho algorithm, this can be done in \(O(m^{1/4})\) time. Finally, sequentially computing \(p^e\)-th roots has total time complexity

\[ O\left(\sum_{e < s}p^{1/2+\varepsilon}\right). \]

Since prime factors \(p\) with \(e < s\) appear at least twice in \(\varphi(m)\), we must have \(p < m^{1/2}\). Hence, the total time complexity is \(O(m^{1/4+\varepsilon})\).

In fact, in this scenario, we can achieve \(O(m^{1/4+\varepsilon})\) time complexity without using Pollard Rho for factorization. Indeed, it suffices to factor \(d\) by trial division up to primes not exceeding \(m^{1/4}\). Let \(z\) be the integer obtained after removing these small prime factors. Then for prime factors \(p > m^{1/4}\) of \(z\), we must have \(\nu_p(\varphi(m)) < 4\), where \(\nu_p(n)\) denotes the exponent of \(p\) in the prime factorization of \(n\). Since we only need to consider cases where

\[ 1 \le e = \nu_p(d) < s = \nu_p(\varphi(m)) < 4, \]

there can be at most one such prime factor \(p\); otherwise, if there were more, each would appear at least twice in \(\varphi(m)\), and their product would exceed \(m\). To isolate this (possibly existing) unique large prime factor, compute

\[ p^\star=\gcd\left(z,\dfrac{\varphi(m)}{z}\right) = \prod_{p : \nu_p(d) < \nu_p(\varphi(m))}p^{\min\{\nu_p(d),\nu_p(\varphi(m))-\nu_p(d)\}}. \]

Enumerating all possibilities for \(\nu_p(d),\nu_p(\varphi(m))\) shows that the exponent of \(p\) in this product is always \(1\), so this yields the unique large prime factor \(p^\star\) (if it exists). As for the remaining part \(z/p^\star\), since it can only contain prime factors with \(e=s\), no further factorization is needed.

A reference implementation for finding \(k\)-th roots modulo a prime is as follows:

Reference Implementation for Problem Library Checker - Kth Root (Mod)
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
#include <algorithm>
#include <cmath>
#include <iostream>
#include <random>
#include <tuple>
#include <unordered_map>
#include <vector>

std::mt19937 rng(std::random_device{}());

// Binary exponentiation.
int pow(int a, int b, int m = 0) {
  int res = 1;
  for (; b; b >>= 1) {
    if (b & 1) res = m ? (long long)res * a % m : res * a;
    a = m ? (long long)a * a % m : a * a;
  }
  return res;
}

// Find a P-th non-residue mod M.
int non_residue(int p, int m, int phi) {
  std::uniform_int_distribution<int> dis(1, m - 1);
  while (true) {
    int c = dis(rng);
    if (pow(c, phi / p, m) != 1) return c;
  }
  return -1;
}

// Euclidean Algorithm.
int gcd(int a, int b) { return b ? gcd(b, a % b) : a; }

// Extended Euclidean Algorithm.
int ex_gcd(int a, int b, int& x, int& y) {
  if (!b) {
    x = 1;
    y = 0;
    return a;
  } else {
    int d = ex_gcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
  }
}

// Returns the modular inverse of A modulo M.
// Assumes that gcd(A, M) = 1, so the inverse exists.
int inv(int a, int m) {
  int x, y;
  ex_gcd(a, m, x, y);
  return (x % m + m) % m;
}

// Subroutine: Find a P^E-th root of A mod M.
int peth_root_mod_m(int p, int e, int a, int m, int phi) {
  int s = 0, r = phi, pe = pow(p, e);
  for (; r % p == 0; r /= p, ++s);
  int q = pe - inv(r, pe);
  int ans = pow(a, ((long long)q * r + 1) / pe % phi, m);
  int eta = non_residue(p, m, phi);
  std::unordered_map<int, int> mp;
  int zeta = pow(eta, r, m);
  int xi = pow(eta, phi / p, m);
  // Precompute powers for BSGS.
  int B = std::sqrt((s - e) * p + 0.25l) + 1;
  int pB = p / B + 1;
  int po0 = pow(xi, pB, m);
  for (int j = 1, po1 = 1; j <= B; ++j) {
    po1 = (long long)po1 * po0 % m;
    mp[po1] = j;
  }
  // Compute p-adic digits of h.
  for (int j = 0; j < s - e; ++j) {
    int err = (long long)pow(ans, pe, m) * inv(a, m) % m;
    int xi_hj = pow(err, pow(p, s - e - j - 1), m);
    long long hj = 0;
    // BSGS query.
    for (int i = 1; i <= pB; ++i) {
      xi_hj = (long long)xi_hj * xi % m;
      if (mp.count(xi_hj)) {
        hj = mp[xi_hj] * pB - i;
        break;
      }
    }
    ans = (long long)ans * pow(zeta, phi - hj * pow(p, j) % phi, m) % m;
  }
  return ans;
}

// Find a K-th root of A modulo prime P.
int kth_root_mod_p(int k, int a, int p) {
  a %= p;
  if (k == 0) return a == 1 ? 0 : -1;
  if (a == 0) return 0;
  int d = gcd(k, p - 1);
  if (pow(a, (p - 1) / d, p) != 1) return -1;
  a = pow(a, inv(k / d, (p - 1) / d), p);
  for (int dp = 2; dp * dp <= d && dp * dp * dp * dp < p; ++dp) {
    if (d % dp == 0) {
      int de = 0;
      for (; d % dp == 0; d /= dp, ++de);
      a = peth_root_mod_m(dp, de, a, p, p - 1);
    }
  }
  if (d != 1) {
    int dp = gcd(d, (p - 1) / d), de = 0;
    if (dp != 1) {
      for (; d % dp == 0; d /= dp, ++de);
      a = peth_root_mod_m(dp, de, a, p, p - 1);
    }
    if (d != 1) a = peth_root_mod_m(d, 1, a, p, p - 1);
  }
  return a;
}

void solve() {
  int k, y, p;
  std::cin >> k >> y >> p;
  std::cout << kth_root_mod_p(k, y, p) << '\n';
}

int main() {
  int t;
  std::cin >> t;
  for (; t; --t) {
    solve();
  }
}

Handling the General Case

Consider the general case where the modulus \(m\) is a prime power \(p^e\), but \(\gcd(a,m)>1\). If \(a\equiv 0\pmod{m}\), then

\[ x = p^{\lceil e/k \rceil}\ell\pmod{p^e},~\ell=0,1,\cdots,p^{e-\lceil e/k\rceil}-1 \]

are all solutions to the original equation. Next, consider the case \(a\not\equiv 0\pmod{m}\). Let \(a = p^sa'\) with \(p\perp a'\). Then, setting \(x=p^zx'\) with \(p\perp x'\), we have

$$