模版

ACM模板

 

用于简化模运算的定理

欧拉定理 a与m互质，a^phi(m)%m==1-->m为质数 a^(m-1)%m==1

威尔逊定理 p为素数 (p-1)!%p==p-1==-1

费马定理 p为素数，a和p互质，a^(p-1)%p==1

费马小定理 p为素数，a^p%p==a-->(a^p/a)%p==a/a-->a^(p-1)%p==1

a和m不互质时若m较小可以直接找a^x%m=1,即直接for循环从1开始找x

A^x % m = A^(x%phi(m)+phi(m)) % m (x >= phi(m))

用于化简(a^(b^(c^...))),先算最上面的次方,然后递归回去

 

换底公式 log(a)(b)=log(c)(b)/log(c)(a)，当结果超过浮点数范围时用log算,负数没有对数

 #pragma comment(linker, "/STACK:102400000,102400000")//手动开大栈区
多次模规律 x=x*x%p;在一定次数后x会不再改变，即x*x%p==x

用scanf(“%[^\n]”,s);时记得初始化s[0]=0;否则输入有空行时会WA,且输入空行时返回值为0

a开n次根:x=pow(a,1.0/n)

处理误差时同时判x*x==a || (x+1)*(x+1)==a (间隔为2x+1,所以基本能处理误差)

 

树上某个点的子树的所有信息为dfs时进入这个点时已经记录的信息和离开这个点时已经记录的信息的差集

//multiset中 s.erase(2)会一次删除所有的 2，若要只删一个应用 s.erase(s.find(2))--易错

多个样例输入 T较大且满足从小到大递推时可以将询问离线出来处理

abs(x)在没有 using namespace std;的情况下超过 long long 会出错

gcd计数相关:莫比乌斯反演,容斥原理,欧拉函数

递归超过 1e5层左右就会爆栈 RE,将递归改成手写栈

 

 

 1 C(n,k+1)=C(n,k)*(n-k)/(k+1)
 2 
 3 //O(n)递推
 4 
 5 C(n,k)=n!/(k!*(n-k!))=n*(n-1)*...*(n-k+1){k个相乘}/k! (k<=n)
 6 
 7 C(n,k)=C(n,n-k)
 8 
 9 C(n,k)=C(n-1,k)+C(n-1,k-1)
10 
11 {杨辉等式，杨辉三角形的对称性}
12 
13 C(n,0)+C(n,1)+...+C(n,n-1)+C(n,n)=2^n {二项式定理}
14 
15 C(n,0)+C(n+1,1)+C(n+2,2)+.+C(n+r,r)=C(n+r+1,r)//把 C(n,0)看出 C(n+1,0)不断往后合并
16 
17 可重复组合（隔板法）
18 
19 从 n个元素中可重复地选 k个，n个盒子放 k个相同的球，选法为
20 
21 C(n+k-1,n-1)=C(n+k-1,k)
22 
23 有重复元素排列
24 
25 n个元素中有 n1个相同，n2个相同。。。
26 
27 n!/(n1!*n2!*...*nk!)
28 
29 圆排列
30 
31 从 n个不同元素中选 k个不分首位围成一个圆圈，方案数为
32 
33 A(n,k)/k 如果 k==n，则有 n!/n=(n-1)!种
34 
35 条件排列
36 
37 将 1,2,3,4,5五个数进行排列，要求 4排在 2前面
38 
39 方法 1：先排 135,3!种,再插 2,4，先插 4，有四种插法，对于每种插法，2有 4,3,2,1种插法，所以种数=3!*10=60
40 
41 方法 2：不考虑 2,4的全排列，共有 5！种，2和 4的位置先后分布情况各半，所以种数=5!/2=60
42 
43 错位排列
44 
45 {a1,a2,...,an}，ai不放在第 i个位置的种数
46 
47 Dn=n!*(1-1/1!+1/2!-1/3!+1/4!-...!(-1)^n/n!)
48 
49 f(n)=(n-1)*(f(n-1)+f(n-2)) 递推式
组合公式
 1 一、球相同，盒子相同，且盒子不能空
 2 
 3 结论 ｎ个相同的球放入ｍ个相同的盒子（n≥m），不能有空盒时的放法种数等于ｎ分解为ｍ个数的和的种数.
 4 
 5 二、球相同，盒子相同，且盒子可以空
 6 
 7 结论 ｎ个相同的球放入ｍ个相同的盒子（n≥m），可以有空盒时的放法种数等于将ｎ分解为ｍ个、(ｍ－1)个、
 8 
 9 (ｍ－2)个、…、2 个、1 个数的和的所有种数之和.
10 
11 三、球相同，盒子不同，且盒子不能空
12 
13 结论 ｎ个相同的球放入ｍ个不同的盒子中（n≥m），不能有空盒的放法数 C(n-1,m-1).
14 
15 四、球相同，盒子不同，且盒子可以空
16 
17 结论 ｎ个相同的球放入ｍ个不同的盒子中（n≥m），可以有空盒的放法数 C(n+m-1,m-1).
18 
19 五、球不同，盒子相同，且盒子不能空
20 
21 .
22 
23 递推式:S(n,m)=S(n-1,m-1)+m*S(n-1,m)
24 
25 结论 ｎ个不同的球放入ｍ个相同的盒子中（n≥m），不能有空盒的放法种数等于ｎ个不同的球分成ｍ堆的种
26 
27 数.
28 
29 六、球不同，盒子相同，且盒子可以空
30 
31 结论 ｎ个不同的球放入ｍ个相同的盒子中（n≥m），可以有空盒的放法种数等于将ｎ个不同的球分成ｍ堆、(ｍ
32 
33 －1)堆、(ｍ－2)堆、…、2 堆、1 堆的所有种数之和.
34 
35 七、球不同，盒子不同，且盒子不能空
36 
37 结论 ｎ个不同的球放入ｍ个不同的盒子中，不能有空盒的放法种数等于ｎ个不同的球分成ｍ堆的种数乘以ｍ！.
38 
39 八、球不同，盒子不同，且盒子可以空
40 
41 结论 ｎ个不同的球放入ｍ个不同的盒子中（n≥m），可以有空盒的放法种数等于ｍ^ｎ种.
球盒问题 

将 n个不同的球和 m个不同球按原来顺序”穿插”起来,即将 n个球看成一样的,m个也是一样的,然后进行排列,有

(n+m)!/(n!m!)中插法,即将 n个球按原来顺序放入(m+1)个盒子

 1 前几项:1 1 2 5 14 42 132 429 1430 4862 16796 58786 2080127429002674440969484535357670129644790
 2 
 3 1.递归公式
 4 
 5 h(n) = h(n-1)*(4*n-2) / (n+1);
 6 
 7 2.递归公式
 8 
 9 f(n)=f(0)*f(n-1)+f(1)*f(n-2)+...+f(n-2)*f(1)+f(n-1)*f(0){累加 f(i)*f(n-i-1),i从 0到 n-1}
10 
11 3.组合公式 1
12 
13 f(n)=C(2*n,n)/(n+1)
14 
15 4.组合公式 2
16 
17 f(n)=C(2*n,n)-C(2*n,n-1)
18 
19 用于
20 
21 1 二叉数形态的计数
22 
23 2 火车进站后的出站顺序种数,乘法加括号的运算顺序种数
24 
25 3 将凸 n边形用 n-3条不相交的对角线分成 n-2个互相没有重叠的三角形的分法
26 
27 4
28 
29 n个 A和 n个 B排成一排，在任何位置 B的个数不超过 A的个数的排列种数
Catalan数(栈相关计数)
 1 假设 x为奇数，y为偶数，则 z为奇数，2z与 2x的最大公因数为 2，2z和 2x可分别写作
 2 
 3 2z = (z + x) + (z - x)
 4 
 5 2x = (z + x) - (z - x)
 6 
 7 那么跟据最大公因数性质，z + x和 z - x的最大公因数也为 2，又因为：
 8 
 9 (z + x)(z - x) = y^2，两边同除以 4得：
10 
11 ((z + x) / 2)((z - x) / 2) = (y / 2)^2
12 
13 故可令：
14 
15 z + x = 2m^2, z - x = 2n^2
16 
17 其中 z = m + n, x = m - n（m与 n互质）
18 
19 则有：
20 
21 y2 = z^2 - x^2 = (2m)^2 (2n)^2 = 4m^2 n^2
22 
23 即 y = 2mn。
24 
25 综上所述，可得到下式：
26 
27 x = m^2 - n^2, y = 2mn, z = m^2 + n^2. (m, n为任意自然数)
//勾股数性质

#ifdef ONLINE_JUDGE

#else

freopen("in","r",stdin);

freopen("out","w",stdout);

#endif

 1 struct node
 2 {
 3 int d,i;
 4 }d[110];
 5 struct cmp
 6 {
 7 bool operator ()(const node &a, const node &b)
 8 {
 9 return a.d<b.d || a.d==b.d && a.i>b.i;
10 }
11 };
12 priority_queue<node, vector<node>, cmp>q;
13 priority_queue<int, vector<int>, less<int> >q;
14 priority_queue<int, vector<int>, greater<int> >q;
//优先队列重载 cmp
 1 template <class T>  
 2 bool scanff(T &ret)
 3 { //Faster Input:short,int,long long,double,float
 4     char c; int sgn; T bit=0.1;  
 5     if(c=getchar(),c==EOF) return 0;  
 6     while(c!='-'&&c!='.'&&(c<'0'||c>'9')) c=getchar();  
 7     sgn=(c=='-')?-1:1;  
 8     ret=(c=='-')?0:(c-'0');  
 9     while(c=getchar(),c>='0'&&c<='9') ret=ret*10+(c-'0');  
10     if(c==' '||c=='\n') { ret*=sgn; return 1; }  
11     while(c=getchar(),c>='0'&&c<='9') ret+=(c-'0')*bit,bit/=10;  
12     ret*=sgn;  
13     return 1;  
14 }
fast input
 1 template <class T>
 2 void printff(T ret)
 3 { //Faster Output:short,int,long long
 4     char s[22];
 5     if(ret<0) { printf("-");ret=-ret; }
 6     int len=0;
 7     while(ret) { s[++len]=ret%10+'0';ret/=10; }
 8     if(!len) s[++len]='0';
 9     while(len) printf("%c",s[len--]);
10 }
fast output 
 1 void phi_table(int *phi,int *f,int *p,int &pn)
 2 {
 3     phi[1]=1;pn=0;
 4     int k;
 5     for(int i=2;i<N;i++)
 6     {
 7         if(!f[i]) { p[pn++]=f[i]=i;phi[i]=i-1; }
 8         for(int j=0;j<pn&&(k=p[j]*i)<N;j++)
 9         {
10             f[k]=p[j];
11             if(f[i]==p[j]) { phi[k]=phi[i]*p[j]; break; }
12             else phi[k]=phi[i]*(p[j]-1);
13         }
14     }
15     //for(int i=2;i<=N;i++) phi[i]+=phi[i-1];//累加
16 }
//欧拉函数线性筛法，phi为欧拉数组，f为素数标记数组，p为素数数组O(n) 
1 void fa(int n,int &dn,int*d)        
2 {
3     dn=0;
4     int m=sqrt(n)+0.5;
5     for(int i=1;i<=m;i++) if(n%i==0) { d[dn++]=i;d[dn++]=n/i; }
6     if(m*m==n) dn--;
7     sort(d,d+dn);
8 }
//找n所有约数并从小到大存到d数组O(n^(1/2))
 1 void fa(int n,int &pn,int *p)
 2 {
 3     pn=0;
 4     int m=sqrt(n)+0.5;
 5     for(int i=2;i<=m;i++) if(n%i==0)
 6     {
 7         p[pn]=1;
 8         while(n%i==0)
 9         {
10             p[pn]*=i;
11             n/=i;
12         }
13         pn++;
14     }
15     if(n!=1) p[pn++]=n;
16 }
//找n所有质因子并将素数从小到大存最高幂次到d数组O(n^(1/2))
 1 void Prin()            //O(nlogn)
 2 {
 3     memset(prin,0,sizeof(prin));
 4     for(int i=2;i<N;i++) if(!prin[i])
 5         for(int j=i;j<N;j+=i)
 6         {
 7             int t=j;
 8             while(t%i==0)
 9             {
10                 prin[j]++;
11                 t/=i;
12             }
13         }
14 }
//筛出1到N中每个数质因子的个数,如24=2*2*2*3,则prin[24]=4,也可以改为筛质因子
 1 template<class T>
 2 T gcd(T a,T b)
 3 {
 4     while(b)
 5     {
 6         T r=b;
 7         b=a%b;
 8         a=r;
 9     }
10     return a;
11 }
//非递归gcd O(logn)
 1 LL slove(LL a,LL b)
 2 {
 3     LL g = gcd(b/a,a);  
 4     while(g!=1)
 5     {
 6         a = a/g;  
 7         g = gcd(b/a, a);  
 8         }
 9     return b/a;
10 }
//返回最小的ans满足lcm(a,ans)=b
 1 LL quick_mod(LL t,LL n)
 2 {
 3     LL ans=1;
 4     while(n)
 5     {
 6         if(n&1) ans=ans*t%mod;
 7         t=t*t%mod;
 8         n>>=1;
 9     }
10     return ans;
11 }
//t^n%mod二分快速幂O(log(n))

//欧拉函数为 phi(n)=小于 n且与 n互质的数的个数,这些数的和为 n*phi(n)/2

证明: 设 共有 m 个与 n 互质的数, m 为 n 的欧拉函数, m=φ(n)；

设 ans=∑i=0,m−1 ai…(1) ,gcd(ai,n)=1&& ai<n;

可知 gcd(ai,n−ai)=1;

则 ans=∑i=0,m−1 (n−ai)…(2) , gcd(ai,n−ai)=1 && n−ai<n;

(1) + (2) 得 2∗ans=n∗m, 即 ans=n∗φ(n)/2;

欧拉函数引理

1.p为素数 phi(p)=p-1

2.p为素数 phi(p^a)=(p-1)*p^(a-1)=p^a*(1-1/p)

3.a,b互质 phi(a*b)=phi(a)*phi(b)

n=p1^a1*p2^a2*...*pk^ak

phi(n)=p1^a1*(1-1/p1)*p2^a2*(1-1/p2)*...*pk^ak*(1-1/pk)=n*(1-1/p1)*(1-1/p2)*...*(1-1/pk)

 1 int euler(int n)
 2 {
 3     int ans = n;
 4     for (int i = 2; i * i <= n; i++)
 5         if (n % i == 0)
 6         {
 7             ans = ans / i * (i - 1);
 8             while (n % i == 0)
 9                 n /= i;
10         }
11     if (n > 1)
12         ans = ans / n * (n - 1);
13     return ans;
14 }
//O(sqrt(n))求单个数的欧拉函数,phi[n]=n/pi*(pi-1) 
 1 void phi_table(int *phi, int *f, int *p, int &pn)
 2 {
 3     phi[1] = 1;
 4     pn = 0;
 5     int k;
 6     for (int i = 2; i < N; i++)
 7     {
 8         if (!f[i])
 9         {
10             p[pn++] = f[i] = i;
11             phi[i] = i - 1;
12         }
13         for (int j = 0; j < pn && (k = p[j] * i) < N; j++)
14         {
15             f[k] = p[j];
16             if (f[i] == p[j])
17             {
18                 phi[k] = phi[i] * p[j];
19                 break;
20             }
21             else
22                 phi[k] = phi[i] * (p[j] - 1);
23         }
24     }
25     //for(int i=2;i<=N;i++) phi[i]+=phi[i-1];//累加
26 }
//欧拉函数线性筛法，phi为欧拉数组，f为素数标记数组，p为素数数组 O(n)
 1 void prime(int &pn, int *p, bool *vis)
 2 {
 3     pn = 0;
 4     for (int i = 2; i < N; i++)
 5     {
 6         if (!vis[i])
 7             p[pn++] = i;
 8         for (int j = 0; j < pn && i * p[j] < N; j++)
 9         {
10             vis[i * p[j]] = 1;
11             if (i % p[j] == 0)
12                 break;
13         }
14     }
15 }
//素数线性筛法，p为素数数组，vis为素数标记数组 O(n)
 1 int vis[N / 31], pn, p[N / 10];
 2 inline bool judge(int x) { return vis[x >> 5] & (1 << (x & 31)); }
 3 inline void setbit(int x) { vis[x >> 5] |= 1 << (x & 31); }
 4 void init()
 5 {
 6     pn = 0;
 7     memset(vis, 0, sizeof(vis));
 8     for (int i = 2; i < M; i++)
 9     {
10         if (!judge(i))
11             p[pn++] = i;
12         for (int j = 0; p[j] * i < M; j++)
13         {
14             setbit(i * p[j]);
15             if (i % p[j] == 0)
16                 break;
17         }
18     }
19 }
//位压素数筛,省内存
 1 int prime(int a, int b)
 2 {
 3     memset(vis, 1, sizeof(vis));
 4     for (int i = 0; i < pn && (LL)p[i] * p[i] <= b; i++)
 5     {
 6         int x = a / p[i] + (a % p[i] != 0);
 7         if (x == 1)
 8             x++;
 9         for (LL j = (LL)x * p[i]; j <= b; j += p[i])
10             vis[j - a] = 0; //,printf("%d\n",j);
11     }
12     int ans = 0;
13     for (LL i = 0; i + a <= b; i++)
14         ans += vis[i];
15     if (a == 1)
16         ans--;
17     return ans;
18 }
//区间素数筛
 1 LL cal(LL p, LL n)
 2 {
 3     p %= mod;
 4     if (p == 1)
 5         return (n + 1) % mod;
 6     LL inv = quick_mod(p - 1, mod - 2);
 7     LL t = quick_mod(p, n);
 8     if (t * p % mod != 1)
 9         return (t * p % mod - 1) * inv % mod;
10     return (t - 1) * inv % mod * t % mod;
11 }
//所有约数的和=(p[i]^(e[i]+1) - 1)/(p[i]-1) ,0<=i<pn,累乘 //计算(p^0+p^1+...+p^n)%mod,p==1时和 p^(n+1)%mod==1时会出错
 1 int log_mod(int a, int b)
 2 {
 3     int m, v, e = 1;
 4     m = (int)sqrt(mod + 0.5);
 5     v = Inv(quick_mod(a, m));
 6     map<int, int> x;
 7     x[1] = 0;
 8     for (int i = 1; i < m; i++)
 9     {
10         e = (LL)e * a % mod;
11         if (!x.count(e))
12             x[e] = i;
13     }
14 
15     for (int i = 0; i < m; i++)
16     {
17         if (x.count(b))
18             return i * m + x[b];
19         b = (LL)b * v % mod;
20     }
21     return -1;
22 }
//离散对数,解 a^x%mod=b,返回最小的 
 1 void gcd(LL a, LL b, LL &d, LL &x, LL &y)
 2 {
 3     if (!b)
 4     {
 5         d = a;
 6         x = 1;
 7         y = 0;
 8     }
 9     else
10     {
11         gcd(b, a % b, d, y, x);
12         y -= a / b * x;
13     }
14 }
//拓展欧几里得求 ax+by=g=gcd(a,b)的解(x,y),通解为(x+kb/g,y+ka/g) //则 ax+by=c(g|c) 的解为(xc/g,yc/g)
 1 LL equation(LL a, LL b, LL r, LL &x, LL &y)
 2 {
 3     LL g;
 4     gcd(a, b, g, x, y); //求出 ax+by==g
 5     if (r % g)
 6         return 0; //无解
 7     x *= r / g;
 8     y *= r / g;
 9     return g > 0 ? g : -g; //返回 gcd
10 }
//解方程 ax+by==r,r==0时按需特殊处理 通解(x+k*b/g,y-k*a/g);
 1 int a, b, c, x0, x2, y0, y2;
 2 LL solve(LL a, LL b, LL c, LL x0, LL x2, LL y0, LL y2)
 3 {
 4     if (a < 0 && b < 0)
 5         a = -a, b = -b, c = -c;
 6     if (!b || a < 0 && b > 0)
 7         swap(a, b), swap(x0, y0), swap(x2, y2);
 8     if (!a && !b)
 9     {
10         if (c)
11             return 0;
12         return (x2 - x0 + 1LL) * (y2 - y0 + 1LL);
13     }
14     if (!a)
15     {
16         if (c % b || y0 > -c / b || y2 < -c / b)
17             return 0;
18         return x2 - x0 + 1;
19     }
20     LL g, x, y;
21     g = equation(a, b, -c, x, y);
22     if (!g)
23         return 0;
24 
25     LL aa = abs(a / g), bb = abs(b / g);
26     x0 += ((x % bb + bb) % bb - (x0 % bb + bb) % bb + bb) % bb;
27     x2 -= ((x2 % bb + bb) % bb - (x % bb + bb) % bb + bb) % bb;
28     y0 += ((y % aa + aa) % aa - (y0 % aa + aa) % aa + aa) % aa;
29     y2 -= ((y2 % aa + aa) % aa - (y % aa + aa) % aa + aa) % aa;
30     if (x0 > x2 || y0 > y2)
31         return 0;
32     x0 = (-c - a * x0) / b;
33     //a*x0注意溢出
34     x2 = (-c - a * x2) / b;
35     if (b > 0)
36         swap(x0, x2);
37     x0 = max(x0, y0);
38     x2 = min(x2, y2);
39     if (x0 > x2)
40         return 0;
41     return (x2 - x0) / aa + 1;
42 }
//ax+by+c=0 在 x0<=x<=x2 y0<=y<=y2范围内的整数解个数
 1 LL multi_equation(int n, int *r, int *m)
 2 {
 3     LL x, y, g, lcm, ans;
 4     lcm = m[0];
 5     ans = r[0];
 6     for (int i = 1; i < n; i++)
 7     {
 8         g = equation(lcm, -m[i], r[i] - ans, x, y);
 9         if (!g)
10             return -1; //无解
11         ans += x * lcm;
12         lcm *= m[i] / g;
13         ans %= lcm;
14         //注意溢出
15     }
16     if (ans % lcm == 0)
17         return lcm; //结果为 0时答案为 lcm
18     return (ans + lcm) % lcm;
19 }
//类似中国剩余定理,可处理 m数组不互质的情况
1 LL Inv(LL a, LL n) //O(logn)
2 {
3     LL g, x, y;
4     gcd(a, n, g, x, y);
5     return g == 1 ? (x + n) % n : -1;
6 }
//逆元求(a^b-c)/d%p时当 a^b-c==0时(a^b-c)*inv(d)==0会出错,需特殊处理
//逆元,a/b%mod=a*b^-1%mod，f(x)=inv(x,mod)或 f(x)=quick_mod(x,mod-2,mod)，mod需是质数
//若 a,m互质,r=a/b%m=a*inv(b)%m,否则 r=a/b%m=a%(m*b)/b 可以算 a%(m*b)但算不了 a/b时用
//mod为素数线性求逆元:inv(x)=-(mod/x)*inv(mod%x)=(mod– mod/x)*inv(mod%x)%mod
 1 LL quick_mod(LL t, LL n)
 2 {
 3     t %= mod;
 4     LL ans = 1;
 5     while (n)
 6     {
 7         if (n & 1)
 8             ans = ans * t % mod;
 9         t = t * t % mod;
10         n >>= 1;
11     }
12     return ans;
13 }
//t^n%mod二分快速幂 O(log(n))
 1 LL fac[N];
 2 LL Inv(LL x)
 3 {
 4     return quick_mod(x, mod - 2);
 5 }
 6 void init()
 7 {
 8     //要先执行
 9     fac[0] = 1;
10     for (int i = 1; i < N; i++)
11         fac[i] = fac[i - 1] * i % mod;
12     //阶乘模
13 }
14 //若 mod不为质数且 n*m较小可以考虑打 n*m的 C[n][m]二维表
15 LL Cnm(int n, int m)
16 {
17     if (n < m || n < 0 || m < 0)
18         return 0;
19     return fac[n] * Inv(fac[m] * fac[n - m] % mod) % mod;
20 }
//求 n,m小于 N的组合数 C(n,m)%mod O(nlogn)
 1 LL Lucas(LL n, LL m)
 2 {
 3     if (n < m || n < 0 || m < 0)
 4         return 0;
 5     LL ans = 1;
 6     while (n && m)
 7     {
 8         ans = ans * Cnm(n % mod, m % mod) % mod; //这儿 Cnm的求逆不可先处理，会超时
 9         n /= mod;
10         m /= mod;
11     }
12     return ans;
13 }
14 LL Lucas(LL n, LL m)
15 //递归
16 {
17     if (m == 0)
18         return 1;
19     else
20         return (Cnm(n % mod, m % mod) * Lucas(n / mod, m / mod)) % mod;
21 }
/*首先给出这个 Lucas定理：
A、B是非负整数，p是质数。AB写成 p进制：A=a[n]a[n-1]...a[0]，B=b[n]b[n-1]...b[0]。
则组合数 C(A,B)与 C(a[n],b[n])*C(a[n-1],b[n-1])*...*C(a[0],b[0]) mod p同余
即：Lucas(n,m,p)=c(n%p,m%p)*Lucas(n/p,m/p,p) */
//求 C(n,m)%mod mod最大为 10^5。a,b可以很大！
  1 //对每个数分解质因子和次幂
  2 int pn, p[N], vis[N], d[N][10], e[N][10], dn[N];
  3 void init()
  4 {
  5     pn = 1;
  6     memset(vis, 0, sizeof(vis));
  7     memset(e, 0, sizeof(e));
  8     memset(dn, 0, sizeof(dn));
  9     for (int i = 2; i < N; i++)
 10         if (!vis[i])
 11         {
 12             p[pn] = i;
 13             for (int j = i; j < N; j += i)
 14             {
 15                 vis[j] = 1;
 16                 int &n = dn[j];
 17                 d[j][n] = i;
 18                 int x = j;
 19                 while (x % i == 0)
 20                     x /= i, e[j][n]++;
 21                 n++;
 22             }
 23             vis[i] = pn++;
 24         }
 25 }
 26 
 27 int c[N];
 28 inline void add(int n, int x)
 29 {
 30     if (!x)
 31         return;
 32     for (int i = 0; i < dn[n]; i++)
 33     {
 34         int id = vis[d[n][i]];
 35         c[id] += x * e[n][i];
 36     }
 37 }
 38 //x个组合数相乘求模,mod可不为质数,n不能太大,复杂度 x*log(n)+n*log(p)
 39 int bit[N];
 40 inline int sum(int x)
 41 {
 42     int res = 0;
 43     while (x > 0)
 44         res += bit[x], x -= x & -x;
 45     return res;
 46 }
 47 inline void update(int x, int y, int d)
 48 {
 49     if (!d)
 50         return;
 51     while (x <= n)
 52         bit[x] += d, x += x & -x;
 53     y++;
 54     while (y <= n)
 55         bit[y] -= d, y += y & -y;
 56 }
 57 int vv[n];
 58 inline void Cnm(int n, int m)
 59 {
 60     vv[n]++;
 61     vv[m]--;
 62     vv[n - m]--;
 63 }
 64 LL cal()
 65 {
 66     LL ans = 1;
 67     for (int i = 2; i <= n; i++)
 68         update(1, i, vv[i]);
 69     for (int i = 2; i <= n; i++)
 70         add(i, sum(i));
 71     for (int i = 1; i < pn && p[i] <= n; i++)
 72         if (c[i])
 73             ans = ans * quick_mod(p[i], c[i]) % mod;
 74     return ans;
 75 }
 76 
 77 //单个组合数取 mod,mod可不为质数 复杂度 n*log(p)
 78 LL Cnm(int n, int m)
 79 {
 80     if (n < m || n < 0 || m < 0)
 81         return 0;
 82     m = min(m, n - m);
 83     LL ans = 1;
 84     memset(c, 0, sizeof(c));
 85     for (int i = n, j = m; j > 0; i--, j--)
 86     {
 87         add(i, 1);
 88         add(j, -1);
 89     }
 90     for (int i = 1; i < pn; i++)
 91         if (c[i])
 92             ans = ans * quick_mod(p[i], c[i], mod) % mod;
 93     return ans;
 94 }
 95 int c[N], vv[N];
 96 inline void add(int n, int x)
 97 {
 98     if (!x)
 99         return;
100     for (int i = 0; i < dn[n]; i++)
101     {
102         int id = vis[d[n][i]];
103         c[id] += x * e[n][i];
104     }
105 }
大组合数求模
 1 int c[N][N], fac[N], part[N][N];
 2 void init()
 3 {
 4     memset(c, 0, sizeof(c));
 5     c[0][0] = fac[0] = 1;
 6     for (int i = 1; i < N; i++) //杨辉三角 Cnm递推
 7     {
 8         fac[i] = (LL)fac[i - 1] * i % mod;
 9         c[i][0] = 1;
10         for (int j = 1; j <= i; j++)
11             c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % mod;
12     }
13     //集合划分,part[i][j]表示 i个不同的元素划分成 j个相同组(i个不同球放到 j个相同盒子,不能有空盒,
14     //可以有空盒的为前缀和),递推式 S(n,m)=S(n-1,m-1)+m*S(n-1,m)
15     memset(part, 0, sizeof(part));
16     part[1][1] = 1;
17     for (int i = 2; i < N; i++)
18         for (int j = 1; j <= i; j++)
19             part[i][j] = (part[i - 1][j] * j + part[i - 1][j - 1]) % mod;
20 }
//集合划分,杨辉三角递推
 1 typedef unsigned int UI;
 2 const int HASH=100007;
 3 char e[N][20],w[N][20];
 4 struct Hash
 5 {
 6     int *head, *next, n;
 7     Hash()
 8     {
 9         head = new int[HASH];
10         next = new int[N];
11         n = 1;
12         memset(head, 0, sizeof(head));
13     }
14     ~Hash()
15     {
16         delete[] head;
17         delete[] next;
18     }
19     UI hash(char *s)
20     {
21         UI ans = 0;
22         for (int i = 0; s[i]; i++)
23             ans = s[i] + ans * 101;
24         return ans % HASH;
25     }
26     int hash(int s) { return s % HASH; }
27     bool cmp(char *a, char *b) { return !strcmp(a, b); }
28     bool cmp(int a, int b) { return a == b; }
29     bool insert(int s)
30     {
31         int h = hash(e[s]);
32         int u = head[h];
33         while (u)
34         {
35             //当插入元素中不会有重复时可不要此循环
36             if (cmp(e[u], e[s]))
37                 return 0;
38             u = next[u];
39         }
40         next[s] = head[h];
41         head[h] = s;
42         return 1;
43     }
44     int find(char *s)
45     {
46         int h = hash(s);
47         int u = head[h];
48         while (u)
49         {
50             if (cmp(e[u], s))
51                 return u;
52             u = next[u];
53         }
54         return 0;
55     }
56 };
//hash模板，当 e数组较大时要开到外面，可用于求四元方程 a+b+c+d=0的解,状态标记 //可用于单词表查找，查找单词表中某单词是否可以由另两个单词合成(1~len-1拆分查找)
 1 const LL base=1e18;
 2 const LL MOD=1e9;
 3 struct int128
 4 {
 5     LL a, b;
 6     //a为低位
 7     int128(LL x = 0, LL y = 0)
 8     {
 9         a = x;
10         b = y;
11     }
12     int128 operator+(const int128 &x) const
13     {
14         return int128(a + x.a, b + x.b);
15     }
16     int128 operator-(const int128 &x) const
17     {
18         return int128(a - x.a, b - x.b);
19     }
20     int128 operator*(const int128 &x) const //相乘结果范围为[-(2^63)*1e18,(2^63-1)*1e18]
21     {
22         LL x0 = a % MOD, x1 = a / MOD;
23         LL y0 = x.a % MOD, y1 = x.a / MOD;
24         int128 c;
25         c.a = x0 * y0;
26         c.b = x1 * y1;
27         c.a += x0 * y1 % MOD * MOD;
28         c.b += x0 * y1 / MOD;
29         c.a += x1 * y0 % MOD * MOD;
30         c.b += x1 * y0 / MOD;
31         return c;
32     }
33     int128 operator/(const LL &x) const
34     {
35         LL ans, r, t = a, e = base / 10;
36         ans = b / x;
37         r = b % x;
38         for (int i = 17; i >= 0; i--, e /= 10)
39         {
40             ans = ans * 10 + (r * 10 + t / e) / x;
41             r = (r * 10 + t / e) % x;
42             t %= e;
43         }
44         return int128(ans);
45     }
46     int128 operator%(const LL &x) const
47     {
48         LL ans, t = a, e = base / 10;
49         ans = b % x;
50         for (int i = 17; i >= 0; i--, e /= 10)
51         {
52             ans = (ans * 10 + t / e) % x;
53             t %= e;
54         }
55         return int128(ans);
56     }
57 };
58 
59  
60 void print(int128 x)
61 {
62     if (x.b)
63         printf("%lld%018lld", x.b, x.a);
64     else
65         printf("%lld", x.a);
66 }
//大整数相关 //手动模拟__int28
 1 const int base=10;
 2 int largediv(char*a,char*b,char*ans)
 3 {
 4     int j = 0, la = strlen(a), lb = strlen(b), l = 0;
 5     while (la - j >= lb)
 6     {
 7         int t = 0;
 8         while ((j > 0 && a[j - 1] != '0') || strncmp(a + j, b, lb) >= 0)
 9         {
10             t++;
11             for (int i = lb - 1; i >= 0; i--)
12             {
13                 a[i + j] -= b[i] - '0';
14                 if (a[i + j] < '0')
15                 {
16                     a[i + j] += base;
17                     a[i + j - 1]--;
18                 }
19             }
20         }
21         j++;
22         ans[l++] = t + '0';
23         if (l == 0 && ans[l] == '0')
24             l--;
25     }
26     if (!l)
27         ans[l++] = '0';
28     ans[l] = 0;
29     int x = 0;
30     while (x < la && a[x] == '0')
31         x++;
32     if (x == la)
33         x--;
34     return x;
35 }
//大整数除法或求模，商为 ans(可能会有一个前导 0,输出要处理，注意结果为 0)，余数为 a+x
 1 struct Fraction
 2 {
 3     LL num;
 4     LL den;
 5     Fraction(LL num = 0, LL den = 1)
 6     {
 7         if (den < 0)
 8         {
 9             num = -num;
10             den = -den;
11         }
12         assert(den != 0);
13         LL g = __gcd(abs(num), den);
14         this->num = num / g;
15         this->den = den / g;
16     }
17     Fraction operator+(const Fraction &x) const
18     {
19         return Fraction(num * x.den + den * x.num, den * x.den);
20     }
21     Fraction operator-(const Fraction &x) const
22     {
23         return Fraction(num * x.den - den * x.num, den * x.den);
24     }
25     Fraction operator*(const Fraction &x) const
26     {
27         return Fraction(num * x.num, den * x.den);
28     }
29     Fraction operator/(const Fraction &x) const
30     {
31         return Fraction(num * x.den, den * x.num);
32     }
33     bool operator<(const Fraction &x) const
34     {
35         return num * x.den < den * x.num;
36     }
37     bool operator==(const Fraction &x) const
38     {
39         return num * x.den == den * x.num;
40     }
41 };
//分数类
  1 struct bigint
  2 {
  3     int *a;
  4     //a[0]存长度
  5     static const int maxlen = 1e4 + 10;
  6     static const int basel = 8;
  7     static const int base = 1e9;
  8     bigint() { a = new int[maxlen]; }
  9     ~bigint() { delete[] a; }
 10     //输出大整数 a,base不一样时记得改 04d
 11     void print()
 12     {
 13         int n = a[0];
 14         printf("%d", a[--n]);
 15         for (--n; n; n--)
 16             printf("%09d", a[n]);
 17         printf("\n");
 18     }
 19     //将正序的字符串 x转为 10^(basel+1)进制逆序大整数 a
 20     bigint &trans(char *x)
 21     {
 22         int lx = strlen(x);
 23         int &la = a[0] = 1;
 24         for (int i = lx - 1; i >= 0;)
 25         {
 26             int t = 0, k = i - (basel < i ? basel : i);
 27             for (int j = k; j <= i; j++)
 28                 t = t * 10 + x[j] - '0';
 29             i = k - 1;
 30             a[la++] = t;
 31         }
 32         return *this;
 33     }
 34     bool equal(int x) { return a[0] == 2 && a[1] == x; }
 35     //<base的整数转为大整数
 36     bigint &intTobig(int x)
 37     {
 38         a[0] = 2;
 39         a[1] = x;
 40         return *this;
 41     }
 42     bigint &copy(bigint & b)
 43     {
 44         memcpy(a, b.a, sizeof(a[0]) * b.a[0]);
 45         return *this;
 46     }
 47     void swap(bigint & b)
 48     {
 49         int *tmp = a;
 50         a = b.a;
 51         b.a = tmp;
 52     }
 53     bigint &add(bigint & tb, bigint & tans)
 54     {
 55         int *b = tb.a, *ans = tans.a;
 56         int l = max(a[0], b[0]);
 57         memset(ans, 0, sizeof(a[0]) * (l + 1));
 58         for (int i = 1; i < l; i++)
 59         {
 60             if (i >= a[0])
 61                 a[i] = 0;
 62             if (i >= b[0])
 63                 b[i] = 0;
 64             ans[i] += a[i] + b[i];
 65             ans[i + 1] += ans[i] / base;
 66             ans[i] %= base;
 67         }
 68         ans[0] = (ans[l] ? l + 1 : l);
 69         return tans;
 70     }
 71     bigint &minus(bigint & tb, bigint & tans)
 72     {
 73         int *b = tb.a, *ans = tans.a;
 74         memset(ans, 0, sizeof(int) * (a[0] + 1));
 75         for (int i = 1; i < a[0]; i++)
 76         {
 77             if (i >= b[0])
 78                 b[i] = 0;
 79             ans[i] += a[i] - b[i];
 80             if (ans[i] < 0)
 81             {
 82                 ans[i] += base;
 83                 ans[i + 1]--;
 84             }
 85         }
 86         ans[0] = a[0];
 87         while (ans[0] > 2 && !ans[ans[0] - 1])
 88             ans[0]--;
 89         return tans;
 90     }
 91     bigint &mul(bigint & tb, bigint & tans)
 92     {
 93         int *b = tb.a, *ans = tans.a;
 94         memset(ans, 0, sizeof(int) * (a[0] + b[0]));
 95         for (int i = 1; i < a[0]; i++)
 96         {
 97             for (int j = 1; j < b[0]; j++)
 98             {
 99                 LL t = ans[i + j - 1] + (LL)a[i] * b[j];
100                 ans[i + j] += t / base;
101                 ans[i + j - 1] = t % base;
102             }
103         }
104         int &l = ans[0];
105         l = a[0] + b[0] - 2;
106 
107         while (l > 1 && !ans[l])
108             l--;
109         l++;
110         return tans;
111     }
112     //大整数除法或求模,被除数<=Base
113     int div(int b, bigint &tans)
114     {
115         int *ans = tans.a;
116         ans[0] = 1;
117         ans[ans[0]++] = 0;
118         LL sum = 0, first = 1;
119         for (int i = a[0] - 1; i > 0; i--)
120         {
121             sum = sum * base + a[i];
122             if (sum < b && first)
123                 continue;
124             if (first)
125             {
126                 first = 0;
127                 ans[0] = i + 1;
128             }
129             ans[i] = sum / b;
130             sum %= b;
131         }
132         return sum;
133     }
134     bool less(bigint & x) //this<x
135     {
136         if (x.a[0] != a[0])
137             return a[0] < x.a[0];
138         for (int i = a[0] - 1; i > 0; i--)
139             if (a[i] != x.a[i])
140                 return a[i] < x.a[i];
141         return 0;
142     }
143 };
//大整数类
 1 int d[2000],dn,c[2000];
 2 void fac(int n)//找出所有约数
 3 {
 4     dn = 0;
 5     for (int i = 1; i * i <= n; i++)
 6         if (n % i == 0)
 7         {
 8             d[dn++] = i;
 9             if (i * i != n)
10                 d[dn++] = n / i;
11         }
12     sort(d, d + dn);
13 }
14 LL polya(int n)
15 {
16     fac(n);
17     memset(c, 0, sizeof(c));
18     for (int i = 0; i < dn; i++) //容斥求所有约数的欧拉函数
19     {
20         c[i] += d[i];
21         for (int j = i + 1; j < dn; j++)
22             if (d[j] % d[i] == 0)
23                 c[j] -= c[i];
24     }
25     LL ans = 0;
26     for (int i = 0; i < dn; i++)
27     {
28         int x = n / d[i];
29         ans = (ans + (LL)c[i] * mar_quick(x, 2, 1)) % mod; //计算 x 个循环节的置换有多少种方案,polya 定理中
30         为 m ^ x
31     }
32     ans = (ans * Inv(n)) % mod;
33     return ans;
34 }
//polya定理
 1 int d[2000],phi[2000],dn,p[20],pn,e[20];
 2 void facdiv(int n) //质因分解成 n=p[0]^e[0] * p[1]^e[1] * ... *p[pn-1]^e[pn-1]
 3 {
 4     pn = 0;
 5     //for(int j=0,i=prime[0];i*i<=n;i=prime[++j]) if(n%i==0)
 6     //筛出质数表加速质因分解
 7     for (int i = 2; i * i <= n; i++)
 8         if (n % i == 0)
 9         {
10             p[pn] = i;
11             e[pn] = 0;
12             while (n % i == 0)
13                 n /= i, e[pn]++;
14             pn++;
15         }
16     if (n != 1)
17         p[pn] = n, e[pn++] = 1;
18 }
//d为约数数组,phi[i]为第 i个约数的欧拉函数,p为 n分解出来的素数数组,e为次幂数组
 1 void GetDivisorPhi() //根据上面的基因分解算出所有约数的欧拉函数(与 d[i]互质的数的个数)
 2 {
 3     dn = 1;
 4     d[0] = 1;
 5     phi[0] = 1;
 6     for (int i = 0; i < pn; i++)
 7     {
 8         int num = dn;
 9         for (int j = 0; j < num; j++)
10         {
11             phi[dn] = phi[j] * (p[i] - 1);
12             d[dn] = d[j] * p[i];
13             dn++;
14             for (int k = 2; k <= e[i]; k++)
15             {
16                 phi[dn] = phi[dn - 1] * p[i];
17                 d[dn] = d[dn - 1] * p[i];
18                 dn++;
19             }
20         }
21     }
22 }
//根据上面的基因分解算出所有约数的欧拉函数(与 d[i]互质的数的个数)
 1 const int S=20;//随机算法判定次数，S越大，判错概率越小
 2 //计算 (a*b)%c. a,b都是 LL的数，直接相乘可能溢出的
 3 // a,b,c <2^63
 4 LL multi_mod( LL x , LL y ,LL mod)
 5 {
 6     return (x * y - (LL)(x / (long double)mod * y + 1e-3) * mod + mod) % mod;
 7 }
 8 //快速幂，计算 x^n %c
 9 LL quick_mod(LL t,LL n,LL mod)
10 {
11     t %= mod;
12     LL ans = 1;
13     while (n)
14     {
15         if (n & 1)
16             ans = multi_mod(ans, t, mod);
17         t = multi_mod(t, t, mod);
18         n >>= 1;
19     }
20     return ans;
21 }
22 
23  
24 //以 a为基,n-1=x*2^t
25 a^(n-1)=1(mod n) 验证 n是不是合数
26 //一定是合数返回 true,不一定返回 false
27 bool check(LL a,LL n,LL x,LL t)
28 {
29     LL ans = quick_mod(a, x, n);
30     LL last = ans;
31     for (int i = 1; i <= t; i++)
32     {
33         ans = multi_mod(ans, ans, n);
34         if (ans == 1 && last != 1 && last != n - 1)
35             return true; //合数
36         last = ans;
37     }
38     if (ans != 1)
39         return true;
40     return false;
41 }
42 // Miller_Rabin()算法素数判定
43 //是素数返回 true.(可能是伪素数，但概率极小)
44 //合数返回 false;
45 bool Miller_Rabin(LL n)
46 {
47     if (n < 2)
48         return false;
49     if (n == 2)
50         return true;
51     if ((n & 1) == 0)
52         return false; //偶数
53     LL x = n - 1;
54     LL t = 0;
55     while ((x & 1) == 0)
56     {
57         x >>= 1;
58         t++;
59     }
60     for (int i = 0; i < S; i++)
61     {
62         LL a = rand() % (n - 1) + 1;
63         if (check(a, n, x, t))
64             return false; //合数
65     }
66     return true;
67 }
// Miller_Rabin 算法进行素数测试 //速度快，而且可以判断 <2^63的数
 1 LL fac[100];//质因数分解结果（刚返回时是无序的）
 2 int tol;//质因数的个数。数组小标从 0开始
 3 LL gcd(LL a,LL b)
 4 {
 5     if (a == 0)
 6         return 1;
 7     if (a < 0)
 8         return gcd(-a, b);
 9     while (b)
10     {
11         LL t = a % b;
12         a = b;
13         b = t;
14     }
15     return a;
16 }
17 
18  
19 LL Pollard_rho(LL x,LL c)
20 {
21     LL i = 1, k = 2;
22     LL x0 = rand() % x;
23     LL y = x0;
24     while (1)
25     {
26         i++;
27         x0 = (multi_mod(x0, x0, x) + c) % x;
28         LL d = gcd(y - x0, x);
29         if (d != 1 && d != x)
30             return d;
31         if (y == x0)
32             return x;
33         if (i == k)
34             y = x0, k += k;
35     }
36 }
37 //对 n进行素因子分解，分解后的因子不是有序的，要 sort
38 void findfac(LL n)
39 {
40     if (n == 1)
41         return;
42     if (Miller_Rabin(n)) //素数
43     {
44         fac[tol++] = n;
45         return;
46     }
47     LL p = n;
48     while (p >= n)
49         p = Pollard_rho(p, rand() % (n - 1) + 1);
50     findfac(p);
51     findfac(n / p);
52 }
//pollard_rho 算法进行质因数分解，复杂度比 sqrt(n)低
 1 LL multi_mod(LL a,LL b,LL mod)//a,b传过来的要是正数
 2 {
 3     LL ans = 0;
 4     a %= mod;
 5     b %= mod;
 6     while (b)
 7     {
 8         if (b & 1)
 9             ans = ans + a;
10         if (ans >= mod)
11             ans -= mod;
12         a <<= 1;
13         if (a >= mod)
14             a -= mod;
15         b >>= 1;
16     }
17     return ans;
18 }
//模乘法，模拟二进制乘法求模,在 long long*long long求模的情况下不溢出 O(64)
1 LL multi_mod( LL x , LL y ,LL mod)
2 {
3     return (x * y - (LL)(x / (long double)mod * y + 1e-3) * mod + mod) % mod;
4 }
//比上面的模乘法快,并不懂原理
 1 //multi_mod 中过多的传参会耗费大量时间,不爆 long long 的话最好写成 a*b%mod
 2 //如果精度超过 double 可以用 NTT 求多个 P,再用中国剩余定理算回去
 3 //f(i)=sigma{a[j]*b[i-j],0<=i<n}
 4 //const LL P = ( 1LL << 55 ) * 5 + 1 ;
 5 //const LL G = 6 ;
 6 const int P =479<<21|1;
 7 //const int P1=119<<23|1;//费马质数,P-1 能整除 2^n
 8 const int G = 3; //原根,最小的数满足 G^(P-1)==1
 9 const int NUM = 30;
