听说明天AU爷ysy要给我们讲数论
赶紧预习一下好了
①高斯消元
嗯这个就不多解释了…
题1 bzoj1013 球形空间产生器
嗯这题就是给你一个n维球体上的n+1个点，求圆心。
所以假设第i个点是(a[i][1],a[i][2]…a[i][n])
那么我们就是要找一个点(x[1]…x[n])使得
(a[i][1]-x[1])^2+(a[i][2]-x[2])^2+…+(a[i][n]-x[n])^2=r^2
那么我们就有这样n+1个式子：
(a[1][1]-x[1])^2+(a[1][2]-x[2])^2+…+(a[1][n]-x[n])^2=r^2
(a[2][1]-x[1])^2+(a[2][2]-x[2])^2+…+(a[2][n]-x[n])^2=r^2
...
(a[n+1][1]-x[1])^2+(a[n+1][2]-x[2])^2+…+(a[n+1][n]-x[n])^2=r^2
那么我们把这n+1个式子两两对减，就会得到一坨类似这样的式子：
$sum_{i=1}^n{(x[i]*2(a[2][i]-a[1][i]))}+sum_{i=1}^n(a[2][i]^2-a[1][i]^2)=0$
然后高斯消元即可.
#include <iostream> #include <stdio.h> #include <stdlib.h> #include <algorithm> #include <string.h> #include <vector> #include <math.h> #include <limits> #include <set> #include <map> using namespace std; typedef double ld; #define SZ 233 namespace Gauss { ld gd[SZ][SZ+1],ans[SZ]; bool isfree[SZ]; void gauss(int n) { int cur=1; for(int i=1;i<=n;i++) { isfree[i]=0; int tg=0; for(int j=cur;j<=n;j++) if(fabs(gd[j][i])>1e-6) {tg=j; break;} if(!tg) {isfree[i]=1; continue;} for(int j=1;j<=n;j++) swap(gd[cur][j],gd[tg][j]); for(int j=1;j<=n;j++) { if(j==cur||fabs(gd[j][i])<1e-6) continue; ld bs=gd[j][i]/gd[cur][i]; for(int k=1;k<=n+1;k++) gd[j][k]-=bs*gd[cur][k]; } ++cur; } for(int i=1;i<=n;i++) { if(isfree[i]) continue; for(int j=1;j<=n;j++) { if(fabs(gd[j][i])>1e-6) ans[i]=gd[j][n+1]/gd[j][i]; } } } } int n; ld as[SZ][SZ]; void pre() { scanf("%d",&n); for(int i=1;i<=n+1;i++) { for(int j=1;j<=n;j++) scanf("%lf",&as[i][j]); if(i==1) continue; for(int j=1;j<=n;j++) { Gauss::gd[i-1][j]=2*(as[i][j]-as[i-1][j]); Gauss::gd[i-1][n+1]+=as[i][j]*as[i][j]-as[i-1][j]*as[i-1][j]; } } } int main() { using namespace Gauss; pre(); gauss(n); for(int i=1;i<=n;i++) printf("%.3lf%c",ans[i],(i!=n)?' ':' '); }
②exgcd
这种傻逼玩意儿不用我说了吧
求一组x,y使得ax+by=1。如果(a,b)≠1就无解。
我们先求出bx'+(a%b)y'=1的解x'、y'。
那么设a=pb+q(0<=q<b)。所以bx'+qy'=1。
所以bx'+(a-pb)y'=1。
所以ay'+b(x'-py')=1。
即ay'+b(x'-(a/b)y')=1。
这样就求出了一组可行解。
代码稍后放出。
③BSGS
正常的BSGS：求一个最小的非负数x使得a^x=b（mod c），其中c为质数。
首先我们肯定知道x<c-1嘛
额我们设x=pg+m(0<=m<g)，这个g是啥等会儿再说。
那么a^(pg+m)=b（mod c）
我们稍微搞一下变成a^m=b*a^(-gp) （mod c）
然后我们发现枚举一下m，然后把a^m map一下，右边枚举一下g，exgcd(费马小定理)一下，在map里面查一下就行了，第一个查到了显然就是最小的解。
复杂度大概是O((g+c/g)logc)或者O(g+c/g)，所以我们令g=sqrt(c)就可以做到科学的O(sqrt(c)logc)或O(sqrt(c))了。
扩展的BSGS：（方老师写的写错了不要打我）
例二 bzoj2242 计算器
搞这种三合一的题啊，excited
似乎特判十分多的样子
//By zzq #include <iostream> #include <stdio.h> #include <stdlib.h> #include <algorithm> #include <string.h> #include <vector> #include <math.h> #include <limits> #include <set> #include <map> using namespace std; typedef long long ll; int T,k; ll qp(ll a,ll b,ll c) { if(b==0) return 1; ll hf=qp(a,b>>1,c); hf=hf*hf%c; if(b&1) hf=hf*a%c; return hf; } ll gcd(ll a,ll b) {return (b==0)?a:gcd(b,a%b);} //ax+by=1 void exgcd(ll a,ll b,ll& x,ll& y) { if(b==0) {x=1; y=0; return;} exgcd(b,a%b,x,y); ll p=x-a/b*y; x=y; y=p; } void q1(ll a,ll b,ll c) {printf("%lld ",qp(a%c,b,c));} #define orzhzwer {puts("Orz, I cannot find x!"); return;} void q2(ll a,ll b,ll c) { ll g=gcd(a,c); if(b%g) orzhzwer a/=g; b/=g; c/=g; ll x,y; exgcd(a,c,x,y); x=x*b%c; x=(x%c+c)%c; printf("%lld ",x); } map<ll,ll> mp; void q3(ll a,ll b,ll c) { a%=c; b%=c; if(!a&&!b) {puts("1"); return;} if(!a) orzhzwer mp.clear(); ll m=sqrt(c); ++m; ll cur=1; mp[1]=m+233; for(int i=1;i<m;i++) { cur=cur*a%c; if(!mp[cur]) mp[cur]=i; } ll tm=b%c,bz=qp(a%c,(c-2)*m%(c-1),c); for(int i=0;i<m;i++) { long long tmp=mp[tm]; if(tmp) { if(tmp==m+233) tmp=0; printf("%lld ",(i*m+tmp)%(c-1)); return; } tm=tm*bz%c; } orzhzwer } int main() { int T,k; scanf("%d%d",&T,&k); while(T--) { int a,b,c; scanf("%d%d%d",&a,&b,&c); if(k==1) q1(a,b,c); else if(k==2) q2(a,b,c); else q3(a,b,c); } }
④矩阵求逆
只要在高斯消元的时候旁边放一个单位矩阵，把原矩阵消元成单位矩阵，同时把单位矩阵也一起拿去变换就行了。注意取模的时候负数要判一下。
例3 bzoj4128 Matrix
BSGS+矩阵求逆的好(裸?)题
#include <iostream> #include <stdio.h> #include <stdlib.h> #include <algorithm> #include <string.h> #include <vector> #include <math.h> #include <limits> #include <set> #include <map> using namespace std; int n; typedef long long ll; ll p; struct mat{ll s[72][72];}A,B; ll qp(ll a,ll b) { if(b==0) return 1; ll hf=qp(a,b>>1); hf=hf*hf%p; if(b&1) hf=hf*a%p; return hf; } void readin(mat& x) { for(int i=1;i<=n;i++) { for(int j=1;j<=n;j++) scanf("%lld",&x.s[i][j]), x.s[i][j]%=p; } } mat dw() { mat a; for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) a.s[i][j]=(i==j); return a; } mat inv(mat x) { mat a=dw(); for(int i=1;i<=n;i++) { int tg=-1; for(int j=i;j<=n;j++) { if(x.s[j][i]) {tg=j; break;} } if(tg==-1) {printf("Fuck You"); exit(-1);} for(int j=1;j<=n;j++) swap(x.s[tg][j],x.s[i][j]), swap(a.s[tg][j],a.s[i][j]); ll cp=qp(x.s[i][i],p-2); for(int j=1;j<=n;j++) x.s[i][j]=(x.s[i][j]*cp)%p, a.s[i][j]=(a.s[i][j]*cp)%p; for(int j=1;j<=n;j++) { if(i==j) continue; ll bs=x.s[j][i]%p; for(int k=1;k<=n;k++) { x.s[j][k]=(x.s[j][k]-bs*x.s[i][k]%p)%p; x.s[j][k]=(x.s[j][k]+p)%p; a.s[j][k]=(a.s[j][k]-bs*a.s[i][k]%p)%p; a.s[j][k]=(a.s[j][k]+p)%p; } } } return a; } mat operator * (mat a,mat b) { mat ans; for(int i=1;i<=n;i++) { for(int j=1;j<=n;j++) { ll ap=0; for(int k=1;k<=n;k++) ap=(ap+a.s[i][k]*b.s[k][j]%p)%p; ans.s[i][j]=ap; } } return ans; } mat qp(mat a,ll cf) { mat ans=dw(); for(;cf;cf>>=1) { if(cf&1) ans=ans*a; a=a*a; } return ans; } const int hn=7; //强行卡常 ll gh(mat x) { ll hash=0; int cc=0; for(int i=1;i<=n;i++) { for(int j=1;j<=n;j++) { hash=hash*199999+x.s[i][j]; if((++cc)>=hn) return hash; } } return hash; } map<ll,int> mp; //这是我见过最难写的bsgs! void bsgs() { int g=sqrt(p); ++g; mp.clear(); mat cur=dw(); mp[gh(cur)]=g+233; for(int i=1;i<g;i++) { cur=cur*A; ll th=gh(cur); if(!mp[th]) mp[th]=i; } cur=B; mat bz=inv(qp(A,g)); for(int i=0;i<g;i++) { ll tmp=mp[gh(cur)]; if(tmp) { if(tmp==g+233) tmp=0; printf("%lld ",i*g+tmp); return; } cur=cur*bz; } } int main() { scanf("%d%lld",&n,&p); readin(A); readin(B); bsgs(); }
至于这样写好不好写…
下面这份是矩阵求逆的…上面这份是比较科学的做法…嗯我们比较一下代码长度…
那科学的做法是什么呢？
我们不妨再回顾一下bsgs的过程…
我们设x=pg+m(0<=m<g)，其中g≈sqrt(p)…
嗯…这样就得求逆元对吧
现在我们不妨换一个思路！我们设x=pg-m(0<=m<g)，其中g≈sqrt(p)+2…
这样原来a^x=b变成了a^(pg)=b*a^m，这样就不用写什么鬼逆元啦！
感觉我bsgs顿时白学了
比较科学的代码：
#include <iostream> #include <stdio.h> #include <stdlib.h> #include <algorithm> #include <string.h> #include <vector> #include <math.h> #include <limits> #include <set> #include <map> using namespace std; int n; typedef long long ll; ll p; struct mat{ll s[72][72];}A,B; ll qp(ll a,ll b) { if(b==0) return 1; ll hf=qp(a,b>>1); hf=hf*hf%p; if(b&1) hf=hf*a%p; return hf; } void readin(mat& x) { for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) scanf("%lld",&x.s[i][j]), x.s[i][j]%=p; } mat dw() { mat a; for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) a.s[i][j]=(i==j); return a; } mat operator * (mat a,mat b) { mat ans; for(int i=1;i<=n;i++) { for(int j=1;j<=n;j++) { ll ap=0; for(int k=1;k<=n;k++) ap=(ap+a.s[i][k]*b.s[k][j]%p)%p; ans.s[i][j]=ap; } } return ans; } mat qp(mat a,ll cf) { mat ans=dw(); for(;cf;cf>>=1) { if(cf&1) ans=ans*a; a=a*a; } return ans; } const int hn=7; //强行卡常 ll gh(mat x) { ll hash=0; int cc=0; for(int i=1;i<=n;i++) { for(int j=1;j<=n;j++) { hash=hash*199999+x.s[i][j]; if((++cc)>=hn) return hash; } } return hash; } map<ll,int> mp; //这是我见过最难写的bsgs! void bsgs() { int g=sqrt(p)+2; mp.clear(); mat cur=B; mp[gh(cur)]=g+233; for(int i=1;i<g;i++) { cur=cur*A; mp[gh(cur)]=i; } cur=dw(); mat bz=qp(A,g); for(int i=0;i<g;i++) { ll tmp=mp[gh(cur)]; if(tmp) { if(tmp==g+233) tmp=0; printf("%lld ",i*g-tmp); return; } cur=cur*bz; } } int main() { scanf("%d%lld",&n,&p); readin(A); readin(B); bsgs(); }