exBSGS算法

BSGS，全称(Baby Step Giant Step)，是用于求解离散对数的一种算法。
就是用来求(A^x equiv B (mod p)) 的x这么一种算法……

理论知识是：在[0,p)之内是一定有解的，因为指数模的周期性。即(A^x)对p的模随x变化有周期性，最大周期不超过p。首先，余数只可能有p个元素，所以x取不超过p个值必定出现相同的余数。根据同余的性质，只要找到两个余数相同的，剩下的全部乘以(A^k,k)为整数，所以(A^n equiv A^{x+n})对于所有x都成立。

朴素算法是枚举0～p-1。
如何改进呢？我们考虑把数分组，每组大小为(n = sqrt{p}),每组(m = frac{p}{n})个数。
对于每组的询问，我们组内的答案可以看作：(A^{im-y} equiv B (mod p))
移项一下就有：(A^{im} equiv BA^y (mod p))
我们枚举(0～m-1),计算所有的(BA^y),存进hash表，之后枚举(1～n)计算(A^{im})的值，在hash表中找答案即可。

时间复杂度(O(sqrt{n}))

然而这个算法只适用于……(gcd(A,p) = 1)，即二者互质的情况。否则此时(A^{-1})不存在，就不能把答案看作(A^{im-y})

解决的方法是预先处理成二者互质。
我们一直把方程同除(A,p)的公约数就行了。如果中途某一次B除不开就是无解的，否则一直除，得到(B_n = frac{A^n}{d_1d_2…d^n}A^{x-n} + p_ny),令(D = frac{A^n}{d_1d_2…d^n}),此时必存在D的逆元，所以式子变为(A^{x-n} equiv B_nD^{-1} (mod p_n))，求解即可。

这有个三合一的板子题。
SDOI2011 计算器

#include<bits/stdc++.h>
#define rep(i,a,n) for(int i = a;i <= n;i++)
#define per(i,n,a) for(int i = n;i >= a;i--)
#define enter putchar('
')
#define I inline

using namespace std;
typedef long long ll;
const int M = 200005;
const int mod = 999979;

int read()
{
   int ans = 0,op = 1;char ch = getchar();
   while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
   while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
   return ans * op;
}

int T,K,y,z,p;

struct Hash
{
   int sta[M],top,head[M<<3],num[M],val[M],nxt[M],ecnt;
   void init(){ecnt = 0;while(top) head[sta[top--]] = 0;}
   void insert(int x,int y)
   {
      int h = x % mod;
      for(int i = head[h];i;i = nxt[i]) if(num[i] == x) {val[i] = y;return;}
      if(!head[h]) sta[++top] = h;
      nxt[++ecnt] = head[h],head[h] = ecnt;
      num[ecnt] = x,val[ecnt] = y;
   }
   int query(int x)
   {
      int h = x % mod;
      for(int i = head[h];i;i = nxt[i]) if(num[i] == x) return val[i];
      return -1;
   }
}H;

int gcd(int a,int b) {return (!b) ?  a : gcd(b,a%b);}
int mul(int a,int b,int t){return 1ll * a * b % t;}
int qpow(int a,int b,int t)
{
   int g = 1;
   while(b)
   {
      if(b&1) g = mul(g,a,t);
      a = mul(a,a,t),b >>= 1;
   }
   return g;
}

int exgcd(int a,int b,int &x,int &y)
{
   if(!b){x = 1,y = 0;return a;}
   int d = exgcd(b,a%b,y,x);
   y -= a / b * x;
   return d;
}

int inv(int a,int b)
{
   int x,y;
   exgcd(a,b,x,y);
   return (x % b + b) % b;
}

int BSGS(int a,int b,int c)
{
   int cnt = 0,G,d = 1;
   while((G = gcd(a,c)) != 1)
   {
      if(b % G) return -1;
      cnt++,b /= G,c /= G,d = mul(d,a/G,c);
   }
   b = mul(b,inv(d,c),c);
   H.init();
   int s = sqrt(c),p = 1;
   rep(i,0,s-1)
   {
      if(p == b) return i + cnt;
      H.insert(mul(p,b,c),i),p = mul(p,a,c);
   }
   int q = p,t;
   for(int i = s;i <= c + s - 2;i += s)
   {
      t = H.query(q);
      if(t != -1) return i - t + cnt;
      q = mul(q,p,c);
   }
   return -1;
}

int main()
{
   T = read(),K = read();
   while(T--)
   {
      y = read(),z = read(),p = read();
      if(K == 1) printf("%d
",qpow(y,z,p));
      if(K == 2)
      {
	 if(z % gcd(y,p)) {printf("Orz, I cannot find x!
");continue;}
	 int X,Y;
	 int d = exgcd(y,p,X,Y);
	 p /= d;
	 printf("%d
",(mul(X,z/d,p) + p) % p);
      }
      if(K == 3)
      {
	 y %= p,z %= p;
	 int ans = BSGS(y,z,p);
	 (ans == -1) ? printf("Orz, I cannot find x!
") : printf("%d
",ans);
      }
   }
   return 0;
}