<|instruction|>
Solve the below problem. The solution should start with an abridged problem statement. Then key observations. Then full solution based on the observations. Then C++ and Python implementations with comments.

511. Fermat's Last Theorem
Time limit per test: 0.75 second(s)
Memory limit: 262144 kilobytes
input: standard
output: standard

Given a positive integer n and a positive prime number p, find x, y and z such that xn+yn=zn modulo p and x, y and z are nonzero modulo p or report that there's no such triple.

Input
The first line of the input file contains the number t of testcases to solve, 1 ≤ t ≤ 1000. Each of the next t lines contains two integers n and p, 3 ≤ n ≤ 106, 2 ≤ p ≤ 106.

Output
For each input testcase, output one line:
when there exists a solution, output three integers x, y and z, 1 ≤ x, y, z ≤ p-1. If there are multiple solutions, output any.
when there's no solution, output one integer -1.


Example(s)
sample input
sample output
2
5 41
3 5
-1
1 2 4

<|response|>
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, output `-1`.

Constraints: 1 <= t <= 1000, 3 <= n <= 10^6, 2 <= p <= 10^6, p is prime.

2. Key Observations

Observation 1: We may assume z = 1

If z is nonzero modulo prime p, then z has an inverse. Dividing by z^n gives (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.

Observation 2: Nonzero residues modulo prime p form a cyclic group

Modulo a prime p, the nonzero residues form a multiplicative group of size p-1. This group is cyclic, so there exists a primitive root rt such that every nonzero residue can be written as rt^k (mod p) for some integer k.

Observation 3: The set of all n-th powers is a subgroup

Let g = gcd(n, p-1). If x = rt^k, then x^n = rt^(kn). The possible exponents are exactly multiples of g, so the set of all nonzero n-th powers is H = {rt^(gk) mod p}. This subgroup has size (p-1)/g and is generated by rt^g.

Observation 4: We only need two elements of H summing to 1

We need a^n + b^n ≡ 1 (mod p). The values a^n and b^n must both belong to H. So we enumerate all values in H. For each current value `cur`, we check whether other ≡ 1 - cur (mod p) has already appeared in H. If yes, we found a pair.

Observation 5: Recovering the actual root

Suppose cur = rt^(gk). We need x such that x^n ≡ rt^(gk) (mod p). Let x = rt^e. Then (rt^e)^n = rt^(en), so en ≡ gk (mod p-1). Dividing 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).

3. Full Solution Approach

For each test case:
1. If p == 2, output -1 (only nonzero residue is 1, but 1+1=0 mod 2).
2. Compute g = gcd(n, p-1).
3. Find a primitive root rt modulo p (by checking all proper divisors of p-1).
4. The subgroup of n-th powers is generated by step = rt^g mod p.
5. Enumerate the subgroup: cur = step^1, step^2, ... For every cur, store its exponent index. For each cur, compute other = (1-cur+p) % p. If other was already seen, recover x and y as n-th roots and output x y 1.
6. If the subgroup cycle returns to 1 before finding a pair, output -1.

Complexity: O(sqrt(p)*log(p) + m*log(p)) per test case where m = (p-1)/g.

4. 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;
}
```

5. 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()
```
