## 1. Abridged problem statement

For each test case, given an integer `n` and a prime `p`, find nonzero residues `x, y, z` modulo `p` such that

\[
x^n + y^n \equiv z^n \pmod p
\]

with `1 <= x, y, z <= p - 1`.

If no such triple exists, print `-1`.

Constraints:

- `1 <= t <= 1000`
- `3 <= n <= 10^6`
- `2 <= p <= 10^6`
- `p` is prime.

---

## 2. Detailed editorial

We need solve

\[
x^n + y^n \equiv z^n \pmod p
\]

where `x`, `y`, `z` are all nonzero modulo prime `p`.

### Key observation: normalize `z = 1`

If we have a solution with some nonzero `z`, then because `p` is prime, `z` has a modular inverse.

Divide the equation by `z^n`:

\[
\left(\frac{x}{z}\right)^n + \left(\frac{y}{z}\right)^n \equiv 1 \pmod p
\]

So it is enough to find nonzero `a`, `b` such that

\[
a^n + b^n \equiv 1 \pmod p
\]

Then `(a, b, 1)` is a valid answer.

So the task becomes: find two nonzero `n`-th power residues modulo `p` whose sum is `1`.

---

### Multiplicative group modulo `p`

For prime `p`, the nonzero residues modulo `p` form a cyclic multiplicative group of size

\[
\varphi(p) = p - 1
\]

Let `rt` be a primitive root modulo `p`.

That means every nonzero residue can be written as

\[
rt^k \pmod p
\]

for some `k`.

Now consider the set of all `n`-th powers:

\[
H = \{x^n \bmod p : x \neq 0\}
\]

If `x = rt^k`, then

\[
x^n = rt^{kn}
\]

So the possible exponents are all multiples of

\[
g = \gcd(n, p - 1)
\]

Therefore:

\[
H = \{rt^{gk} \bmod p\}
\]

This is a subgroup of the nonzero residues, generated by

\[
rt^g
\]

and its size is

\[
\frac{p - 1}{g}
\]

---

### Searching for two elements of `H` summing to `1`

We enumerate all elements of `H`:

\[
rt^g, rt^{2g}, rt^{3g}, \dots
\]

For each current element `cur`, we need another element:

\[
other \equiv 1 - cur \pmod p
\]

If `other` was already seen in the subgroup enumeration, then we found

\[
cur + other \equiv 1 \pmod p
\]

So both are `n`-th powers, and we can recover their `n`-th roots.

If we return to `1` without finding such a pair, then we have traversed the entire subgroup, so no solution exists.

---

### Recovering the actual roots

Suppose during enumeration we found

\[
cur = rt^{gk}
\]

We need some `x` such that

\[
x^n \equiv rt^{gk} \pmod p
\]

Let

\[
x = rt^e
\]

Then:

\[
x^n = rt^{en}
\]

So we need

\[
en \equiv gk \pmod {p - 1}
\]

Let:

\[
g = \gcd(n, p - 1)
\]

Divide by `g`:

\[
e \cdot \frac{n}{g} \equiv k \pmod {\frac{p - 1}{g}}
\]

Because

\[
\gcd\left(\frac{n}{g}, \frac{p - 1}{g}\right) = 1
\]

the modular inverse exists.

Thus:

\[
e \equiv k \cdot \left(\frac{n}{g}\right)^{-1}
\pmod {\frac{p - 1}{g}}
\]

Then

\[
x = rt^e \pmod p
\]

The same is done for `y`.

---

### Edge case: `p = 2`

Modulo `2`, the only nonzero residue is `1`.

Then:

\[
1^n + 1^n \equiv 1 + 1 \equiv 0 \pmod 2
\]

but

\[
1^n \equiv 1 \pmod 2
\]

So no solution exists for `p = 2`.

---

### Complexity

Let

\[
m = \frac{p - 1}{\gcd(n, p - 1)}
\]

be the size of the subgroup of `n`-th powers.

The algorithm enumerates at most `m` subgroup elements.

Finding a primitive root modulo `p` can be done by factoring `p - 1`.

Overall per test case:

\[
O(\sqrt p \log p + m \log p)
\]

Since `p <= 10^6`, this is practical.

---

## 3. Provided C++ solution with detailed comments

