#include <bits/stdc++.h>

using i64 = long long;

constexpr int N = 1E5;

bool isprime[N];
std::vector<int> primes;

i64 mod(i64 a, i64 p) {
    if (a < 2 * p) {
        return a;
    }
    i64 v = (a - p) / p;
    a -= v * p;
    return a;
}

i64 mulmod(i64 a, i64 b, i64 p) {
    if (b == 0 || a <= 2 * p / b) {
        return mod(a * b, p);
    }
    i64 res = a * b - i64(1.0L * a * b / p) * p;
    res %= p;
    if (res < 0) {
        res += p;
    }
    res += p;
    return res;
}

i64 power(i64 a, i64 b, i64 p) {
    i64 res = 1 % p;
    for (; b; b /= 2, a = mulmod(a, a, p)) {
        if (b % 2) {
            res = mulmod(res, a, p);
        }
    }
    return res;
}

using Matrix = std::array<std::array<i64, 2>, 2>;

Matrix mul(const Matrix &a, const Matrix &b, i64 p) {
    Matrix c{};
    for (int i = 0; i < 2; i++) {
        for (int j = 0; j < 2; j++) {
            for (int k = 0; k < 2; k++) {
                c[i][j] = mod(c[i][j] + mulmod(a[i][k], b[k][j], p), p);
            }
        }
    }
    return c;
}

i64 fib(i64 n, i64 p, i64 f0 = 0, i64 f1 = 1, i64 x = 1, i64 y = 1) {
    Matrix m{x, 1, y, 0};
    Matrix a{f1, f0, 0, 0};
    for (; n; n /= 2, m = mul(m, m, p)) {
        if (n % 2) {
            a = mul(a, m, p);
        }
    }
    return a[0][1];
}

i64 exgcd(i64 a, i64 b, i64 &x, i64 &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    i64 g = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return g;
}

i64 inverse(i64 a, i64 m) {
    i64 x, y;
    exgcd(a, m, x, y);
    x %= m;
    if (x < 0) {
        x += m;
    }
    return x;
}

int get(int x, int y, int g, int a, int p, int e) {
    i64 m = 1;
    for (int i = 0; i < e; i++) {
        m *= p;
    }
    i64 phim = m / p * (p - 1);
    a %= p;
    
    i64 Fx = fib(x, phim);
    i64 Fy = fib(y, phim);
    i64 Fg = fib(g, phim);
    
    i64 ax = (power(a, Fx, m) - 1) % m;
    i64 ay = (power(a, Fy, m) - 1) % m;
    i64 ag = (power(a, Fg, m) - 1) % m;
    if (ag != 0) {
        i64 m1 = 1;
        while (ag % p == 0) {
            ag /= p;
            m1 *= p;
        }
        Fx = fib(x, phim * m1);
        Fy = fib(y, phim * m1);
        Fg = fib(g, phim * m1);
        ax = (power(a, Fx, m * m1) - 1) % (m * m1) / m1;
        ay = (power(a, Fy, m * m1) - 1) % (m * m1) / m1;
        ag = (power(a, Fg, m * m1) - 1) % (m * m1) / m1;
        i64 inv = inverse(ag, m);
        return ax * ay % m * inv % m * inv % m;
    }
    i64 t = fib(g, m, m - 2, 1);
    i64 v = g % 2 ? 1 : m - 1;
    // std::cerr << "x : " << x << ", y : " << y << ", a : " << a << ", p : " << p << ", e : " << e << "\n";
    // std::cerr << Fx << " " << Fy << " " << Fg << "\n";
    // std::cerr << ax << " " << ay << " " << ag << "\n";
    // std::cerr << "t : " << t << "\n";
    return fib(x / g, m, 0, 1, t, v) * fib(y / g, m, 0, 1, t, v) % m;
}

void solve() {
    int x, y, a, m;
    std::cin >> x >> y >> a >> m;
    
    std::vector<std::pair<int, int>> fac;
    int tmp = m;
    for (auto p : primes) {
        if (p * p > m) {
            break;
        }
        if (tmp % p == 0) {
            int e = 0;
            while (tmp % p == 0) {
                tmp /= p;
                e++;
            }
            fac.emplace_back(p, e);
        }
    }
    if (tmp > 1) {
        fac.emplace_back(tmp, 1);
    }
    
    int n = fac.size();
    
    std::vector<int> v(n);
    int g = std::gcd(x, y);
    for (int i = 0; i < n; i++) {
        auto [p, e] = fac[i];
        v[i] = get(x, y, g, a, p, e);
    }
    int ans = 0;
    for (int i = 0; i < n; i++) {
        auto [p, e] = fac[i];
        int pe = 1;
        for (int j = 0; j < e; j++) {
            pe *= p;
        }
        ans = (ans + 1LL * (m / pe) * inverse(m / pe, pe) % m * v[i]) % m;
    }
    std::cout << ans << "\n";
}

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    
    std::fill(isprime + 2, isprime + N, true);
    for (int i = 2; i < N; i++) {
        if (isprime[i]) {
            primes.push_back(i);
        }
        for (auto p : primes) {
            if (i * p >= N) {
                break;
            }
            isprime[i * p] = false;
            if (i % p == 0) {
                break;
            }
        }
    }
    
    int t;
    std::cin >> t;
    
    while (t--) {
        solve();
    }
    
    return 0;
}