<|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 \equiv z^n \pmod p
\]

with

\[
1 \le x,y,z \le 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 needed to solve the problem

### Observation 1: We may assume `z = 1`

If `z` is nonzero modulo prime `p`, then `z` has an inverse modulo `p`.

Starting from

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

divide 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.

---

### 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 \pmod 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} \pmod p\}
\]

This subgroup has size

\[
\frac{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 \equiv 1 \pmod 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 \equiv 1 - cur \pmod p
\]

has already appeared in `H`.

If yes, then

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

and both are `n`-th power residues.

---

### Observation 5: Recovering the actual root

Suppose we found

\[
cur = rt^{gk}
\]

We need an `x` such that

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

Let

\[
x = rt^e
\]

Then

\[
(rt^e)^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}}
\]

Since

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

the modular inverse exists:

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

Then

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

---

## 3. Full solution approach based on the observations

For each test case:

1. If `p == 2`, output `-1`.

   The only nonzero residue modulo `2` is `1`.

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

   but

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

   so no solution exists.

2. Compute

   ```cpp
   g = gcd(n, p - 1)
   ```

3. Find a primitive root `rt` modulo `p`.

4. The subgroup of all nonzero `n`-th powers is generated by

   ```cpp
   step = rt^g mod p
   ```

5. Enumerate the subgroup:

   ```text
   cur = step^1, step^2, step^3, ...
   ```

   For every `cur`, store which exponent index `k` produced it:

   \[
   cur = rt^{gk}
   \]

6. For each `cur`, compute:

   ```cpp
   other = (1 - cur + p) % p
   ```

   If `other` was already seen, then we found

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

7. Recover roots for both values using:

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

   Then

   ```cpp
   x = rt^e mod p
   ```

8. Output:

   ```text
   x y 1
   ```

9. If the subgroup cycle returns to `1` before finding a pair, output `-1`.

### Complexity

Let

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

be the number of nonzero `n`-th power residues.

For each test case:

- Subgroup enumeration: `O(m)`
- Root recovery: `O(log p)`
- Primitive root search is fast for `p <= 10^6`.

With precomputed smallest prime factors, factoring `p - 1` is also fast.

---

## 4. C++ implementation with detailed comments