```cpp
#include <bits/stdc++.h> // Includes almost all standard C++ headers.

using namespace std;

// Output operator for pairs.
// Allows writing: cout << pair_value;
template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

// Input operator for pairs.
// Allows writing: cin >> pair_value;
template<typename T1, typename T2>
istream& operator>>(istream& in, pair<T1, T2>& x) {
    return in >> x.first >> x.second;
}

// Input operator for vectors.
// Reads every element of the vector.
template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for(auto& x: a) {
        in >> x;
    }
    return in;
};

// Output operator for vectors.
// Prints elements separated by spaces.
template<typename T>
ostream& operator<<(ostream& out, const vector<T>& a) {
    for(auto x: a) {
        out << x << ' ';
    }
    return out;
};

// Global variables for the current test case.
// n is the exponent, p is the prime modulus.
int n, p;

// Reads one test case.
void read() {
    cin >> n >> p;
}

// Fast modular exponentiation.
// Computes a^e mod m in O(log e).
int64_t modpow(int64_t a, int64_t e, int64_t m) {
    int64_t r = 1; // Current result.

    // Binary exponentiation loop.
    for(; e; e >>= 1) {
        // If the current bit of e is set, multiply the result by a.
        if(e & 1) {
            r = r * a % m;
        }

        // Square the base for the next bit.
        a = a * a % m;
    }

    // Return a^e mod m.
    return r;
}

// Finds a primitive root modulo prime p.
//
// A primitive root g modulo p has order p - 1.
// This implementation enumerates all divisors d of p - 1,
// and checks that g^d != 1 for every proper divisor d.
int primitive_root(int p) {
    vector<int> divs; // Divisors of phi = p - 1.

    int phi = p - 1; // Size of the multiplicative group modulo p.

    // Enumerate all divisors of phi.
    for(int i = 1; (int64_t)i * i <= phi; i++) {
        if(phi % i == 0) {
            divs.push_back(i);

            // Add the paired divisor phi / i if it is different.
            if(i * i < phi) {
                divs.push_back(phi / i);
            }
        }
    }

    // Sort divisors increasingly.
    sort(divs.begin(), divs.end());

    // Try candidates g = 2, 3, ..., p - 1.
    for(int g = 2; g < p; g++) {
        bool ok = true; // Whether g is primitive.

        // Check all proper divisors of phi.
        // The last divisor is phi itself, so skip it.
        for(int i = 0; i + 1 < (int)divs.size() && ok; i++) {
            // If g^d == 1 for a proper divisor d,
            // then the order of g is not phi.
            if(modpow(g, divs[i], p) == 1) {
                ok = false;
            }
        }

        // If no proper divisor produced 1, g has order phi.
        if(ok) {
            return g;
        }
    }

    // Should not happen for odd prime p.
    return -1;
}

// Computes modular inverse of a modulo m using extended Euclid.
//
// Assumes gcd(a, m) = 1.
// Returns x such that a * x == 1 mod m.
int64_t mod_inverse(int64_t a, int64_t m) {
    // This is an iterative extended Euclidean algorithm.
    // g and r are current remainders.
    // x and y track coefficients.
    int64_t g = m, x = 0, y = 1;

    // Continue until the current remainder becomes zero.
    for(int64_t r = a; r != 0;) {
        // Euclidean quotient.
        int64_t q = g / r;

        // Update remainders.
        g -= q * r;
        swap(g, r);

        // Update coefficients.
        x -= q * y;
        swap(x, y);
    }

    // Normalize inverse into [0, m - 1].
    return (x % m + m) % m;
}

// Recovers an n-th root of rt^(g*k).
//
// We want x such that:
//
//     x^n = rt^(g*k)
//
// Let x = rt^e. Then:
//
//     rt^(e*n) = rt^(g*k)
//
// Therefore:
//
//     e*n = g*k mod p-1
//
// Divide by g:
//
//     e*(n/g) = k mod (p-1)/g
//
// Since gcd(n/g, (p-1)/g) = 1,
// inverse(n/g) exists modulo (p-1)/g.
int nth_root(
    int k,            // We need the root of rt^(g*k).
    int rt,           // Primitive root modulo p.
    int g,            // gcd(n, p - 1). This parameter is not directly used here.
    int n_div_g,      // n / g.
    int phi_div_g,    // (p - 1) / g.
    int p             // Prime modulus.
) {
    // Compute e = k * inverse(n/g) modulo (p-1)/g.
    // The code reduces modulo p - 1, which is also safe for the exponent.
    int64_t e = mod_inverse(n_div_g, phi_div_g) * k % (p - 1);

    // Return rt^e mod p.
    return modpow(rt, e, p);
}

// Solves one test case.
void solve() {
    // Important edge case for the stated constraints:
    // If p == 2, the only nonzero residue is 1.
    // Then 1^n + 1^n = 0 mod 2, but z^n = 1,
    // so no solution exists.
    //
    // The originally provided source did not include this guard,
    // but it is necessary if p = 2 is present in tests.
    if(p == 2) {
        cout << -1 << "\n";
        return;
    }

    // g = gcd(n, p - 1).
    // The set of n-th powers modulo p is generated by rt^g.
    int g = gcd(n, p - 1);

    // Find a primitive root modulo p.
    int rt = primitive_root(p);

    // rt_g is the generator of the subgroup of n-th powers.
    int rt_g = modpow(rt, g, p);

    // seen[value] stores which exponent count produced this subgroup value.
    // -1 means not seen yet.
    vector<int> seen(p, -1);

    // Current subgroup element.
    // We start from 1, then repeatedly multiply by rt_g.
    int64_t cur = 1;

    // cnt means cur = rt_g^cnt = rt^(g*cnt).
    int cnt = 0;

    // Enumerate the subgroup of n-th powers.
    while(true) {
        // Move to the next element of the subgroup.
        cur = cur * rt_g % p;

        // Store the exponent index for this subgroup element.
        seen[cur] = ++cnt;

        // If we returned to 1, the full subgroup cycle is complete.
        // Since 1 + something = 1 would require something = 0,
        // and 0 is not allowed, no solution was found.
        if(cur == 1) {
            cout << -1 << "\n";
            return;
        }

        // We need another subgroup element other such that:
        //
        //     cur + other = 1 mod p
        int other = (1 - cur % p + p) % p;

        // If that complementary value was already seen,
        // then both cur and other are n-th powers.
        if(seen[other] != -1) {
            // cur = rt^(g*cnt), so recover x with x^n = cur.
            int x = nth_root(cnt, rt, g, n / g, (p - 1) / g, p);

            // other = rt^(g*seen[other]), so recover y with y^n = other.
            int y = nth_root(seen[other], rt, g, n / g, (p - 1) / g, p);

            // We normalized z to 1.
            cout << x << " " << y << " " << 1 << "\n";
            return;
        }
    }
}

int main() {
    // Fast I/O.
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    // Number of test cases.
    int T = 1;
    cin >> T;

    // Process all test cases.
    for(int test = 1; test <= T; test++) {
        // Read n and p.
        read();

        // Solve current test case.
        solve();
    }

    return 0;
}
```

