[BZOJ3122]-[Sdoi2013]随机数生成器-BSGS+exgcd

说在前面

me发现me连等比数列和这种东西都不会算了
是脱离文化课太久，还是之前就掌握不扎实啊=A=

---

题目

BZOJ3122传送门

题目大意

给出质数 p $p$，常数 a,b,X1,t $a,b,X_1,t$。求一个最小的 n $n$满足 Xn=t $X_n=t$。其中，关于 X $X$的递推式为 Xi=(aXi−1+b)%p $X_i=(a{X}_{i-1}+b)\mathrm{\%}p$
数据范围： 0≤a,b,X1,t≤p−1 $0\le a,b,X_1,t\le p-1$， p≤1e9 $p\le 1e9$

输入输出格式

输入格式：
第一行一个整数T，表示数据组数
接下来每行五个数字，p，a，b，X1，t，表示一组数据

输出格式：
对于每组数据，如果n存在则输出n，否则输出-1

---

解法

直接把递推式拆开，用等比数列求和公式，得到通项公式：
（下面为了偷懒，直接把等式换成同余式了）
Xn≡X1an−1+b(an−1−1)a−1(modp) $X_n\equiv X_1{a}^{n-1}+\frac{b(a^{n-1}-1)}{a-1}(\mathrm{mod}p)$
把有公因式 an−1 $a^{n-1}$的提出来，把 Xn $X_n$换成 t $t$，得到
t≡an−1(X1+ba−1)−ba−1(modp) $t\equiv a^{n-1}(X_1+\frac{b}{a-1})-\frac{b}{a-1}(\mathrm{mod}p)$
继续移项，得到
(t+ba−1)/(X1+ba−1)≡an−1(modp) $(t+\frac{b}{a-1})/(X_1+\frac{b}{a-1})\equiv a^{n-1}(\mathrm{mod}p)$
然后到这里就用BSGS解决了。注意特判a=1，a=0，b=0等各种特殊情况即可

貌似对这种题还有一种通法：
设 Xi+1+c=a(Xi+c) $X_{i+1}+c=a(X_i+c)$
然后把 Xi+1 $X_{i+1}$用 aXi+b $a{X}_i+b$换掉，可以解出 c $c$，然后可以得到和上面一样的式子
这样做的好处在哪里呢？移项得到： Xi+1=a(Xi+c)−c $X_{i+1}=a(X_i+c)-c$，于是就可以把 Xi $X_i$用带有 −c $-c$的式子替换掉，进而和括号里原本的 +c $+c$抵消，即可瞬间得出通项公式

---

下面是自带小常数的代码

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std ;

int T ;
long long mmod , a , b , X1 , t ;

long long s_pow( long long x , int b ){
    long long rt = 1 ;
    while( b ){
        if( b&1 ) rt = rt * x %mmod ;
        x = x * x %mmod ; b >>= 1 ;
    } return rt ;
}
long long inv( long long x ){ return s_pow( x , mmod - 2 ) ; }

long long exgcd( long long a , long long b , long long &x , long long &y ){
    if( !b ){
        x = 1 , y = 0 ;
        return a ;
    } else {
        long long xx , yy , tmp = exgcd( b , a%b , xx , yy ) ;
        x = yy % mmod , y = ( xx - a / b * yy ) %mmod ;
        return tmp ;
    }
}

struct hashTable{
    int head[10007] , pre[40000] , num[40000] , ex[40000] , tp ;
    void clear(){
        tp = 0 ;
        memset( head , 0 , sizeof( head ) ) ;
    }
    void Insert( int x , int k ){
        int id = x%10007 ;
        for( int i = head[id] ; i ; i = pre[i] )
            if( num[i] == x ) return ;
        pre[++tp] = head[id] ;
        head[id] = tp ;
        num[tp] = x , ex[tp] = k ;
    }
    int Query( int x ){
        int id = x%10007 ;
        for( int i = head[id] ; i ; i = pre[i] )
            if( num[i] == x ) return ex[i] ;
        return -1 ;
    }
} Hs ;

int solve(){
    Hs.clear() ;
    long long sq = sqrt( mmod - 1 ) + 1 , base = inv( s_pow( a , sq ) ) , now = 1 ;
    long long tmp = b * inv(a-1)%mmod , aim = ( t + tmp )%mmod * inv( X1 + tmp )%mmod ;
    for( int i = 0 ; i <= sq ; i ++ ){
        Hs.Insert( now , i ) ;
        if( now == aim ) return i + 1 ;
        else now = now * a %mmod ;
    }
    for( int i = 0 ; i <= sq ; i ++ ){
        int tmp = Hs.Query( aim ) ;
        if( tmp != -1 ) return i * sq + tmp + 1 ;
        else aim = aim * base %mmod ;
    } return -1 ;
}

int main(){
    scanf( "%d" , &T ) ;
    while( T -- ){
        scanf( "%lld%lld%lld%lld%lld" , &mmod , &a , &b , &X1 , &t ) ;
        if( X1 == t ){
            puts( "1" ) ;
        } else if( a == 0 ){
            if( t == b ) puts( "2" ) ;
            else puts( "-1" ) ;
        } else if( a == 1 ){
            if( b == 0 ) puts( "-1" ) ;
            else{
                long long x , y , C = t - X1 , d = exgcd( mmod , b , x , y ) ;
                // gcd( mmod , b ) == 1 
                y = ( y * C%mmod + mmod )%mmod ;
                printf( "%lld\n" , y + 1 ) ;
            }
        } else printf( "%d\n" , solve() ) ;
    }
}