[Baby steps giant steps] BSGS EXT-BSGS

BSGS

给定y,z,p，计算满足Y^x ≡ Z ( mod P)的最小非负整数。P为质数 或者更确切的 gcd(Y,P)|Z

还是引用黄学长的吧

用hash显然比map优

EXT_BSGS

当模不一定为质数时

我们要对算法进行一些修改，使得其可以进行上述的处理。具体证明去看AekdyCoin犇的题解，我就是简单讲一下我的理解。

一开始的方程等价于A^x+C*b=B(a,b是整数)，现在gcd(A，C)!=1,不妨令t=gcd(A,C)

我们考虑方程(A/t)A^x’%(C/t)=(B/t)的解x’与原来的解x有什么关系。

显然现在的方程等价于(A/t)A^x’+(C/t)*b’=B/t.

两端均乘以t得到：A^(x’+1)+C*b’=B

由系数相等有:x=x’+1,b=b’.

那么我们就有一种方法算出x了。先算出x’,再加上1就行了。

考虑如何求出x’,现在的方程如果A,C依旧不互质，就继续迭代将系数除以最大公约数，同时前面的常数增大。

然后求出解之后再加上总共除的次数就好了。

BZOJ3239

裸的BSGS 同poj2417

#include<cstdio>#include<cstdlib>#include<algorithm>#include<map>#include<cmath>using namespace std;typedef long long ll;ll A,B,P;inline ll Pow(ll a,ll b){if (!b) return 1LL;ll ret=Pow(a,b>>1);if (b&1) return ret*ret%P*a%P;return ret*ret%P;}map<ll,int> M;inline void solve(){M.clear();A%=P;if (!A && !B) { printf("1\n"); return; }if (!A) { printf("no solution\n"); return; }ll m=ceil(sqrt((double)P)),tmp=1;M[1]=m+1;for (int i=1;i<m;i++){(tmp*=A)%=P;if (!M[tmp]) M[tmp]=i;}ll inv=1,T=Pow(A,P-1-m);for (int k=0;k<m;k++){int i=M[B*inv%P];if (i){if (i==m+1) i=0;printf("%lld\n",k*m+i);return;}(inv*=T)%=P;}printf("no solution\n");}int main(){freopen("t.in","r",stdin);freopen("t.out","w",stdout);while (~scanf("%lld%lld%lld",&P,&A,&B))solve();return 0;}

BZOJ2242

一堆乱七八糟的操作放在一起

#include<cstdio>#include<cstdlib>#include<algorithm>#include<map>#include<cmath>using namespace std;typedef long long ll;inline char nc(){static char buf[100000],*p1=buf,*p2=buf;if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }return *p1++;}inline void read(ll &x){char c=nc(),b=1;for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;}inline void write(ll x,char c){if (x==0) { putchar('0'); putchar(c); return; }if (x<0) putchar('-'),x=-x;char s[55]={0},len=0;while (x) s[++len]=x%10+'0',x/=10;for (int i=len;i;i--) putchar(s[i]); putchar(c);}namespace Work1{ll A,B,P;inline ll Pow(ll a,ll b){    if (!b) return 1LL;    ll ret=Pow(a,b>>1);    if (b&1) return ret*ret%P*a%P;    return ret*ret%P;}inline void work(){read(A); read(B); read(P);write(Pow(A,B),'\n');}}namespace Work2{#define X first#define Y secondtypedef pair<ll,ll> data;ll A,B,P;ll D;inline data EXGCD(ll a,ll b){if (!b) { D=a; return data(1,0); }data ret=EXGCD(b,a%b);return data(ret.Y,ret.X-a/b*ret.Y);}inline void solve(){data ret=EXGCD(A,P);if (B%D==0){ll ans=((ret.X*B/D)%P+P)%P;write(ans,'\n');}elseprintf("Orz, I cannot find x!\n");}inline void work(){read(A); read(B); read(P);solve();}}namespace Work3{ll A,B,P;inline ll Pow(ll a,ll b){    if (!b) return 1LL;    ll ret=Pow(a,b>>1);    if (b&1) return ret*ret%P*a%P;    return ret*ret%P;}map<ll,int> M;inline void solve(){    M.clear();    A%=P;    if (!A && !B) { printf("1\n"); return; }    if (!A) { printf("Orz, I cannot find x!\n"); return; }    ll m=ceil(sqrt((double)P)),tmp=1;    M[1]=m+1;    for (int i=1;i<m;i++)    {        (tmp*=A)%=P;        if (!M[tmp]) M[tmp]=i;    }    ll inv=1,T=Pow(A,P-1-m);    for (int k=0;k<m;k++)    {        int i=M[B*inv%P];        if (i)        {            if (i==m+1) i=0;            write(k*m+i,'\n');            return;        }        (inv*=T)%=P;    }    printf("Orz, I cannot find x!\n");}inline void work(){read(A); read(B); read(P);solve();}} ll Q,K; int main(){freopen("t.in","r",stdin);freopen("t.out","w",stdout);read(Q); read(K);while (Q--){if (K==1)Work1::work();else if (K==2)Work2::work();else if (K==3)Work3::work();}    return 0;}

