数学模板
1 #include<algorithm> 2 #include<iostream> 3 #include<cstdlib> 4 #include<cstring> 5 #include<cstdio> 6 #include<string> 7 #include<cmath> 8 #include<ctime> 9 #include<queue> 10 #include<stack> 11 #include<map> 12 #include<set> 13 #define rre(i,r,l) for(int i=(r);i>=(l);i--) 14 #define re(i,l,r) for(int i=(l);i<=(r);i++) 15 #define Clear(a,b) memset(a,b,sizeof(a)) 16 #define inout(x) printf("%d",(x)) 17 #define douin(x) scanf("%lf",&x) 18 #define strin(x) scanf("%s",(x)) 19 #define LLin(x) scanf("%lld",&x) 20 #define op operator 21 #define CSC main 22 typedef unsigned long long ULL; 23 typedef const int cint; 24 typedef long long LL; 25 using namespace std; 26 void inin(int &ret) 27 { 28 ret=0;int f=0;char ch=getchar(); 29 while(ch<'0'||ch>'9'){if(ch=='-')f=1;ch=getchar();} 30 while(ch>='0'&&ch<='9')ret*=10,ret+=ch-'0',ch=getchar(); 31 ret=f?-ret:ret; 32 } 33 namespace C 34 { 35 int c[1010][1010]; 36 void js1(int n)//????¨¨y?? 37 { 38 c[1][1]=1; 39 re(i,1,n) 40 { 41 c[i][1]=i; 42 re(j,2,i)c[i][j]=c[i-1][j-1]+c[i-1][j]; 43 } 44 } 45 void js2(int n,int k)//O(n)¦ÌY¨ª? 46 { 47 c[n][1]=n; 48 re(i,2,n)c[n][i]=c[n][i-1]*(n-i+1)/i; 49 } 50 } 51 namespace prime 52 { 53 int prime[100010],prime_tot; 54 bool bo[1000010]; 55 void shai(int limit)//??D?¨¦? 56 { 57 re(i,2,limit) 58 { 59 if(!bo[i])prime[++prime_tot]=i; 60 for(int j=1;j<=prime_tot&&i*prime[j]<=limit;j++) 61 { 62 bo[i*prime[j]]=1; 63 if(i%prime[j]==0)break; 64 } 65 } 66 } 67 bool pd(int x)//??¨ºy?D?¡§ 68 { 69 shai(sqrt(x)+1); 70 re(i,1,prime_tot)if(x%prime[i]==0)return 0; 71 return 1; 72 } 73 } 74 namespace nt 75 { 76 LL gcd(LL a,LL b){LL c;while(a%b)c=a%b,a=b,b=c;return b;} 77 void exgcd(LL a,LL b,LL &d,LL &x,LL &y)//¨¤??1?¡¤??¨¤?¦Ì? 78 { 79 if(!b){d=a;x=1,y=0;} 80 else exgcd(b,a%b,d,y,x),y-=x*(a/b); 81 } 82 LL qpow(LL a,LL b,LL mod)//?¨¬?¨´?Y 83 { 84 a%=mod;LL ret=1; 85 while(b) 86 { 87 if(b&1)ret=ret*a%mod; 88 b>>=1; 89 a=a*a%mod; 90 } 91 return ret; 92 } 93 LL PHI(LL x) 94 { 95 LL ret=x; 96 for(LL i=2;i*i<=x;i++)if(x%i==0) 97 { 98 ret=ret-ret/i; 99 while(x%i==0)x/=i; 100 } 101 if(x>1) 102 ret=ret-ret/x; 103 return ret; 104 } 105 int phi[1000010],prime[1000010],tot; 106 bool a[1000010]; 107 void shai(int limit) 108 { 109 phi[1]=1; 110 re(i,2,limit) 111 { 112 if(!a[i])prime[++tot]=i; 113 for(int j=1;j<=tot&&i*prime[j]<=limit;j++) 114 { 115 a[i*prime[j]]=1; 116 if(i%prime[j]==0) 117 { 118 phi[i*prime[j]]=phi[i]*prime[j]; 119 break; 120 } 121 else phi[i*prime[j]]=phi[i]*(prime[j]-1); 122 } 123 } 124 } 125 LL inv1(LL a,LL n)//???a(exgcd) 126 { 127 LL d,x,y; 128 exgcd(a,n,d,x,y); 129 return d==1?(x+n)%n:-1; 130 } 131 LL inv2(LL a,LL n)//???a(¡¤??¨ªD??¡§¨¤¨ª) 132 { 133 if(gcd(a,n)!=1)return -1; 134 LL zhi=PHI(n); 135 return qpow(a,zhi-1,n); 136 } 137 LL inv[1000010]; 138 void getinv(LL mod) 139 { 140 inv[1]=1; 141 re(i,2,1000000)inv[i]=(mod-mod/i)*inv[mod%i]%mod; 142 } 143 //x=a[i](mod m[i])(1<=i<=n) 144 LL china(int n,int *a,int *m)//?D1¨²¨º¡ê¨®¨¤?¡§¨¤¨ª 145 { 146 LL M=1,d,y,x=0; 147 re(i,1,n)M*=m[i]; 148 re(i,1,n) 149 { 150 LL w=M/m[i]; 151 exgcd(m[i],w,d,d,y); 152 x=(x+y*w*a[i])%M; 153 } 154 return (x+M)%M; 155 } 156 //?a?¡ê¡¤?3¨¬a^x=b(mod n),n?a?¨º¨ºy 157 int log_mod(int a,int b,int n)//BSGS 158 { 159 int m,v,e=1; 160 m=(int)sqrt(n+1); 161 v=inv2(qpow(a,m,n),n); 162 map<int,int>x; 163 x[1]=0; 164 re(i,1,m-1) 165 { 166 e=e*a%n; 167 if(!x.count(e))x[e]=i; 168 } 169 re(i,0,m-1) 170 { 171 if(x.count(b))return i*m+x[b]; 172 b=b*v%n; 173 } 174 return -1; 175 } 176 int prime_tot,mu[100010]; 177 bool bo[100010]; 178 void get_mu(int limit)//?a¡À¨¨?¨²?1o¡¥¨ºy 179 { 180 mu[1]=1; 181 re(i,2,limit) 182 { 183 if(!bo[i])prime[++prime_tot]=i,mu[i]=-1; 184 for(int j=1;j<=prime_tot&&i*prime[j]<=limit;j++) 185 { 186 bo[i*prime[j]]=1; 187 if(i%prime[j])mu[i*prime[j]]=-mu[i]; 188 else {mu[i*prime[j]]=0;break;} 189 } 190 } 191 } 192 } 193 namespace Matrix 194 { 195 struct mat 196 { 197 int a[111][111],r,c; 198 void clear(int rr,int cc,int x) 199 { 200 r=rr,c=cc; 201 re(i,1,r)re(j,1,c)if(i==j)a[i][j]=x; 202 else a[i][j]=0; 203 } 204 int* op [] (int x){return x[a];} 205 mat op + (mat &rhs) 206 { 207 mat ret=*this; 208 re(i,1,r)re(j,1,c)ret[i][j]+=rhs[i][j]; 209 return ret; 210 } 211 mat op - (mat &rhs) 212 { 213 mat ret=*this; 214 re(i,1,r)re(j,1,c)ret[i][j]-=rhs[i][j]; 215 return ret; 216 } 217 mat op * (int &rhs) 218 { 219 mat ret=*this; 220 re(i,1,r)re(j,1,c)ret[i][j]*=rhs; 221 return ret; 222 } 223 mat op / (int &rhs) 224 { 225 mat ret=*this; 226 re(i,1,r)re(j,1,c)ret[i][j]/=rhs; 227 return ret; 228 } 229 mat op * (mat &rhs) 230 { 231 mat ret;ret.clear(r,rhs.c,0); 232 re(i,1,r)re(j,1,rhs.c)re(k,1,c) 233 ret[i][j]+=a[i][k]*rhs[k][j]; 234 return ret; 235 } 236 mat qpow(int n) 237 { 238 mat ret,temp=*this;ret.clear(r,c,1); 239 while(n) 240 { 241 if(n&1)ret=ret*temp; 242 n>>=1; 243 temp=temp*temp; 244 } 245 return ret; 246 } 247 }; 248 } 249 namespace gauss 250 { 251 const double eps=1e-8; 252 int n; 253 double a[111][111],x[111]; 254 void xy1() 255 { 256 re(i,1,n) 257 { 258 if(abs(a[i][i])<eps) 259 re(j,i+1,n)if(abs(a[j][i])>eps) 260 { 261 re(k,i,n+1)swap(a[j][k],a[i][k]); 262 break; 263 } 264 re(j,i+1,n) 265 { 266 if(abs(a[j][i])<eps)continue; 267 double bi=a[i][i]/a[j][i]; 268 re(k,i,n+1)a[j][k]*=bi; 269 re(k,i,n+1)a[j][k]-=a[i][k]; 270 } 271 } 272 rre(i,n,1) 273 { 274 x[i]=a[i][n+1]/a[i][i]; 275 rre(j,i-1,1)a[j][n+1]-=a[j][i]*x[i],a[j][i]=0; 276 } 277 } 278 LL fu[111],ans[111],aa[111][111],mod; 279 void xy2() 280 { 281 re(i,0,n-1) 282 re(j,i+1,n-1) 283 { 284 LL xia=aa[j][i],shang=aa[i][i]; 285 re(k,i,n) 286 { 287 fu[k]=aa[i][k]*xia%mod; 288 aa[j][k]=aa[j][k]*shang%mod; 289 aa[j][k]=((aa[j][k]-fu[k])%mod+mod)%mod; 290 } 291 } 292 rre(i,n-1,0) 293 { 294 LL d,x,y; 295 using namespace nt; 296 exgcd(aa[i][i],mod,d,x,y); 297 aa[i][n]=aa[i][n]*x%mod,aa[i][i]=1; 298 rre(j,i-1,0)(aa[j][n]-=aa[j][i]*aa[i][n]%mod)%=mod,aa[j][i]=0; 299 } 300 } 301 } 302 int main() 303 { 304 return 0; 305 }