```cpp
#include <bits/stdc++.h>
using namespace std;

using int64 = long long;

/*
    Fast modular exponentiation.

    Computes:
        a^e mod mod

    in O(log e).
*/
int64 modpow(int64 a, int64 e, int64 mod) {
    int64 result = 1;

    while (e > 0) {
        if (e & 1) {
            result = result * a % mod;
        }

        a = a * a % mod;
        e >>= 1;
    }

    return result;
}

/*
    Extended Euclidean algorithm.

    Returns x such that:
        a * x == 1 mod m

    Assumes gcd(a, m) = 1.

    If m == 1, all numbers are congruent modulo 1,
    so returning 0 is safe.
*/
int64 mod_inverse(int64 a, int64 m) {
    if (m == 1) {
        return 0;
    }

    int64 b = m;
    int64 u = 1, v = 0;

    while (b != 0) {
        int64 t = a / b;

        a -= t * b;
        swap(a, b);

        u -= t * v;
        swap(u, v);
    }

    u %= m;
    if (u < 0) {
        u += m;
    }

    return u;
}

/*
    Build smallest prime factor array up to maxN.

    spf[x] = smallest prime divisor of x.
*/
vector<int> build_spf(int maxN) {
    vector<int> spf(maxN + 1, 0);

    for (int i = 2; i <= maxN; i++) {
        if (spf[i] == 0) {
            spf[i] = i;

            if (1LL * i * i <= maxN) {
                for (int64 j = 1LL * i * i; j <= maxN; j += i) {
                    if (spf[j] == 0) {
                        spf[j] = i;
                    }
                }
            }
        }
    }

    return spf;
}

/*
    Factor x using the precomputed smallest prime factor array.

    Returns the distinct prime factors of x.
*/
vector<int> distinct_prime_factors(int x, const vector<int>& spf) {
    vector<int> factors;

    while (x > 1) {
        int q = spf[x];
        factors.push_back(q);

        while (x % q == 0) {
            x /= q;
        }
    }

    return factors;
}

/*
    Find a primitive root modulo prime p.

    A candidate r is a primitive root modulo p if for every prime factor q
    of phi = p - 1:

        r^(phi / q) != 1 mod p
*/
int primitive_root(
    int p,
    const vector<int>& spf,
    vector<int>& root_cache
) {
    if (p == 2) {
        return 1;
    }

    if (root_cache[p] != 0) {
        return root_cache[p];
    }

    int phi = p - 1;
    vector<int> factors = distinct_prime_factors(phi, spf);

    for (int r = 2; r < p; r++) {
        bool ok = true;

        for (int q : factors) {
            if (modpow(r, phi / q, p) == 1) {
                ok = false;
                break;
            }
        }

        if (ok) {
            root_cache[p] = r;
            return r;
        }
    }

    // This should never happen for prime p.
    return -1;
}

/*
    Given that:
        value = rt^(g * k)

    recover x such that:
        x^n = value mod p

    We need:
        x = rt^e

    Then:
        e * n == g * k mod (p - 1)

    Divide by g:
        e * (n / g) == k mod ((p - 1) / g)

    Therefore:
        e = k * inverse(n / g) mod ((p - 1) / g)
*/
int recover_root(
    int k,
    int n,
    int p,
    int g,
    int rt
) {
    int n_div_g = n / g;
    int phi_div_g = (p - 1) / g;

    int64 inv = mod_inverse(n_div_g, phi_div_g);
    int64 e = 1LL * (k % phi_div_g) * inv % phi_div_g;

    return (int)modpow(rt, e, p);
}

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

    int T;
    cin >> T;

    vector<pair<int, int>> tests(T);
    int maxP = 2;

    for (int i = 0; i < T; i++) {
        int n, p;
        cin >> n >> p;

        tests[i] = {n, p};
        maxP = max(maxP, p);
    }

    /*
        Precompute smallest prime factors once.

        Since p <= 10^6, this is cheap.
    */
    vector<int> spf = build_spf(maxP);

    /*
        Cache primitive roots for repeated primes.
    */
    vector<int> root_cache(maxP + 1, 0);

    /*
        Instead of clearing a seen array of size p for every test case,
        we use timestamps.

        seen_stamp[value] == current_stamp means value was seen
        in the current test case.

        seen_k[value] stores the exponent index k such that:
            value = rt^(g * k)
    */
    vector<int> seen_stamp(maxP + 1, 0);
    vector<int> seen_k(maxP + 1, 0);
    int current_stamp = 0;

    for (auto [n, p] : tests) {
        current_stamp++;

        /*
            Special case p = 2.

            Only nonzero residue is 1:
                1^n + 1^n = 0 mod 2
            but:
                1^n = 1 mod 2

            Therefore no solution exists.
        */
        if (p == 2) {
            cout << -1 << '\n';
            continue;
        }

        int g = std::gcd(n, p - 1);

        /*
            rt is a primitive root modulo p.
        */
        int rt = primitive_root(p, spf, root_cache);

        /*
            The subgroup of n-th powers is generated by:
                rt^g mod p
        */
        int step = (int)modpow(rt, g, p);

        /*
            Enumerate:
                cur = step^1, step^2, ...

            If cur = rt^(g * cnt), then cnt is stored.
        */
        int64 cur = 1;
        int cnt = 0;

        bool found = false;

        while (true) {
            cur = cur * step % p;
            cnt++;

            /*
                Store the current subgroup element.
            */
            seen_stamp[cur] = current_stamp;
            seen_k[cur] = cnt;

            /*
                If we returned to 1, we finished the entire subgroup cycle.

                The only possible complement of 1 is 0, but 0 is not
                an allowed nonzero n-th power residue.
            */
            if (cur == 1) {
                break;
            }

            /*
                We need:
                    cur + other == 1 mod p

                Hence:
                    other = 1 - cur mod p
            */
            int other = (1 - (int)cur + p) % p;

            /*
                If other was already seen in this subgroup,
                we found two n-th power residues summing to 1.
            */
            if (seen_stamp[other] == current_stamp) {
                int kx = cnt;
                int ky = seen_k[other];

                int x = recover_root(kx, n, p, g, rt);
                int y = recover_root(ky, n, p, g, rt);

                /*
                    We normalized z to 1.
                */
                cout << x << ' ' << y << ' ' << 1 << '\n';

                found = true;
                break;
            }
        }

        if (!found) {
            cout << -1 << '\n';
        }
    }

    return 0;
}
```