10 LL wn[NUM];
11 LL x[N*4],y[N*4];
12 void GetWn()
13 {
14     for (int i = 0; i < NUM; i++)
15     {
16         int t = 1 << i;
17         wn[i] = quick_mod(G, (P - 1) / t);
18     }
19 }
20 void NTT(LL *a,int len,bool rev)
21 {
22     for (int i = 1, j, t, k; i < len; i++)
23     {
24         for (k = len >> 1, t = i, j = 0; k; k >>= 1, t >>= 1)
25             j = j << 1 | t & 1;
26         if (i < j)
27             swap(a[i], a[j]);
28     }
29     for (int s = 2, ds = 1, id = 1; s <= len; ds = s, s <<= 1, id++)
30     {
31         for (int k = 0; k < len; k += s)
32         {
33             LL w = 1, t;
34             for (int i = k; i < k + ds; i++)
35             {
36                 a[i + ds] = (a[i] - (t = multi_mod(a[i + ds], w)) + P) % P;
37                 a[i] = (a[i] + t) % P;
38                 w = multi_mod(w, wn[id]);
39             }
40         }
41     }
42     if (rev)
43     {
44         for (int i = 1; i<len>> 1; i++)
45             swap(a[i], a[len - i]);
46         LL inv = quick_mod(len, P - 2);
47         for (int i = 0; i < len; i++)
48             a[i] = multi_mod(a[i], inv);
49     }
50 }
51 void Conv(LL *a, LL *b, int n)
52 {
53     NTT(a, n, 0);
54     NTT(b, n, 0);
55     for (int i = 0; i < n; i++)
56         a[i] = multi_mod(a[i], b[i]);
57     NTT(a, n, 1);
58 }
//快速数论变换 NTT,得到的多项式的每一位的结果已经模了 P //不比 FFT 快,因为求模花的时间较多
 1 const double pi=acos(-1);
 2 struct cplx
 3 {
 4     //复数类/多项式/大整数
 5     double r, i;
 6     cplx(double r = 0, double i = 0) : r(r), i(i) {}
 7     cplx operator+(const cplx &a) const
 8     {
 9         return cplx(r + a.r, i + a.i);
10     }
11     cplx operator-(const cplx &a) const
12     {
13         return cplx(r - a.r, i - a.i);
14     }
15     cplx operator*(const cplx &a) const
16     {
17         return cplx(r * a.r - i * a.i, r * a.i + i * a.r);
18     }
19 };
20 //快速傅里叶变换 FFT(nlogn)
21 //len = 2^M,reverse F[i] with F[j] j为 i二进制反转
22 void rader(cplx *f,int len)
23 {
24     int j = len >> 1;
25     for (int i = 1; i < len - 1; ++i)
26     {
27         if (i < j)
28             swap(f[i], f[j]); // reverse
29         int k = len >> 1;
30         while (j >= k)
31         {
32             j -= k;
33             k >>= 1;
34         }
35         if (j < k)
36             j += k;
37     }
38 }
39 //做 FFT,len必须为 2^k形式，on==1时是 DFT，on==-1时是 IDFT
40 void FFT(cplx *f,int len,int t)
41 {
42     rader(f, len);
43     for (int h = 2; h <= len; h <<= 1)
44     {
45         cplx wn(cos(-t * 2 * pi / h), sin(-t * 2 * pi / h));
46         for (int j = 0; j < len; j += h)
47         {
48             cplx e(1, 0); //旋转因子
49             for (int k = j; k < j + h / 2; ++k)
50             {
51                 cplx u = f[k];
52 
53                 cplx v = e * f[k + h / 2];
54                 f[k] = u + v;
55                 f[k + h / 2] = u - v;
56                 e = e * wn;
57             }
58         }
59     }
60     if (t == -1)
61         for (int i = 0; i < len; ++i)
62             f[i].r /= len; //IDFT
63 }
64 //传入 a，b之前需将 a，b的长度处理为 len=2^k,a，b数组需开到原长度的 3-4倍,变换结果为 a数组
65 void Conv(cplx *a,cplx *b,int len) //求卷积
66 {
67     FFT(a, len, 1);
68     FFT(b, len, 1);
69     for (int i = 0; i < len; ++i)
70         a[i] = a[i] * b[i];
71     FFT(a, len, -1);
72 }
//快速傅里叶变换 FFT(nlogn)
  1 const int N = 5e6+2;
  2 bool np[N];
  3 int p[N],pi[N];
  4 int getprime()
  5 {
  6     //pi[i]表示不超过 i的质数个数，p存素数
  7     int cnt = 0;
  8     np[0] = np[1] = true;
  9     pi[0] = pi[1] = 0;
 10     for (int i = 2; i < N; ++i)
 11     {
 12         if (!np[i])
 13             p[++cnt] = i;
 14         for (int j = 1; j <= cnt && i * p[j] < N; ++j)
 15         {
 16             np[i * p[j]] = true;
 17         }
 18         pi[i] = cnt;
 19     }
 20     return cnt;
 21 }
 22 const int M = 7;
 23 const int PM = 2*3*5*7*11*13*17;
 24 int phi[PM+1][M+1],sz[M+1];
 25 void init()
 26 {
 27     getprime();
 28     sz[0] = 1;
 29     for (int i = 0; i <= PM; ++i)
 30         phi[i][0] = i;
 31     for (int i = 1; i <= M; ++i)
 32     {
 33         sz[i] = p[i] * sz[i - 1];
 34         for (int j = 1; j <= PM; ++j)
 35         {
 36             phi[j][i] = phi[j][i - 1] - phi[j / p[i]][i - 1];
 37         }
 38     }
 39 }
 40 
 41  
 42 int sqrt2(LL x)
 43 {
 44     LL r = (LL)sqrt(x - 0.1);
 45     while (r * r <= x)
 46         ++r;
 47     return int(r - 1);
 48 }
 49 int sqrt3(LL x)
 50 {
 51     LL r = (LL)cbrt(x - 0.1);
 52     while (r * r * r <= x)
 53         ++r;
 54     return int(r - 1);
 55 }
 56 LL getphi(LL x,int s)
 57 {
 58     if (s == 0)
 59         return x;
 60     if (s <= M)
 61         return phi[x % sz[s]][s] + (x / sz[s]) * phi[sz[s]][s];
 62     if (x <= p[s] * p[s])
 63         return pi[x] - s + 1;
 64     if (x <= p[s] * p[s] * p[s] && x < N)
 65     {
 66         int s2x = pi[sqrt2(x)];
 67         LL ans = pi[x] - (s2x + s - 2) * (s2x - s + 1) / 2;
 68         for (int i = s + 1; i <= s2x; ++i)
 69         {
 70             ans += pi[x / p[i]];
 71         }
 72         return ans;
 73     }
 74     return getphi(x, s - 1) - getphi(x / p[s], s - 1);
 75 }
 76 LL getpi(LL x)
 77 {
 78     //返回小于 x 的素数个数
 79     if (x < N)
 80         return pi[x];
 81     LL ans = getphi(x, pi[sqrt3(x)]) + pi[sqrt3(x)] - 1;
 82     for (int i = pi[sqrt3(x)] + 1, ed = pi[sqrt2(x)]; i <= ed; ++i)
 83         ans -= getpi(x / p[i]) - i + 1;
 84     return ans;
 85 }
 86 LL lehmer_pi(LL x)
 87 {
 88     if (x < N)
 89         return pi[x];
 90     int a = (int)lehmer_pi(sqrt2(sqrt2(x)));
 91     int b = (int)lehmer_pi(sqrt2(x));
 92     int c = (int)lehmer_pi(sqrt3(x));
 93     LL sum = getphi(x, a) + LL(b + a - 2) * (b - a + 1) / 2;
 94     for (int i = a + 1; i <= b; i++)
 95     {
 96         LL w = x / p[i];
 97         sum -= lehmer_pi(w);
 98         if (i > c)
 99             continue;
100         LL lim = lehmer_pi(sqrt2(w));
101         for (int j = i; j <= lim; j++)
102         {
103             sum -= lehmer_pi(w / p[j]) - (j - 1);
104         }
105     }
106     return sum;
107 
108  
109 }
110 LL getans(LL x)
111 { // x < 1e11
112     LL ans = 0;
113     //返回相乘不大于 x的质数对数
114     for (int i = 1, ed = pi[sqrt2(x - 1)]; i <= ed; ++i)
115         ans += lehmer_pi(x / p[i]) - i;
116     return ans;
117 }
//算 1e11内质数个数，copy的，先要执行 init()函数
1 void swap(int &a,int &b)
2 {
3     a ^= b;
4     b ^= a;
5     a ^= b;
6 }
//交换两个数
1 struct node
2 {
3     int i, w;
4     bool operator<(const node &a) const
5     {
6         return w > a.w;
7     }
8 };
//优先队列结构体
 1 int prim()
 2 {
 3     //O(nlogn)
 4     memset(d, 0x3f, sizeof(d));
 5     memset(vis, 0, sizeof(vis));
 6     d[1] = 0;
 7     priority_queue<node> q;
 8     q.push(node{1, 0});
 9     while (!q.empty())
10     {
11         node u = q.top();
12         q.pop();
13         if (vis[u.i])
14             continue;
15         vis[u.i] = 1;
16         for (vector<int>::iterator i = v[u.i].begin(), j = w[u.i].begin(); i != v[u.i].end(); i++, j++)
17         {
18             if (vis[*i] || d[*i] <= *j)
19                 continue;
20             d[*i] = *j;
21             q.push(node{*i, d[*i]});
22         }
23     }
24     int sum = 0;
25     for (int i = 1; i <= n; i++)
26         sum += d[i];
27     return sum;
28 }
29 
30  
31 int prim()
32 {
33     //O(n^2)
34     memset(d, 0x3f, sizeof(d));
35     memset(vis, 0, sizeof(vis));
36     d[1] = 0;
37     for (int i = 2; i <= n; i++)
38     {
39         int x = -1;
40         for (int j = 1; j <= n; j++)
41             if (!vis[j] && (x == -1 || d[j] < d[x]))
42                 x = j; //j要未作为过起点
43         vis[x] = 1;
44         for (vector<int>::iterator j = v[x].begin(), k = w[x].begin(); j != v[x].end(); j++, k++)
45             if (!vis[*j] && *k < d[*j])
46                 d[*j] = *k;
47     }
48     int sum = 0;
49     for (int i = 1; i <= n; i++)
50         sum += d[i];
51     return sum;
52 }
//Prim加点法最小生成树
 1 struct maxMacth
 2 {
 3     int n;
 4     int vm[N], um[N];
 5     bool vis[N];
 6     vector<int> g[N];
 7     int dx[N], dy[N], dis;
 8     void init(int num)
 9     {
10         n = num;
11         memset(vm, -1, sizeof(vm));
12         memset(um, -1, sizeof(um));
13         for (int i = 0; i <= n; i++)
14             g[i].clear();
15     }
16     void inserts(int u, int v)
17     {
18         g[u].push_back(v);
19     }
20     bool searchP()
21     {
22         queue<int> q;
23         dis = INF;
24         memset(dx, -1, sizeof(dx));
25         memset(dy, -1, sizeof(dy));
26         for (int i = 1; i <= n; i++)
27             if (um[i] == -1)
28             {
29                 q.push(i);
30                 dx[i] = 0;
31             }
32         while (!q.empty())
33 
34         {
35             int u = q.front();
36             q.pop();
37             if (dx[u] > dis)
38                 break;
39             for (int i = 0; i < g[u].size(); i++)
40             {
41                 int v = g[u][i];
42                 if (dy[v] == -1)
43                 {
44                     dy[v] = dx[u] + 1;
45                     if (vm[v] == -1)
46                         dis = dy[v];
47                     else
48                     {
49                         dx[vm[v]] = dy[v] + 1;
50                         q.push(vm[v]);
51                     }
52                 }
53             }
54         }
55         return dis != INF;
56     }
57     bool dfs(int u)
58     {
59         for (int i = 0; i < g[u].size(); i++)
60         {
61             int v = g[u][i];
62             if (!vis[v] && dy[v] == dx[u] + 1)
63             {
64                 vis[v] = 1;
65                 if (vm[v] != -1 && dy[v] == dis)
66                     continue;
67                 if (vm[v] == -1 || dfs(vm[v]))
68                 {
69                     vm[v] = u;
70                     um[u] = v;
71                     return 1;
72                 }
73             }
74         }
75         return 0;
76     }
77     int maxMatch()
78     {
79         int res = 0;
80         while (searchP())
81         {
82             memset(vis, 0, sizeof(vis));
83             for (int i = 1; i <= n; i++)
84                 if (um[i] == -1 && dfs(i))
85                     res++;
86         }
87         return res;
88     }
89 }MM;
//二分图最大匹配
 1 const int N=1e6+10;
 2 const int M=26;
 3 int ch[N][M];
 4 int val[N];
 5 int sz;
 6 void init()
 7 {
 8     sz = 1;
 9     memset(ch[0], 0, sizeof(ch[0]));
10 }
11 inline int idx(char c)
12 {
13     return c - 'a'; 
14 }
15 void insert(char *s,int v)
16 {
17     int u = 0, i;
18     for (i = 0; s[i]; i++)
19     {
20         int c = idx(s[i]);
21         if (!ch[u][c])
22         {
23             memset(ch[sz], 0, sizeof(ch[sz]));
24             val[sz] = 0;
25             ch[u][c] = sz++;
26         }
27         u = ch[u][c];
28     }
29     val[u] = v;
30     //val[u]=max(val[u],i);//val记录长度
31 }
32 int f[N],last[N];//last数组连接着所有相同后缀结束的模板串
33 void getFail()
34 {
35     queue<int> q;
36     f[0] = 0;
37     for (int c = 0; c < M; c++)
38     {
39         int u = ch[0][c];
40         if (u)
41         {
42             f[u] = 0;
43             q.push(u);
44             last[u] = 0;
45         }
46     }
47     while (!q.empty())
48     {
49         int r = q.front();
50         q.pop();
51         for (int c = 0; c < M; c++)
52         {
53             int u = ch[r][c];
54             if (!u)
55                 continue;
56             q.push(u);
57 
58             int v = f[r];
59             while (v && !ch[v][c])
60                 v = f[v];
61             f[u] = ch[v][c];
62             last[u] = val[f[u]] ? f[u] : last[f[u]];
63         }
64     }
65 }
66 void print(int i,int j)//输出所有匹配位置
67 {
68     if (j)
69     {
70         //找到后的处理,i是结束位置,字符串长度或值为 val[i]
71         print(i, last[j]);
72     }
73 }
74 int find(char *T)
75 {
76     int j = 0, i = 0;
77     for (; T[i]; i++)
78     {
79         int c = idx(T[i]);
80         while (j && !ch[j][c])
81             j = f[j];
82         j = ch[j][c];
83         if (val[j])
84             print(i, j);
85         else if (last[j])
86             print(i, last[j]);
87     }
88     return i;
89 }
//AC自动机 //字符串字典树
 1 const int maxc=2;
 2 const int BIT=31;
 3 int ch[N<<5][maxc],si[N<<5];
 4 void init()
 5 {
 6     //每次建一棵新的字典树都要初始化
 7     memset(ch[0], 0, sizeof(ch[0]));
 8     si[0] = 1;
 9 }
