[bzoj 3122] [Sdoi2013]随机数生成器：数论，同余，分类讨论，BSGS

题意：给四个参数p, a, b, x1，生成一个数列，以x1为首项，x[i+1] = (ax[i]+b) mod p (i>=1)，问是否存在x[i]=t，如果存在，求出最小的i。0<=a, b, x1, t < p，2<=p<=10^9，p为素数。

以前看过这道题，知道思路。然而AC还是不太容易啊……

如果这是数学试卷上的一道数列大题，如果递推式不mod p，作为高中生的我们会求一求通项，解一个方程。

现在，这是一道OI题，每次mod p——但这并不妨碍我们求通项，因为同余号和等号具有许多一致之处。

求一求特征根，先不急着做除法，得到

(a−1)x≡−b(modp)

即便p是素数，也要考察a−1≡0的情形。在模算术中，0没有逆元，正如0不能作分母。

当a=1，原递推式为

xi+1=(xi+b)modp

容易写出通项

xn=x1+(n−1)bmodp

当a≠1，原递推式化为

xi+1+ba−1≡a(xi+ba−1)(modp)

通项公式

xn=(an−1(x1+ba−1)−ba−1)modp

通项求完了，让我们来解方程。

令

ai=t

且i是满足此方程的最小正整数。

当a=1，

(i−1)b≡t−x1

若b=0，{an}是常数数列。若x1=t，i=1；否则，无解。
若b≠0，i≡(t−x1b+1)，注意i是正整数。

当a=0，
数列首项为x1，其余项为b。若x1=t，i=1；否则，若t=b，i=2；否则，无解。

当a≠1且a≠0，

ai−1(x1+ba−1)≡ba−1+t

若x1+ba−1≡0，返回通项公式看一看，发现{xn}是常数数列。若x1=t，i=1；否则，无解。也可以直接看同余式右边是否为0，这两个条件是等价的。
若x1+ba−1≢0，则ai+1≡ba−1+tx1+ba−1。离散对数问题。我们用大步小步算法求出此方程的最小正整数解，或者说明无解。

关于大步小步算法，请听下回分解~

综上所述，我们要考察a是否为1或0。若a=1，考察b是否为0。若a!=1且a!=0，考察(x1+b*inv(a-1))%p是否为0。

不单独考察a=0也可以，这要求在BSGS过程中特判a^0=1。但是0^0感觉有点怪怪的……

#include <cstdio>#include <cmath>#include <map>using namespace std;typedef long long ll;ll p;inline ll fast_exp(ll x, ll n){    ll z = 1;    for (ll y = x; n; n >>= 1, y = y*y%p)        if (n & 1)            z = z*y%p;    return z;   }inline ll inv(ll x){    return fast_exp(x, p-2);}ll bsgs(ll a, ll y) // minimum x such that a^x mod p = y{    ll m = sqrt(p);    map<ll, ll> M;    for (ll i = 1, t = a*y%p; i <= m; ++i, t = t*a%p)        M[t] = i;    for (ll i = m, b = fast_exp(a, m), t = b; i-m < p-1; i += m, t = t*b % p)        if (M.count(t))            return i-M[t];    return -2;}int main(){    int T;    ll a, b, x1, t;    for (scanf("%d", &T); T--; ) {        scanf("%lld %lld %lld %lld %lld", &p, &a, &b, &x1, &t);        ll ans = -1;        if (a == 0)            ans = x1 == t ? 1 : (b == t ? 2 : -1);        else if (a == 1) {            if (b) {                ans = ((t-x1+p)%p*inv(b)+1)%p;                ans = ans ? ans : p;                            } else                ans = x1 == t ? 1 : -1;        } else {            ll c = b*inv(a-1)%p, d = (x1+c)%p;            if (d)                ans = bsgs(a, (t+c)*inv((x1+c)%p)%p) + 1;            else                ans = (t+c)%p ? -1 : 1;        }        printf("%lld\n", ans);    }    return 0;}