BZOJ3122

求数列通项 再用BSGS解

#include<cstdio>#include<cstdlib>#include<algorithm>#include<cmath>#include<map>using namespace std;typedef long long ll;inline char nc(){static char buf[100000],*p1=buf,*p2=buf;if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }return *p1++;}inline void read(ll &x){char c=nc(),b=1;for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;}inline void write(ll x,char c){if (x==0) { putchar('0'); putchar(c); return; }if (x<0) putchar('-'),x=-x;char s[55]={0},len=0;while (x) s[++len]=x%10+'0',x/=10;for (int i=len;i;i--) putchar(s[i]); putchar(c);}namespace MOD{#define X first#define Y secondtypedef pair<ll,ll> data;ll D;inline data EXGCD(ll a,ll b){if (a<b) { data ret=EXGCD(b,a); return data(ret.Y,ret.X); }if (!b) { D=a; return data(1,0); }data ret=EXGCD(b,a%b);return data(ret.Y,ret.X-a/b*ret.Y);}inline ll Inv(ll A,ll P){data ret=EXGCD(A,P);return (ret.X%P+P)%P;}inline ll Solve(ll A,ll B,ll P){data ret=EXGCD(A,P);if (B%D==0){ll ans=((ret.X*B/D)%P+P)%P;return ans;}elsereturn -1;}}namespace BSGS{inline ll Pow(ll a,ll b,ll p){ll ret=1;for (;b;b>>=1,a=a*a%p)if (b&1)(ret*=a)%=p;return ret;}map<ll,int> M;inline ll Solve(ll A,ll B,ll P){M.clear();A%=P;if (!A && !B) return 1LL;if (!A) return -1;ll m=sqrt((double)P)+1,tmp=1;M[1]=m+1;for (int i=1;i<m;i++){(tmp*=A)%=P;if (!M[tmp]) M[tmp]=i;}ll inv=1,T=Pow(A,P-1-m,P);for (int k=0;k<m;k++){int i=M[B*inv%P];if (i){if (i==m+1) i=0;return k*m+i;}(inv*=T)%=P;}return -1;}}ll P,A,B,X,T;ll N;int main(){ll Q;freopen("t.in","r",stdin);freopen("t.out","w",stdout);read(Q);while (Q--){read(P); read(A); read(B); read(X); read(T);if (X==T) { printf("1\n"); continue; }if (A==0){if (B==T) { printf("2\n"); continue; }else { printf("-1\n"); continue; }}if (A==1){N=MOD::Solve(B,((T-X+B)%P+P)%P,P);if (N==0) N=P;}else {ll C=MOD::Inv(A-1,P);ll AN=MOD::Solve((X+B*C)%P,(B*C+T)%P,P);if (AN==-1)N=-1;else{N=BSGS::Solve(A,AN,P);if (N!=-1) N++;}}write(N,'\n');}return 0;}

poj3243

扩展BSGS

