解高次同余方程 (A^x=B(mod C),0<=x


先給出我所參考的兩個鏈接:

http://hi.baidu.com/aekdycoin/item/236937318413c680c2cf29d4 (AC神,數論帝  擴展Baby Step Giant Step解決離散對數問題)

http://blog.csdn.net/a601025382s/article/details/11747747

Baby Step Giant Step算法:復雜度O( sqrt(C) )

我是綜合上面兩個博客,才差不多懂得了該算法。 

 

先給出AC神的方法:

 

原創帖!轉載請注明作者 AekdyCoin !

【普通Baby Step Giant Step】

【問題模型】
求解
A^x = B (mod C) 中 0 <= x < C 的解,C 為素數

【思路】
我們可以做一個等價
x = i * m + j ( 0 <= i < m, 0 <=j < m) m = Ceil ( sqrt( C) )
而這么分解的目的無非是為了轉化為:
(A^i)^m * A^j = B ( mod C)

之后做少許暴力的工作就可以解決問題:
(1) for i = 0 -> m, 插入Hash (i, A^i mod C)
(2) 枚舉 i ,對於每一個枚舉到的i,令 AA = (A^m)^i mod C
我們有
AA * A^j = B (mod C)
顯然AA,B,C均已知,而由於C為素數,那么(AA,C)無條件為1
於是對於這個模方程解的個數唯一(可以利用擴展歐幾里得或 歐拉定理來求解)
那么對於得到的唯一解X,在Hash表中尋找,如果找到,則返回 i * m + j
注意:由於i從小到大的枚舉,而Hash表中存在的j必然是對於某個剩余系內的元素X 是最小的(就是指標)
所以顯然此時就可以得到最小解

如果需要得到 x > 0的解,那么只需要在上面的步驟中判斷 當 i * m + j > 0 的時候才返回


到目前為止,以上的算法都不存在爭議,大家實現的代碼均相差不大。可見當C為素數的時候,此類離散對數的問題可以變得十分容易實現。


【擴展Baby Step Giant Step】

【問題模型】
求解
A^x = B (mod C) 中 0 <= x < C 的解,C 無限制(當然大小有限制……)

【寫在前面】
這個問題比較麻煩,目前網絡上流傳許多版本的做法,不過大部分已近被證明是完全錯誤的!

這里就不再累述這些做法,下面是我的做法(有問題歡迎提出)

下面先給出算法框架,稍后給出詳細證明:

(0) for i = 0 -> 50 if(A^i mod C == B) return i    O(50)
(1)  d=0                D=1 mod C
while((tmp=gcd(A,C))!=1)
{
if(B%tmp)return -1; // 無解!
++d;
C/=tmp;
B/=tmp;
D=D*A/tmp%C;        //不理解處1
}
(2) m = Ceil ( sqrt(C) ) //Ceil是必要的     O(1)
(3) for i = 0 -> m 插入Hash表(i, A^i mod C)  O( m)
(4) K=pow_mod(A,m,C)
for i = 0 -> m
解 D * X = B (mod C) 的唯一解 (如果存在解,必然唯一!)
之后Hash表中查詢,若查到(假設是 j),則 return i * m + j + d    //不理解處2
否則
D=D*K%C,繼續循環
(5) 無條件返回 -1 ;//無解!

 

上面有兩處我一開始不理解的地方,稍后我會解釋,先把AC神的方法給copy完

 

下面是證明:
推論1:
A^x = B(mod C)
等價為
A^a * A^b = B ( mod C) (a+b) == x a,b >= 0

證明:
A^x = K * C + B (模的定義)
A^a * A^b = K*C + B( a,b >=0, a + b == x)
所以有
A^a * A^b = B(mod C)
推論 2:
令 AA * A^b = B(mod C)
那么解存在的必要條件為: 可以得到至少一個可行解 A^b = X (mod C)
使上式成立
推論3
AA * A^b = B(mod C)
中解的個數為 (AA,C)
由推論3 不難想到對原始Baby Step Giant Step的改進
For I = 0 -> m
For any solution that AA * X = B (mod C)
If find X
Return I * m + j+d(這里原本少了d
而根據推論3,以上算法的復雜度實際在 (AA,C)很大的時候會退化到幾乎O(C)
歸結原因,是因為(AA,C)過大,而就是(A,C)過大
於是我們需要找到一中做法,可以將(A,C)減少,並不影響解

下面介紹一種“消因子”的做法

一開始D = 1 mod C
進行若干論的消因子,對於每次消因子
令 G = (A,C[i]) // C[i]表示經過i輪消因子以后的C的值
如果不存在 G | B[i] //B[i]表示經過i輪消因子以后的B的值
直接返回無解
否則
B[i+1] = B[i] / G
C[i+1] = C[i] / G
D = D * A / G

具體實現只需要用若干變量,細節參考代碼

假設我們消了a'輪(假設最后得到的B,C分別為B',C')
那么有
D * A^b = B' (mod C')

於是可以得到算法

for i = 0 -> m
解 ( D* (A^m) ^i ) * X = B'(mod C')
由於 ( D* (A^m) ^i , C') = 1 (想想為什么?)
於是我們可以得到唯一解
之后的做法就是對於這個唯一解在Hash中查找

這樣我們可以得到b的值,那么最小解就是a' + b !!

現在問題大約已近解決了,可是細心看來,其實還是有BUG的,那就是
對於
A^x = B(mod C)
如果x的最小解< a',那么會出錯
而考慮到每次消因子最小消 2
故a'最大值為log(C)
於是我們可以暴力枚舉0->log(C)的解,若得到了一個解,直接返回
否則必然有 解x > log(C)

PS.以上算法基於Hash 表,如果使用map等平衡樹維護,那么復雜度會更大

 

接下來我將講解下面兩處我一開始不理解的地方:

(0) for i = 0 -> 50 if(A^i mod C == B) return i    O(50)
(1)  d=0                D=1 mod C
while((tmp=gcd(A,C))!=1)
{
if(B%tmp)return -1; // 無解!
++d;
C/=tmp;
B/=tmp;
D=D*A/tmp%C;        //不理解處1
}
(2) m = Ceil ( sqrt(C) ) //Ceil是必要的     O(1)
(3) for i = 0 -> m 插入Hash表(i, A^i mod C)  O( m)
(4) K=pow_mod(A,m,C)
for i = 0 -> m
解 D * X = B (mod C) 的唯一解 (如果存在解,必然唯一!)
之后Hash表中查詢,若查到(假設是 j),則 return i * m + j + d    //不理解處2
否則
D=D*K%C,繼續循環
(5) 無條件返回 -1 ;//無解!

 

后來看了上面給的第二個博客中的解釋(被作者附在代碼的后面),瞬間明白了啊!!!

下面給出那個博客里寫的內容,做些適當改動,看着舒服點:


求解模方程A^x=B(mod C),C不為素數。拓展Baby Step Giant Step
(該作者的方法一開始沒有暴力先枚舉答案,即沒有AC神的第(0)步,直接從第(1)步開始了。)
初始D=1,d=0,i=0;
1.令tmp=gcd(A,C),若tmp==1則執行下一步。否則由於A^x=k*C+B;(k為某一整數),則(A/tmp)*A^(x-1)=k*(C/tmp)+B/tmp,(B/tmp為整除,若不成立則無解)
令C=C/tmp,D=D*A/tmp,B=B/tmp,d++。則D*A^(x-d)=B(mod C),接着重復1步驟。
看到沒,while循環的原因瞬間明白了不?d的含義瞬間懂了有木有!!!也就是后面求DX=B(mod C)的解其實是x-d!!!所以最后才要i*m+j+d,加上個d!!!
2.通過1步驟后,保證了A和D都與C互質,方程轉換為D*A^(x-d)=B(mod C)。由於A和C互質,所以由歐拉定理A^phi(C)==1(mod C),(A,C互質)
可知,phi(C)<=C,A^0==1(mod C),所以構成循環,且循環節不大於C。從而推出如果存在解,則必定1<=x<C。(在此基礎上我們就可以用
Baby Step Giant Step方法了)
3.令m=ceil(sqrt(C)),則m*m>=C。用哈希表存儲A^0,A^1,..,A^(m-1),接着判斷1~m*m-1中是否存在解。
4.為了減少時間,所以用哈希表縮減復雜度。分成m次遍歷,每次遍歷A^m長度。由於A和D都與C互質,所以gcd(D,C)=1。
接下來實際求DX=B(mod C)的解。用拓展的歐幾里德定理求得DX=1(mod C),即D*x+C*y=gcd(D,C)=1的一組整數解(x,y)。
則D*x+C*y=1-->D*x%C=(D*x+C*y)%C=1-->D*(x*B)%C=((D*x)%C*B%C)%C=B。
所以最終DX=B(mod C)的解X=x*B。若X在哈希表中存在,值為j,則D*A^j=B(mod C),最后我們要求解的答案就是ans=i*m+j+d。
如果不存在,則令D=D*A^m%C,i++后遍歷下一個A^m,直到遍歷(A^m)^0到(A^m)^(m-1)還未找到,則說明不解並退出。

本人注明:該博客中i循環是0~m-1,在AC神的方法中是0~m,本人覺得應該前者比較好。因為解x是1<=x<C的,而m=ceil(sqrt(C)),如果i=m的話,那m*m>C,顯然不正確。

 

下面是AC神的兩道題的代碼:

HDU 2815

#include<iostream>
#include<map>
#include<cmath>
using namespace std;
typedef long long LL;
int gcd(int a,int b){return b?gcd(b,a%b):a;}
int ext_gcd(int a,int b,int& x,int& y){
int t,ret;
if (!b){x=1,y=0;return a;}
ret=ext_gcd(b,a%b,x,y);
t=x,x=y,y=t-a/b*y;
return ret;
}
int Inval(int a,int b,int n){
int x,y,e;
ext_gcd(a,n,x,y);
e=(LL)x*b%n;
return e<0?e+n:e;
}
int pow_mod(LL a,int b,int c){LL ret=1%c;a%=c;while(b){if(b&1)ret=ret*a%c;a=a*a%c;b>>=1;}return ret;}
int BabyStep(int A,int B,int C){
map<int,int> Hash;
LL buf=1%C,D=buf,K;
int i,d=0,tmp;
for(i=0;i<=100;buf=buf*A%C,++i)if(buf==B)return i;
while((tmp=gcd(A,C))!=1)
{
if(B%tmp)return -1;
++d;
C/=tmp;
B/=tmp;
D=D*A/tmp%C;
}
Hash.clear();
int M=(int)ceil(sqrt((double)C));
for(buf=1%C,i=0;i<=M;buf=buf*A%C,++i)if(Hash.find((int)buf)==Hash.end())Hash[(int)buf]=i;
for(i=0,K=pow_mod((LL)A,M,C);i<=M;D=D*K%C,++i)
{
tmp=Inval((int)D,B,C);
if(tmp>=0&&Hash.find(tmp)!=Hash.end())return i*M+Hash[tmp]+d;
}
return -1;
}
int main()
{
int A,B,C;
while(scanf("%d%d%d",&A,&C,&B)!=EOF)
{
if(B>=C){puts("Orz,I can’t find D!");continue;}
int tmp=BabyStep(A,B,C);
if(tmp<0)puts("Orz,I can’t find D!");else printf("%d\n",tmp);
}
return 0;
}
View Code

 

POJ 3243 Clever Y

#include<iostream>
#include<map>
#include<cmath>
using namespace std;
typedef long long LL;
const int maxn = 65535;
struct hash{
int a,b,next;
}Hash[maxn << 1];
int flg[maxn + 66];
int top,idx;
void ins(int a,int b){
int k = b & maxn;
if(flg[k] != idx){
flg[k] = idx;
Hash[k].next = -1;
Hash[k].a = a;
Hash[k].b = b;
return ;
}
while(Hash[k].next != -1){
if(Hash[k].b == b) return ;
k = Hash[k].next;
}
Hash[k].next = ++ top;
Hash[top].next = -1;
Hash[top].a = a;
Hash[top].b = b;
}
int find(int b){
int k = b & maxn;
if(flg[k] != idx) return -1;
while(k != -1){
if(Hash[k].b == b) return Hash[k].a;
k = Hash[k].next;
}
return -1;
}
int gcd(int a,int b){return b?gcd(b,a%b):a;}
int ext_gcd(int a,int b,int& x,int& y){
int t,ret;
if (!b){x=1,y=0;return a;}
ret=ext_gcd(b,a%b,x,y);
t=x,x=y,y=t-a/b*y;
return ret;
}
int Inval(int a,int b,int n){
int x,y,e;
ext_gcd(a,n,x,y);
e=(LL)x*b%n;
return e<0?e+n:e;
}
int pow_mod(LL a,int b,int c){LL ret=1%c;a%=c;while(b){if(b&1)ret=ret*a%c;a=a*a%c;b>>=1;}return ret;}
int BabyStep(int A,int B,int C){
top = maxn; ++ idx; 
LL buf=1%C,D=buf,K;
int i,d=0,tmp;
for(i=0;i<=100;buf=buf*A%C,++i)if(buf==B)return i;
while((tmp=gcd(A,C))!=1){
if(B%tmp)return -1;
++d;
C/=tmp;
B/=tmp;
D=D*A/tmp%C;
}
int M=(int)ceil(sqrt((double)C));
for(buf=1%C,i=0;i<=M;buf=buf*A%C,++i)ins(i,buf);
for(i=0,K=pow_mod((LL)A,M,C);i<=M;D=D*K%C,++i){
tmp=Inval((int)D,B,C);int w ;
if(tmp>=0&&(w = find(tmp)) != -1)return i*M+w+d;
}
return -1;
}
int main(){
int A,B,C;
while(scanf("%d%d%d",&A,&C,&B)!=EOF,A || B || C){
B %= C;
int tmp=BabyStep(A,B,C);
if(tmp<0)puts("No Solution");else printf("%d\n",tmp);
}
return 0;
}
View Code

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM