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

261. Discrete Roots
time limit per test: 1 sec.
memory limit per test: 65536 KB
input: standard
output: standard



There are a lot of mysteries and legends around computer science. One of the stories tells us about three Russian hackers who know the secret of breaking down widely used cryptographic algorithm. The fact itself threatens security of economics of many countries. Until recent time nobody knew anything about these hackers but now federal security bureau knows their names (Roman, Sergey and Andrew) and they also know that their hack method somehow uses discrete roots finding algorithm. And of course nobody knows this algorithm. We suggest you to try to solve much simpler task.
Given two prime numbers P and K (2 <= P <= 10^9, 2 <= K <= 100000) and integer number A (0 <= A < P) you are to find all the roots of the equation x^K = A mod P.

Input
Integer numbers P, K, A.

Output
The first line of output should contain number of roots of the equation. On the second line all the roots should be listed in ascending order.
Note: all the roots should be in the range [0..P-1].

Sample test(s)

Input
11 3 8

Output
1
2
Author:	Igor A. Kulkin
Resource:	Saratov SU Contest: Golden Fall 2004
Date:	October 2, 2004

<|response|>
1. Abridged Problem Statement
Given a prime P, integers K and A (0 ≤ A < P), find all x in [0..P−1] satisfying
 x^K ≡ A  (mod P).
Output the count of solutions and the solutions in ascending order.

2. Key Observations
- If A = 0, then x^K ≡ 0 mod P ⇒ x ≡ 0 is the only solution (since P is prime).
- Let φ = P−1. In the multiplicative group modulo P (which is cyclic of order φ), any nonzero x can be written as g^y for a primitive root g.
- x^K ≡ A  ⇔ (g^y)^K ≡ A ⇔ g^{K·y} ≡ A.
  Taking discrete logarithm base g: if g^T ≡ A, then we need K·y ≡ T (mod φ).
- The linear congruence K·y ≡ T (mod φ) has solutions iff d = gcd(K, φ) divides T. If so, there are exactly d solutions mod φ, spaced by φ/d.
- Each solution y gives x = g^y (mod P).

3. Full Solution Approach
Step 1: Handle A=0.
  • If A==0, output 1 and “0”. Done.

Step 2: Compute d = gcd(K, P−1).
  • A necessary condition for any solution is A^{(P−1)/d} ≡ 1 (mod P).
  • If this fails, output 0.

Step 3: Find a primitive root g modulo P.
  • Factor φ = P−1 into distinct prime factors.
  • Test candidates g=2,3,…: g is primitive iff for every prime factor q of φ,
    g^{φ/q} mod P ≠ 1.

Step 4: Compute discrete logarithm T = log_g(A) (mod P) via Baby-Step Giant-Step.
  • Let m = ⌈√φ⌉. Precompute baby steps: for j=0..m−1 store (g^j mod P) → j in a hash.
  • Compute factor = g^{−m} mod P = g^{φ−m}. Then for i=0..m:
      check if (A·factor^i mod P) is in the baby‐step table; if yes, recover T = i·m + j.

Step 5: Solve linear congruence K·y ≡ T (mod φ).
  • If T mod d ≠ 0, no solutions → output 0.
  • Else let K' = K/d, φ' = φ/d, T' = T/d. Compute inv = (K')^{−1} mod φ' via extended GCD.
  • One solution y0 = T'·inv mod φ'. All solutions are y_i = y0 + i·φ' for i=0..d−1.

Step 6: Convert y_i to x_i = g^{y_i} mod P.
  • Collect all x_i, sort them, remove duplicates (though they should already be distinct), and output.

Overall complexity is dominated by factoring φ (up to √φ) and BSGS (O(√φ)), which is acceptable for P up to 10^9.

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

int64_t p, k, a;

void read() { cin >> p >> k >> a; }

int64_t mul_mod(int64_t x, int64_t y, int64_t mod) {
    return (int64_t)((__int128)x * y % mod);
}