#include<cstdio>#include<cstdlib>#include<cmath>#include<map>#include<algorithm>#include<cstring>using namespace std;typedef long long ll;inline char nc(){static char buf[100000],*p1=buf,*p2=buf;if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }return *p1++;}inline void read(ll &x){char c=nc(),b=1;for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;}inline void write(ll x,char c){if (x==0) { putchar('0'); putchar(c); return; }if (x<0) putchar('-'),x=-x;char s[55]={0},len=0;while (x) s[++len]=x%10+'0',x/=10;for (int i=len;i;i--) putchar(s[i]); putchar(c);}#define mod 1313131struct Hashset {int head[mod], next[35010], f[35010], v[35010], ind;void reset() {ind = 0;memset(head, -1, sizeof head);}void insert(int x, int _v) {int ins = x % mod;for(int j = head[ins]; j != -1; j = next[j])if (f[j] == x) {v[j] = min(v[j], _v);return;}f[ind] = x, v[ind] = _v;next[ind] = head[ins], head[ins] = ind++;}int operator [] (const int &x) const {int ins = x % mod;for(int j = head[ins]; j != -1; j = next[j])if (f[j] == x)return v[j];return -1;}}M;#define X first#define Y secondtypedef pair<ll,ll> data;inline data EXGCD(ll a,ll b){if (a<b) { data ret=EXGCD(b,a); return data(ret.Y,ret.X); }if (!b) { return data(1,0); }data ret=EXGCD(b,a%b);return data(ret.Y,ret.X-a/b*ret.Y);}inline ll Inv(ll A,ll P){data ret=EXGCD(A,P);return (ret.X%P+P)%P;}inline ll Pow(ll a,ll b,ll p){ll ret=1;for (;b;b>>=1,a=a*a%p)if (b&1)(ret*=a)%=p;return ret;}//map<ll,int> M;inline ll BSGS(ll A,ll B,ll P){//M.clear();M.reset();A%=P;if (!A && !B) return 1LL;if (!A) return -1;ll m=sqrt((double)P)+1,tmp=1;M.insert(1,m+1);for (int i=1;i<m;i++){(tmp*=A)%=P;if (M[tmp]==-1) M.insert(tmp,i);}ll inv=1,T=Inv(Pow(A,m,P),P);for (int k=0;k<m;k++){int i=M[B*inv%P];if (i!=-1){if (i==m+1) i=0;return k*m+i;}(inv*=T)%=P;}return -1;}ll A,B,P,ans;inline ll Calc(ll A,ll B,ll P){ll tmp=1;for (int i=0;i<100;(tmp*=A)%=P,i++)if (tmp==B)return i;ll C=1,D,num=0,ret;data G;while (1){G=EXGCD(A,P);D=A*G.first+P*G.second;if (D==1){(B*=Inv(C,P))%=P;ret=BSGS(A,B,P);if (ret==-1) return -1;return ret+num;}if (B%D) return -1;B/=D; P/=D;(C*=A/D)%=P;num++;}}int main(){freopen("t.in","r",stdin);freopen("t.out","w",stdout);while (1){read(A); read(P); read(B);if (!A && !B && !P) break;ans=Calc(A,B,P);if (ans==-1)printf("No Solution\n");elsewrite(ans,'\n');}return 0;}

BZOJ3283

又是一堆乱七八糟的，组合数取模参见另一篇博文

