[poj 2417]Discrete Logging 数论 BSGS

Time Limit: 5000MS		Memory Limit: 65536K
Total Submissions: 3516		Accepted: 1651

Description

Given a prime P, 2 <= P < 231, an integer B, 2 <= B < P, and an integer N, 1 <= N < P, compute the discrete logarithm of N, base B, modulo P. That is, find an integer L such that

BL == N (mod P)

Input

Read several lines of input, each containing P,B,N separated by a space.
Output

For each line print the logarithm on a separate line. If there are several, print the smallest; if there is none, print "no solution".
Sample Input

5 2 1
5 2 2
5 2 3
5 2 4
5 3 1
5 3 2
5 3 3
5 3 4
5 4 1
5 4 2
5 4 3
5 4 4
12345701 2 1111111
1111111121 65537 1111111111

Sample Output

0
1
3
2
0
3
1
2
0
no solution
no solution
1
9584351
462803587

Hint

The solution to this problem requires a well known result in number theory that is probably expected of you for Putnam but not ACM competitions. It is Fermat's theorem that states

B(P-1) == 1 (mod P)

for any prime P and some other (fairly rare) numbers known as base-B pseudoprimes. A rarer subset of the base-B pseudoprimes, known as Carmichael numbers, are pseudoprimes for every base between 2 and P-1. A corollary to Fermat's theorem is that for any m

B(-m) == B(P-1-m) (mod P) .

Source
Waterloo Local 2002.01.26
题目大意
给出B,N,P 求出 B^L == N (mod P)中L的最小值，无解输出no solution

解题思路
直接按照poj 1395的思路模拟铁定超时，我们这个时候考虑：
P是质数那么B^(P-1)=1(mod P) [据费马小定理，题中B<P]
可知剩余系的循环节为p-1
根据同余定理，我们可以把L拆解成两部分
此时 B^A * B^(L-A) = N (MOD P)
考虑A为何值时我们可以简化计算
BSGS就是利用这样的思想 将剩余系的余数分拆成floor(sqrt(P))组
再设每一组元素为c个
每次我们的A都可以取得A=c*i ,同时满足c*i<(p-1)
则L-A<c
c约为sqrt(P)，我们可以直接打出前sqrt(p)项的表存到map之中待处理

那么我们再枚举i，来改变A的值
通过求逆元来推测一个B^(L-A) mod P可能的值，查map，看这种方案是否有对应解。

#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <map>
#define LL long long 
using namespace std;
typedef int Ma[10][10];
map <LL , LL>inv;
void egcd(LL a,LL b,LL &x,LL &y)
{
  if (b==0)
  {
    x=1;
    y=0;
    return;
  }
  egcd(b,a%b,x,y);
  LL t=x;
  x=y;y=t-a/b*y;
}
int main()
{
  LL p,b,n;
  while (~scanf("%I64d%I64d%I64d",&p,&b,&n))
  {
    inv.clear();
    LL c=sqrt((double) p);
    // printf("%I64d\n",c);
    LL cur=1;
    for (LL i=0;i<c;i++)
    {
      inv[cur]=i+1;
      cur=cur*b%p;
    }
    LL t=cur;
    LL ans=p+1;
    cur=1;
    for (LL i=0;i<=(p/c+1);i++)
    {
      LL x,y;
      egcd(cur,p,x,y);
      x=(x+p)%p;
      // printf("%I64d %I64d\n",cur,i);
      if (inv[n*x%p]) {
        ans=min(ans,inv[n*x%p]+i*c-1);
        // printf("%I64d\n",ans);
      }
      cur=cur*t%p;
    }
    if (ans>p) puts("no solution");
    else printf("%I64d\n",ans);
  }
  return 0;
}