int64_t ext_gcd(int64_t a, int64_t b, int64_t& x, int64_t& y) {
    if(b == 0) {
        x = 1;
        y = 0;
        return a;
    }

    int64_t x1, y1;
    int64_t g = ext_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

int64_t inv_mod(int64_t a, int64_t m) {
    int64_t x, y;
    int64_t g = ext_gcd(((a % m) + m) % m, m, x, y);
    if(g != 1) {
        return -1;
    }
    return ((x % m) + m) % m;
}

int64_t mod_exp(int64_t base, int64_t exp, int64_t mod) {
    int64_t result = 1;
    base %= mod;
    if(base < 0) {
        base += mod;
    }

    while(exp > 0) {
        if(exp & 1) {
            result = mul_mod(result, base, mod);
        }
        base = mul_mod(base, base, mod);
        exp >>= 1;
    }
    return result;
}

vector<int64_t> factorize(int64_t n) {
    vector<int64_t> factors;
    for(int64_t i = 2; i * i <= n; i++) {
        if(n % i == 0) {
            factors.push_back(i);
            while(n % i == 0) {
                n /= i;
            }
        }
    }
    if(n > 1) {
        factors.push_back(n);
    }
    return factors;
}

int64_t find_primitive_root(int64_t mod) {
    if(mod == 2) {
        return 1;
    }

    int64_t phi = mod - 1;
    vector<int64_t> prime_factors = factorize(phi);
    for(int64_t g = 2; g < 200000; g++) {
        bool ok = true;
        for(int64_t f: prime_factors) {
            if(mod_exp(g, phi / f, mod) == 1) {
                ok = false;
                break;
            }
        }
        if(ok) {
            return g;
        }
    }
    return -1;
}

int64_t baby_step_giant_step(int64_t base, int64_t target, int64_t mod) {
    if(target % mod == 1) {
        return 0;
    }

    base %= mod;
    target %= mod;

    int64_t m = (int64_t)sqrtl((long double)(mod - 1)) + 2;

    unordered_map<int64_t, int64_t> baby;
    int64_t cur = 1;
    for(int64_t j = 0; j < m; j++) {
        if(!baby.count(cur)) {
            baby[cur] = j;
        }
        cur = mul_mod(cur, base, mod);
    }

    int64_t a_inv_m = mod_exp(base, (mod - 1) - m, mod);

    int64_t giant = target;
    for(int64_t i = 0; i < m; i++) {
        auto it = baby.find(giant);
        if(it != baby.end()) {
            return i * m + it->second;
        }
        giant = mul_mod(giant, a_inv_m, mod);
    }
    return -1;
}

void solve() {
    // Solve x^K = A (mod P) for prime P.
    //
    // - If A == 0 the only root is 0.
    //
    // - Let d = gcd(K, P-1). A solution exists iff A^((P-1)/d) == 1 (mod P).
    //
    // - Take a primitive root alpha of P (found by checking candidates against
    //   the prime factors of P-1). Solve the discrete log T with alpha^T == A
    //   via baby-step giant-step, so x = alpha^L means we need K*L == T
    //   (mod P-1).
    //
    // - That linear congruence is solvable iff d | T. Reduce to K'*L == T'
    //   (mod M) with T'=T/d, M=(P-1)/d, K'=K/d, giving L0 = T'*inv(K') (mod M).
    //   The d distinct roots are alpha^(L0 + i*M) for i in [0, d), which we
    //   sort and deduplicate.

    if(a == 0) {
        cout << 1 << "\n";
        cout << 0 << "\n";
        return;
    }

    int64_t gx, gy;
    int64_t d = ext_gcd(k, p - 1, gx, gy);
    if(mod_exp(a, (p - 1) / d, p) != 1) {
        cout << 0 << "\n";
        return;
    }

    int64_t alpha = find_primitive_root(p);
    if(alpha == -1) {
        cout << 0 << "\n";
        return;
    }

    int64_t t = baby_step_giant_step(alpha, a, p);
    if(t == -1 || t % d != 0) {
        cout << 0 << "\n";
        return;
    }

    int64_t t_prime = t / d;
    int64_t m = (p - 1) / d;
    int64_t k_prime = ((k / d) % m + m) % m;

    int64_t inv_k_prime = (m == 1) ? 0 : inv_mod(k_prime, m);
    if(inv_k_prime == -1) {
        cout << 0 << "\n";
        return;
    }

    int64_t l0 = (m == 1) ? 0 : mul_mod(t_prime % m, inv_k_prime, m);

    set<int64_t> roots;
    for(int64_t i = 0; i < d; i++) {
        int64_t exponent = (l0 + i * m) % (p - 1);
        roots.insert(mod_exp(alpha, exponent, p));
    }

    cout << roots.size() << "\n";
    bool first = true;
    for(int64_t x: roots) {
        if(!first) {
            cout << ' ';
        }
        cout << x;
        first = false;
    }
    cout << "\n";
}

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 Implementation
```python
import sys
from math import gcd, isqrt

# Fast modular exponentiation
def modexp(base, exp, mod):
    result = 1
    base %= mod
    while exp > 0:
        if exp & 1:
            result = (result * base) % mod
        base = (base * base) % mod
        exp >>= 1
    return result

# Extended GCD: returns (g, x, y) s.t. a*x + b*y = g = gcd(a,b)
def extgcd(a, b):
    if b == 0:
        return (a, 1, 0)
    g, x1, y1 = extgcd(b, a % b)
    # back-substitute
    x = y1
    y = x1 - (a // b) * y1
    return (g, x, y)

# Modular inverse of a mod m
def invmod(a, m):
    g, x, _ = extgcd(a, m)
    if g != 1:
        return None
    return x % m

# Factor n into distinct prime factors
def factorize(n):
    fac = []
    i = 2
    while i * i <= n:
        if n % i == 0:
            fac.append(i)
            while n % i == 0:
                n //= i
        i += 1
    if n > 1:
        fac.append(n)
    return fac

# Find a primitive root modulo prime p
def find_primitive(p):
    if p == 2:
        return 1
    phi = p - 1
    primes = factorize(phi)
    # test candidates
    for g in range(2, p):
        ok = True
        for q in primes:
            # if g^(phi/q) ≡ 1 => not primitive
            if modexp(g, phi // q, p) == 1:
                ok = False
                break
        if ok:
            return g
    return None

# Baby-Step Giant-Step to solve g^x ≡ a (mod p)
def discrete_log(g, a, p):
    g %= p
    a %= p
    phi = p - 1
    m = isqrt(phi) + 1
    # baby steps
    baby = {}
    cur = 1
    for j in range(m):
        if cur not in baby:
            baby[cur] = j
        cur = cur * g % p
    # factor = g^{-m} = g^(phi - m)
    factor = modexp(g, phi - m, p)
    giant = a
    # giant steps
    for i in range(m):
        if giant in baby:
            return i * m + baby[giant]
        giant = giant * factor % p
    return None

def solve(P, K, A):
    # Step 1: handle A == 0
    if A == 0:
        return [0]
    phi = P - 1
    # Step 2: solvability check
    d = gcd(K, phi)
    if modexp(A, phi // d, P) != 1:
        return []
    # Step 3: primitive root
    g = find_primitive(P)
    # Step 4: discrete log T
    T = discrete_log(g, A, P)
    if T is None or T % d != 0:
        return []
    # Step 5: solve K·y ≡ T (mod phi)
    Kp, phip, Tp = K // d, phi // d, T // d
    invKp = invmod(Kp, phip)
    y0 = (Tp * invKp) % phip
    # Step 6: build solutions
    roots = []
    for i in range(d):
        exp = (y0 + i * phip) % phi
        roots.append(modexp(g, exp, P))
    return sorted(set(roots))

def main():
    data = sys.stdin.read().split()
    P, K, A = map(int, data)
    ans = solve(P, K, A)
    print(len(ans))
    if ans:
        print(*ans)

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