每日小结2 baby-step_giant_step

维基百科：baby-step_giant_step

名字已经暗示了算法的精髓。感觉和桶计数很像。就是把小单位合并。用大单位和小单位一起精确查找。和杭州网络赛的中值滤波很像呀。

an==b(mod c) 给定你n的范围。可能很大。但是由数论的知识知道n<=c.

简单论述一下啊：

{an},n~1,2,3... 可以看成a的生成群。其实也不是严格的。严格的a必须与c互质才行。因为不互质的话，a是没有逆元的。再次对上句话进行证明。（ps：为啥我老是证明呢。因为我刚学。很多东西都不能像加减乘除那样一眼看出来。还有可能是反应慢吧。反正就是喜欢证明。不证明的话用起来没底，也不想用。）我们都知道ax+by=c有解的充分必要条件是gcd(a,b)|c;这个的证明可以化为基本的加减乘除了。可以方程两边同时除以gcd(a,b).这样左边是整数，右边如果gcd（a，b）不能整除c的话就是分数，显然无解。 所以如果a和c不互质的话，a*x+c*y==1无解的。所以a无逆元。没有逆元就不能叫群了。以上可以不用了解。我们目的是证明n<=c.直观来看，a n在0到c-1变化（是显然的）。这样数列an就肯定是个周期数列（抽屉原理）。通过数学归纳法，可以严格证明周期是存在的。且一定小于等于n。

知道了大概n的范围了。但是n有可能在10的10此方左右。一般程序108方是可以忍受的。所以直接枚举肯定不行。但sqrt（1010）就很小了。baby-step_giant_step就是利用了这个复杂度的算法。an=am*i-j,则am*i-j==b(mod c)<==>am*i==b*aj(mod c);i和j的范围都是m，m=sqrt（n）；

在判断是否存在j的时候，用二分，可以在logm时间内查找出。这样的话，总的时间复杂度sqrt（n）*log（sqrt（n））.维基百科说有O（1）的hash。应该是很牛的hash吧。有待研究呀。

写了个测试小程序。输入a，b，c，d（n的范围）; 输出最小的非负n。

#include <iostream>
#include <cmath>
#include <algorithm>
using namespace std;
#define maxn 100000
#define ll long long
struct baby{
ll r;
int m;
bool operator <(const baby &a)const {
if(r==a.r)
return m<a.m;
return r<a.r;
}
void set(ll a,int b){
r=a,m=b;
}
} baby_step[40000];

ll gaint_step;
ll a,b,c;
int find(ll s,int m){
int h=0, t=m, mid;
while(h<=t){
mid=(h+t)/2;
if(baby_step[mid].r<s){
h=mid+1;
}
else if(baby_step[mid].r>s)
t=mid-1;
else
return mid;
}
return -1;
}
int main(){
ll t=1;
int j,ans,d,lim;
while(cin>>a>>b>>c>>d){
lim=int(sqrt(1.0*d));
t=1;
ans=-1;
for(int i=0; i<lim; ++i){
baby_step[i].set((t*b)%c,i);
t*=a;
t%=c;
}
sort(baby_step,baby_step+lim);
gaint_step=t;

for(int i=1; i<=lim+1; ++i){
if((j=find(t,lim-1))!=-1){
ans=i*lim-baby_step[j].m;
break;
}
t*=gaint_step;
t%=c;
}
if(ans==-1)
cout<<"NO"<<endl;
else
cout<<ans<<endl;
}
return 0;
}

////////////////////////////////////////////////////////////////////////////////////更新

http://hi.baidu.com/pojoflcc/blog/item/78bbd24b7d5d732908f7effa.html在这个blog上看到了这句话：

经过AekdyCoin教导 发现这个算法当a和n不互质的时候会死 因为没有逆元 i*m-j不能随便减

上面我讨论过群的概念。但是这一点却没有想到。当a和c有公约数时，我的测试程序会出错的。这个是为什么呢。我来举个例子。

a=4，c=12，b=1；容易知道an==4（mod 12）那么如果42==4，我的程序会认为4==1.那么前面的推导am*i-j==b(mod c)<==>am*i==b*aj(mod c)是错的。am*i-j==b(mod c)==>am*i==b*aj(mod c)是可以的。但am*i-j==b(mod c)<==am*i==b*aj(mod c)是缺乏根据的。这些也可以用基本的加减乘除运算推演。如果抽象一点的话，可以说。ａ和ｃ不互质，则ａ不是Ｚｃ×里的元素，不存在除法运算。

看了上面blog后面的分析，没怎么看懂。。。想了下，可以这样的，就是把a和c的gcd找到，然后除一下。如果b不能整除，那么肯定无解。然后就可以继续按照上文做了。

////////////////////////////////////更新中 修改了前面代码的bug

下面代码在poj上实验可以的。

#include <iostream>
#include <cmath>
#include <algorithm>
using namespace std;
#define maxn 100000
#define ll long long
int lim;
int gcd(int a,int b){
if(!b)
return a;
else
return gcd(b,a%b);
}

struct baby{
ll r;
int m;
bool operator <(const baby &a)const {
if(r==a.r)
return m<a.m;
return r<a.r;
}
void set(ll a,int b){
r=a,m=b;
}
} baby_step[40000];

ll gaint_step;
ll a,b,c;
int find(ll s,int m){
int h=0, t=m, mid;
while(h<=t){
mid=(h+t)/2;
if(baby_step[mid].r<s){
h=mid+1;
}
else if(baby_step[mid].r>s)
t=mid-1;
else
{
int i;
for( i=mid+1; i<lim&&baby_step[i].r==s; ++i);
return i-1;
}
}
return -1;
}
int main(){
ll t=1;
int j,ans,d;
while(cin>>a>>c>>b){
d=c;
if(a==0&&b==0&&c==0)
break;
if(c==1) {if(b==0) ans=0; goto L;}
if(a==0) {if(c==0) ans=0; goto L;}
if(b==1) { ans=0; goto L;}
lim=int(sqrt(1.0*d));
t=1;
ans=-1;
if(b%gcd(a,c)) goto L;
for(int i=0; i<lim; ++i){
baby_step[i].set((t*b)%c,i);
t*=a;
t%=c;
}
sort(baby_step,baby_step+lim);
gaint_step=t;

for(int i=1; i<=lim+1; ++i){
if((j=find(t,lim-1))!=-1){
ans=i*lim-baby_step[j].m;
break;
}
t*=gaint_step;
t%=c;
}
L:    if(ans==-1)
cout<<"No Solution"<<endl;
else
cout<<ans<<endl;
}
return 0;
}