高次不定方程BSGS算法

学习数学真是一件赛艇的事.

BSGS名字听起来非常有意思,力拔山兮气盖世,北上广深,小步大步…算法其实更有意思,它是用来求解一个方程的

AxBmodP A x ≡ B m o d P

是不是特别眼熟,有几个式子长的特别像,先观察一下:

一:快速幂: 求 ABmodP A B m o d P 的值

二:乘法逆元    Ax1(modP) A ∗ x ≡ 1 ( m o d P )

或者 AxB(modP) A ∗ x ≡ B ( m o d P )

三:欧拉定理 Aφ(P)1(modP) A φ ( P ) ≡ 1 ( m o d P ) (A,P互质)

四:费马小定理 AP11(modP) A P − 1 ≡ 1 ( m o d P ) (P是质数)

先说下这四者关系:快速幂可以快速求后三个,费马小定理是欧拉定理的特殊情况,逆元可以通过费马小定理和快速幂解

可如果有这样的一个式子:

AxB(modP) A x ≡ B ( m o d P ) 我们先假设A,P互质

好像和这四个式子都很像,所以呢?所以呢?

当年我们证明费马小定理的时候发现这个x的范围是在 [0,p1] [ 0 , p − 1 ] 之间

那么我们就可以枚举x从0到p-1,复杂度为O(P);

应该能拿上30分

接下来就是一波骚操作:我们令 m=P m = ⌈ P ⌉ ,然后就可以设 x=im+j x = i ∗ m + j ,其中 i=xm,j=x i = ⌊ x m ⌋ , j = x ,把x代入原来的式子可以得到 Aim+jB(modP) A i ∗ m + j ≡ B ( m o d P ) 两边乘上一个 Aim A − i ∗ m ,就可以得到 AjBAim(modP) A j ≡ B ∗ A − i ∗ m ( m o d P )

所以呢?

所以就可以求了啊

我们只要枚举左边的 j j ,把左边的答案和j存起来 left l e f t _ ans[j]=Aj a n s [ j ] = A j (可是存不下怎么办,哈希蛤一下就存下了)然后再枚举右边的 i,计算右边的值,看看我们右边的值是否在数组里出现过,如果出现过那么我们通过i和j找到的 i*m+j 就是一个答案了

然后就会发现复杂度被我们开了一个方

冷静分析:

这个算法先枚举j需要 P P 的时间,再枚举i需要 P P 的时间,不过枚举i是要算下逆元需要 log2P l o g 2 P 的时间,看起来复杂度=O( P P + P P log2(P))=O( P P log2P l o g 2 P ), 不过我们再看看右边的式子:

BAimmodP=B(Am)imodP B ∗ A − i ∗ m m o d P = B ∗ ( A − m ) i m o d P

然后我们就得到右边的递推式,只要先求出 AmmodP A − m m o d P 就可以O(1)计算右边的式子了,其中 AmAP1mmodP A − m ≡ A P − 1 − m m o d P ,因为费马小定理…所以复杂度被我们降到了O( P P + log2P l o g 2 P + P P )=O( P P )

灼热分析:

算法的思想其实就是分块,把x分成 P P * P P 的块,会到设x的式子,x=i*m+j,我们先Baby_Step枚举小的j,再Giant_Step枚举大的i,名字听起来很形象哈哈哈哈哈哈.因为先枚举小的,后枚举大的,所以当出现i满足条件时,可以保证此时答案是最小的正整数解,这时直接return i*m+j

科学分析:

哈希好用呐~之前懒得用哈希,总觉得用STL的map能省很多事,然后就很尴尬的调了两天…一直TLE,最后绝望的手写了哈希表,然后居然就p+的A掉了,千万别用map,千万别用map,千万别用map,STL里面的玄学操作看起来很好用,我们最好还是乖乖学一学正常操作,老老实实手写哈希….

然后看看代码

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
#define ll long long
const int mod=1048573;
int hcnt=0,head[mod+10];
struct Haha{
    int val,id,next;
}hash[mod+10];
void insert(int x,int pos){
    int k=x%mod;
    hash[++hcnt].val=x;
    hash[hcnt].id=pos;
    hash[hcnt].next=head[k];
    head[k]=hcnt;
}
int find(int x){
    int k=x%mod;
    for(int i=head[k];i;i=hash[i].next){
        if(hash[i].val==x) return hash[i].id;
    }
    return -1;
}
ll ksm(int a,int b,int p){
    int x=a;
    ll ret=1;
    if(b<0) return -1;
    while(b){
        if(b&1) ret=1ll*(ret*x)%p;
        b>>=1;
        x=1ll*x*x%p;
    }
    return ret;
}
int BSGS(int a,int b,int p){
    int m=(int)(sqrt(p)+0.999999);
    if(b==1) return 0;
    if(a==b) return 1;
    if(!b){
        if(!a) return 1;
        return -1;
    }
    ll x=1;
    for(int i=1;i<=m;++i){
        x=x*a%p;
        insert(x,i);
    }
    ll inv=1;
    int inv2=ksm(a,p-m-1,p)%p; 
    for(int i=0;i<m;++i){
        int k=i*m;
        if(inv==-1) return -1;
        int ans=inv*b%p;
        int jgy=find(ans);
        if(~jgy){
            return k+jgy;
        }
        inv=1ll*inv*inv2%p;
    }
    return -1;
}
void init(){
    for(int i=0;i<mod;++i){
        hash[i].val=-1;
        hash[i].next=0;
        hash[i].id=0;
    }
    memset(head,0,sizeof(head));
    hcnt=0;
}
int main(){
    int a,b,p;
    while(~scanf("%d%d%d",&p,&a,&b)){
        init();
        int dove=BSGS(a,b,p);
        if(~dove) printf("%d\n",dove);
        else printf("no solution\n");
    }
    return 0;
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值