所谓高次同余不等式,其实就是形如 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 }