题意：求$a^x \equiv b(\mod p^e)$ 的最小整数解。
$3 < p \leq 50,p^e \leq 10^18,a\nmid p,b \nmid p$，询问不超过 $1000$ 组。

暴力分跑 bsgs 其实就可以了...因为保证了 $(a,p^e) = 1$。
但是这样显然会 T，于是我们就不要往 bsgs 的方面去想了。
我们考虑 $p$ 很小并且 $>2$，这就强烈暗示我们使用原根解离散对数了。
首先关于原根有一点小性质：

1. $2,4,p^e,2p^e$ 都有原根
2. $p$ 和 $p^e$ 有相同的原根。

所以我们就可以通过暴力找出原根了，接下来我们对两边取 $log$。
$$(log_ga)x \equiv log_gb (\mod \varphi(p^e))$$

考虑我们如果能求出 $(log_ga)$ 和 $(log_gb)$ 的话，我们就代入解一个同余方程就好了。我们考虑求 $(log_ga)$：设 $x = log_ga$,将其化为指数形式：
$$g^x \equiv a (\mod p^e)$$
（我们需要注意如果对指数取模就要用欧拉定理加上 $\varphi$，如果是对整个数取模就要用原来的模数）
首先一个直观的想法就是使用 bsgs （那和暴力有什么区别）所以我们考虑如何利用质数比较小的这一特性来计算。
考虑我们如果已经计算出下式:
$$g^x \equiv a(\mod p^i)$$
尝试是否能快速推到 $\mod p^{i+1}$ 的形式。我们考虑用模数上的指数给指数方程编号（例如上式的编号是 $i$），那么我们考虑第 $i$ 个质数方程解出的解满足 $x \equiv x_i (\mod \varphi(p^i))$。而我们现在想要得到新解 $x \equiv x_{i+1} (\mod \varphi(p^{i+1}))$。注意观察这两个解的规律，不难发现 $x_{i+1} \equiv x_i (\mod \varphi(p^i))$ 即 $x_{i+1} = x_i+k\varphi(p^i)$。
我们枚举这个 $k$ 然后检验答案即可。不难发现 $k$ 的数量只有 $O(p)$ 个。
所以我们这样一层层迭代下去就可以了。

<h1>代码</h1>

/*
 * Author: RainAir
 * Time: 2019-07-21 12:39:06
 */
#include <algorithm>
#include <iostream>
#include <cstring>
#include <climits>
#include <cstdlib>
#include <cstdio>
#include <bitset>
#include <vector>
#include <cmath>
#include <ctime>
#include <queue>
#include <stack>
#include <map>
#include <set>

#define fi first
#define se second
#define U unsigned
#define P std::pair
#define Re register
#define LL long long
#define pb push_back
#define MP std::make_pair
#define all(x) x.begin(),x.end()
#define CLR(i,a) memset(i,a,sizeof(i))
#define FOR(i,a,b) for(Re int i = a;i <= b;++i)
#define ROF(i,a,b) for(Re int i = a;i >= b;--i)
#define DEBUG(x) std::cerr << #x << '=' << x << std::endl
#define int LL

inline int qmul(int x,int y,int p){
    LL k = (LL)((1.0Lxy)/(1.0Lp)),t = xy-k*p;
    t -= p;while(t < 0) t += p;
    return t;
}

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

inline int qpow(int a,int n,int ha){
    int res = 1;
    while(n){
        if(n & 1) res = qmul(res,a,ha);
        a = qmul(a,a,ha);
        n >>= 1;
    }
    return res;
}


inline int uqpow(int a,int n){
    int res = 1;
    while(n){
        if(n & 1) res = 1llresa;
        a = 1llaa;
        n >>= 1;
    }
    return res;
}
const int MAXN = 10000+5;
bool vis[MAXN];

inline bool pdg(int x,int p){
    int t = 1;CLR(vis,0);
    FOR(i,0,p-2){
        if(vis[t]) return false;
        vis[t] = true;
        t = qmul(t,x,p);
    }
    return true;
}

int g;

inline int getg(int p){
    FOR(i,2,p-1){
        if(pdg(i,p)) return i;
    }
}

inline int calc(int p,int e,int a){
    if(e == 1){
        FOR(i,0,p-1){
            //DEBUG(qpow(g,i,p));DEBUG(a%(p));
            if(qpow(g,i,p) == a%p) return i;
        }
        return -1;
    }
    int x0 = calc(p,e-1,a);
    if(x0 == -1) return -1;
    int base = uqpow(p,e-2)*(p-1);
    int ha = uqpow(p,e);
    FOR(i,0,p-1){
        int t = x0+i*base;
        if(qpow(g,t,ha) == a%ha) return t;
    }
    return -1;
}

inline void Solve(){
    int a,b,p,e;scanf("%lld%lld%lld%lld",&a,&b,&p,&e);
    g = getg(p);//DEBUG(g);
    int aa = calc(p,e,a);
    int bb = calc(p,e,b); // aa*ax = bb(mod pp)
    int pp = (p-1)*uqpow(p,e-1);int x=0,y=0;// aax-ppy = bb
    int _gcd = exgcd(aa,pp,x,y);
    if(bb%_gcd){
        puts("-1");return;
    }
    x = qmul(x,bb/_gcd,pp);
    pp /= _gcd;
    x = (x%pp+pp)%pp;
    printf("%lldn",x);
}

signed main(){
    int T;scanf("%lld",&T);
    while(T--) Solve();
    return 0;
}