10 void insert(int s)
11 {
12     int u = 0;
13     for (int i = BIT - 1; i >= 0; i--)
14     {
15         int c = (s & (1 << i) ? 1 : 0);
16         if (!ch[u][c])
17         {
18             memset(ch[si[0]], 0, sizeof(ch[0]));
19             si[si[0]] = 0;
20             ch[u][c] = si[0]++;
21         }
22         u = ch[u][c];
23         si[u]++;
24     }
25 }
26 int find(int s)
27 {
28     int u = 0;
29     for (int i = BIT - 1; i >= 0; i--)
30     {
31         int c = (s & (1 << i) ? 1 : 0);
32         if (!ch[u][c] || si[ch[u][c]] <= 0)
33             return 0;
34         u = ch[u][c];
35     }
36     return 1;
37 }
38 bool erase(int s)
39 {
40     int u = 0;
41     for (int i = BIT - 1; i >= 0; i--)
42     {
43         int c = (s & (1 << i) ? 1 : 0);
44         if (!ch[u][c] || si[ch[u][c]] <= 0)
45             return 0;
46         u = ch[u][c];
47         si[u]--;
48     }
49     return 1;
50 }
//二进制数字典树
 1 int ch[N*20][2],si[N*20];
 2 //1 bro 0 son
 3 char val[N*20];
 4 void init()
 5 {
 6     ch[0][0] = ch[0][1] = 0;
 7     si[0] = 1;
 8     val[0] = 0;
 9 }
10 void insert(char *s)
11 {
12     int u = 0;
13     for (int i = 0; s[i]; i++)
14     {
15         int v = ch[u][0];
16         while (v && val[v] != s[i])
17             v = ch[v][1];
18         if (!v)
19         {
20             //链表头插
21             int &t = si[0];
22             ch[t][1] = ch[u][0];
23 
24             si[t] = ch[t][0] = 0;
25             val[t] = s[i];
26             v = ch[u][0] = t++;
27         }
28         u = v;
29     }
30     si[u] = 1;
31 }
32 bool find(char *s)
33 {
34     int u = 0;
35     for (int i = 0; s[i]; i++)
36     {
37         for (u = ch[u][0]; u && val[u] != s[i]; u = ch[u][1])
38             ;
39         if (!u)
40             return 0;
41     }
42     if (si[u])
43         return 1;
44     return 0;
45 }
46 //字符串哈希 H(i)=H(i+1)x+s[i],Hash(i,L)=H(i)-H(i+L)*x^L;
47 int bfind()
48 {
49     int l = 1, r = n;
50     while (l < r)
51     {
52         int mid = r - (r - l) / 2;
53         if (hash_l(mid))
54             l = mid;
55         else
56             r = mid - 1;
57     }
58     return l;
59 }
} //左兄弟右儿子字符串字典树(二叉，省空间，时间换空间)
 1 const int N=1e6+10;
 2 int bit=8;
 3 int base=1<<bit;
 4 int mask=base-1;
 5 int maxbit=20;
 6 int n,x[N],h,next[N],m;
 7 int head[N],tial[N];
 8 void radixsort(int &h,int *x)
 9 {
10     for (int i = 0; i < maxbit; i += bit)
11     {
12         memset(head, 0, sizeof(head[0]) * base);
13         for (int j = h, tmp; j; j = tmp)
14         //分桶
15 
16         {
17             tmp = next[j];
18             int c = (x[j] >> i) & mask;
19             if (!head[c])
20             {
21                 head[c] = tial[c] = j;
22                 next[j] = 0;
23             }
24             else
25             {
26                 next[tial[c]] = j;
27                 next[j] = 0;
28                 tial[c] = j;
29             }
30         }
31         h = 0;
32         int ti;
33         //for(int j=0;j<base;j++) if(head[j]) //合并,升序
34         for (int j = mask; j >= 0; j--)
35             if (head[j]) //合并，降序
36             {
37                 if (!h)
38                     h = head[j];
39                 else
40                     next[ti] = head[j];
41                 ti = tial[j];
42             }
43     }
44 }
45 void init()
46 {
47     bit = 0;
48     while (1 << bit < n)
49         bit++;
50     bit -= 3;
51     bit = std::min(20, bit); //最大的比 maxbit小一点
52     bit = std::max(4, bit);
53     base = 1 << bit;
54     mask = base - 1;
55 }
56 int main()
57 {
58     while (scanf("%d%d", &n, &m) == 2)
59     {
60         for (int i = 1; i <= n; i++)
61         {
62             scanf("%d", x + i);
63             x[i] += 500000;
64             next[i] = i + 1;
65         }
66         init();
67         h = 1;
68         next[n] = 0;
69         radixsort(h, x);
70 
71         for (int i = h; i && m; i = next[i])
72             printf("%d%c", x[i] - 500000, (--m) ? ' ' : '\n');
73         return 0;
74     }
75 }
//多关键基数排序
 1 template <class T>
 2 void qsort(int l,int r,T *p)
 3 {
 4     if (l >= r)
 5         return;
 6     T t = p[l];
 7     int i = l, j = r;
 8     while (i < j)
 9     {
10         //以第一个为基准排序，大的放左边，小的放右边
11         while (i < j && p[j]->score <= t->score)
12             j--; //找到一个小的
13         p[i] = p[j];
14         while (i < j && p[i]->score > t->score)
15             i++; //找到一个大的
16         p[j] = p[i];
17     }
18     p[i] = t;
19     qsort(l, i - 1, p);
20     qsort(i + 1, r, p);
21 }
//快排
 1 template<class T>
 2 void mergesort(int l,int r,T *p,T *tmp)
 3 {
 4     if (l == r)
 5         return;
 6     int mid = r - (r - l) / 2;
 7     mergesort(l, mid - 1, p, tmp);
 8     mergesort(mid, r, p, tmp);
 9     int i = l, j = mid, k = l;
10     while (i < mid || j <= r)
11     {
12         //归并
13         if (j > r || i < mid && p[i]->score >= p[j]->score)
14             tmp[k++] = p[i++];
15         else
16             tmp[k++] = p[j++];
17     }
18     while (l <= r)
19         p[l] = tmp[l++];
20     //复制回原数组
21 }
//归并排序
 1 template <class T>
 2 void HeapAdjust(T *p,int s,int m)
 3 {
 4     T rc = p[s];
 5     //沿 key较大的孩子结点向下筛选
 6     for (int j = s << 1; j <= m; j <<= 1)
 7     {
 8         if (j < m && p[j] < p[j + 1])
 9             ++j;
10         if (rc >= p[j])
11             break;
12         p[s] = p[j];
13         s = j;
14         //j为 key较大的记录的下标
15         //rc应插入在位置 s上
16     }
17     p[s] = rc;
18     //插入
19 }
20 template <class T>
21 void HeapSort(T *p)
22 {
23     //数组下标从 1开始
24     for (int i = n / 2; i > 0; --i)
25         HeapAdjust(p, i, n); //建堆,不断将子二叉树调整为堆，然后整个二叉树就满足
26     了堆性质
27     for (int i = n, j = 1; j <= m; i--, j++)
28     //每次将当前最大(小)(堆顶)元素,放到尾部，调整堆
29     {
30         swap(p[1], p[i]);
31         HeapAdjust(p, 1, i - 1);
32     }
33 }
//堆排序 //将某个子二叉树调整为堆,满足 p[i]>=k[i*2]&&p[i]>=k[i*2+1](或小于)
 1 template <class T>
 2 void HeapUpAdjust(T *p,int s)//从 s向上到 1调整堆
 3 {
 4     int x = p[s];
 5     while (s > 1)
 6     {
 7         if (x < p[s >> 1])
 8             p[s] = p[s >> 1];
 9         else
10             break;
11         s >>= 1;
12     }
13     p[s] = x;
14 }
15 template <class T>
16 void HeapDownAdjust(T *p,int s,int n)//从 s向下到 n调整堆
17 {
18     int x = p[s];
19     s <<= 1;
20     while (s <= n)
21     {
22         if (s < n && p[s | 1] < p[s])
23             s |= 1;
24         if (p[s] >= x)
25             break;
26         p[s >> 1] = p[s];
27 
28         s <<= 1;
29         p[s >> 1] = x;
30     }
31 }
//堆调整

优先队列入队：加元素到尾部，然后向上调整，size++
优先队列出队：将最后一个元素放到顶部，向下调整堆，size--

 1 int n,a[N],d[N][20];
 2 void RMQ_init()
 3 {
 4     for (int i = 1; i <= n; i++)
 5         d[i][0] = a[i];
 6     for (int j = 1; (1 << j) <= n; j++)
 7         for (int i = 1; i + (1 << j) - 1 <= n; i++)
 8             d[i][j] = min(d[i][j - 1], d[i + (1 << (j - 1))][j - 1]);
 9 }
10 int RMQ(int L,int R)
11 {
12     int k = 0;
13     while ((1 << (k + 1)) <= R - L + 1)
14         k++;
15     return min(d[L][k], d[R - (1 << k) + 1][k]);
16 }
//RMQ
公平游戏:一个游戏者可以把状态 A变成状态 B,另一个也可以.象棋是不公平游戏
规则 1:一个状态是必败状态当且仅当它的所有后继都是必胜状态
规则 2:一个状态是必胜状态当且仅当它至少有一个后继是必败状态
dp时,对于每个状态,查找它的后继状态有没有必败状态,有就是必胜状态,没有就是必败状态
当数据较大,dp打表会超时时就打出一个比较小的表找数(0出现)的规律
  1 struct Point
  2 {
  3     double x, y;
  4     Point(double x = 0, double y = 0) : x(x), y(y) {}
  5 };
  6 typedef Point Vector;
  7 Vector operator + (Vector A,Vector B) {
  8     return Vector(A.x + B.x, A.y + B.y); }
  9 Vector operator - (Point A,Point B) {
 10     return Vector(A.x - B.x, A.y - B.y); }
 11 Vector operator * (Vector A,double p) {
 12     return Vector(A.x * p, A.y * p); }
 13 Vector operator / (Vector A,double p) {
 14     return Vector(A.x / p, A.y / p); }
 15 bool operator < (const Point & a,const Point &b) {
 16     return a.x < b.x || (a.x == b.x && a.y < b.y); }
 17 const double eps=1e-9;
 18 const double pi=acos(-1);
 19 int dcmp(double x) {
 20     if (fabs(x) < eps)
 21         return 0;
 22     else
 23         return x < 0 ? -1 : 1; }
 24 bool operator == (const Point& a,const Point &b) {
 25     return dcmp(a.x - b.x) == 0 && dcmp(a.y - b.y) == 0; }
 26 
 27  
 28 //基本运算
 29 double Dot(Vector A,Vector B) //向量的点积
 30 {
 31     return A.x * B.x + A.y * B.y;
 32 }
 33 double Length(Vector A) //向量的长度
 34 {
 35     return sqrt(Dot(A, A));
 36 }
 37 double Angle(Vector A,Vector B) //两个向量的夹角
 38 {
 39     return acos(Dot(A, B) / Length(A) / Length(B));
 40 }
 41 double Cross(Vector A,Vector B) //向量叉积
 42 {
 43     return A.x * B.y - A.y * B.x;
 44 }
 45 double Area2(Point A,Point B,Point C)//三角形面积
 46 {
 47     return Cross(B - A, C - A);
 48 }
 49 Vector Rotate(Vector A,double rad)//向量逆时针旋转
 50 {
 51     return Vector(A.x * cos(rad) - A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad));
 52 }
 53 Vector Normal(Vector A)
 54 {//计算向量的单位法向量(左转 90度以后把长度归一化),调用前确保 A不是 0向量
 55     double L = Length(A);
 56     return Vector(-A.y / L, A.x / L);
 57 }
 58 //点和直线
 59 Point GetLineIntersection(Point P,Vector v,Point Q,Vector w)
 60 {//调用前确保两条直线 P+tv和 Q+tw有唯一交点,即 Cross(v,w)非 0
 61     Vector u = P - Q;
 62     double t = Cross(w, u) / Cross(v, w);
 63     return P + v * t;
 64 }
 65 double DistanceToLine(Point P,Point A,Point B)
 66 {//点到直线的距离
 67     Vector v1 = B - A, v2 = P - A;
 68     return fabs(Cross(v1, v2)) / Length(v1); //不取绝对值得到有向距离
 69 }
 70 
 71  
 72 double DistanceToSegment(Point P,Point A,Point B)
 73 {//点到线段的距离
 74     if (A == B)
 75         return Length(P - A);
 76     Vector v1 = B - A, v2 = P - A, v3 = P - B;
 77     if (dcmp(Dot(v1, v2)) < 0)
 78         return Length(v2);
 79     else if (dcmp(Dot(v1, v3)) > 0)
 80         return Length(v3);
 81     else
 82         return fabs(Cross(v1, v2)) / Length(v1);
 83 }
 84 Point GetLineProjection(Point P,Point A,Point B)
 85 {//点在直线上的投影
 86     Vector v = B - A;
 87     return A + v * (Dot(v, P - A) / Dot(v, v));
 88 }
 89 bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2)
 90 {//判线段相交,两线段恰好有一个公共点且不再任何一条线段的端点
 91     double c1 = Cross(a2 - a1, b1 - a1), c2 = Cross(a2 - a1, b2 - a1);
 92     double c3 = Cross(b2 - b1, a1 - b1), c4 = Cross(b2 - b1, a2 - b1);
 93     return dcmp(c1) * dcmp(c2) < 0 && dcmp(c3) * dcmp(c4) < 0;
 94 }
 95 bool PointOnSegment(Point p,Point a1,Point a2)
 96 {//判断点是否在线段上(不含端点,端点直接判 p和 a1或 a2是否相等)
 97     //if(p==a1||p==a2) return 1;//判点相同
 98     return dcmp(Cross(a1 - p, a2 - p)) == 0 && dcmp(Dot(a1 - p, a2 - p)) < 0;
 99 }
