Higher-Order Residues & Roots of Unity

Prerequisite knowledge: Discrete Logarithm

This article discusses higher-order residues and roots of unity in modular arithmetic, and introduces algorithms for extracting roots modulo a modulus.

Higher-Order Residues

A higher-order residue in modular arithmetic can be viewed as a question of whether it is possible to extract a higher-order root modulo a modulus. It is a generalization of a quadratic residue.

In other words, $a$ modulo $m$ has a $k$-th root if and only if $a$ modulo $m$ is a $k$-th residue.

Properties

As with quadratic residues, one can discuss deciding $k$-th residues, their number, and the number of $k$-th residue classes. As with other congruence equation problems, the Chinese remainder theorem can reduce them to the case of prime-power moduli. Depending on whether a primitive root exists, this further splits into odd prime-power moduli and powers of $2$.

The odd prime-power case is relatively simple. In fact, whenever a primitive root exists, we have the following result:

The case of powers of $2$ is more special. To handle it, we use a result about the structure of the reduced residue system modulo $2^e$: every odd $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}$. This gives the following result:

This completely resolves the problem of deciding $k$-th residues for different moduli. The Legendre symbol and quadratic reciprocity law for quadratic residues can also be generalized to higher-order residues, but this is difficult and requires concepts such as cyclotomic fields. In algebraic number theory, quadratic reciprocity is ultimately generalized by Artin reciprocity.

Roots of Unity

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

As in the complex case, we have the following definition:

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

When, modulo $m$, a primitive $k$-th root of unity exists, its algebraic properties agree with those of the primitive complex $k$-th root of unity $\omega_k$, so it can replace $\omega_k$ in various computations. For example, applying it to the fast Fourier transform gives the number-theoretic transform over the finite field¹.

Properties

Every (primitive) root of unity of every order exists in the complex field. In number theory, however, (primitive) roots of unity do not always exist.

These properties show that, relative to the case where primitive roots exist, modulo $m$ the $\lambda$-primitive root plays a similar foundational role. Unlike primitive roots, powers of the $\lambda$-primitive root cannot generate all roots of unity modulo $m$. Nevertheless, because $\lambda$-primitive roots have a non-negligible density², if we need to find a primitive $k$-th root of unity, we can first find a $\lambda$-primitive root randomly and then obtain a primitive $k$-th root of unity by taking a power.

If $a$ modulo $m$ has a root of degree $k$, one can use, modulo $m$, all roots of unity of degree $k$ to generate all roots of $a$ modulo $m$ of degree $k$.

Generating all $k$-th roots using the $k$-th roots of unity is analogous to using the solutions of a homogeneous linear system to generate the general solution of a nonhomogeneous linear system.

The discussion so far has been general. When a primitive root exists, the structure of the roots of unity is simpler:

This theorem shows that, when a primitive root exists, all $k$-th roots of unity form a cyclic group, and the primitive $k$-th roots of unity are generators of this cyclic group. We will see later that the Tonelli–Shanks algorithm exploits this fact to accelerate the discrete-logarithm part of root extraction.

Extracting Roots Modulo a Modulus

Finally, this article discusses how to find $k$-th roots. For $k=2$, there are many efficient algorithms for extracting square roots modulo a modulus. For general $k$, however, no polynomial-time algorithm is known. This section introduces two common algorithms that find a root in $O(m^{1/2})$ and $O(m^{1/4+\varepsilon})$ time, respectively, for degree $k$. The Chinese remainder theorem can always reduce the problem to prime-power moduli, so this section focuses on the prime-power case.

Naive Algorithm

The preceding analysis of the properties of $k$-th residues has already indicated a method for finding $k$-th roots modulo a prime power. Strictly speaking, the case solved above is when the radicand $a$ is coprime to the modulus $m$. The algorithm can be summarized as follows:

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

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

  Here, $\operatorname{ind}_g a$ can be found using the BSGS algorithm, and the solutions of the linear congruence equation yield all $a$'s $k$-th roots $x\equiv g^y\pmod m$.

  There is also a similar approach. Again let $x\equiv g^y\pmod m$; by rewriting,

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

  This becomes finding the discrete logarithm with base $g^k$ of $a$. A set of particular solutions can likewise be found using BSGS. The general solution can be obtained from the solution expression above: multiply the particular solution by every root of unity of degree $k$.

  Whichever approach is used, when the primitive root is known, the complexity of finding one solution is $O(m^{1/2})$. Since a primitive root can be found in $o(m^{1/2})$ time, the total complexity remains $O(m^{1/2})$.