#include<cstdio>#include<cstdlib>#include<algorithm>#include<cmath>#include<map>#include<cstring>#define cl(x) memset(x,0,sizeof(x))using namespace std;typedef long long ll;inline char nc(){static char buf[100000],*p1=buf,*p2=buf;if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }return *p1++;}inline void read(ll &x){char c=nc(),b=1;for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;}inline void write(ll x,char c){if (x==0) { putchar('0'); putchar(c); return; }if (x<0) putchar('-'),x=-x;char s[55]={0},len=0;while (x) s[++len]=x%10+'0',x/=10;for (int i=len;i;i--) putchar(s[i]); putchar(c);}inline ll Mul(ll A,ll B,ll P){A%=P,B%=P;ll ret=A*B-(ll)((double)A/P*B+0.5)*P;return ret<0?ret+P:ret;}inline ll Pow(ll A,ll B,ll P){ll ret=1;for (;B;B>>=1,A=Mul(A,A,P))if (B&1)ret=Mul(ret,A,P);return ret;}#define X first#define Y secondtypedef pair<ll,ll> data;inline data EXGCD(ll a,ll b){if (a<b) { data ret=EXGCD(b,a); return data(ret.Y,ret.X); }if (!b) { return data(1,0); }data ret=EXGCD(b,a%b);return data(ret.Y,ret.X-a/b*ret.Y);}inline ll Inv(ll A,ll P){data ret=EXGCD(A,P);return (ret.X%P+P)%P;}namespace Work1{ll A,B,P;inline void Solve(){read(A); read(B); read(P);write(Pow(A,B,P),'\n');}}namespace Work2{#define mod 1313131struct Hashset {int head[mod], next[35010], f[35010], v[35010], ind;void reset() {ind = 0;memset(head, -1, sizeof head);}void insert(int x, int _v) {int ins = x % mod;for(int j = head[ins]; j != -1; j = next[j])if (f[j] == x) {v[j] = min(v[j], _v);return;}f[ind] = x, v[ind] = _v;next[ind] = head[ins], head[ins] = ind++;}int operator [] (const int &x) const {int ins = x % mod;for(int j = head[ins]; j != -1; j = next[j])if (f[j] == x)return v[j];return -1;}}M;inline ll BSGS(ll A,ll B,ll P){M.reset();A%=P;if (!A && !B) return 1LL;if (!A) return -1;ll m=sqrt((double)P)+1,tmp=1;M.insert(1,0);for (int i=1;i<m;i++){(tmp*=A)%=P;if (M[tmp]==-1) M.insert(tmp,i);}ll inv=1,T=Inv(Pow(A,m,P),P);for (int k=0;k<m;k++){int i=M[B*inv%P];if (i!=-1)return k*m+i;(inv*=T)%=P;}return -1;}ll A,B,P,ans;inline ll Calc(ll A,ll B,ll P){for (ll tmp=1,i=0;i<100;(tmp*=A)%=P,i++)if (tmp==B)return i;ll C=1,D,num=0,ret;data G;while (1){G=EXGCD(A,P);D=A*G.first+P*G.second;if (D==1){(B*=Inv(C,P))%=P;ret=BSGS(A,B,P);if (ret==-1) return -1;return ret+num;}if (B%D) return -1;B/=D; P/=D;(C*=A/D)%=P;num++;}}inline void Solve(){read(A); read(B); read(P);ans=Calc(A,B,P);if (ans==-1)printf("Math Error\n");elsewrite(ans,'\n');}}namespace Work3{namespace China{ll A[100005],N[100005],M[100005];int tot=0;ll P=1;inline void clear(){P=1; tot=0;}inline void add(ll a,ll n){A[++tot]=a; N[tot]=n; P*=n;}inline ll Solve(){ll ret=0,inv;for (int i=1;i<=tot;i++){M[i]=P/N[i];inv=Inv(M[i],N[i]);(ret+=Mul(A[i],Mul(M[i],inv,P),P))%=P;}return ret;}}int vst[100005],prime[100005],num;  inline void Pre()   {      const int maxn=100000;      for (int i=2; i<=maxn; i++)          if (!vst[i])           {              prime[++num]=i;              for (int j=i+i; j<=maxn; j+=i) vst[j]=1;          }  } ll M,ans;  ll p[50005],a[50005],p_a[50005],cnt;ll n,m;inline void Machine(ll M)   {      for (int i=1; i<=num; i++)          if (M%prime[i]==0) {              p[++cnt]=prime[i];              while (M%prime[i]==0) a[cnt]++,M/=prime[i];              p_a[cnt]=Pow(p[cnt],a[cnt],1000000);          } } inline data Calc(ll x,ll p,ll a,ll p_a){ll A=1,B=1;for (ll i=1;i<p_a;i++)if (i%p)(A*=i)%=p_a;A=Pow(A,x/p_a,p_a);for (ll i=x-x%p_a+1;i<=x;i++)if (i%p)(A*=i)%=p_a;B=x/p;if (B){data tmp=Calc(B,p,a,p_a);(A*=tmp.first)%=p_a;B+=tmp.second;}return data(A,B);}inline void Solve(){China::clear();cnt=0; cl(p); cl(a); cl(p_a);data ret,tmp;read(m); read(n); read(M);Machine(M);for (int i=1;i<=cnt;i++){ret=Calc(n,p[i],a[i],p_a[i]);tmp=Calc(m,p[i],a[i],p_a[i]);ret.second-=tmp.second;(ret.first*=Inv(tmp.first,p_a[i]))%=p_a[i];tmp=Calc(n-m,p[i],a[i],p_a[i]);ret.second-=tmp.second;(ret.first*=Inv(tmp.first,p_a[i]))%=p_a[i];(ret.first*=Pow(p[i],ret.second,p_a[i]))%=p_a[i];China::add(ret.first,p_a[i]);}ans=China::Solve();write(ans,'\n');}}int main(){ll Q,order;freopen("t.in","r",stdin);freopen("t.out","w",stdout);read(Q);Work3::Pre();while (Q--){read(order);if (order==1) Work1::Solve();if (order==2) Work2::Solve();if (order==3) Work3::Solve();}return 0;}