标签： 数学
[模板]ntt和fft
FFT：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 |
struct Complex { double real, image; Complex( double _real, double _image) { real = _real; image = _image; } Complex(){} }; Complex operator + (const Complex &c1, const Complex &c2) { return Complex(c1.real + c2.real, c1.image + c2.image); } Complex operator - (const Complex &c1, const Complex &c2) { return Complex(c1.real - c2.real, c1.image - c2.image); } Complex operator * (const Complex &c1, const Complex &c2) { return Complex(c1.real*c2.real - c1.image*c2.image, c1.real*c2.image + c1.image*c2.real); } int rev(int id, int len) { int ret = 0; for(int i = 0; (1 << i) < len; i++) { ret <<= 1; if(id & (1 << i)) ret |= 1; } return ret; } Complex A[1 << 19]; void FFT(Complex *a, int len, int DFT) { for(int i = 0; i < len; i++) A[rev(i, len)] = a[i]; for(int s = 1; (1 << s) <= len; s++) { int m = (1 << s); Complex wm = Complex(cos(DFT*2*PI/m), sin(DFT*2*PI/m)); for(int k = 0; k < len; k += m) { Complex w = Complex(1, 0); for(int j = 0; j < (m >> 1); j++) { Complex t = w*A[k + j + (m >> 1)]; Complex u = A[k + j]; A[k + j] = u + t; A[k + j + (m >> 1)] = u - t; w = w*wm; } } } if(DFT == -1) for(int i = 0; i < len; i++) A[i].real /= len, A[i].image /= len; for(int i = 0; i < len; i++) a[i] = A[i]; return; } |
NTT
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 |
#include <iostream> #include <cstdio> #include <cstring> #include <algorithm> #include <cmath> using namespace std; const int maxn=300020; const int P=880803841; const int G=26; int wn[25]; long long mul(long long x,long long y) { return (x*y-(long long)(x/(long double)P*y+1e-3)*P+P) % P; } int qpow(int x,int k,int p) { int ret=1; while(k) { if (k&1) ret=1LL*ret*x%p; k>>=1; x=1LL*x*x%p; } return ret; } void getwn() { for (int i=1;i<=21;i++) { int t=1<<i; wn[i]=qpow(G,(P-1)/t,P); } } void change(int *y,int len) { for (int i=1,j=len/2;i<len-1;i++) { if (i<j) swap(y[i],y[j]); int k=len/2; while(j>=k) { j-=k; k/=2; } j+=k; } } void NTT(int *y,int len,int on) { change(y,len); int id=0; for (int h=2;h<=len;h<<=1) { id++; for (int j=0;j<len;j+=h) { int w=1; for (int k=j;k<j+h/2;k++) { int u=y[k]; int t=1LL*y[k+h/2]*w%P; y[k]=u+t; if (y[k]>=P) y[k]-=P; y[k+h/2]=u-t+P; if (y[k+h/2]>=P) y[k+h/2]-=P; w=1LL*w*wn[id]%P; } } } if (on==-1) { for (int i=1;i<len/2;i++) { swap(y[i],y[len-i]); } int inv=qpow(len,P-2,P); for (int i=1;i<len;i++) { y[i]=1LL*y[i]*inv%P; } } } int main() { getwn(); return 0; } |
插值？插值！
这是一个神奇的东西，mark！
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 |
#include <cstdio> #include <cstring> #include <cmath> #include <algorithm> #include <vector> #include <string> #include <map> #include <set> #include <cassert> using namespace std; #define rep(i,a,n) for (int i=a;i<n;i++) #define per(i,a,n) for (int i=n-1;i>=a;i--) #define pb push_back #define mp make_pair #define all(x) (x).begin(),(x).end() #define fi first #define se second #define SZ(x) ((int)(x).size()) typedef vector<int> VI; typedef long long ll; typedef pair<int,int> PII; const ll mod=1000000007; ll powmod(ll a,ll b) {ll res=1;a%=mod; assert(b>=0); for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;} // head int _; ll n; namespace linear_seq { const int N=10010; ll res[N],base[N],_c[N],_md[N]; vector<ll> Md; void mul(ll *a,ll *b,ll k) { rep(i,0,k+k) _c[i]=0; rep(i,0,k) if (a[i]) rep(j,0,k) _c[i+j]=(_c[i+j]+a[i]*b[j])%mod; for (int i=k+k-1;i>=k;i--) if (_c[i]) rep(j,0,SZ(Md)) _c[i-k+Md[j]]=(_c[i-k+Md[j]]-_c[i]*_md[Md[j]])%mod; rep(i,0,k) a[i]=_c[i]; } int solve(ll n,VI a,VI b) { ll ans=0,pnt=0; ll k=SZ(a); assert(SZ(a)==SZ(b)); rep(i,0,k) _md[k-1-i]=-a[i];_md[k]=1; Md.clear(); rep(i,0,k) if (_md[i]!=0) Md.push_back(i); rep(i,0,k) res[i]=base[i]=0; res[0]=1; while ((1ll<<pnt)<=n) pnt++; for (int p=pnt;p>=0;p--) { mul(res,res,k); if ((n>>p)&1) { for (int i=k-1;i>=0;i--) res[i+1]=res[i];res[0]=0; rep(j,0,SZ(Md)) res[Md[j]]=(res[Md[j]]-res[k]*_md[Md[j]])%mod; } } rep(i,0,k) ans=(ans+res[i]*b[i])%mod; if (ans<0) ans+=mod; return ans; } VI BM(VI s) { VI C(1,1),B(1,1); int L=0,m=1,b=1; rep(n,0,SZ(s)) { ll d=0; rep(i,0,L+1) d=(d+(ll)C[i]*s[n-i])%mod; if (d==0) ++m; else if (2*L<=n) { VI T=C; ll c=mod-d*powmod(b,mod-2)%mod; while (SZ(C)<SZ(B)+m) C.pb(0); rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod; L=n+1-L; B=T; b=d; m=1; } else { ll c=mod-d*powmod(b,mod-2)%mod; while (SZ(C)<SZ(B)+m) C.pb(0); rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod; ++m; } } return C; } int gao(VI a,ll n) { VI c=BM(a); c.erase(c.begin()); rep(i,0,SZ(c)) c[i]=(mod-c[i])%mod; return solve(n,c,VI(a.begin(),a.begin()+SZ(c))); } }; int main() { while(~scanf("%lld",&n)){ printf("%d\n",linear_seq::gao(VI{1,1, 5, 11, 36, 95, 281, 781, 2245, 6336, 18061, 51205, 145601, 413351, 1174500, 3335651, 9475901, 26915305, 76455961, 217172736},n)); } } |
Read More »
插值求多项式系数
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 |
#include <iostream> #include <cstring> #include <algorithm> #include <cstring> #include <map> #include <vector> #include <cmath> using namespace std; const int maxm=100; const int mod=1e9+7; const double EPS=1e-8; typedef vector<double> vec; typedef vector<vec> mat; vec gauss_jordan(const mat&A,const vec&b) { int n=A.size(); mat B(n,vec(n+1)); for (int i=0;i<n;i++) { for (int j=0;j<n;j++) B[i][j]=A[i][j]; } for (int i=0;i<n;i++) B[i][n]=b[i]; for (int i=0;i<n;i++) { int pivot=i; for (int j=i;j<n;j++) { if (abs(B[j][i])>abs(B[pivot][i])) pivot=j; } swap(B[i],B[pivot]); if (abs(B[i][i])<EPS) return vec(); for (int j=i+1;j<=n;j++) B[i][j]/=B[i][i]; for (int j=0;j<n;j++) { if (i!=j) { for (int k=i+1;k<=n;k++) B[j][k]-=B[j][i]*B[i][k]; } } } vec x(n); for (int i=0;i<n;i++) { x[i]=B[i][n]; } return x; } mat cheng(mat m,mat mm) { int n=m.size(); mat res(n,vec(n)); for (int i=0;i<n;i++) { for (int j=0;j<n;j++) { for (int k=0;k<n;k++) { res[i][k]=res[i][k]+m[i][j]*mm[j][k]; } } } return res; } mat qpow(mat m,long long b) { int n=m.size(); mat ret(n,vec(n)); for (int i=0;i<n;i++) { ret[i][i]=1; } while(b) { if (b&1) ret=cheng(ret,m); m=cheng(m,m); b/=2; } return ret; } long long a[maxm]; long long sum[maxm];//插值的公式值 map<int,int> f; long long get(int n) { int crt=1; while(f[crt]) { crt++; } return crt; } const int maxn=50;//项数 mat A(maxn,vec(maxn)); vec b(maxn); int main() { a[1]=1; a[2]=2; f[1]=1; sum[1]=1; sum[2]=3; for (int i=3;i<maxm;i++) { if (i%2) { a[i]=2*a[i-1]; } else { int crt=get(i); a[i]=a[i-1]+crt; } for (int j=1;j<i;j++) { f[a[i]-a[j]]=1; } a[i]%=mod; sum[i]=(a[i]+sum[i-1]) % mod; } for (int i=maxn;i<=maxn+maxn-1;i++) { b[i-maxn]=sum[i]; } for (int i=0;i<maxn;i++) { for (int j=0;j<maxn-1;j++) { A[i][j]=sum[i+j+1]; } A[i][maxn-1]=1; } vec x=gauss_jordan(A,b); for (int i=0;i<maxn;i++) cout<<x[i]<<endl; mat B(maxn,vec(maxn)); for (int i=0;i<maxn-1;i++) { B[0][i]=x[maxn-2-i]; } B[0][maxn-1]=x[maxn-1]; for (int i=1;i<maxn;i++) { for (int j=0;j<maxn;j++) { B[i][j]=0; } if (i<maxn-1) B[i][i-1]=1; else B[i][i]=1; } for (int i=0;i<maxn;i++) { for (int j=0;j<maxn;j++) { cout<<B[i][j]<<" "; } cout<<endl; } for (int i=maxn;i<100;i++) { mat ret=qpow(B,i-maxn+1); double ans=0; for (int j=0;j<maxn-1;j++) { ans+=sum[maxn-1-j]*ret[0][j]; } ans+=ret[0][maxn-1]; cout<<ans<<" "<<sum[i]<<endl; } return 0; } |
数论模板（一）
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 |
//最大公约数 long long gcd(long long a,long long b) { return b==0?a:gcd(b,a%b); } //大整数相乘取模 long long multi(long long a,long long b,long long mod) { long long ret=0; while(b) { if (b&1) ret=(ret+a) %mod; a=(a<<1)%mod; b=b>>1; } return ret; } //快速幂1 long long quickpow(long long a,long long b,long long mod) { long long ret=1; while(b) { if (b&1) ret=multi(ret,a,mod); a=multi(a,a,mod); b>>=1; } return ret; } //快速幂2 long long quick_pow(long long a,long long b) { long long res=1; while(b>0) { if (b&1) res=(res*a)%p; a=(a*a)%p; b>>=1; } return res; } //欧拉函数 long long eular(long long n) { long long ret=1,i; for (i=2;i*i<=n;i++) { if (n%i==0) { n=n/i; ret*=i-1; while(n%i==0) { n=n/i; ret*=i; } } } if (n>1) ret*=n-1; return ret; } //预处理maxn范围内的欧拉函数 void euler_pre() { for (int i = 0; i <= maxn; ++i) euler[i] = i; for (int i = 2; i <= maxn; ++i) { if (euler[i] == i) { for (int j = i; j <= maxn; j += i) euler[j] = euler[j] / i * (i - 1); } } } //求一个数的所有约数，返回约数个数，约数存在c数组 int getdivision(int n) { int i; int length=0; for(i=1;i*i<n;i++){ if(n%i == 0){ c[length ++] = i; c[length ++] = n/i; } } if(i*i == n) c[length ++] = i; return length; } //矩阵（n×n）,p 为模数 struct Mat { long long mat[maxn][maxn]; long long n; Mat operator * (Mat b) { Mat c; c.n=n; memset(c.mat, 0, sizeof(c.mat)); int i, j, k; for(k = 0; k < n; ++k) { for(i = 0; i < n; ++i) { for(j = 0; j < n; ++j) { c.mat[i][j] =(c.mat[i][j]+(mat[i][k] * b.mat[k][j]) %p)%p; } } } return c; } }; //矩阵快速幂 Mat quick_mat(Mat a,long long k，long long p;) { Mat c; c.n=a.n; int i, j; for(i = 0; i < a.n; ++i) for(j = 0; j < a.n; ++j) c.mat[i][j] = (i == j); for(; k; k >>= 1) { if(k&1) c = c*a; a = a*a; } return c; } //高斯消元 typedef vector<double> vec; typedef vector<vec> mat; mat A; vec gauss_jordan(int n) { for (int c=0,r=0;c<n;++c,++r) { int pivot =r; for (int i=r+1;i<=n;++i) { if (abs(A[i][c])>abs(A[pivot][c])) pivot = i; } swap(A[r],A[pivot]); for (int j=n;j>=c;--j) A[r][j]/=A[r][c]; for (int i=0;i<=n;i++) { if (i!=r) { for (int k=c+1;k<=n;k++) { A[i][k]-=A[i][c]*A[r][k]; } } } } vec x(n); for (int i=0;i<n;i++) x[i]=A[i][n]; return x; } //组合数取模，Lucas定理，p为模数 LL C(LL n, LL m) { if(m > n) return 0; LL ans = 1; for(int i=1; i<=m; i++) { LL a = (n + i - m) % p; LL b = i % p; ans = ans * (a * quick_mod(b, p-2) % p) % p; } return ans; } LL Lucas(LL n, LL m) { if(m == 0) return 1; return C(n % p, m % p) * Lucas(n / p, m / p) % p; } //莫比乌斯函数 const int maxn=1000000; bool check[maxn+10]; int prime[maxn+10]; int mu[maxn+10]; void moblus() { memset(check,false,sizeof check); mu[1]=1; int tot=0; for (int i=2;i<=maxn;i++) { if (!check[i]) { prime[tot++]=i; mu[i]=-1; } for (int j=0;j<tot;j++) { if (i*prime[j]>maxn) break; check[i*prime[j]]=true; if (i%prime[j]==0) { mu[i*prime[j]]=0; break; } else { mu[i*prime[j]]=-mu[i]; } } } } //扩展欧几里德定理 typedef long long LL; void extend_Euclid(LL a, LL b, LL &x, LL &y) { if(b == 0) { x = 1; y = 0; return; } extend_Euclid(b, a % b,x, y); LL tmp = x; x = y; y = tmp - a / b * y; } //中国剩余定理 typedef long long LL; LL RemindChina(LL a[],LL m[],int k) { LL M = 1; LL ans = 0; for(int i=0; i<k; i++) M *= m[i]; for(int i=0; i<k; i++) { LL x, y; LL Mi = M / m[i]; extend_Euclid(Mi, m[i], x, y); ans = (ans + Mi * x * a[i]) % M; } if(ans < 0) ans += M; return ans; } |
Homework of Politics
题目： Homework of Politics After the Politics class, the teacher gives the student a homework.But Nobita and Takeshi don’t like this homework at all.So they decide to play a game, and the person who los…
Read More »