100 //多边形
101 double PolygonArea(Point *p,int n)
102 {//多边形有向面积
103     double area = 0;
104     for (int i = 1; i < n - 1; i++)
105         area += Cross(p[i] - p[0], p[i + 1] - p[0]);
106     return area / 2;
107 }
108 int isPointInPolygon(Point p,Point*poly,int n)
109 {//点在多边形判定
110     int wn = 0;
111     for (int i = 0; i < n; i++)
112     {
113         if (PointOnSegment(p, poly[i], poly[(i + 1) % n]))
114             return -1; //在边界上
115         int k = dcmp(Cross(poly[(i + 1) % n] - poly[i], p - poly[i]));
116         int d1 = dcmp(poly[i].y - p.y);
117         int d2 = dcmp(poly[(i + 1) % n].y - p.y);
118         if (k > 0 && d1 <= 0 && d2 > 0)
119             wn++;
120         if (k < 0 && d2 <= 0 && d1 > 0)
121             wn--;
122     }
123     if (wn != 0)
124         return 1; //内部
125     return 0;
126     //外部
127 }
128 
129  
130 //计算凸包,输入点数组 p,个数为 p,输出点数组 ch.函数返回凸包顶点数
131 //输入不能有重复点.注意:函数执行完之后输入点的顺序被破坏
132 //如果不希望在凸包的边上有输入点,把两个<=改成<
133 //在精度要求高时建议用 dcmp比较
134 int ConvexHull(Point* p,int n,Point *ch)
135 {
136     sort(p, p + n);
137     int m = 0;
138     for (int i = 0; i < n; i++)
139     {
140         while (m > 1 && Cross(ch[m - 1] - ch[m - 2], p[i] - ch[m - 2]) <= 0)
141             m--;
142         ch[m++] = p[i];
143     }
144     int k = m;
145     for (int i = n - 2; i >= 0; i--)
146     {
147         while (m > k && Cross(ch[m - 1] - ch[m - 2], p[i] - ch[m - 2]) <= 0)
148             m--;
149         ch[m++] = p[i];
150     }
151     if (n > 1)
152         m--;
153     return m; //凸包顶点数
154 }
155 double PointoRecDist(Point p,Point ra,Point rb)
156 {//点到长方形的距离
157     if (ra.x > rb.x)
158         swap(ra.x, rb.x); //变成左上和右下角
159     if (ra.y > rb.y)
160         swap(ra.y, rb.y);
161     double dx = 0, dy = 0;
162     if (p.x < ra.x)
163         dx = ra.x - p.x;
164     else if (p.x > rb.x)
165         dx = p.x - rb.x;
166     if (p.y < ra.y)
167         dy = ra.y - p.y;
168     else if (p.y > rb.y)
169         dy = p.y - rb.y;
170     return sqrt(dx * dx + dy * dy);
171 }
172 double getcycy(double x,double neg,Point cyc,double ri)
173 {//由 x得到圆上点的 y坐标,neg表示上 1下-1
174     return cyc.y + neg * sqrt(ri * ri - (x - cyc.x) * (x - cyc.x));
175 }
//几何模板
 1 const double eps=1e-9;
 2 int gauss(double a[N][N],bool *l,double *ans,int n)
 3 {
 4     int res = 0, r = 0;
 5     for (int i = 0; i < n; i++)
 6         l[i] = false; //全部初始化为自由元
 7     for (int i = 0; i < n; i++)
 8     {
 9         for (int j = r; j < n; j++)
10             if (fabs(a[j][i] > eps))
11             {
12                 for (int k = i; k <= n; k++)
13                     swap(a[j][k], a[r][k]);
14                 break;
15             }
16         if (fabs(a[r][i]) < eps)
17         {
18             res++;
19             continue;
20         }
21         for (int j = 0; j < n; j++)
22             if (j != r && fabs(a[j][i]) > eps)
23             {
24                 double tmp = a[j][i] / a[r][i];
25                 for (int k = i; k <= n; k++)
26                     a[j][k] -= tmp * a[r][k];
27             }
28         l[i] = true;
29         r++;
30     }
31     for (int i = 0; i < n; i++)
32         if (l[i])
33             for (int j = 0; j < n; j++)
34                 if (fabs(a[j][i]) > 0)
35                     ans[i] = a[j][n] / a[j][i];
36     return res;
37 }
//gauss消元模板,返回矩阵的秩,l数组标记自由元,ans为结果,n为未知数个数
 1 //若只是记录某个数是否存在可用 bitset优化，原来已存在的数记为 bitset x ；
 2 //则 x中所有标记位为 1的数+a[i]之后变成 x=x<<a[i]|x
 3 //如果数最大小于 64可用 long long代替 bitset,更快,最大小于 128时用__int128
 4 //无限背包有排列种数
 5 //(1,4,9,...,17^2 中可重复选能有多少种方法拼成 1~300,114,141,411,111111 4种拼成 6的方法)
 6 memset(dp,0,sizeof(dp));
 7 dp[0]=1;
 8 for(int i=0;i<=300;i++) //每个数都可以往后面加 1,4,...,17^2
 9 for(int j=1;j*j+i<=300;j++)