• When $m=2^e$ and $e\in\mathbf N_+$, first find the exponents in $a\equiv (-1)^s5^r\pmod m$, namely $s,r$. Of these two 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}$$

  The value $r=\operatorname{ind}_5((-1)^sa)$ can be found using BSGS 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 linear system is easy to obtain, and $x=(-1)^z5^y$ is the desired root. The complexity of finding one solution with this algorithm is still $O(m^{1/2})$.

Of course, in the no-solution case, the criterion described above can be used to decide this quickly in $O(\log m)$ time, without checking during the solving process.

The following is a reference implementation for finding a $k$-th root modulo a prime (the code is only illustrative; its time complexity is too high to pass the problem):

// 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

A suitable generalization of the Tonelli–Shanks algorithm for extracting square roots modulo a modulus can solve root extraction modulo prime powers. One relatively direct generalization is the Adleman–Manders–Miller algorithm³, but its complexity is still not good enough⁴. This section introduces an improved Tonelli–Shanks algorithm proposed by sugarknri, Min_25, 37zigen, and others. It can find, in $O(m^{1/4+\varepsilon})$ time, a $k$-th root.

The central idea of the Tonelli–Shanks algorithm is to move the discrete-logarithm computation into a group of order $2^e$, thereby reducing the time complexity. Similarly, a discrete logarithm in a group of prime-power order $p^e$ can be computed relatively efficiently, but the complexity is $\Omega\left(\sqrt{p}\right)$. The Adleman–Manders–Miller algorithm decomposes the computation of a $k$-th root into discrete logarithms in several groups of prime-power order, but, limited by $k$'s largest prime factor $p_\text{max}(k)$, its complexity is still $\Omega\left(\sqrt{p_\text{max}(k)}\right)$. The algorithm in this section improves the process further, avoiding discrete logarithms for large prime factors and controlling the overall complexity at $O(m^{1/4+\varepsilon})$.

Process

Consider, for a prime power $m$, computing the root of $a$ of degree $k$, that is, solving the congruence equation:

For the case $m=2^e$, we additionally need $a\equiv 1\pmod{4}$, so that $a$ can be written as a power of $g=5$. As discussed earlier, modulo $2^e$, extracting a $k$-th root can always be reduced to this case. When handling this case, namely $2^e$, every occurrence of $\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 degree divides $\varphi(m)$. Let $d=\gcd(k,\varphi(m))$. By the properties of $k$-th residues, when $a$ modulo $m$ is a $k$-th residue, then $a$ is always, modulo $m$, a $\dfrac{\varphi(m)}{d}$-th root of unity. By the properties of roots of unity, for any $\ell\perp\dfrac{\varphi(m)}{d}$, the map $x\mapsto x^{\ell}$ is a bijection among the $\dfrac{\varphi(m)}{d}$-th roots of unity. Therefore, we can take

Raising both sides of the original congruence equation to the power $\ell$ gives

The leftmost congruence uses Euler's theorem and the following congruence relation (where $c\in\mathbf Z$):

For the transformed problem, consider the prime factorization of $d$:

Starting from $b=a^\ell$, for each $p^e\neq 1$, successively extract a $p^e$-th root; the result is the root of $b$ of degree $d$, which is the root of $a$ of degree $k$.

Finally, the problem becomes solving the following equation:

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

is a root of $b$ modulo $m$ of degree $p^e$.

To compute $x$, find, modulo $m$, a $p$-th nonresidue $\eta$. By the properties above, it is enough to choose $\eta\perp m$ at random and check that $\eta^{\varphi(m)/p}\bmod{m}\neq 1$. The density of such numbers is

