离散对数&&大步小步算法及扩展

2018/05/29 22:13
阅读数 14

bsgs algorithm

ax≡b(mod n)

大步小步算法,这个算法有一定的局限性,只有当gcd(a,m)=1

 

原理

此处讨论n为素数的时候。

ax≡b(mod n)(n为素数)

由费马小定理可知,只需要验证0,1,2...n-1是不是解即可,因为an-1 = 1mod(n)

算法过程

1、首先求出a0,a1,a2,...,am-1 模上n的值是否为b,存储在e[i]中,求出am的逆a-m

2、下面考虑am,am+1,...,a2m-1 模上n的值是否为b

此时不用一一检查,如果当中有解,相当于存在e[i],使得e[i] * am = b mod(n)

两边乘上a-m,e[i] = b * a-m mod(n),只需要检查存不存在这样的e[i]即可

3、同理,可以递推检查出a2m - a3m-1中解的情况

为了方便,把e[i]存储在map<int, int>x中,x[j]表示满足ei =j 的最小下标(因为可能有多个值相同)

 

 1 map<ll, int>x;
 2 ll log_mod(ll a, ll b, ll n)//n为素数
 3 {
 4     //if(b >= n)return -1;
 5     a %= n;b %= n;//注意题目,如果b >= n是不存在解的
 6     if(a == 0)
 7     {
 8         if(b == 0)return 1;if(b == 1)return 0;
 9         else return -1;
10     }
11     ll m = ceil(sqrt(n + 0.5));
12     ll v = pow(a, n - 1 - m, n);
13     x.clear();
14     x[1] = m;
15     ll e = 1;
16     for(int i = 1; i < m; i++)//计算a的i次方mod n,并存下来
17     {
18         e = e * a % n;
19         if(!x[e])x[e] = i;
20     }
21     for(int i = 0; i < m; i++)//计算a^(im), a^(im+1),...,a^(im+m-1)
22     {
23         int num = x[b];
24         if(num)return i * m + (num == m ? 0 : num);
25         b = b * v % n;
26     }
27     return -1;
28 }

 

扩展:

当n不为素数的时候,如何求解ax≡b(mod n) ?

转化成gcd(a, n) = 1即可

如何转化呢?

利用公式:axb(mod n)ax/d ≡ b/(mod n/d),d=gcd(a,n)(这里是a乘上x,上述求解方程为ax

 

每次用一个a和m和b消去gcd(a, n),消δ次,每次n /= g, b /= g,迭代更新下一个g,直到g=1

最后ax变成了ax - δ *a', 方程变成:ax - δ *a' = b' (mod n')     这里a' b' n'都是消因子之后的

满足:a' * g = aδ  b' * g = b  n' * g = n    g = gcd(aδ, n)

如果某一步g不整除b',直接返回-1

如果某一步b' = a',假设此时消因子次数达到cnt次,那就返回cnt

  因为消因子次数cnt次等价于

  原方程:ax-cntacnt = b (mod n)

  消因子后:ax-cnta' = b' (mod n')

  若此时b' = a' 那么等式就可以直接消去a' b',得到:ax-cnt = 1 (mod n),解就是cnt

 

下面求解:ax - δ *a' = b' (mod n')  

可以求出a'的逆元,然后乘过去,用上述的大步小步算法求解

下面介绍一个技巧不用求逆元。

先把上述式子写成ax * tmp = b (mod n)   求出x之后加上δ就是解。

设m = sqrt(n +0.5), x = k * m - q (1 <= k <= m, q <= m)

上式可写成tmp * a k*m-q = b (mod n)

等价于  tmp * a k*m = b * aq(mod n)

可以从1-m枚举k,每次求出tmp * a k*m 判断是不是存在b*aq

最开始的大步小步算法也可以这样写

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 typedef long long ll;
 4 ll pow(ll a, ll b, ll m)
 5 {
 6     ll ans = 1;
 7     a %= m;
 8     while(b)
 9     {
10         if(b & 1)ans = (ans % m) * (a % m) % m;
11         b /= 2;
12         a = (a % m) * (a % m) % m;
13     }
14     ans %= m;
15     return ans;
16 }
17 ll ext_log_mod(ll a, ll b, ll n)
18 {
19     if(b >= n)return -1;//一些特殊情况的判断
20     a %= n;
21     if(b == 1)return 0;
22     //if(n == 1)return -1;
23     ll cnt = 0;//记录消因子次数
24     ll tmp = 1;//存当前a'的值
25     for(ll g = __gcd(a, n); g != 1; g = __gcd(a, n))
26     {
27         if(b % g)return -1;//不能整除
28         b /= g; n /= g; tmp = tmp * a / g % n;
29         cnt++;
30         if(b == tmp)return cnt;
31     }
32 
33     ll m = sqrt(n + 0.5);
34     ll t = b;
35     map<ll, ll>Map;//记录b * a ^ i, i
36     Map[b] = 0;
37     for(int i = 1; i <= m; i++)
38     {
39         b = b * a % n;
40         Map[b] = i;
41     }
42     a = pow(a, m, n);
43     for(int k = 1; k <= m; k++)//枚举k
44     {
45         tmp = tmp * a % n;//求出tmp*a^(k*m)
46         if(Map.count(tmp))return k * m - Map[tmp] + cnt;
47     }
48     return -1;
49 }
50 int main()
51 {
52     ll a, b, p;
53     while(scanf("%lld%lld%lld", &a, &p, &b) != EOF)
54     {
55         ll ans = ext_log_mod(a, b, p);
56         if(ans == -1)printf("Orz,I can’t find D!\n");
57         else printf("%lld\n", ans);
58     }
59     return 0;
60 }

 

展开阅读全文
amp
打赏
0
0 收藏
分享
加载中
更多评论
打赏
0 评论
0 收藏
0
分享
返回顶部
顶部