10 dp[j*j+i]+=dp[i];//取多少个 i
11 //无限背包组合种数
12 //(1,4,9,...,17^2 中可重复选能有多少种方法拼成 1~300,114,111111 2种拼成 6的方法)
13 memset(dp,0,sizeof(dp));
14 dp[0]=1;
15 for(int i=1;i*i<=300;i++)//取的数递增,即取了 1后只能取 4再 9再...,不能继续取 1了
16 for(int j=0;j+i*i<=300;j++)
17 dp[j+i*i]+=dp[j];
18 将 n拆分为 m个数(n个相同的球放入 m个盒子)的组合方法数<=>将 n拆分为最大数为 m相加的组合方法数
19 推论: 设 m<=n,n拆分成最多 m个数的和的拆分数等于将 n拆分成最大数不超过 m的拆分数
无限背包计数
  1 #include <bits/stdc++.h>
  2 typedef long long LL;
  3 const int N=1e6+10;
  4 int cas=1,T;
  5 int m,n;
  6 char s[40][40];
  7 LL dp[2][N];
  8 const int HASH=30007;
  9 int head[HASH],next[N],state[2][N],hn[2];
 10 inline int insert(int s,int *state,LL *dp)
 11 {
 12     int h = state[s] * 7 % HASH; //必要时乘个质数会散一些,如果还有一个数要变成(state[s]+x*1007)%HASH
 13     int u = head[h];
 14     while (u)
 15     {
 16         if (state[u] == state[s])
 17         {
 18             dp[u] += dp[s];
 19             return 0;
 20         }
 21         u = next[u];
 22     }
 23     next[s] = head[h];
 24     head[h] = s;
 25     return 1;
 26 }
 27 inline void ha(int x,int k,int d,int &h)
 28 {
 29     state[d ^ 1][h] = x;
 30     dp[d ^ 1][h] = dp[d][k];
 31     h += insert(h, state[d ^ 1], dp[d ^ 1]);
 32 }
 33 LL ans[40][40];
 34 int main()
 35 {
 36     while (scanf("%d%d", &n, &m) == 2)
 37     {
 38         for (int i = 0; i < n; i++)
 39             scanf("%s", s[i]);
 40         for (int i = 0; i < n; i++)
 41             for (int j = 0; j < m; j++)
 42                 s[i][j] = (s[i][j] == '.');
 43         int d = 0;
 44         int M = m << 1;
 45         memset(ans, 0, sizeof(ans));
 46         state[0][1] = 0;
 47         dp[0][1] = 1;
 48         hn[0] = 2;
 49         for (int i = 0; i < n; i++)
 50         {
 51 
 52             for (int j = 0; j < m; j++, d ^= 1)
 53             {
 54                 hn[d ^ 1] = 1;
 55                 if (!s[i][j]) //障碍
 56                 {
 57                     for (int k = 1, &h = hn[d ^ 1]; k < hn[d]; k++)
 58                     {
 59                         int t = state[d][k];
 60                         if ((t & (3 << M)) || (t & 3))
 61                             continue;
 62                         state[d ^ 1][h] = t << 2;
 63                         dp[d ^ 1][h] = dp[d][k];
 64                         h++;
 65                     }
 66                     continue;
 67                 }
 68                 memset(head, 0, sizeof(head));
 69                 for (int k = 1, &h = hn[d ^ 1]; k < hn[d]; k++) //每个非障碍格子都可以走
 70                 {
 71                     int t = state[d][k];
 72                     if ((t & (3 << M)) && (t & 3))
 73                     {
 74                         int x = t & 3;
 75                         int y = t & (3 << M);
 76                         t = t ^ x ^ y;
 77                         y >>= M;
 78                         if (x == 1 && y == 2)
 79                         {
 80                             if (t == 0)                //没有没有其它插头时可以是一条(i,j)结束的回路
 81                                 ans[i][j] += dp[d][k]; //()
 82                         }
 83                         else if (x == 2 && y == 1)
 84                             ha(t << 2, k, d, h);   //)(
 85                         else if (x == 2 && y == 2) //))
 86                         {
 87                             for (int l = 2, c = 0; l < M; l += 2) //往前找到匹配的(并改为)
 88                             {
 89                                 if (((t >> l) & 3) == 2)
 90                                     c--;
 91                                 else if (((t >> l) & 3) == 1)
 92                                     c++;
 93                                 if (c == 1)
 94                                 {
 95                                     t = t ^ (3 << l);
 96                                     break;
 97                                 }
 98                             }
 99                             ha(t << 2, k, d, h);
100                         }
101                         else
102                         {
103                             //((
104                             for (int l = M - 2, c = 0; l > 0; l -= 2) //往后找到匹配的)并改为(
105                             {
106                                 if (((t >> l) & 3) == 1)
107                                     c--;
108 
109                                 else if (((t >> l) & 3) == 2)
110                                     c++;
111                                 if (c == 1)
112                                 {
113                                     t = t ^ (3 << l);
114                                     break;
115                                 }
116                             }
117                             ha(t << 2, k, d, h);
118                         }
119                     }
120                     else if (!(t & (3 << M)) && !(t & 3)) //##
121                     {
122                         //if(s[i][j]) ha(t<<2,k,d,h);//不走
123                         if (i + 1 < n && s[i + 1][j] && j + 1 < m && s[i][j + 1])
124                             ha(t << 2 | 6, k, d, h);
125                     }
126                     else
127                     {
128                         //#( #) (# )#
129                         int x = t & 3;
130                         int y = t & (3 << M);
131                         if (x)
132                             t ^= x;
133                         else if (y)
134                             t ^= y, x = y >> M;
135                         if (i + 1 < n && s[i + 1][j])
136                             ha(t << 2 | (x << 2), k, d, h);
137                         if (j + 1 < m && s[i][j + 1])
138                             ha(t << 2 | x, k, d, h);
139                     }
140                 }
141             }
142         }
143         for (int i = n - 1; i >= 0; i--)
144         {
145             for (int j = m - 1; j >= 0; j--)
146                 if (s[i][j])
147                 {
148                     printf("%lld\n", ans[i][j]);
149                     i = -1;
150                     break;
151                 }
152         }
153     }
154     return 0;
155 }
//轮廓线 dp求只有一个闭合回路方法数(括号匹配判断是否形成回路)
 1 int n,a[N],idx[N];
 2 int lbig[N],lsmall[N],rsmall[N],rbig[N],bit[N],id[N];
 3 //lbig左边比当前数大,rbig右边比当前数大(不等于)
 4 LL sum(int x)
 5 {
 6     LL rec = 0;
 7     while (x)
 8     {
 9         rec += bit[x];
10         x -= x & -x;
11     }
12     return rec;
13 }
14 LL add(int x)
15 {
16     while (x <= n)
17     {
18         bit[x]++;
19         x += x & -x;
20     }
21 }
22 void cal()
23 {
24     for (int i = 1; i <= n; i++)
25         idx[i] = a[i];
26     sort(idx + 1, idx + n + 1);
27     memset(bit, 0, sizeof(bit));
28     map<int, int> mp; //用于处理相同数
29     for (int i = 1; i <= n; i++)
30     {
31         id[i] = lower_bound(idx + 1, idx + 1 + n, a[i]) - idx;
32         lsmall[i] = sum(id[i] - 1);
33         add(id[i]);
34         int &x = mp[a[i]];
35         lbig[i] = i - 1 - lsmall[i] - x;
36         x++;
37     }
38     memset(bit, 0, sizeof(bit));
39     mp.clear();
40     for (int i = n; i > 0; i--)
41     {
42         rsmall[i] = sum(id[i] - 1);
43         add(id[i]);
44         int &x = mp[a[i]];
45         rbig[i] = n - i - rsmall[i] - x;
46         x++;
47     }
48 }
//BIT统计逆序对
 1 LL count(LL n,LL num)
 2 {
 3     LL e = 1, ans = 0;
 4     for (int i = 0; n >= e; e *= 10, i++)
 5     {
 6         LL pre = n / e / 10;
 7         LL dig = n / e % 10;
 8         LL suf = n % e;
 9         if (dig > num)
10             pre++;
11         if (!num)
12             pre--; //num==0时去掉前导 0的情况
13         ans += pre * e;
14         if (dig == num)
15             ans += suf;
16     }
17     if (n && !num)
18         ans++; //num==0且 n>0时 0也要算上
19     return ans;
20 }
//数位 dp,数 1~n-1中的数位总共出现了多少个 num(0<=num<=9)
  1 #include <bits/stdc++.h>
  2 using namespace std;
  3 typedef long long LL;
  4 const int N=1e5+10;
  5 const int INF=0x3f3f3f3f;
  6 const int bit=7;
  7 int cas=1,T;
  8 int dp[10][3];
  9 void init()
 10 {
 11     memset(dp, 0, sizeof(dp));
 12     dp[0][2] = 1;
 13     for (int i = 1; i <= bit; i++)
 14     {
 15         dp[i][0] = dp[i - 1][0] * 10 + 2 * dp[i - 1][1] + dp[i - 1][2]; //62,4
 16         dp[i][1] = dp[i - 1][1] + dp[i - 1][2];                         //no62,4 head=2
 17         dp[i][2] = 7 * dp[i - 1][1] + 8 * dp[i - 1][2];                 //no62,4 head!=2
 18     }
 19     //for(int i=1;i<=bit;i++) printf("%d %d %d %d\n",i,dp[i][0],dp[i][1],dp[i][2]);
 20 }
 21 int count(int n)
 22 {
 23     int ans = 0;
 24     for (int i = bit, e = 1e6 + 0.5; i >= 0; i--, e /= 10)
 25     {
 26         if (e > n)
 27             continue;
 28         int pre = n / e % 10;  //当前位数字
 29         ans += pre * dp[i][0]; //0~pre-1 *e
 30         if (pre > 6 || pre == 6 && n * 10 / e % 10 > 2)
 31             ans += dp[i][1];
 32         //6+ no62,4 head=2
 33         当前位 > 6 或者 当前位 = 6且下一位 > 2
 34 
 35                                  if (pre > 4) ans += dp[i][2] + dp[i][1]; //4+ no62,4 head=2 + no62,4, head!=2
 36         if (pre == 4 || n / e % 100 == 62)                                //+ 后边所有
 37         {
 38             ans += n % e;
 39             break;
 40         }
 41     }
 42     return ans;
 43 }
 44 int main()
 45 {
 46     init();
 47     int l, r;
 48     //while(scanf("%d%d",&l,&r)==2 && (l&&r)) printf("%d %d\n",count(r+1),count(l));
 49     while (scanf("%d%d", &l, &r) == 2 && (l && r))
 50         printf("%d\n", r - l + 1 - count(r + 1) + count(l));
 51     return 0;
 52 }
 53 //数位 dp记忆化搜索,从高位往低位搜
 54 //这里是[l,r]中数位中连续奇数数字有偶数个,连续偶数数字有奇数个
 55 //第三维表示数位的数字,还没有数字为 2,奇数为 1,偶数为 0
 56 //第四维表示数位的连续个数,还没有数字为 2,奇数个为 1,偶数个为 0
 57 //第二维表示当前数位为 dig时有多少种
 58 #include <bits/stdc++.h>
 59 using namespace std;
 60 typedef long long LL;
 61 const int N=1e5+10;
 62 const int INF=0x3f3f3f3f;
 63 int cas=1,T;
 64 LL dp[20][10][3][3];//dp[a][b][c][d]记忆化前面数位为 d个 c,当前位为 b,后面有 a个 9时的种数
 65 LL dfs(int n,LL e,LL x,int odd,int even)
 66 {
 67     if (n < 0)
 68         return odd != even;
 69     if (e > x && odd == 2)
 70         return dfs(n - 1, e / 10, x, odd, even);
 71     int dig = x / e % 10;
 72     if ((x + 1) % e == 0 && dp[n][dig][odd][even] != -1)
 73         return dp[n][dig][odd][even]; //如果已经搜过直
 74     接返回
 75     if (odd == 2)
 76     {
 77         //还没有数位时不同转移
 78         LL ans = 0;
 79         ans += dfs(n - 1, e / 10, e - 1, 2, 2); //前导 0,还是没有数位
 80         for (int i = 1; i < dig; i++)
 81             ans += dfs(n - 1, e / 10, (i + 1) * e - 1, i & 1, 1);
 82         ans += dfs(n - 1, e / 10, x % e, dig & 1, 1);
 83         if ((x + 1) % e == 0)
 84             dp[n][dig][odd][even] = ans;
 85         return ans;
 86     }
 87     LL ans = 0;
 88     for (int i = 0; i < dig; i++)
 89         if ((i & 1) == odd || (i & 1) != odd && odd != even)
 90             ans += dfs(n - 1, e / 10, (i + 1) * e - 1, i & 1, ((i & 1) == odd) ? (even ^ 1) : 1);
 91 
 92     if ((dig & 1) == odd || (dig & 1) != odd && odd != even)
 93         ans += dfs(n - 1, e / 10, x % e, dig & 1, ((dig & 1) == odd) ? (even ^ 1) : 1);
 94     if ((x + 1) / 10 == e)
 95         dp[n][dig][odd][even] = ans; //x=a99…99 时记忆化
 96     return ans;
 97 }
 98 int main()
 99 {
100     memset(dp, -1, sizeof(dp));
101     scanf("%d", &T);
102     while (T--)
103     {
104         LL l, r;
105         scanf("%lld%lld", &l, &r);
106         LL lans = dfs(18, 1e18, l - 1, 2, 2);
107         //printf("lans:%lld\n",lans);
108         LL rans = dfs(18, 1e18, r, 2, 2);
109         //printf("rans:%lld\n",rans);
110         printf("Case #%d: %lld\n", cas++, rans - lans);
111     }
112     return 0;
113 }
//数位 dp模板(递推过程),考虑前导 0是注意 0~9999…9999 这种的处理 //这里是算[l,r]中没有 62也没有 4的个数,从低位往高位推,每次在最高位加一个数字进行转移
TSP问题:有 n个城市,两两之间均有道路相连,距离为 dis(i,j),求最短的路径经过每个点一次并回到原点
做法:状压 dp,O(n*2^n)
如果给定起点和终点 s,e,则可以用 O(n^2)dp 解决,每次在当前最小路径中插入下一个点获得加一个点进来后的最
小路径
//圆周涂色
把一个圆分成 n个扇形,每个扇形可以涂 m种颜色,要求相邻的颜色互不相同,有多少种涂法?
把 n个扇形的涂法记为 f(n),则 f(n)+f(n-1)=m*(m-1)^(n-1),包含第一个和第 n个扇形涂色相同的情况,相当于
将第 1个扇形和第 n个扇形合并成一个扇形的涂法,化简得 f(n)=(m-1)^n+(n&1?-1:1)*(m-1);
 1 const int M=75;
 2 int
 3 ro[3][M]={{0,10,11,12,13,14,15,16,17,18,46,47,48,49,50,51,52,53,54,25,22,19,26,23,20,27,24,
 4 21,9,8,7,6,5,4,3,2,1,39,42,45,38,41,44,37,40,43,36,35,34,33,32,31,30,29,28,//3*3*6=54个方格
 5 63,59,55,62,64,56,58,66,65,60,57,61,
 6 //十二条边
 7 71,72,68,67,74,73,69,70},//8个顶点
 8 {0,43,40,37,44,41,38,45,42,39,16,13,10,17,14,11,18,15,12,7,4,1,8,5,2,9,6,3,30,33,36,29,32,3
 9 5,28,31,34,52,49,46,53,50,47,54,51,48,25,22,19,26,23,20,27,24,21,
10 62,58,61,66,55,57,65,63,59,56,60,64,
11 71,67,70,74,72,68,69,73},
12 {0,7,4,1,8,5,2,9,6,3,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,4
13 2,43,44,45,10,11,12,13,14,15,16,17,18,48,51,54,47,50,53,46,49,52,
14 56,57,58,55,60,61,62,59,64,65,66,63,
15 68,69,70,67,72,73,74,71}};
16 
17  
18 int pu[M],t[M],pemu[100][M],pn;
19 void count(int *v)
20 {
21     int vis[M];
22     memset(vis, 0, sizeof(vis));
23     for (int i = 1; i < M; i++)
24         if (!vis[i])
25         {
26             int c = 0;
27             int u = i;
28             while (!vis[u])
29             {
30                 vis[u] = 1;
31                 u = pu[u];
32                 c++;
33             }
34             v[c]++;
35         }
36 }
37 void init()
38 {
39     memset(pemu, 0, sizeof(pemu));
40     pn = 0;
41     set<LL> s;
42     map<int, int> cc;
43     for (int i = 0; i < M; i++)
44         pu[i] = i;
45     for (int i = 1; i <= 4; i++)
46     {
47         for (int l = 1; l < M; l++)
48             t[l] = pu[ro[0][l]];
49         memcpy(pu, t, sizeof(pu));
50         for (int j = 1; j <= 4; j++)
51         {
52             for (int l = 1; l < M; l++)
53                 t[l] = pu[ro[1][l]];
54             memcpy(pu, t, sizeof(pu));
55             for (int k = 1; k <= 4; k++)
56             {
57                 for (int l = 1; l < M; l++)
58                     t[l] = pu[ro[2][l]];
59                 memcpy(pu, t, sizeof(pu));
60                 //printf("i:%d j:%d k:%d\n",i,j,k);
61                 LL ha = 0;
62                 for (int l = 1; l < M; l++)
63                     ha = ha * 1000007 + pu[l];
64                 //for(int l=1;l<=12;l++) printf("%d ",pu[l]);printf("\n");
65                 if (s.count(ha))
66                     continue;
67                 s.insert(ha);
68                 count(pemu[pn]);
69                 int sum = 0;
70                 for (int l = 1; l < M; l++)
71                     printf("%d ", pemu[pn][l]), sum += pemu[pn][l];
72                 printf("\n");
73                 cc[sum]++;
74                 pn++;
75             }
76         }
77     }
78     for (map<int, int>::iterator it = cc.begin(); it != cc.end(); it++)
79         printf("%d %d\n", it->first, it->second);
80     printf("%d\n", s.size());
81 }
//置换打表
 1 //f(n)=a*x^n+b*x^n,得到 f(n)=c1*f(n-1)+c2*f(n-2),a*x^n=f(n)-b*x^n
 2 //求 a^n%mod用矩阵快速幂求解 f(n)时 f(n)%mod的循环节为(mod*mod-1),最小循环节是(mod*mod-1)的约数
 3 //即 f(n)%mod=f(n%(mod*mod-1))%mod
 4 //当循环节为 p时 f(0)==f(p) && f(1)==f(p+1)
 5 //可以暴力找循环节
 6 const int M=2;
 7 struct mar
 8 {
 9     LL a[M][M];
10     mar operator*(const mar &b)
11     {
12         mar ans;
13         memset(&ans, 0, sizeof(ans));
14         for (int i = 0; i < M; i++)
15             for (int j = 0; j < M; j++)
16                 for (int k = 0; k < M; k++)
17                     ans.a[i][j] = (ans.a[i][j] + a[i][k] * b.a[k][j]) % mod;
18         return ans;
19     }
20 };
21 LL a[M]={2,1};//递推式 f(n)=a[0]f(n-1)+a[1]f(n-2)+...+a[M-1]f(n-M)
22 LL f[M]={1,0};//f[0]=f(M-1),f[1]=f(M-2),...,f(M-2)=f(1),f(M-1)=f(0) 倒过来的
23 LL mar_pow(LL n)
24 {
25     mar ans, t;
26     memset(&ans, 0, sizeof(ans));
27     memset(&t, 0, sizeof(t));
28     for (int i = 0; i < M; i++)
29         ans.a[i][i] = 1; //单位矩阵,相当于快速幂中初始化为 1
30     memcpy(t.a[0], a, sizeof(a));
31     for (int i = 1; i < M; i++)
32         t.a[i][i - 1] = 1;
33     while (n)
34     {
35         if (n & 1)
36             ans = ans * t;
37         t = t * t;
38         n >>= 1;
39     }
40     LL res = 0;
41     for (int i = 0; i < M; i++)
42         res = (res + ans.a[M - 1][i] * f[i]) % mod;
43     //f(n)=f(n-1)*ans.a[M-1][i]+...+f(0)*ans.a[M-1][0];
44     return res;
45 }
//通用矩阵快速幂
 1 void dfs(int x)
 2 {
 3     int pos = upper_bound(b + 1, b + 1 + n, k / a[x]) - b - 1;
 4     ans -= sum(pos);
 5     for (int i = 0; i < u[x].size(); i++)
 6         dfs(u[x][i]);
 7     ans += sum(pos);
 8     add(upper_bound(b + 1, b + 1 + n, a[x]) - b - 1);
 9 }
10 struct node
11 {
12     int i, x, pos;
13 };
14 void slove(int x)//由上面的 dfs改成的手写栈
15 {
16     stack<node> s;
17     int pos = upper_bound(b + 1, b + 1 + n, k / a[x]) - b - 1;
18     s.push(node{0, x, pos});
19     while (!s.empty())
20     {
21         node v = s.top();
22         s.pop();
23         if (v.i == u[v.x].size()) //dfs结束时
24         {
25             ans += sum(v.pos);
26             add(upper_bound(b + 1, b + 1 + n, a[v.x]) - b - 1);
27             continue; //return
28         }
29         x = u[v.x][v.i++]; //子树
30         s.push(v);
31         pos = upper_bound(b + 1, b + 1 + n, k / a[x]) - b - 1;
32         ans -= sum(pos);
33         s.push(node{0, x, pos}); //子树节点入栈
34     }
35 }
//dfs改手写栈
求<=n中除 p外所有质数的乘积%p
要解决这个问题首先要理解一个动态规划的定义，定义 S(v,p)为 1−v 中经过小于等于 p 的素数筛过之后的数的个数，
即此时整个数集里只包含<=p 的素数以及素因子>p 的数。
相应的有状态转移 S(v,p)=S(v,p−1)–(S(v/p,p−1)−S(p−1,p−1))
相应的可以定义状态 F(v,p)为 1−v 中经过小于等于 p 的素数筛过之后的数的乘积。相应的转移方程为：
F(v,p)=F(v,p−1)/(F(v/p,p−1)/F(p−1,p−1)∗pS(v/p,p−1)−S(p−1,p−1))
原题中要模的是 M，所以在定义 S(v,p)和相应的 F(v,p)时应当改为 1−v 中（除了 M,2M…）经过小于等于 p 的素数筛过
之后的数的个数，以及这些数的乘积，这个实际上影响到的是初始化的阶段。
实现上由于 v 的取值只需要用到 n/1,n/2…的值，所以分别用 l,h 的两种数组记录，其中 cl[x]表示的 S(x,p),ch[x]表示
S(n/x,p)，pl[x]表示 F(x,p)，ph[x]表示 F(n/x,p)。
 1 syntax on
 2 "高亮
 3 set nu
 4 "显示行号
 5 set softtabstop=4
 6 set shiftwidth=4
 7 set smartindent
 8 set scrolloff=7
 9 set mouse=a