---

## 4. Python solution with detailed comments

```python
import sys
from math import gcd


def mod_inverse(a, m):
    """
    Return modular inverse of a modulo m.

    That is, return x such that:

        a * x == 1 mod m

    Assumes gcd(a, m) = 1.

    If m == 1, every integer is congruent to 0 modulo 1.
    In this problem, that case is never needed for an actual root recovery,
    but returning 0 keeps the function safe.
    """
    if m == 1:
        return 0

    # Extended Euclidean algorithm.
    old_r, r = a, m
    old_s, s = 1, 0

    while r:
        q = old_r // r

        old_r, r = r, old_r - q * r
        old_s, s = s, old_s - q * s

    # old_s is the coefficient of a.
    return old_s % m


def primitive_root(p):
    """
    Find a primitive root modulo prime p.

    The multiplicative group modulo p has size phi = p - 1.

    A number g is a primitive root iff for every prime divisor q of phi:

        g^(phi / q) != 1 mod p
    """

    # For p = 2, 1 is the only nonzero residue.
    if p == 2:
        return 1

    phi = p - 1
    x = phi
    factors = []

    # Factor phi and collect distinct prime factors.
    d = 2
    while d * d <= x:
        if x % d == 0:
            factors.append(d)

            while x % d == 0:
                x //= d

        d += 1

    if x > 1:
        factors.append(x)

    # Try candidates starting from 2.
    g = 2
    while True:
        ok = True

        # If g^(phi/q) == 1 for some prime factor q,
        # then the order of g divides phi/q, so it is not primitive.
        for q in factors:
            if pow(g, phi // q, p) == 1:
                ok = False
                break

        if ok:
            return g

        g += 1


def solve_case(n, p):
    """
    Solve one test case.
    Return either "-1" or "x y z".
    """

    # Edge case: modulo 2, only nonzero residue is 1.
    # 1^n + 1^n = 0 mod 2, but z^n = 1,
    # so no solution exists.
    if p == 2:
        return "-1"

    # Let g = gcd(n, p - 1).
    # The set of n-th powers modulo p is generated by rt^g,
    # where rt is a primitive root modulo p.
    g = gcd(n, p - 1)

    # Primitive root modulo p.
    rt = primitive_root(p)

    # Generator of the subgroup of n-th powers.
    step = pow(rt, g, p)

    # Size-related values.
    n_div_g = n // g
    phi_div_g = (p - 1) // g

    # Inverse of n/g modulo (p-1)/g.
    #
    # This exists because after dividing by g,
    # the two numbers become coprime.
    inv_n_div_g = mod_inverse(n_div_g, phi_div_g)

    # seen[value] = k means:
    #
    #     value = rt^(g*k)
    #
    # So value is the k-th enumerated element of the subgroup.
    seen = {}

    # Start from identity and repeatedly multiply by step.
    cur = 1
    cnt = 0

    while True:
        # Move to next subgroup element:
        #
        #     cur = step^cnt = rt^(g*cnt)
        cur = (cur * step) % p
        cnt += 1

        # Mark this subgroup element as seen.
        seen[cur] = cnt

        # If cur becomes 1, we completed the whole subgroup cycle.
        # No pair summing to 1 was found.
        if cur == 1:
            return "-1"

        # We need:
        #
        #     cur + other == 1 mod p
        #
        # Therefore:
        #
        #     other == 1 - cur mod p
        other = (1 - cur) % p

        # If other is already in the subgroup, we found:
        #
        #     cur + other == 1 mod p
        if other in seen:
            kx = cnt
            ky = seen[other]

            # Recover x from:
            #
            #     x^n = rt^(g*kx)
            #
            # Let x = rt^e.
            # Need:
            #
            #     e * (n/g) == kx mod ((p-1)/g)
            #
            # Therefore:
            #
            #     e = kx * inverse(n/g) mod ((p-1)/g)
            ex = (kx * inv_n_div_g) % (p - 1)
            ey = (ky * inv_n_div_g) % (p - 1)

            x = pow(rt, ex, p)
            y = pow(rt, ey, p)

            # z is normalized to 1.
            return f"{x} {y} 1"


def main():
    data = sys.stdin.read().strip().split()

    if not data:
        return

    t = int(data[0])
    idx = 1

    ans = []

    for _ in range(t):
        n = int(data[idx])
        p = int(data[idx + 1])
        idx += 2

        ans.append(solve_case(n, p))

    sys.stdout.write("\n".join(ans))


if __name__ == "__main__":
    main()
```