Thus, in expectation, no more than $4$ random integers are needed to find one. Note that $\eta^{rp^{s-1}}\not\equiv 1\pmod m$ and $\eta^{rp^s}\equiv 1\pmod m$. Therefore, if we set $\zeta=\eta^r\bmod m$ and $\xi=\eta^{rp^{s-1}}\bmod m$, they are respectively primitive $p^s$-th and $p$-th roots of unity.

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

Compute these digits one by one. Once the first $j$ digits have been computed, we necessarily have

Thus, $h_j$ can be found by computing a discrete logarithm to the base $\xi$. To obtain better time complexity, use BSGS. In total, $(s-e)$ discrete logarithms must be computed. If $B$ powers of $\xi$ are precomputed, the time complexity of one discrete-logarithm computation is $O(p/B)$, and the total time complexity is

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

Time Complexity

The time complexity of this algorithm is $O(m^{1/4+\varepsilon})$. In this complexity discussion, we always assume that one multiplication takes $O(1)$ time, and that exponentiation always uses Euler's theorem to reduce the exponent, so computing any individual power takes $O(\log m)$ time.

First consider computing a single $p^e$-th root. Finding a $p$-th nonresidue requires checking an expected $O(1)$ numbers, for a total time complexity of $O(\log m)$. Computing each of $s,r,\zeta,\eta,b^{qr}$ takes only $O(\log m)$ time. When computing $h$, each digit requires $O(\log m)$ time for exponentiation, and there are $(s-e)$ digits, so the total time complexity is $O((s-e)\log m)$. As stated above, preprocessing for the discrete logarithms and the $(s-e)$ queries together take $O\left(\sqrt{(s-e)p}\right)$ time. Since $s-e\in O(\log m)$, computing a single $p^e$-th root takes $O(p^{1/2+\varepsilon})$ time. In particular, when $s=e$, the time complexity can be reduced further to $O(\log m)$.

We can then consider the total time complexity of the algorithm. Computing $\varphi(m),d,\ell$ all takes $O(\log m)$ time. Next, factor $d=\prod_p{p^e}$; this can be done in $O(m^{1/4})$ time using Pollard Rho. Finally, the total time for successively computing $p^e$-th roots is

Since $e < s$, every prime factor $p$ occurs in $\varphi(m)$ at least $2$ times, so we necessarily have $p < m^{1/2}$. Therefore, the total time complexity is $O(m^{1/4+\varepsilon})$.

In fact, in this situation we can still achieve $O(m^{1/4+\varepsilon})$ time without using Pollard Rho to factor the integer. It is enough to factor $d$ by trial division, enumerating only prime factors no greater than $m^{1/4}$. Let the integer obtained after removing these small prime factors be $z$. For $z$'s prime factor $p > m^{1/4}$, we necessarily have $\nu_p(\varphi(m)) < 4$, where $\nu_p(n)$ denotes, for $n$, the exponent of $p$ in the prime factorization of the integer. Since we only need to consider

the case described above: at most one prime factor $p$ can satisfy this condition; otherwise, each would occur in $\varphi(m)$ at least $2$ times, making the total product exceed $m$. To isolate this one (possibly existing) large prime factor, compute

Enumerating all possible values of $\nu_p(d),\nu_p(\varphi(m))$ shows that the exponent of $p$ in the product is always $1$, so the result is exactly the unique large prime factor $p^\star$ (if one exists). The remaining part $z/p^\star$ can contain only prime factors satisfying $e=s$, so it does not need to be factored further.

The following is a reference implementation for finding a $k$-th root modulo a prime:

#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, still with modulus $m$, a prime power $p^e$, but now $\gcd(a,m)>1$. If $a\equiv 0\pmod{m}$, then

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

Since $(x')^k\perp p$, this holds if and only if $kz = s$ and $(x')^k\equiv a'\pmod{p^{e-s}}$. The first equation has a solution if and only if $k\mid s$, in which case $z=\dfrac{s}{k}$; the second equation has already been solved. Note that the modulus of the general solution of the second equation differs from that of the general solution of the original equation, so each solution $x'$ of the second equation corresponds to several solutions of the original equation:

The following is a reference implementation for finding all $k$-th roots modulo an arbitrary modulus:

