Codeforces 919E 中国剩余定理 + 线性时间求逆元 + 折半搜索

       计算有多少 $n(1\ge n\ge x)$ 满足

$$n\cdot a^n\equiv b(modp)$$

       $a,b,p,x(2\le p\le 10^6+3)(1\le a,b<p)(1\le x\le 10^{12})$， $p$ 保证为素数。

       看到 $p$ 的数据范围很容易联想，要枚举 $p$ 范围内的 $n$ ，再根据 $x$ 范围内的 $n$ 和当前 $n$ 的关系计算贡献。

       我们注意到当 $n$ 大于 $p$ 的时候会出现 $n\cdot a^n\equiv (n\%p)\cdot a^{n\%(p-1)}(modp)$。

       我们令 $x=n\%(p-1),y=n\%p$，则有 $y\cdot a^x\equiv b(modp)$，则有 $a^x\equiv b\cdot inv(y)$，其中 $inv(y)$ 是关于 $p$ 的乘法逆元。

       于是我们线性时间算出 $1$ 到 $min(x,p-1)$ 的逆元（这里应该用快速幂暴力算时间应该也允许）

       记录下所有的 $b\cdot inv(y)$ 和它对应的 $y$。

       然后我们枚举 $x$，找到最小的 $n$ 使得 $n\%(p-1)=x,n\%p=y$，也就是中国剩余定理。

       然后对于一个小于等于 $x$ 的数，它所能产生的贡献即为 $(x-n)/(p\cdot (p-1))+1$。

#include<bits/stdc++.h>
using namespace std;
#define endl '\n'
#define IOS ios::sync_with_stdio(false); cin.tie(0); cout.tie(0)
typedef long long LL;
const int maxn = 1e6 + 5;
LL num[maxn];
unordered_map<int, int> id;

LL a[3], b[3];
int n = 2;

LL qmod(LL a, LL b, LL mod){		
    LL ans = 0;
    while(b > 0){
        if(b & 1) ans = (ans + a) % mod;
        a = (a + a) % mod;
        b >>= 1;
    }
    return ans;
}

LL exgcd(LL a, LL b, LL &x, LL &y){
    if(b == 0){
        x = 1, y = 0;
        return a;
    }else{
        LL gcd = exgcd(b, a % b, x, y);
        LL t = y;
        y = x - (a / b) * y;
        x = t;
        return gcd;
    }
}

LL china(){
    LL m = a[1], ans = b[1], x, y;
    for(int i = 2; i <= n; i++){
        LL c = (b[i] - ans % a[i] + a[i]) % a[i];
        LL gcd = exgcd(m, a[i], x, y);
        if(c % gcd) return -1;		// 无解
        LL ag = a[i] / gcd;
        x = qmod(x, c / gcd, ag);
        ans += x * m;
        m *= ag;
        ans = (ans % m + m) % m;
    }
    return (ans % m + m) % m;
}


int main(){
    IOS;
    int aa, bb, p;
    LL x, pp;
    cin >> aa >> bb >> p >> x;
    a[1] = p;
    a[2] = p - 1;
    pp = (LL)p * (p - 1);
    int limit = min(x, p - 1LL);
    num[1] = 1;
    for(int i = 2; i <= limit; i++)
        num[i] = ((p - p / i) * num[p % i]) % p;
    for(int i = 1; i <= limit; i++)
        id[num[i] * bb % p] = i;
    LL ans = 0, now = aa, minone;
    for(int i = 1; i <= limit; i++){
        if(id.count(now)){
            b[1] = id[now];
            b[2] = i;
            minone = china();
            if(x >= minone) ans += (x - minone) / pp + 1;
        }
        now = now * aa % p;
    }
    cout << ans << endl;
}