---

## 5. Compressed editorial

Normalize the equation by dividing by `z^n`. Since `z` is nonzero modulo prime `p`, this is valid. Therefore it is enough to find nonzero `x, y` such that:

\[
x^n + y^n \equiv 1 \pmod p
\]

Let `rt` be a primitive root modulo `p`. Every nonzero residue is `rt^k`.

The set of all `n`-th powers is:

\[
H = \{rt^{nk}\}
\]

Let:

\[
g = \gcd(n, p - 1)
\]

Then:

\[
H = \{rt^{gk}\}
\]

So `H` is generated by `rt^g`.

Enumerate all elements of `H`. For each element `cur`, check whether:

\[
other = 1 - cur \pmod p
\]

has already appeared. If yes, then `cur + other = 1`, and both are `n`-th powers.

To recover the actual root of `rt^(gk)`, solve:

\[
(rt^e)^n = rt^{gk}
\]

so:

\[
en \equiv gk \pmod {p - 1}
\]

Divide by `g`:

\[
e \cdot \frac{n}{g} \equiv k \pmod {\frac{p - 1}{g}}
\]

Thus:

\[
e = k \cdot \left(\frac{n}{g}\right)^{-1}
\]

and the root is:

\[
rt^e \pmod p
\]

If the subgroup cycle ends without finding a pair, output `-1`.

For `p = 2`, output `-1` immediately.