#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 odd prime power.
int primitive_root(PrimePower pp) {
  std::vector<int> exp;
  int phi = pp.pe / pp.p * (pp.p - 1);
  for (auto factor : factorize(pp.p - 1)) {
    exp.push_back(phi / factor.p);
  }
  if (pp.e != 1) exp.push_back(phi / pp.p);
  int ans = 0;
  bool succ = false;
  while (!succ) {
    ++ans;
    succ = true;
    for (int b : exp) {
      if (pow(ans, b, pp.pe) == 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;
  }
}

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

// 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 all the K-th roots with a primitive root G known.
std::vector<int> calc(int g, int k, int a, int p, int pe) {
  int ind = log(g, a, pe);
  if (ind == -1) return {};
  int mm = p == 2 ? pe / 4 : pe / p * (p - 1);
  int y0, d;
  std::tie(y0, d) = solve_linear(k, ind, mm);
  if (y0 == -1) return {};
  int ans = pow(g, y0, pe), po = pow(g, d, pe);
  std::vector<int> res(mm / d);
  for (auto& x : res) {
    x = ans;
    ans = (long long)ans * po % pe;
  }
  return res;
}

// Find all the K-th roots of A modulo prime power P^E.
std::vector<int> kth_roots_mod_pe(int k, int a, PrimePower pp) {
  int p = pp.p, e = pp.e, pe = pp.pe;
  a %= pe;
  if (a == 0) {
    int d = pow(p, (e - 1) / k + 1);
    std::vector<int> res(pe / d);
    for (int i = 0; i < pe / d; ++i) {
      res[i] = i * d;
    }
    return res;
  }
  int s = 0;
  for (; a % p == 0; a /= p, ++s);
  if (s % k) return {};
  int psk = pow(p, s / k), pss = pow(p, s - s / k), pes = pow(p, e - s);
  std::vector<int> res;
  if (p != 2) {
    int g = primitive_root(PrimePower(p, e - s, pes));
    res = calc(g, k, a, p, pes);
  } else if (pes == 2) {
    res.push_back(a);
  } else if (k & 1) {
    int z = a % 4 == 3;
    a = z ? pes - a : a;
    res = calc(5, k, a, p, pes);
    if (z) {
      for (auto& x : res) x = pes - x;
    }
  } else {
    if (a % 4 == 3) return {};
    res = calc(5, k, a, p, pes);
    int m = res.size();
    res.reserve(m * 2);
    for (int i = 0; i < m; ++i) {
      res.push_back(pes - res[i]);
    }
  }
  int m = res.size();
  res.reserve(m * pss);
  for (int j = 1; j < pss; ++j) {
    for (int i = 0; i < m; ++i) {
      res.push_back(res.end()[-m] + pes);
    }
  }
  for (auto& x : res) x *= psk;
  return res;
}

// Find all the K-th roots of A modulo positive integer M.
std::vector<int> kth_roots_mod_m(int k, int a, int m) {
  auto factors = factorize(m);
  int m0 = 0;
  std::vector<std::vector<int>> sols;
  for (const auto& pp : factors) {
    sols.push_back(kth_roots_mod_pe(k, a, pp));
    if (sols.back().empty()) return {};
  }
  std::vector<int> ans;
  for (int i = 0; i < (int)factors.size(); ++i) {
    auto pp = factors[i];
    if (!i) {
      m0 = pp.pe;
      ans = sols[i];
    } else {
      long long m1 = pp.pe * inv(pp.pe, m0);
      long long m2 = m0 * inv(m0, pp.pe);
      m0 *= pp.pe;
      std::vector<int> _ans;
      _ans.reserve(ans.size() * sols[i].size());
      for (auto x : ans) {
        for (auto y : sols[i]) {
          _ans.push_back((m1 * x + m2 * y) % m0);
        }
      }
      ans.swap(_ans);
    }
  }
  return ans;
}

void solve() {
  int n, m, k;
  std::cin >> n >> m >> k;
  auto ans = kth_roots_mod_m(n, k, m);
  if (ans.empty()) {
    std::cout << 0 << '\n';
    return;
  }
  std::cout << ans.size() << '\n';
  std::sort(ans.begin(), ans.end());
  for (auto x : ans) std::cout << x << ' ';
  std::cout << '\n';
}

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