luogu2485 [SDOI2011]计算器 poj3243 Clever Y BSGS算法

BSGS 算法,即 Baby Step,Giant Step 算法、拔山盖世算法。
计算 \(a^x \equiv b \pmod p\)

\(p\)为质数时

特判掉 \(a,p\) 不互质的情况。
由于费马小定理 \(x^{p-1} \equiv 1 \pmod p\)\(p\) 为质数,则要是暴力的话只需要枚举到 \(p-1\) 即可。
假设 \(x=it-j\),其中 \(t= \lceil \sqrt p \rceil,j \in [0,t]\),方程变为 \(a^{it-j} \equiv b \pmod p\),即 \(a^{it} \equiv ba^j \pmod p\)。我们惊喜地发现,左右最多也就 \(t\) 个左右种可能的取值(这就是 \(t\) 为什么取那个值的原因),那我们枚举 \(j\),把 \(ba^j\) 所对应的 \(j\) 都存起来,然后枚举 \(i\) 找有无对应即可。
保存的话用手写 hash 比 map 快很多的。
\(ba^j\) 冲突怎么办?答:存 \(j\) 大的。因为要想 \(x\) 尽量小,就要让 \(j\) 尽量大。
第一个找到的 \(i\) 和它对应的 \(j\) 就是答案。为什么?答:因为 \(i\) 变化 \(1\)\(x\) 变化的幅度是 \(t\)
时间复杂度 \(\mathrm{O}(\sqrt{p})\)

计算器:

#include <iostream>
#include <cstdio>
#include <cmath>
#include <map>
using namespace std;
typedef long long ll;
int T, k;
ll y, z, p;
map<int,int> d;
ll work1(ll a, ll b, ll p){
    ll re=1;
    while(b){
        if(b&1) re = (re * a) % p;
        a = (a * a) % p;
        b >>= 1;
    }
    return re;
}
ll exgcd(ll a, ll b, ll &x, ll &y){
    if(!b){
        x = 1;
        y = 0;
        return a;
    }
    ll re=exgcd(b, a%b, x, y);
    ll qwq=x;
    x = y;
    y = qwq - a / b * y;
    return re;
}
void work2(){
    ll u, v;
    ll gcd=exgcd(y, p, u, v);
    if(z%gcd)   printf("Orz, I cannot find x!\n");
    else
        printf("%lld\n", ((u*z/gcd)%(p/gcd)+(p/gcd))%(p/gcd));
}
void work3(){
    d.clear();
    int m=sqrt(p);
    if(m*m!=p)  m++;
    for(int i=0; i<=m; i++)
        d[int((z*work1(y, i, p))%p)] = i;//当然这里也可以换掉快速幂,因为递推顺便就能求幂了
    if(y%p==0){
        if(z%p==1)  printf("0\n");
        else if(z%p==0) printf("1\n");
        else    printf("Orz, I cannot find x!\n");
        return ;
    }
    for(int i=1; i<=m; i++){
        int tmp=work1(y, i*m, p);
        if(d.find(tmp)!=d.end()){
            printf("%lld\n", (ll)i*m-d[tmp]);
            return ;
        }
    }
    printf("Orz, I cannot find x!\n");
}
int main(){
    cin>>T>>k;
    while(T--){
        scanf("%lld %lld %lld", &y, &z, &p);
        if(k==1)    printf("%lld\n", work1(y, z, p));
        if(k==2)    work2();
        if(k==3)    work3();
    }
    return 0;
}

测试数据:

2 3
39 26 13
6 11 5
1
0

\(p\)无限制时

计算 \(a^x \equiv b \pmod p\),最难过的事情就是 \((a,p) \not = 1\)。那我们就想办法消掉一些 \(p\) 的因子让 \(a,p\) 互素。记 \(d=(a,p)\)
倘若 \(d \nmid b\)(不整除的语法是 \nmid),则 \(b=1\) 时方程有解,解为 \(x=0\)
倘若 \(d \mid b\)\(d=1\),则就是普通的 bsgs。
否则,
\[a^x \equiv b \pmod p \Rightarrow \frac{a}{d} \times a^{x-1} \equiv \frac{b}{d} \pmod {\frac{p}{d}}\]
当然,还有更难过的是 \((a,p/d)\) 可能还不为 \(1\)……那就继续搞下去,直到
\[\frac{a^k}{\prod_{i=1}^kd_i} \times a^{x-k} \equiv \frac{b}{\prod_{i=1}^kd_i} \pmod {\frac{p}{\prod_{i=1}^kd_i}}\]
\((a,\frac{p}{\prod_{i=1}^kd_i})=1\)
当然还有更更难过的是还没到互素的时候就有 \(\frac{p}{\prod_{i=1}^\kappa d_i} \nmid b\)……那唯一的可能就是 \(x=0\) 了。
然后发现如果真正的解 \(x<k\) 的话比较难过,那我们在累积公约数的时候可以顺便判断 \(x \in [0,k)\) 时是否是解。
对于 \(x \geq k\) 的情况,我们为了赏心悦目,把方程换元为
\[a' \times a^{x'} \equiv b' \pmod {p'}\]
其中 \((a',p')=1\)
这样就是一个普通的 bsgs。最终的解是 \(x=x'+k\)

Clever Y:

#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
using namespace std;
typedef long long ll;
int a, b, p;
const int mod=32767;//这不是质数,不好
int gcd(int a, int b){
    return !b?a:gcd(b,a%b);
}
int ksm(int a, int b, int p){
    int re=1;
    while(b){
        if(b&1) re = ((ll)re * a) % p;
        a = ((ll)a * a) % p;
        b >>= 1;
    }
    return re;
}
struct HASHTABLE{
    int nxt[200005], dai[200005], ret[200005], hea[200005], hmn;
    void clear(){
        memset(nxt, 0, sizeof(nxt));
        memset(hea, 0, sizeof(hea));
        hmn = 0;
    }
    void insert(int x, int y){
        int tmp=x%mod;
        for(int i=hea[tmp]; i; i=nxt[i])
            if(dai[i]==x){
                ret[i] = y;
                return ;
            }
        nxt[++hmn] = hea[tmp];
        dai[hmn] = x; ret[hmn] = y;
        hea[tmp] = hmn;
    }
    int query(int x){
        int tmp=x%mod;
        for(int i=hea[tmp]; i; i=nxt[i])
            if(dai[i]==x)
                return ret[i];
        return -1;
    }
}hashTable;
int exbsgs(int a, int b, int p){
    hashTable.clear();
    a %= p; b %= p;
    if(b==1)    return 0;
    int cnt=0, d=1, g;
    do{
        g = gcd(a, p);
        if(b%g) return -1;
        p /= g; b /= g;
        d = ((ll)d * a / g) % p;
        cnt++;
        if(b==d)    return cnt;
    }while(g!=1);
    int m=ceil(sqrt(p)), t=b;
    for(int i=0; i<=m; i++){
        hashTable.insert(t, i);
        t = ((ll)t * a) % p;
    } 
    g = ksm(a, m, p);
    t = ((ll)d * g) % p;
    for(int i=1; i<=m; i++){
        int re=hashTable.query(t);
        if(re>=0)   return i*m-re+cnt;
        t = ((ll)t * g) % p;
    }
    return -1;
}
int main(){
    while(scanf("%d %d %d", &a, &p, &b)!=EOF){
        if(a+b+p==0)    break;
        int re=exbsgs(a, b, p);
        if(re==-1)  printf("No Solution\n");
        else    printf("%d\n", re);
    }
    return 0;
}

转载于:https://www.cnblogs.com/poorpool/p/8511570.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值