BSGS算法逆元

【写在前面】：
    其实很早之前就已经知道这个东西了
    最近做了一道BSGS的题目，有感而发
    决定扯一下淡，顺便加深一下记忆
    为啥要把BSGS和逆元一起讲呢？因为BSGS里面也牵涉到了逆元的相关知识（虽然并没有什么关系）
【参考】：

【线性求逆元】：
    我们称 A^(-1) mod p为A在mod p意义下的逆元
    显然在p为质数且(A,p)=1时可以用欧拉定理搞定
    然A p没有任何联系时，我们可以用扩展欧几里德求出，这样求解[1,n]的逆元复杂度为O(n log n)
    我们考虑一个SB的事实p mod p=0
    然后我们令p=k*i+r (1<=i<p r<i)
    显然k*i+r mod p=0
    两边同乘i^(-1)*r^(-1)得到 k*r^(-1)+i^(-1) mod p=0
    i^(-1) mod p=-k*r^(-1)
    显然k=p/i，r=p%i代入原式可以得到i^(-1) mod p=-p/i*(p%i)^(-1)
    因此我们令rev[i]表示i的逆元，那么rev[i]=-(p/i)*rev[p%i]
    递推处理即可，复杂度O(N)
    同时这篇blog还提供了O(log p)复杂度求出单个逆元
    p%i<=p/2，然后就可以递归处理了
【简介】：
    首先BSGS的全称是：Baby Step Giant Step
    它主要用来求解这样的问题：A^x mod p=B的最小非负整数解
    BSGS分为质数版和extend版，我们分开总结一下
【质数版BSGS】：
    所谓质数版指的是p为指数的情况
    我们考虑p为指数时，根据费马小定理可以知道 0<=x<p
    显然暴力枚举的复杂度是O(p)的，当p较大时显然不适用
    我们可以这样考虑：令x=C*i+j（其中C为常数）
    代入原式可得(A^C)^i*A^j mod p=B
    显然我们可以用0<=i<=p/C
    于是我们可以用O(C)的时间求出A^C 然后再用O(p/C)的时间枚举i
    于是问题就转化为是否存在一个j满足D*A^j mod p=B，D为确定的值
    A^j mod p=B/D，其中0<=j<C
    我们可以用O(C)的时间预处理出A^j mod p的值，存入map里面
    然后可以在log的时间复杂度内查询最小的j满足条件
    显然我们取C=sqrt(p)，那么整个算法的复杂度就是O(sqrt(p) log sqrt(p))
    B/D可以利用扩展欧几里德或者欧拉定理算出来
【extend版BSGS】
    如果p不是质数，那么可能会出现gcd(A,p)!=1的时候
    我们考虑在求解D*A^j mod p=B时候如果gcd(D,p)!=1时的情况
    由于D=A^C，因此就是gcd(A,p)!=1的情况，此时不能用扩展欧几里德求出A^j
    我们考虑一个重要的性质：Ad mod Bd = Cd等价于A mod B = C
    因此我们考虑对于gcd(A,p)!=1的情况，我们可以暴力消因子
    我们考虑A^x mod p=B
    我们不断的求出gcd(A,p)直到(A,p)=1为之，令t=gcd(A,p)
    如果B mod t!=0显然当前方程无解，那么对应的原方程也无解
    我们令总共消掉的数为D，那么消完后的方程等价于(A/D)^x mod (p/D) = B/D
    然而这样子做对于A/D还需要求解逆元然后就不好做了
    我们这样做：（参考：https://oi.abcdabcd987.com/bsgs/）

    D = 1, cnt = 0
    while gcd(A, p) != 1:
    if (B % gcd(A, p) != 0)
        No Solution
    B /= gcd(A, p)
    p /= gcd(A, p)
    D *= A / gcd(A, p)
    cnt += 1

    那么之后的方程就变成D*A^(x-cnt) mod p=B（注意此时的B,p对应消过因子的后的数）
    此时gcd(D,p)=1就可以用BSGS求出x-cnt，然后加上cnt即可

    但是上述的方法其实是存在bug的
    我们求出的x-cnt>=0，因此x>=cnt，那么x<cnt的解怎么办？
    我们知道cnt最多为p中2的个数，因此cnt是O(log p)级别的，我们暴力枚举这一部分的情况即可
【一个很实用的trick】：
    我们考虑质数版BSGS
    我们不妨令x=C*i-j (1<=i<=p/C)
    那么(A^C)^i mod p=B*A^j，我们可以预处理出B*A^j，然后map即可
    然后这样似乎没什么大的改进？真的吗？你会求逆矩阵吗，反正我是不会
    我们看一下这题：bzoj 4128: Matrix
    我们利用这个trick，既不用逆矩阵，也不用hash，直接裸上map即可
    不过这么实用的技巧为啥很少人用捏？
【总结】：
    简言之，BSGS在求解一类高次同余问题中是非常实用的
    而它也可以与各种其他的东西相结合
    我们需要针对不同的题目实用不同的做法
    当然extend版的BSGS是非常实用的东西
    就这样了...
    （下面是一些例题）

