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 ≡ z^n (mod 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 to solve x^n + y^n ≡ z^n (mod 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:

(x/z)^n + (y/z)^n ≡ 1 (mod p)

So it is enough to find nonzero a, b such that a^n + b^n ≡ 1 (mod p). Then (a, b, 1) is a valid answer.

Multiplicative group modulo p

For prime p, the nonzero residues modulo p form a cyclic multiplicative group of size p - 1. Let rt be a primitive root modulo p; every nonzero residue can be written as rt^k (mod p) for some k.

The set of all n-th powers is:
H = {x^n mod p : x ≠ 0}

If x = rt^k, then x^n = rt^(kn). So the possible exponents are all multiples of g = gcd(n, p-1), and H = {rt^(gk) mod p}. This is a subgroup generated by rt^g, with size (p-1)/g.

Searching for two elements of H summing to 1

Enumerate all elements of H: rt^g, rt^(2g), rt^(3g), ... For each current element `cur`, we need another element other ≡ 1 - cur (mod p). If `other` was already seen in the subgroup enumeration, then cur + other ≡ 1 (mod p), and we found our pair.

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 ≡ rt^(gk) (mod p). Let x = rt^e. Then x^n = rt^(en), so we need en ≡ gk (mod p-1). Let g = gcd(n, p-1). Divide by g: e*(n/g) ≡ k (mod (p-1)/g). Since gcd(n/g, (p-1)/g) = 1, the modular inverse exists:

e ≡ k * (n/g)^(-1) (mod (p-1)/g)

Then x = rt^e (mod p). The same is done for y.

Edge case: p = 2

Modulo 2, the only nonzero residue is 1. Then 1^n + 1^n ≡ 0 (mod 2) but 1^n ≡ 1 (mod 2), so no solution exists for p = 2.

Complexity per test case

Let m = (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 with trial division. Overall: O(sqrt(p) * log(p) + m * log(p)) per test case.

3. C++ Solution
```cpp
#include <bits/stdc++.h>

using namespace std;

template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

template<typename T1, typename T2>
istream& operator>>(istream& in, pair<T1, T2>& x) {
    return in >> x.first >> x.second;
}

template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for(auto& x: a) {
        in >> x;
    }
    return in;
};

template<typename T>
ostream& operator<<(ostream& out, const vector<T>& a) {
    for(auto x: a) {
        out << x << ' ';
    }
    return out;
};

int n, p;

void read() { cin >> n >> p; }

int64_t modpow(int64_t a, int64_t e, int64_t m) {
    int64_t r = 1;
    for(; e; e >>= 1) {
        if(e & 1) {
            r = r * a % m;
        }
        a = a * a % m;
    }
    return r;
}

int primitive_root(int p) {
    vector<int> divs;
    int phi = p - 1;
    for(int i = 1; (int64_t)i * i <= phi; i++) {
        if(phi % i == 0) {
            divs.push_back(i);
            if(i * i < phi) {
                divs.push_back(phi / i);
            }
        }
    }
    sort(divs.begin(), divs.end());
    for(int g = 2; g < p; g++) {
        bool ok = true;
        for(int i = 0; i + 1 < (int)divs.size() && ok; i++) {
            if(modpow(g, divs[i], p) == 1) {
                ok = false;
            }
        }
        if(ok) {
            return g;
        }
    }
    return -1;
}

int64_t mod_inverse(int64_t a, int64_t m) {
    int64_t g = m, x = 0, y = 1;
    for(int64_t r = a; r != 0;) {
        int64_t q = g / r;
        g -= q * r;
        swap(g, r);
        x -= q * y;
        swap(x, y);
    }
    return (x % m + m) % m;
}

int nth_root(int k, int rt, int g, int n_div_g, int phi_div_g, int p) {
    // Find x such that x^n = rt^(g*k), i.e., x = rt^(k * inverse(n/g, phi/g))
    int64_t e = mod_inverse(n_div_g, phi_div_g) * k % (p - 1);
    return modpow(rt, e, p);
}

void solve() {
    // Starting with a tiny bit of history, Fermat's last theorem was actually a
    // conjecture from 1637 saying that a^n + b^n = c^n has no integer solutions
    // for n > 2. It was finally proven by Andrew Wiles, but here the modular
    // setting makes things much easier. In particular, Schur proved in 1916
    // that for every n, there is some p_0, such that the above has a solution
    // for p >= p_0. That proof was not constructive, but later works show how
    // to give more insight. An example is:
    //
    //     https://www.scirp.org/pdf/apm20241410_35302479.pdf.
    //
    // Although not directly given, an algorithmic way of finding this would be:
    //
    //     1) The n-th powers mod p form a subgroup of (Z/pZ)*.
    //        Let rt be a primitive root of p. For any prime there is a
    //        primitive root, and also if we look at the smallest ones, they
    //        aren't huge - for example, under 10^6 the largest primitive root
    //        is 73 at p=760321. Checking if rt is a primitive root can be done
    //        by making sure there is no smaller cycle that phi(p) = p - 1. The
    //        complexity of finding this primitive root is then O(sqrt(p) *
    //        max_rt), which for the given constraints is quick. As rt is the
    //        primitive root, every nonzero element is rt^k for some k in [0,
    //        p-2]. Then (rt^k)^n = rt^(kn), and the image of the map x
    //        -> x^n is {rt^(kn) : k in [0, p-2]} = {rt^m : g | m} where g =
    //        gcd(n, p-1). This is a classic result about what values can be
    //        achieved from ax mod q, for a = n, q = phi(p), and x is the k.
    //        Clearly, this subgroup has size (p-1)/g and is generated by rt^g.
    //
    //     2) We can search for two elements a, b in this subgroup with a + b
    //        = 1 (mod p). If we find such a pair, then a = x^n and b = y^n for
    //        some x, y, and we have x^n + y^n = 1 = 1^n (mod p), giving
    //        solution (x, y, 1). Why is it enough to search for a + b = 1 (mod
    //        p)? Say we only had a + b = c. Then (a*c^-1) + (b*c^-1) = 1 mod p,
    //        so we know there is also some a' + b' = 1 (mod p). Note that c^-1
    //        (mod p) will always exist a p is a prime.
    //
    //     3) To find such a pair efficiently, let us iterate through powers of
    //        rt^g. Let st = (rt^g)^cnt for cnt = 1, 2, ... Store each st in a
    //        dictionary. For each st, check if (1 - st) mod p is already in the
    //        table. If so, we found a + b = 1. If we complete the full
    //        cycle (st = 1) without finding a pair, no solution exists.
    //
    //     4) To recover x from x^n = rt^(g*k), we need to find the n-th root.
    //        We want x = rt^e such that e*n = g*k (mod p-1). Dividing by g:
    //        e*(n/g) = k (mod (p-1)/g). Since gcd(n/g, (p-1)/g) = 1, the
    //        inverse exists: e = k * inverse(n/g, (p-1)/g), giving x = rt^e.
    //
    //     5) When does no solution exist? Essentially, when the subgroup of
    //        n-th powers has a small size. In particular, for any x^n, there is
    //        only one y^n that satisfies x^n + y^n = 1 mod p, which means the
    //        numbers [1; p) are partitioned into pairs. To not have a solution,
    //        the subgroup has to not have two numbers from the same pair. The
    //        subgroup is not quite random, but a way of thinking about this is
    //        that the process is similar to the birthday paradox, or in
    //        O(sqrt(N)) time we will either find a solution (match), or the
    //        cycle will be too small and we will repeat meaning we terminate.
    //
    //        For some empirical evidence, consider the following: we have a bag
    //        with N balls of N/2 colours (2 each). What's the expected number
    //        of draws until we get the same colour twice? Turns out E[T] =
    //        O(sqrt(N)), which we can either show analytically (compute P(T >
    //        k) and sum over k), or convince ourselves through simulations.
    //        This means that in O(sqrt(N)) time we will either find a solution,
    //        or the cycle will be too small and we will repeat meaning we
    //        terminate. As a reference, below is a table of how simulations
    //        look:
    //
    //                 m        Exact       Approx    Simulated     Err%
    //      ------------------------------------------------------------
    //                10       5.6755       5.6050       5.6709    1.24%
    //                50      12.5645      12.5331      12.6321    0.25%
    //               100      17.7467      17.7245      17.7046    0.12%
    //               500      39.6432      39.6333      39.9551    0.02%
    //              1000      56.0569      56.0499      56.1090    0.01%
    //              5000     125.3345     125.3314     124.9028    0.00%
    //             10000     177.2476     177.2454     179.4765    0.00%
    //
    //    *In the above, the names are:
    //        Exact:     E[T] = 2 + sum_{j=1}^{m-1} 2^j * C(m-1,j) / C(2m-1,j)
    //        Approx:    E[T] ~= sqrt(m*PI)
    //        Simulated: We actually simulate the process.

    //
    // Combining all of the above, we have a solution with time complexity of
    // O(sqrt(N+P) * log(N+P)).

    int g = gcd(n, p - 1);
    int rt = primitive_root(p);
    int rt_g = modpow(rt, g, p);

    vector<int> seen(p, -1);
    int64_t cur = 1;
    int cnt = 0;

    while(true) {
        cur = cur * rt_g % p;
        seen[cur] = ++cnt;
        if(cur == 1) {
            cout << -1 << "\n";
            return;
        }
        int other = (1 - cur % p + p) % p;
        if(seen[other] != -1) {
            int x = nth_root(cnt, rt, g, n / g, (p - 1) / g, p);
            int y = nth_root(seen[other], rt, g, n / g, (p - 1) / g, p);
            cout << x << " " << y << " " << 1 << "\n";
            return;
        }
    }
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    int T = 1;
    cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        // cout << "Case #" << test << ": ";
        solve();
    }

    return 0;
}
```

4. Python Solution
```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 ≡ 1 (mod p).

Let rt be a primitive root modulo p. The set of all n-th powers is H = {rt^(gk)} where g = gcd(n, p-1), and H is generated by rt^g.

Enumerate all elements of H. For each element `cur`, check whether other = 1 - cur (mod p) has already appeared. If yes, then cur + other = 1, and both are n-th power residues.

To recover the actual root of rt^(gk), solve (rt^e)^n = rt^(gk), so en ≡ gk (mod p-1). Divide by g: e*(n/g) ≡ k (mod (p-1)/g), giving e = k * (n/g)^(-1).

If the subgroup cycle ends without finding a pair, output -1. For p = 2, output -1 immediately.
