高次同余方程求解

所谓高次同余不等式，其实就是形如 a^x = b (mod c)的方程，求解 x

一般来说如果 a，b，c很小的直接枚举就行了

但是，当 a,b,c很大的时候，就必须用更优的算法。。

这里，参考AekdyCoin大神的博客。如果你英文好的话也可以直接看wiki

这里谈谈自己的理解，如有错误，欢迎指出（假设你已经看了该博客）：

1) 对于拓展的baby-step gaint-step为什么可以用消因子：

根据同余性质 ac = bc (mod m),且 c ！= 0，那么

a = b (mod (m / (c, m)))

所以，消因子只不过用了此定理，把 c = gcd(c, m)也就是把共同因子消掉而已。。

2）为什么最终答案是 i * m + j + d

因为最终等于解 D*a ^ (i * m) * x = b' (mod c')

d对应的就是 D的次数， j对应的就是 x的次数，所以答案很显然。。

code：

  1 /*
  2    Time:2013.4.29
  3    State:Accepted
  4 */ 
  5 #include <iostream>
  6 #include <cstring>
  7 #include <cstdlib>
  8 #include <cstdio>
  9 #include <cmath>
 10 #define MXN 131071
 11 #define LL __int64
 12 using namespace std;
 13 
 14 struct oo{
 15     int  p , z, next;
 16 } hash[MXN << 1];
 17 bool bo[MXN << 1];
 18 int  a, b, c, top;
 19 
 20 void insert_h(LL p, LL z){ //往hash表里添加元素，采用链式哈希表 
 21     LL k = z & MXN; 
 22     if (bo[k] == false){
 23         bo[k] = true;
 24         hash[k].next = -1;
 25        // printf("p = %lld z = %lld  k = %lld\n",p, z, k);
 26         hash[k].p = p;
 27         hash[k].z = z;    
 28         return;   
 29     }
 30     while (hash[k].next != -1){
 31          if (hash[k].z == z) return;
 32          k = hash[k].next;  
 33     }
 34     if (hash[k].z == z) return;
 35     hash[k].next = ++top;  
 36   //  printf("p = %lld  z = %lld  k = %lld\n",p , z,  top);
 37     
 38     hash[top].next = -1;
 39     hash[top].z = z;
 40     hash[top].p = p;
 41 }
 42 
 43 LL find(LL z){ //在哈希表里查找元素 
 44     LL k = z & MXN;
 45     if (bo[k] == false) return -1;
 46     while (k != -1){
 47         if( hash[k].z == z ) return hash[k].p;
 48         k = hash[k].next;
 49     }
 50     return -1;
 51 }
 52 
 53 LL extend_gcd(LL a, LL b, LL &x, LL &y){ //拓展欧几里德求出解 
 54      if (b == 0){
 55         x = 1;
 56         y = 0;
 57         return a;
 58      }
 59      LL ret = extend_gcd(b , a % b , x, y),  t = x;
 60      x = y;
 61      y = t - a / b * y;
 62 }
 63 
 64 LL gcd(LL a, LL b){
 65     return b == 0 ? a : gcd(b , a % b);
 66 }
 67 
 68 LL pow_mod(LL a, LL n, LL c){ //非递归形式求出 a^n mod c 的结果，
 69      LL ret = 1 % c;  
 70      a %= c;
 71      while (n){
 72         if (n & 1)  ret = ret * a % c;
 73         a = a * a % c;
 74         n = n >> 1;
 75      }
 76      return ret % c;
 77 }
 78 
 79 
 80 
 81 LL baby_step(LL a, LL b, LL c){ //主要程序  
 82      LL tmp = 1 % c, i;
 83      top = MXN; 
 84      for (i = 0; i <= 100; ++i, tmp = tmp * a % c)  //答案可能小于D对应的次数是枚举一下 
 85           if (tmp == b) return i;
 86      int cnt = 0; 
 87      LL d = 1;
 88      while ((tmp = gcd(a, c)) != 1){ //消因子 
 89           if (b % tmp) return -1; //如果 b不是tmp倍数显然无解 
 90           b /= tmp;
 91           c /= tmp;
 92           d =  d * a / tmp % c;
 93           ++cnt;  
 94      }
 95      LL m = (LL)ceil(sqrt(c + 0.00));
 96    //  printf("m = %lld\n", m);
 97      for (tmp = 1 % c, i = 0; i <= m;  ++i, tmp = tmp *a % c)
 98             insert_h(i, tmp); 
 99      LL x, y, k = pow_mod(a, m ,c); //求出 a^m mod c 
100  //    printf("d = %lld\n", d );
101      for (i = 0; i <= m; ++i){
102          extend_gcd(d, c, x, y); //求出 dx = 1 （mod c） 
103          tmp = ((b * x) % c + c) % c;
104         // printf("tmp = %lld\n", tmp);
105          if ((y = find(tmp)) != -1){
106            //  printf("y = %lld tmp = %lld\n", y, tmp);
107              return i * m + y + cnt;
108          }
109          d = d * k % c;
110      }
111      return -1;
112 }
113 
114 int main(){
115     freopen("poj3243.in", "r", stdin);
116     freopen("poj3243.out","w", stdout);
117     while (scanf("%d%d%d", &a, &c, &b) != EOF){
118          if (!a && !b && !c) break;
119          b %= c;
120          memset(bo, 0, sizeof(bo));
121          LL tmp = baby_step(a, b, c);  
122          if (tmp < 0) puts("No Solution");
123          else printf("%lld\n", tmp); 
124     }
125     fclose(stdin);  fclose(stdout);
126 }