bzoj 4128: Matrix

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <map>
using namespace std;
const int N=75;
int n,maxmod,ans,m;
struct Matrix
{
    int a[N][N];
    bool operator < (const Matrix &b) const
    {
        for (int i=1; i<=n; i++)
            for (int j=1; j<=n; j++)
            {
                if (a[i][j]<b.a[i][j]) return true;
                if (a[i][j]>b.a[i][j]) return false;
            }
        return false;
    }
    bool operator == (const Matrix &b) const
    {
        for (int i=1; i<=n; i++)
            for (int j=1; j<=n; j++)
            if (a[i][j]!=b.a[i][j]) return false;
        return true;
    }
    inline void init(int flag)
    {
        for (int i=1; i<=n; i++)
            for (int j=1; j<=n; j++)
            a[i][j]=min(flag,(int)(i==j));
    }
} a,b,now,Now,A,B;
map<Matrix,int> mp;
inline Matrix operator * (const Matrix &a,const Matrix &b)
{
    Matrix c; c.init(0);
    for (int i=1; i<=n; i++)
        for (int j=1; j<=n; j++)
            for (int k=0; k<=n; k++)
            c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%maxmod)%maxmod;
    return c;
}
inline long long getint()
{
    long long x=0; char c=getchar(); bool flag=false;
    while ((c!='-')&&((c<'0')||(c>'9'))) c=getchar();
    if (c=='-') flag=true,c=getchar();
    while ((c>='0')&&(c<='9')) x=x*10+(long long)(c-'0'),c=getchar();
    if (flag) return -x; else return x;
}
void init()
{
    n=getint(); maxmod=getint(); m=ceil(sqrt((double)maxmod));
    for (int i=1; i<=n; i++) for (int j=1; j<=n; j++) a.a[i][j]=getint();
    for (int i=1; i<=n; i++) for (int j=1; j<=n; j++) b.a[i][j]=getint();
}
void build()
{
    now.init(1); now=now*b; mp[now]=0;
    for (int i=1; i<m; i++) now=now*a,mp[now]=i;
    now.init(1); for (int i=1; i<=m; i++) now=now*a;
}
void solve()
{
    Now.init(1);
    for (int i=1; i<=m; i++)
    {
        Now=Now*now;
        if (mp.count(Now)) {ans=i*m-mp[Now]; break;}
    }
    printf("%d\n",ans);
}
int main()
{
    init();
    build();
    solve();
    return 0;
}

bzoj 2242

type
  point=record
    id:longint; cnt:int64;
  end;
var
  baby:array [0..100001] of point;
  n,m,maxn,tot:longint; a,b,p,x,y:int64;
procedure init;
begin
  readln(n,m);
end;
function quickmul(a,b,p:int64):int64;
var
  t:int64;
begin
  t:=1; a:=a mod p;
  while b<>0 do
  begin
    if b mod 2=1 then t:=(t*a) mod p;
    a:=(a*a) mod p; b:=b div 2;
  end;
  exit(t);
end;
function ex_gcd(a,b:int64):int64;
var
  tt,dd:int64;
begin
  if b=0 then begin x:=1; y:=0; exit(a); end;
  dd:=ex_gcd(b,a mod b);
  tt:=x; x:=y; y:=tt-(a div b)*y; exit(dd);
end;
function solve(a,b,p:int64):int64;
var
  t,d,plus:int64;
begin
  d:=ex_gcd(a,p);
  if b mod d<>0 then exit(-1);
  t:=b div d; x:=x*t; x:=(x mod p+p) mod p; exit(x);
end;
procedure qsort(x,y:longint);
var
  p,q:longint; mid:point;
begin
  p:=x; q:=y; mid:=baby[(x+y) div 2];
  repeat
    while (baby[x].cnt<mid.cnt) or ((baby[x].cnt=mid.cnt) and (baby[x].id<mid.id)) do inc(x);
    while (baby[y].cnt>mid.cnt) or ((baby[y].cnt=mid.cnt) and (baby[y].id>mid.id)) do dec(y);
    if x<=y then
    begin
      baby[0]:=baby[x];
      baby[x]:=baby[y];
      baby[y]:=baby[0];
      inc(x); dec(y);
    end;
  until x>y;
  if p<y then qsort(p,y);
  if x<q then qsort(x,q);
end;
function find(x:int64):longint;
var
  l,r,mid,k:longint;
begin
  k:=0; l:=1; r:=tot;
  while l<=r do
  begin
    mid:=(l+r) div 2;
    if baby[mid].cnt=x then exit(mid);
    if baby[mid].cnt<x then l:=mid+1 else r:=mid-1;
  end;
  exit(0);
end;
procedure baby_giant_step(a,b,p:int64);
var
  i,j,t:longint; k:int64;
begin
  baby[1].id:=0; baby[1].cnt:=1; maxn:=trunc(sqrt(p-1));
  if sqrt(p-1)<>trunc(sqrt(p-1)) then inc(maxn); inc(maxn);
  for i:=2 to maxn do
  begin
    baby[i].id:=i-1;
    baby[i].cnt:=(baby[i-1].cnt*a) mod p;
  end;
  qsort(1,maxn); tot:=1;
  for i:=2 to maxn do
  begin
    if baby[i].cnt=baby[i-1].cnt then continue;
    inc(tot); baby[tot]:=baby[i];
  end;
  for i:=0 to maxn do
  begin
    k:=quickmul(a,maxn,p);
    k:=quickmul(k,i,p);
    k:=solve(k,b,p);
    t:=find(k);
    if t=0 then continue;
    writeln((i*maxn+baby[t].id) mod p);
    exit;
  end;
  writeln('Orz, I cannot find x!');
end;
procedure main;
var
  i:longint; ans:int64;
begin
  for i:=1 to n do
  begin
    readln(a,b,p);
    if m=1 then writeln(quickmul(a,b,p));
    if m=2 then begin ans:=solve(a,b,p); if ans=-1 then writeln('Orz, I cannot find x!') else writeln(ans); end;
    if m=3 then baby_giant_step(a,b,p);
  end;
end;
begin
  init;
  main;
end.

zhouzixuan的博客

博客

评论

发表评论

zhouzixuan的博客

博客

BSGS算法 逆元

评论

发表评论

BSGS算法逆元