---

## 5. Python implementation with detailed comments

```python
import sys
from math import gcd


def mod_inverse(a, m):
    """
    Return x such that:

        a * x == 1 mod m

    Assumes gcd(a, m) = 1.

    If m == 1, returning 0 is safe because every integer
    is congruent to 0 modulo 1.
    """

    if m == 1:
        return 0

    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

    return old_s % m


def build_spf(max_n):
    """
    Build smallest prime factor array.

    spf[x] = smallest prime divisor of x.
    """

    spf = [0] * (max_n + 1)

    for i in range(2, max_n + 1):
        if spf[i] == 0:
            spf[i] = i

            if i * i <= max_n:
                for j in range(i * i, max_n + 1, i):
                    if spf[j] == 0:
                        spf[j] = i

    return spf


def distinct_prime_factors(x, spf):
    """
    Factor x using the precomputed smallest prime factor array.

    Return only distinct prime factors.
    """

    factors = []

    while x > 1:
        q = spf[x]
        factors.append(q)

        while x % q == 0:
            x //= q

    return factors


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

    A candidate r is primitive if for every prime factor q of p - 1:

        r^((p - 1) / q) != 1 mod p
    """

    if p == 2:
        return 1

    if p in root_cache:
        return root_cache[p]

    phi = p - 1
    factors = distinct_prime_factors(phi, spf)

    r = 2

    while True:
        ok = True

        for q in factors:
            if pow(r, phi // q, p) == 1:
                ok = False
                break

        if ok:
            root_cache[p] = r
            return r

        r += 1


def recover_root(k, n, p, g, rt):
    """
    We know:

        value = rt^(g * k)

    We want x such that:

        x^n = value mod p

    Let:

        x = rt^e

    Then:

        e * n == g * k mod (p - 1)

    Divide by g:

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

    Therefore:

        e = k * inverse(n / g) mod ((p - 1) / g)
    """

    n_div_g = n // g
    phi_div_g = (p - 1) // g

    inv = mod_inverse(n_div_g, phi_div_g)
    e = (k % phi_div_g) * inv % phi_div_g

    return pow(rt, e, p)


def solve_case(
    n,
    p,
    spf,
    root_cache,
    seen_stamp,
    seen_k,
    current_stamp
):
    """
    Solve one test case.

    Returns a string:
        "-1"
    or:
        "x y z"
    """

    # Modulo 2, the only nonzero residue is 1.
    # 1^n + 1^n = 0 mod 2, but z^n = 1 mod 2.
    if p == 2:
        return "-1"

    g = gcd(n, p - 1)

    # Primitive root modulo p.
    rt = primitive_root(p, spf, root_cache)

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

    cur = 1
    cnt = 0

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

        # cur = rt^(g * cnt)
        seen_stamp[cur] = current_stamp
        seen_k[cur] = cnt

        # Completed full subgroup cycle.
        if cur == 1:
            return "-1"

        # Need cur + other == 1 mod p.
        other = (1 - cur) % p

        # If other was already seen, then both are n-th powers.
        if seen_stamp[other] == current_stamp:
            kx = cnt
            ky = seen_k[other]

            x = recover_root(kx, n, p, g, rt)
            y = recover_root(ky, n, p, g, rt)

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


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

    if not data:
        return

    t = int(data[0])
    tests = []

    max_p = 2
    idx = 1

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

        tests.append((n, p))
        max_p = max(max_p, p)

    # Precompute smallest prime factors up to the largest p.
    spf = build_spf(max_p)

    # Cache primitive roots for repeated primes.
    root_cache = {}

    # Timestamp arrays avoid clearing a large seen array every test case.
    seen_stamp = [0] * (max_p + 1)
    seen_k = [0] * (max_p + 1)

    answers = []

    for current_stamp, (n, p) in enumerate(tests, start=1):
        answers.append(
            solve_case(
                n,
                p,
                spf,
                root_cache,
                seen_stamp,
                seen_k,
                current_stamp
            )
        )

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


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