10 "按 tab缩进 4个空格大小
11 "默认缩进 4个空格大小（indent的缩进）
12 "自动缩进
13 "光标移动到 buffer的顶部和底部时保持 7行距离
14 "鼠标所有功能可用
15 set cursorline
16 "当前行下划线划出
17 "出现 swap文件时用这个
18 set nobackup
19 set noswapfile
20 winpos -4 -4
21 "设置窗口位置,使用 vim时自动最大化窗口
22 set lines=50 columns=200 "设定窗口大小 "按 F5编译运行
23 set guifont=UbuntuMono\ 14
24 colorscheme desert
25 "Ubuntu下 GVIM设置英文字体,字体名和字体大小
26 "F5编译运行
27 "w保存文本 ! g++ % -o %< 编译代码 ! ./%< linux执行 ! %<.exe windows执行
28 map <F5> :call Compile()<CR>
29 func Compile()
30 :w
31 :! g++ -std=c++11 % -o %<
32 :! ./%<
33 endfunc
34 "noremap insert模式（打字时的模式）下映射
35 "vmap/vnoremap visaul模式（选择时的模式）下映射
36 "map/nnoremap normal模式（按下 Esc后的模式）下映射
37 map <C-a> ggVG
38 vmap <C-c> "+y
39 map <C-v> "+p
40 map zz I//<Esc>
41 "Ctrl+a全选
42 "Ctrl+c选择后复制，可粘贴到外编辑器，如果只在 vim中粘贴按 y复制即可
43 "Ctrl+v normal模式下粘贴外部复制的文本
44 "zz行首注释
45 map tt I<Right><Right><Backspace><Backspace><Esc> "tt取消注释
46 "cpp文件自动输出，F2手动输出
47 auto BufNewFile *.cpp exec ":call Preset()"
48 map<F2> :call Preset()<CR>
49 func Preset()
50 "setline将当前行设为字符串，append在当前行后追加
51 call setline(line("."),"#include <bits/stdc++.h>")
52 let l=line(".")-1
53 let l=l+1 | call append(l,"using namespace std;")
54 let l=l+1 | call append(l,"typedef long long LL;")
55 let l=l+1 | call append(l,"const int N=1e5+10;")
56 let l=l+1 | call append(l,"const int INF=0x3f3f3f3f;")
57 
58  
59 let l=l+1 | call append(l,"const int mod=1e9+7;")
60 let l=l+1 | call append(l,"int cas=1,T;")
61 let l=l+1 | call append(l,"int main()")
62 let l=l+1 | call append(l,"{")
63 let l=l+1 | call append(l," //freopen(\"in\",\"r\",stdin);")
64 let l=l+1 | call append(l," //scanf(\"%d\",&T);")
65 let l=l+1 | call append(l," return 0;")
66 let l=l+1 | call append(l,"}")
67 endfunc
68 """""""""""下面的不需要用到就不要打了
69 "F3输出头文件,bits/stdc++.h可用时不需要
70 map <F3> :call AddHeader()<CR>
71 func! AddHeader()
72 call setline(line("."),"#include <cstdio>")
73 let l=line(".")-1
74 let l=l+1 | call append(l,"#include <queue>")
75 let l=l+1 | call append(l,"#include <cstring>")
76 let l=l+1 | call append(l,"#include <string>")
77 let l=l+1 | call append(l,"#include <iostream>")
78 let l=l+1 | call append(l,"#include <cstdlib>")
79 let l=l+1 | call append(l,"#include <algorithm>")
80 let l=l+1 | call append(l,"#include <vector>")
81 let l=l+1 | call append(l,"#include <map>")
82 let l=l+1 | call append(l,"#include <set>")
83 let l=l+1 | call append(l,"#include <ctime>")
84 let l=l+1 | call append(l,"#include <cmath>")
85 let l=l+1 | call append(l,"#include <cctype>")
86 endfunc
87 "set guifont=Consolas:h14:cANSI
88 "windows下设置英文字体,字体名和字体大小
"Vim配置,引号开头的是注释，不用打
 1 1）、单点增减+区间求和
 2 思路：C[x]表示该点的元素:sum(x)=C[1]+C[2]+……C[x]
 3 int arr[MAXN];
 4 inline int sum(int x){
 5     int res = 0;
 6     while (x)
 7         res += arr[x], x -= lowbit(x);
 8     return res;}
 9 inline void add(int x,int n){
10     while (x < MAXN)
11         arr[x] += n, x += lowbit(x);}
12 inline int query(int x,int y){
13     return sum(y) - sum(x - 1);}
14 （2）、区间增减+单点查询
15 思路：C[x]表示该点元素与左边元素的差值:num[x]=C[1]+C[2]+……C[x]
16 int arr[MAXN]
17 inline int sum(int x){
18     int res = 0;
19     while (x)
20         res += arr[x], x -= lowbit(x);
21     return res;}
22 inline void add(int x,int n){
23     while (x < MAXN)
24         arr[x] += n, x += lowbit(x);}
25 inline int update(int x,int y,int n){
26     add(x, n);
27     add(y + 1, -n);}
28 （3）、区间增减+区间查询
29 思路：C1[x]表示该点元素与左边的差值，C2[x]表示的是 x*C[x]
30 sum(sum(C[j],j<=i)i<=x)
31 = x*C[1]+(x-1)*C[2]+……+C[x]
32 =(x+1)*sum(C[i],i<=x)-sum(i*C[i],i<=x);
33 则可以想到用 C1[x]维护 C[x]的值,C2[x]维护 x*C[X]的值
34 template <typename X>
35 struct tree_array{
36     struct tree_array_single
37     {
38         X arr[MAXN];
39 
40         void add(int x, X n)
41         {
42             while (x <= N)
43                 arr[x] += n, x += lowbit(x);
44         }
45         X sum(int x)
46         {
47             X sum = 0;
48             while (x)
49                 sum += arr[x], x -= lowbit(x);
50             return sum;
51         }
52     } T1, T2;
53     void reset()
54     {
55         CLR(T1.arr, 0);
56         CLR(T2.arr, 0);
57     }
58     void add(int x, X n)
59     {
60         T1.add(x, n);
61         T2.add(x, x * n);
62     }
63     void update(int L, int R, int n)
64     {
65         add(L, n);
66         add(R + 1, -n);
67     }
68     X sum(int x) { return (x + 1) * T1.sum(x) - T2.sum(x); }
69     X query(int L, int R) { return sum(R) - sum(L - 1); }
70 };
14、树状数组
  1 ①单点增减（add） + 矩形求和（query）
  2 ②矩形增减（update） + 单点求值（sum）
  3 int arr[MAXN][MAXN]
  4 inline void add(int x,int y,int n) {
  5     for (int i = x; i < MAXN; i += lowbit(i))
  6         for (int j = y; j < MAXN; j += lowbit(j))
  7             arr[i][j] += n;
  8 }
  9 inline int sum(int x,int y){
 10     int res = 0;
 11     for (int i = x; i; i -= lowbit(i))
 12         for (int j = y; j; j -= lowbit(j))
 13             res += arr[i][j];
 14     return res;
 15 }
 16 inline int query(int L,int B,int R,int T) {
 17     return sum(R, T) + sum(L - 1, B - 1) - sum(R, B - 1) - sum(L - 1, T);
 18 }
 19 inline void update(int L,int B,int R,int T,int n){
 20     update(L, B, n);
 21     update(L, T + 1, n);
 22     update(R + 1, B, n);
 23     update(R + 1, T + 1, n);
 24 }
 25 ③矩形增减（update）+ 矩形求和（query）
 26 template<typename X>
 27 class tree_array{
 28     struct tree_array_single
 29     {
 30         X arr[MAXN][MAXN];
 31         void add(int x, int y, X n)
 32         {
 33             for (int i = x; i < MAXN; i += lowbit(i))
 34                 for (int j = y; j < MAXN; j += lowbit(j))
 35                     arr[i][j] += n;
 36         }
 37         X sum(int x, int y)
 38         {
 39             X res = 0;
 40             for (int i = x; i; i -= lowbit(i))
 41                 for (int j = y; j; j -= lowbit(j))
 42 
 43                     res += arr[i][j];
 44             return res;
 45         }
 46     } T1, T2, T3, T4;
 47     void add(int x, int y, int n)
 48     {
 49         T1.add(x, y, n);
 50         T2.add(x, y, y * n);
 51         T3.add(x, y, x * n);
 52         T4.add(x, y, x * y * n);
 53     }
 54     X sum(int x, int y)
 55     {
 56         return (x + 1) * (y + 1) * T1.sum(x, y) - (x + 1) * T2.sum(x, y) - (y + 1) * T3.sum(x, y) + T4.sum(x, y);
 57     }
 58 
 59   public:
 60     void init()
 61     {
 62         CLR(T1.arr, 0);
 63         CLR(T2.arr, 0);
 64         CLR(T3.arr, 0);
 65         CLR(T4.arr, 0);
 66     }
 67     void update(int L, int B, int R, int T, int n)
 68     {
 69         add(L, B, n);
 70         add(L, T + 1, -n);
 71         add(R + 1, B, -n);
 72         add(R + 1, T + 1, n);
 73     }
 74     X query(int L, int B, int R, int T)
 75     {
 76         return sum(R, T) - sum(L - 1, T) - sum(R, B - 1) + sum(L - 1, B - 1);
 77     }
 78 };
 79 ④单点增减（add） + 立方体求和（query）
 80 ⑤立方体增减（update） + 单点求值（sum）
 81 int arr[MAXN][MAXN][MAXN];
 82 inline int sum(int x,int y,int z){
 83     int res = 0;
 84     for (int i = x; i; i -= lowbit(i))
 85         for (int j = y; j; j -= lowbit(j))
 86             for (int k = z; k; k -= lowbit(k))
 87                 res ^= arr[i][j][k];
 88     return res;
 89 }
 90 inline void add(int x,int y,int z,int n){
 91     for (int i = x; i < MAXN; i += lowbit(i))
 92         for (int j = y; j < MAXN; j += lowbit(j))
 93             for (int k = z; k < MAXN; k += lowbit(k))
 94                 arr[i][j][k] += n;
 95 }
 96 inline void update(int x1,int y1,int z1,int x2,int y2,int z2,int n){
 97     add(x1, y1, z1, n);
 98     add(x2 + 1, y1, z1, -n);
 99     add(x1, y2 + 1, z1, -n);
100     add(x1, y1, z2 + 1, -n);
101     add(x2 + 1, y2 + 1, z1, n);
102     add(x2 + 1, y1, z2 + 1, n);
103     add(x1, y2 + 1, z2 + 1, n);
104     add(x2 + 1, y2 + 1, z2 + 1, -n);
105 }
106 inline int query(int x1,int y1,int z1,int x2,int y2,int z2){
107     return sum(x2, y2, z2) - sum(x2, y2, z1 - 1) - sum(x2, y1 - 1, z2) - sum(x1 - 1, y2, z2) + sum(x2, y1 - 1, z1 - 1) + sum(x1 - 1, y2, z1 - 1) + sum(x1 - 1, y1 - 1, z2) - sum(x1 - 1, y1 - 1, z1 - 1);
108 }
109 
110  
111 ⑥立方体增减（update） + 立方体求和（query）///随便写写……复杂度较高
112 template<typename X>
113 class tree_array_Cube{
114     struct tree_array_single
115     {
116         X arr[MAXN][MAXN][MAXN];
117         X sum(int x, int y, int z)
118         {
119             X res = 0;
120             for (int i = x; i; i -= lowbit(i))
121                 for (int j = y; j; j -= lowbit(j))
122                     for (int k = z; k; k -= lowbit(k))
123                         res += arr[i][j][k];
124             return res;
125         }
126         void add(int x, int y, int z, X n)
127         {
128             for (int i = x; i < MAXN; i += lowbit(i))
129                 for (int j = y; j < MAXN; j += lowbit(j))
130                     for (int k = z; k < MAXN; k += lowbit(k))
131                         arr[i][j][k] += n;
132         }
133     } T1, T2, T3, T4, T5, T6, T7, T8;
134     void add(int x, int y, int z, X n)
135     {
136         T1.add(x, y, z, n);
137         T2.add(x, y, z, x * n);
138         T3.add(x, y, z, y * n);
139         T4.add(x, y, z, z * n);
140         T5.add(x, y, z, x * y * n);
141         T6.add(x, y, z, y * z * n);
142         T7.add(x, y, z, x * z * n);
143         T8.add(x, y, z, x * y * z * n);
144     }
145     X sum(int x, int y, int z)
146     {
147         return (x + 1) * (y + 1) * (z + 1) * T1.sum(x, y, z) - (y + 1) * (z + 1) * T2.sum(x, y, z) - (x + 1) * (z + 1) * T3.sum(x, y, z) - (x + 1) * (y + 1) * T4.sum(x, y, z) + (z + 1) * T5.sum(x, y, z) + (x + 1) * T6.sum(x, y, z) + (y + 1) * T7.sum(x, y, z) - T8.sum(x, y, z);
148     }
149 
150   public:
151     void init()
152     {
153         CLR(T1.arr, 0);
154         CLR(T2.arr, 0);
155         CLR(T3.arr, 0);
156         CLR(T4.arr, 0);
157         CLR(T5.arr, 0);
158         CLR(T6.arr, 0);
159         CLR(T7.arr, 0);
160         CLR(T8.arr, 0);
161     }
162     void update(int x1, int y1, int z1, int x2, int y2, int z2, X n)
163     {
164         add(x1, y1, z1, n);
165         add(x2 + 1, y1, z1, -n);
166         add(x1, y2 + 1, z1, -n);
167         add(x1, y1, z2 + 1, -n);
168         add(x2 + 1, y2 + 1, z1, n);
169         add(x2 + 1, y1, z2 + 1, n);
170         add(x1, y2 + 1, z2 + 1, n);
171         add(x2 + 1, y2 + 1, z2 + 1, -n);
172     }
173     X query(int x1, int y1, int z1, int x2, int y2, int z2)
174     {
175         return sum(x2, y2, z2) - sum(x2, y2, z1 - 1) - sum(x2, y1 - 1, z2) - sum(x1 - 1, y2, z2) + sum(x2, y1 - 1, z1 - 1) + sum(x1 - 1, y2, z1 - 1) + sum(x1 - 1, y1 - 1, z2) - sum(x1 - 1, y1 - 1, z1 - 1);
176     }
177 };
15、多维树状数组
 1 inline void init()
 2 {
 3     CLR(arr, 0);
 4     for (int i = 1; i <= N; ++i)
 5         for (int j = i; j <= N && arr[j] < num[i]; j += lowbit(j))
 6             arr[j] = num[i];
 7 }
 8 inline int query(int L,int R)
 9 {
10     int res = 0;
11     for (--L; L < R;)
12     {
13         if (R - lowbit(R) >= L)
14         {
15             res = max(res, arr[R]);
16             R -= lowbit(R);
17         }
18         else
19         {
20             res = max(res, num[R]);
21             --R;
22         }
23     }
24     return res;
25 }
26 inline void update(int x,int val)
27 {
28     int ori = num[x];
29     num[x] = val;
30     if (val >= ori)
31         for (int i = x; i <= N && arr[i] < val; i += lowbit(i))
32             arr[i] = val;
33     else
34     {
35         for (int i = x; i <= N && arr[i] == ori; i += lowbit(i))
36         {
37             arr[i] = val;
38             for (int j = lowbit(i) >> 1; j; j >>= 1)
39                 arr[i] = max(arr[i], arr[i - j]);
40         }
41     }
42 }
16、树状数组—区间最大值

 

posted @ 2017-05-08 23:21  cdongyang  阅读(289)  评论(0)  编辑  收藏  举报