CSP2019前夕整理一下模板，顺便供之后使用
0. 非算法内容
0.1. 读入优化
描述:
使用getchar()实现的读入优化。
代码:
inline int read()
{
int x=0; bool f=1; char c=getchar();
for(;!isdigit(c);c=getchar()) if(c=='-') f=0;
for(; isdigit(c);c=getchar()) x=x*10+(c^'0');
if(f) return x;
return -x;
}
0.2. 高级读入优化
描述:
使用fread()实现的读入优化。
代码:
暂无
1. 数据结构
1.1. 虚树
描述:
给定树上的\(k\)个关键点，构建出一棵虚树，只有关键点和任意两个关键点的LCA会被保留，且原树上的祖先关系和虚树上祖先关系保持一致。可以证明虚树最多有\((2k-1)\)个点。
把所有关键点按DFS序排序，用一个栈维护动态加点连边即可。时间复杂度\(O(k\log n)\).
注意事项:
一定要处理好最先加入栈中的点(所有点的LCA).
代码:
addedge0_(u,v)表示在\(u\)和\(v\)之间连边，边权根据实际情况确定。
bool cmp(int x,int y) {return dfn[x]<dfn[y];}
void build()
{
sort(ky+1,ky+kyn+1,cmp); rt = ky[1]; for(int i=2; i<=kyn; i++) rt = LCA(rt,ky[i]);
tp = 1; stk[tp] = rt; n0++; kid[n0] = rt; isky[ky[1]] = true;
for(int i=(rt==ky[1]?2:1); i<=kyn; i++)
{
int u = ky[i],lca = LCA(stk[tp],u);
if(lca==stk[tp]) {tp++; stk[tp] = u; n0++; kid[n0] = u; isky[u] = true;}
else
{
while(tp>1 && dep[stk[tp-1]]>dep[lca]) {addedge0_(stk[tp],stk[tp-1]); tp--;}
addedge0_(stk[tp],lca); tp--;
if(stk[tp]!=lca)
{
tp++; stk[tp] = lca; n0++; kid[n0] = lca;
}
tp++; stk[tp] = u; n0++; kid[n0] = u; isky[u] = true;
}
}
while(tp>1)
{
addedge0_(stk[tp],stk[tp-1]);
tp--;
}
}
1.2. 左偏树
描述:
可并堆的一种实现，支持插入、删除、合并、取最小值等堆操作，单个操作时间复杂度均不超过\(O(\log n)\).
代码:
int merge(int u,int v)
{
if(u==0||v==0) return u+v;
if(val[u]>val[v]) swap(u,v);
son[u][1] = merge(son[u][1],v);
if(dis[son[u][1]]>dis[son[u][0]]) {swap(son[u][1],son[u][0]);}
dis[u] = dis[son[u][1]]+1;
return u;
}
void insert(int u,int x)
{
siz++; val[siz] = x;
rtn[u] = merge(rtn[u],siz);
}
int popnode(int u)
{
return merge(son[u][0],son[u][1]);
}
1.3. KD树
暂无
1.4. 李超线段树
描述: 维护线段树支持如下操作: 在区间内插入一条直线，询问区间内所有直线的最高点。时间复杂度\(O(n\log^2n)\).
代码:
暂无
2. 字符串
2.1. 后缀数组
描述:
倍增法求后缀数组，时间复杂度\(O(n\log n)\).
sa[i]: 排在第\(i\)位的后缀
rk[i]: 后缀\(i\)的排名
h[i]: 后缀\(i\)与其前一名的后缀的LCP。
height[i]: 第\(i\)名与其前一名的后缀的LCP.
注意:
(1) w1处不可令h[1]=1然后从\(2\)开始循环。
代码:
namespace SA
{
int height[N+3],h[N+3],tmp[N+3],st[lgN+2][N+3],buc[N+3],rk[N+3],sa[N+3];
void buildSA()
{
int *x = rk,*y = tmp;
for(int i=1; i<=S; i++) buc[i] = 0;
for(int i=1; i<=n; i++) buc[x[i]=a[i]]++;
for(int i=1; i<=S; i++) buc[i] += buc[i-1];
for(int i=n; i>=1; i--) sa[buc[x[i]]--] = i;
int p = 0,s = S;
for(int j=1; p<n; j<<=1)
{
p = 0;
for(int i=n-j+1; i<=n; i++) y[++p] = i;
for(int i=1; i<=n; i++) {if(sa[i]>j) y[++p] = sa[i]-j;}
for(int i=1; i<=s; i++) buc[i] = 0;
for(int i=1; i<=n; i++) buc[x[y[i]]]++;
for(int i=1; i<=s; i++) buc[i] += buc[i-1];
for(int i=n; i>=1; i--) sa[buc[x[y[i]]]--] = y[i];
p = 1; swap(x,y); x[sa[1]] = 1;
for(int i=2; i<=n; i++) x[sa[i]] = y[sa[i]]==y[sa[i-1]]&&y[sa[i]+j]==y[sa[i-1]+j]?p:++p;
s = p;
}
for(int i=1; i<=n; i++) rk[sa[i]] = i;
for(int i=1; i<=n; i++) //w1
{
h[i] = h[i-1]==0?0:h[i-1]-1;
while(i+h[i]<=n && sa[rk[i]-1]+h[i]<=n && a[i+h[i]]==a[sa[rk[i]-1]+h[i]]) {h[i]++;}
}
for(int i=1; i<=n; i++) height[i] = h[sa[i]];
for(int i=1; i<=n; i++) st[0][i] = height[i];
for(int j=1; j<lgN; j++)
{
for(int i=1; i+(1<<j)-1<=n; i++) {st[j][i] = min(st[j-1][i],st[j-1][i+(1<<j-1)]);}
}
}
int querymin(int l,int r)
{
int g = lg2[r-l+1]; int *adr = st[g];
return min(adr[l],adr[r-(1<<g)+1]);
}
int LCP(int x,int y)
{
if(x==y) return n-x+1; if(rk[x]>rk[y]) swap(x,y);
return querymin(rk[x]+1,rk[y]);
}
}
using SA::buildSA;
using SA::LCP;
2.2. 扩展KMP
描述:
求出一个串的每个后缀与整个串的LCP. 与Manacher算法极其类似，时间复杂度\(O(n)\).
z[i]: 后缀\(i\)与母串的LCP长度。
代码:
void Z_box()
{
int mx = 0,mxid = 0; z[1] = m;
for(int i=2; i<=m; i++)
{
if(a[i]!=a[1]) {z[i] = 0;}
else if(i>mx) {z[i] = 1;}
else {z[i] = min(z[i-mxid+1],mx-i+1);}
while(i+z[i]<=m && a[i+z[i]]==a[z[i]+1]) {z[i]++;}
if(i+z[i]-1>mx) {mx = i+z[i]-1; mxid = i;}
}
}
2.3. Manacher算法
描述:
求出以每个位置为中心的最长回文子串，时间复杂度\(O(n)\).
\(p_i\)表示以\(i\)为中心的最长回文串长度。
注意事项:
(1) 不要把\(2n+1\)写成\(n\).
代码:
void manacher()
{
for(int i=n; i>=1; i--) a[2*i] = a[i];
for(int i=1; i<=2*n+1; i+=2) a[i] = '#';
int mxid = 1,mx = 1; p[1] = 1;
for(int i=2; i<=2*n+1; i++)
{
if(i>mx) {p[i] = 1;}
else {p[i] = min(p[2*mxid-i],mx-i+1);}
while(i-p[i]>0 && i+p[i]<=2*n+1 && a[i+p[i]]==a[i-p[i]]) {p[i]++;}
if(i+p[i]-1>mx) {mx = i+p[i]-1; mxid = i;}
}
}
2.4. 后缀自动机
描述:
增量法构造后缀自动机，时间复杂度\(O(n)\).
注意广义后缀自动机必须对Trie树进行BFS建立才能保证复杂度为节点数，否则复杂度退化为所有叶子节点总深度。
注意事项:
(1) 不要忘记初始化三个变量。
代码:
void init()
{
siz = lstpos = rtn = 1;
}
void insertchar(char ch)
{
int p = lstpos,np; siz++; np = lstpos = siz; len[np] = len[p]+1; sz[np]++;
for(; p && son[p][ch]==0; p=fa[p]) {son[p][ch] = np;}
if(p==0) {fa[np] = rtn;}
else
{
int q = son[p][ch];
if(len[q]==len[p]+1) {fa[np] = q;}
else
{
siz++; int nq = siz; len[nq] = len[p]+1;
memcpy(son[nq],son[q],sizeof(son[q]));
fa[nq] = fa[q]; fa[q] = fa[np] = nq;
for(; p && son[p][ch]==q; p=fa[p]) {son[p][ch] = nq;}
}
}
}
2.5. 回文自动机
描述:
一个字符串本质不同的回文子串个数不超过\(n\). 回文自动机上一个节点代表一个回文子串。
增量法构造回文自动机，时间复杂度\(O(n)\).
代码:
void initPAM()
{
siz = 1; fail[0] = fail[1] = 1; len[0] = 0; len[1] = -1; lstpos = 1;
}
void insertchar(int id)
{
int p = lstpos;
while(a[id-1-len[p]]!=a[id]) {p = fail[p];}
if(!son[p][a[id]])
{
siz++; int u = siz,v = fail[p];
while(a[id-1-len[v]]!=a[id]) {v = fail[v];}
fail[u] = son[v][a[id]]; len[u] = len[p]+2; son[p][a[id]] = u;
}
lstpos = son[p][a[id]];
}
2.6. Lyndon分解
描述:
若一个串的最小循环表示为它本身，则称作Lyndon串。
将一个串划分为若干字典序不增的Lyndon串，称作Lyndon划分。一个串的Lyndon划分方案唯一。
Lyndon划分的Duval算法，时间复杂度\(O(n)\).
代码:
void Duval()
{
for(int i=1; i<=n;)
{
int j = i,k = i+1;
while(k<=n && a[j]<=a[k])
{
if(a[j]<a[k]) {j = i-1;}
j++; k++;
}
while(i<=j)
{
r[++cnt] = i+k-j-1;
i+=k-j;
}
}
}
2.7. 最小循环表示
暂无
2.8. 双模 Hash
仅用于线上比赛省下敲这几行代码的时间。
const pll P = mkpr(193000403ll,880230293ll);
pll E,pwE[mxN+3];
llong randw() {return (rand()<<15)|rand();}
pll operator +(const pll &x,const pll &y) {pll ret(x.first+y.first-P.first,x.second+y.second-P.second); ret.first+=(ret.first>>31)&P.first,ret.second+=(ret.second>>31)&P.second; return ret;}
pll operator -(const pll &x,const pll &y) {pll ret(x.first-y.first,x.second-y.second); ret.first+=(ret.first>>31)&P.first,ret.second+=(ret.second>>31)&P.second;}
pll operator *(const pll &x,const pll &y) {return mkpr(x.first*y.first%P.first,x.second*y.second%P.second);}
3. 图论
3.1. Tarjan算法
3.1.1. 强连通分量
描述:
Tarjan求有向图的强连通分量。构建DFS树，统计每个点的DFS时间戳(dfn[u])以及其所到达的时间戳最小的点(low[u]), 若后者为该点本身，则求出一个极大强连通分量。时间复杂度\(O(n+m)\).
注意事项:
注意要从每个未遍历的点出发进行遍历。
代码:
void tarjan(int u)
{
cnt++; dfn[u] = cnt; low[u] = cnt; ins[u] = true;
tp++; sta[tp] = u;
for(int i=fe0[u]; i; i=e0[i].nxt)
{
if(!dfn[e0[i].v]) {tarjan(e0[i].v); low[u] = min(low[u],low[e0[i].v]);}
else if(ins[e0[i].v]) low[u] = min(low[u],dfn[e0[i].v]);
}
if(low[u]==dfn[u])
{
num++; ca[num] = a[u];
while(sta[tp]!=u)
{
ins[sta[tp]] = false;
clr[sta[tp]] = num;
ca[num] += a[sta[tp]];
tp--;
}
ins[u] = false; clr[u] = num; tp--;
}
}
3.1.2. 点双连通分量
描述:
Tarjan算法求无向图的点双连通分量。构建DFS树，low[u]的定义改为该点经过至多一条非树边到达的最小的时间戳，若某个点\(u\)存在至少一个儿子\(v\)满足\(low[v]\ge dfn[u]\), 则\(u\)是割点。时间复杂度\(O(n+m)\).
点双连通分量一定是边双连通分量，割边的两个端点一定是割点。
代码:
(建圆方树)
namespace Graph
{
const int N = 1e6;
const int M = 4e6;
struct Edge
{
int v,nxt;
} e[(M<<1)+3];
int fe[N+3];
int fa[N+3];
int dfn[N+3],low[N+3],stk[N+3];
int n,en,cnt,tp;
void addedge(int u,int v)
{
en++; e[en].v = v;
e[en].nxt = fe[u]; fe[u] = en;
}
void Tarjan(int u)
{
cnt++; dfn[u] = low[u] = cnt;
tp++; stk[tp] = u;
for(int i=fe[u]; i; i=e[i].nxt)
{
int v = e[i].v;
if(v==fa[u]) continue;
if(!dfn[v])
{
fa[v] = u; Tarjan(v);
low[u] = min(low[u],low[v]);
if(low[v]>=dfn[u])
{
Tree::n++; Tree::addedge(u,Tree::n); Tree::addedge(Tree::n,u);
while(tp>0)
{
Tree::addedge(Tree::n,stk[tp]);
Tree::addedge(stk[tp],Tree::n);
tp--;
if(stk[tp+1]==v) {break;}
}
}
}
else {low[u] = min(low[u],dfn[v]);}
}
}
void buildTree()
{
Tree::n = n;
for(int i=1; i<=n; i++) Tree::a[i] = 1;
for(int i=1; i<=n; i++)
{
if(!dfn[i]) {Tarjan(i);}
}
}
}
3.1.3. 边双连通分量
描述:
把点双连通分量的\(\ge\)改成\(\gt\)即可。
代码:
略。
3.2. 欧拉回路
3.2.1. 无向图欧拉回路
描述:
若无向图连通且度数为奇数的点个数不超过\(2\), 则存在欧拉回路。DFS时记录回溯的路径，即可构造出一条欧拉回路，时间复杂度\(O(m)\).
注意事项:
(1) 一定要使用自杀式遍历，否则时间复杂度退化为平方级。
(2) 注意按此写法自杀式遍历不可以用取地址实现。
(3) 注意遍历的时候一定是for(int i=fe[u]; i; i=fe[u])而不是for(int i=fe[u]; i; i=e[i].nxt).
代码:
namespace Undirected
{
struct Edge{int v,nxt;} e[M+2];
int fe[N+2];
int dgr[N+2];
int ans[(N<<1)+2];
bool vis[M+2];
int n,m,en,tp;
void addedge(int u,int v)
{
en++; e[en].v = v;
e[en].nxt = fe[u]; fe[u] = en;
dgr[u]++;
}
void dfs(int u)
{
for(int i=fe[u]; i; i=fe[u])
{
fe[u] = e[i].nxt;
if(vis[(i+1)>>1]) continue;
vis[(i+1)>>1] = true;
dfs(e[i].v);
tp++; ans[tp] = (i&1) ? ((i+1)>>1) : -((i+1)>>1);
}
}
void solve()
{
scanf("%d%d",&n,&m);
for(int i=1; i<=m; i++){int x,y;scanf("%d%d",&x,&y);addedge(x,y);addedge(y,x);}
for(int i=1; i<=n; i++)
{
if(dgr[i]&1){printf("NO");return;}
}
tp = 0; dfs(e[1].v);
if(tp<m) {printf("NO"); return;}
printf("YES\n");
for(int i=tp; i>=1; i--)
{
printf("%d ",ans[i]);
}
}
}
3.2.2. 有向图欧拉回路
描述:
有向图中每个点入度等于出度是存在欧拉回路的必要条件。依然可以通过记录回溯路径的方法构造欧拉回路。时间复杂度\(O(m)\).
代码:
namespace Directed
{
struct Edge{int v,nxt;} e[M+2];
int fe[N+2];
int ind[N+2];
int oud[N+2];
int ans[(N<<1)+2];
bool vis[M+2];
int n,m,en,tp;
void addedge(int u,int v)
{
en++; e[en].v = v;
e[en].nxt = fe[u]; fe[u] = en;
oud[u]++; ind[v]++;
}
void dfs(int u)
{
for(int i=fe[u]; i; i=fe[u])
{
fe[u] = e[i].nxt;
if(vis[i]) continue;
vis[i] = true;
dfs(e[i].v);
tp++; ans[tp] = i;
}
}
void solve()
{
scanf("%d%d",&n,&m); tp = 0;
for(int i=1; i<=m; i++)
{
int x,y; scanf("%d%d",&x,&y);
addedge(x,y);
}
for(int i=1; i<=n; i++)
{
if(ind[i]!=oud[i]) {printf("NO"); return;}
}
dfs(e[1].v);
if(tp<m) {printf("NO"); return;}
printf("YES\n");
for(int i=tp; i>=1; i--) printf("%d ",ans[i]);
}
}
3.3. 斯坦纳树
描述: 给定带边权无向图中\(k\)个关键点，要求将它们联通，求最小代价。设\(dp[i][S]\)表示以\(i\)点为根，关键点的\(S\)集合已经联通的最小花费。则有两种转移: (1) 通过两个子集合并，不改变根进行转移; (2) 改变根，\(dp[u][S]=dp[v][S]+w(u,v)\), 使用SPFA进行转移。时间复杂度\(O(n3^k+SPFA(n,m)2^k)\).
代码:
struct Edge
{
int v,w,nxt;
} e[(M<<1)+3];
int fe[N+3];
int ky[NN+3];
int dp[N+3][(1<<NN)+3];
int ans[(1<<NN)+3];
bool inq[M+3];
int que[M+3];
void SPFA(int sta)
{
int head = 1,tail = 1;
for(int i=1; i<=n; i++)
{
if(dp[i][sta]<INF)
{
que[tail] = i; tail++; if(tail>n+1) tail = 1;
inq[i] = true;
}
}
while(head!=tail)
{
int u = que[head]; head++; if(head>n+1) head = 1;
for(int i=fe[u]; i; i=e[i].nxt)
{
int v = e[i].v;
if(dp[u][sta]+e[i].w<dp[v][sta])
{
dp[v][sta] = dp[u][sta]+e[i].w;
if(!inq[v])
{
que[tail] = v; tail++; if(tail>n+1) tail = 1;
inq[v] = true;
}
}
}
inq[u] = false;
}
}
void StainerTree()
{
memset(dp,42,sizeof(dp));
for(int i=0; i<nn; i++) dp[ky[i]][(1<<i)] = 0;
for(int i=1; i<(1<<nn); i++)
{
for(int j=(i-1)&i; j; j=(j-1)&i)
{
for(int k=1; k<=n; k++)
{
dp[k][i] = min(dp[k][i],dp[k][i^j]+dp[k][j]);
}
}
SPFA(i);
}
}
4. 数论及其他数学知识
4.1. 快速幂、预处理阶乘和逆元
代码:
const int P = 1e9+7;
llong fact[10000003],facti[10000003],inv[10000003];
llong quickpow(llong x,llong y)
{
llong cur = x,ret = 1ll;
for(int i=0; y; i++)
{
if(y&(1ll<<i)) {y-=(1ll<<i); ret = ret*cur%P;}
cur = cur*cur%P;
}
return ret;
}
void initfact(int n)
{
fact[0] = 1ll; for(int i=1; i<=n; i++) fact[i] = fact[i-1]*i%P;
facti[n] = quickpow(fact[n],P-2); for(int i=n-1; i>=0; i--) facti[i] = facti[i+1]*(i+1ll)%P;
for(int i=1; i<=n; i++) inv[i] = facti[i]*fact[i-1]%P;
}
4.2. 快速乘
描述:
骆可强论文中提到的快速乘方法，利用整数溢出实现，时间复杂度\(O(1)\).
注意事项:
\(x,y\) 必须小于 \(mod\).
代码:
llong quickmul(llong x,llong y,llong mod)
{
llong tmp=(x*y-(llong)((ldouble)x/mod*y+0.5)*mod);
return tmp<0?tmp+mod:tmp;
}
4.3. 扩展欧几里得算法
描述:
求解不定方程\(ax+by=\gcd(a,b)\). 时间复杂度\(O(\log(a+b))\).
代码:
llong exgcd(llong a,llong b,llong &x,llong &y)
{
if(b==0) {x = 1,y = 0; return a;}
llong nx,ny; llong ret = exgcd(b,a%b,nx,ny);
x = ny; y = nx-a/b*ny;
return ret;
}
4.4. 扩展中国剩余定理
描述:
求解同余方程组，模数可以不互质。具体做法是将模数不互质的方程进行合并。
注意事项:
(1) 注意取模溢出问题。
代码:
llong quickmul(llong x,llong y,llong mod)
{
llong tmp = x*y-(llong)((ldouble)x/mod*y)*mod;
return tmp<0?tmp+mod:tmp;
}
llong exgcd(llong a,llong b,llong &x,llong &y)
{
if(b==0) {x=1,y=0; return a;}
llong nx,ny,ret = exgcd(b,a%b,nx,ny);
x = ny; y = nx-a/b*ny; return ret;
}
llong excrt(int n,llong a[],llong p[])
{
llong mod = p[1],ans = a[1];
for(int i=2; i<=n; i++)
{
llong c = (a[i]-ans%p[i]+p[i])%p[i];
llong x,y; llong g = exgcd(mod,p[i],x,y);
if(c%g) return -1;
x = quickmul(x,c/g,p[i]/g);
ans += mod*x;
mod *= p[i]/g;
ans = (ans%mod+mod)%mod;
}
return ans;
}
4.5. 扩展BSGS
暂无
4.6. 二次剩余之Cipolla算法
暂无
4.7. 类欧几里得算法
描述: 求\(\sum^n_{x=0}x^{k_1}\lfloor\frac{ax+b}{c}\rfloor^{k_2}\), 基本思路是若\(a\ge c,b\ge c\)则提出其中的常数项，否则考虑“旋转坐标系”，从枚举横坐标转为枚举纵坐标，求出对每个\(y\)其权值乘上覆盖其的个数之和。时间复杂度\(O(\text{poly}(k_1k_2)\log (a+b+c))\).
代码:
struct Element
{
llong a[N+3][N+3];
Element() {for(int i=0; i<=N; i++) for(int j=0; j<=N; j++) a[i][j] = 0ll;}
};
Element f(llong n,llong a,llong b,llong c)
{
Element ret;
if(a==0)
{
for(int k1=0; k1<=N; k1++)
{
for(int k2=0; k2+k1<=N; k2++)
{
ret.a[k1][k2] = quickpow(b/c,k2)*PowSum::calc(n,k1)%P;
}
}
return ret;
}
if(a>=c||b>=c)
{
llong a1 = a/c,a2 = a%c,b1 = b/c,b2 = b%c;
Element tmp = f(n,a2,b2,c);
for(int k1=0; k1<=N; k1++)
{
for(int k2=0; k2+k1<=N; k2++)
{
if(k2==0)
{
ret.a[k1][k2] = PowSum::calc(n,k1);
}
else
{
for(int i=0; i<=k2; i++)
{
for(int j=0; j+i<=k2; j++)
{
ret.a[k1][k2] = (ret.a[k1][k2]+fact[k2]*finv[j]%P*finv[i]%P*finv[k2-j-i]%P*quickpow(a1,i)%P*quickpow(b1,k2-j-i)%P*tmp.a[k1+i][j])%P;
}
}
}
}
}
return ret;
}
llong m = (a*n+b)/c;
Element tmp = f(m-1,c,c-b-1,a);
for(int k1=0; k1<=N; k1++)
{
for(int k2=0; k2+k1<=N; k2++)
{
if(k2==0)
{
ret.a[k1][k2] = PowSum::calc(n,k1);
}
else
{
ret.a[k1][k2] = PowSum::calc(n,k1)*quickpow(m,k2)%P;
for(int i=0; i<=k2-1; i++)
{
for(int j=0; j<=k1+1; j++)
{
ret.a[k1][k2] = (ret.a[k1][k2]-comb(k2,i)*PowSum::coe[k1][j]%P*tmp.a[i][j]%P+P)%P;
}
}
}
}
}
return ret;
}
4.8. Nim 积
描述: 时间复杂度 \(O(\log^2n)\)，记忆化 \(x,y\le 255\). 注意初始时要把记忆化数组重置为 \(-1\).
代码:
ullong nimmul(ullong x,ullong y,int len=32)
{
if(x==0||y==0) return 0ll; if((x|y)==1ll) {return 1ll;}
if(len<=4 && nim_prod[x][y]!=-1) {return nim_prod[x][y];}
ullong xa = x>>len,xb = x^(xa<<len),ya = y>>len,yb = y^(ya<<len);
ullong z1 = nimmul(xb,yb,len>>1),z2 = nimmul(xa^xb,ya^yb,len>>1),z3 = nimmul(xa,ya,len>>1),z4 = nimmul(z3,(1llu<<len-1),len>>1);
ullong ret = z1^((z2^z1)<<len)^z4;
if(len<=4) {nim_prod[x][y] = ret;}
return ret;
}
5. 多项式
5.1. FFT及其应用
5.1.1. FFT多项式乘法
描述: 复数意义下多项式乘法，时间复杂度\(O(n\log n)\).
代码:
struct Complex
{
double x,y;
Complex() {}
Complex(double _x,double _y) {x = _x,y = _y;}
};
Complex operator +(Complex x,Complex y) {return Complex(x.x+y.x,x.y+y.y);}
Complex operator -(Complex x,Complex y) {return Complex(x.x-y.x,x.y-y.y);}
Complex operator *(Complex x,Complex y) {return Complex(x.x*y.x-x.y*y.y,x.y*y.x+x.x*y.y);}
namespace FFT
{
const int N = 1<<18;
const int lgN = 18;
const double PI = acos(-1);
int fftid[N+3];
Complex sexp[N+3];
Complex tmp1[N+3],tmp2[N+3];
int getdgr(int n) {int ret = 1; while(ret<=n) ret<<=1; return ret;}
void init_fftid(int dgr)
{
int len = 0; for(int i=1; i<=lgN; i++) {if(dgr==(1<<i)) {len = i; break;}}
fftid[0] = 0; for(int i=1; i<dgr; i++) fftid[i] = (fftid[i>>1]>>1)|((i&1)<<(len-1));
}
void fft(int dgr,int coe,Complex poly[],Complex ret[])
{
init_fftid(dgr);
if(poly!=ret) {for(int i=0; i<dgr; i++) ret[fftid[i]] = poly[i];}
else {for(int i=0; i<dgr; i++) {if(i<fftid[i]) swap(ret[i],ret[fftid[i]]);}}
for(int i=1; i<=(dgr>>1); i<<=1)
{
Complex tmp(cos(PI/i),coe*sin(PI/i));
sexp[0] = Complex(1,0); for(int j=1; j<i; j++) {sexp[j] = sexp[j-1]*tmp;}
for(int j=0; j<dgr; j+=(i<<1))
{
Complex *expn=sexp;
for(Complex *k=ret+j; k<ret+i+j; k++,expn++)
{
Complex x = (*k),y = (k[i])*(*expn);
(*k) = x+y; k[i] = x-y;
}
}
}
if(coe==-1) {for(int i=0; i<dgr; i++) ret[i].x/=dgr,ret[i].y/=dgr;}
}
}
5.1.2. NTT多项式乘法、求逆、对数函数、指数函数
描述: 模意义下多项式乘法。
代码:
namespace FFT
{
llong tmp1[N+3],tmp2[N+3],tmp3[N+3],tmp4[N+3],tmp5[N+3],tmp6[N+3],tmp7[N+3],tmp8[N+3],tmp9[N+3],tmp10[N+3];
llong tst1[N+3],tst2[N+3],tst3[N+3];
llong sexp[N+3];
int fftid[N+3];
int getdgr(int n) {int ret = 1; while(ret<=n) ret<<=1; return ret;}
void init_fftid(int dgr)
{
int len = 0; for(int i=1; i<=LGN; i++) {if((1<<i)==dgr) {len = i; break;}}
fftid[0] = 0; for(int i=1; i<dgr; i++) fftid[i] = (fftid[i>>1]>>1)|((i&1)<<(len-1));
}
void ntt(int dgr,int coe,llong poly[],llong ret[])
{
init_fftid(dgr);
if(poly==ret) {for(int i=0; i<dgr; i++) {if(i<fftid[i]) swap(ret[i],ret[fftid[i]]);}}
else {for(int i=0; i<dgr; i++) ret[i] = poly[fftid[i]];}
for(int i=1; i<=(dgr>>1); i<<=1)
{
llong tmp = quickpow(G,(P-1)/(i<<1));
if(coe==-1) {tmp = mulinv(tmp);}
sexp[0] = 1ll; for(int j=1; j<i; j++) sexp[j] = sexp[j-1]*tmp%P;
for(int j=0; j<dgr; j+=(i<<1))
{
for(llong *k=ret+j,*kk=sexp; k<ret+i+j; k++,kk++)
{
llong y = k[i]*(*kk)%P;
k[i] = (*k)-y<0 ? (*k)-y+P : (*k)-y;
(*k) = (*k)+y>=P ? (*k)+y-P : (*k)+y;
}
}
}
if(coe==-1)
{
llong tmp = mulinv(dgr);
for(int i=0; i<dgr; i++) ret[i] = ret[i]*tmp%P;
}
}
void polymul(int dgr,llong poly1[],llong poly2[],llong ret[])
{
ntt((dgr<<1),1,poly1,tmp1); ntt((dgr<<1),1,poly2,tmp2);
for(int i=0; i<(dgr<<1); i++) ret[i] = tmp1[i]*tmp2[i]%P;
ntt((dgr<<1),-1,ret,ret);
}
void polyinv(int dgr,llong poly[],llong ret[])
{
for(int i=0; i<(dgr<<1); i++) ret[i] = tmp3[i] = tmp4[i] = tmp5[i] = 0ll;
ret[0] = mulinv(poly[0]);
for(int i=1; i<=(dgr>>1); i<<=1)
{
for(int j=0; j<(i<<1); j++) tmp3[j] = poly[j];
ntt((i<<2),1,tmp3,tmp4); ntt((i<<2),1,ret,tmp5);
for(int j=0; j<(i<<2); j++) tmp3[j] = tmp4[j]*tmp5[j]%P*tmp5[j]%P;
ntt((i<<2),-1,tmp3,tmp4);
for(int j=0; j<(i<<1); j++) ret[j] = (ret[j]+ret[j]-tmp4[j]+P)%P;
}
}
void polyder(int dgr,llong poly[],llong ret[])
{
for(int i=0; i<dgr-1; i++) ret[i] = poly[i+1]*(i+1)%P;
ret[dgr-1] = 0ll;
}
void polyint(int dgr,llong poly[],llong ret[])
{
for(int i=1; i<dgr; i++) ret[i] = poly[i-1]*mulinv(i)%P;
ret[0] = 0ll;
}
void polyln(int dgr,llong poly[],llong ret[])
{
for(int i=0; i<(dgr<<1); i++) ret[i] = tmp6[i] = tmp7[i] = tmp8[i] = 0ll;
polyder(dgr,poly,tmp6);
polyinv(dgr,poly,tmp7);
polymul(dgr,tmp6,tmp7,tmp8);
polyint(dgr,tmp8,ret);
}
void polyexp(int dgr,llong poly[],llong ret[])
{
for(int i=0; i<(dgr<<1); i++) ret[i] = tmp9[i] = tmp10[i] = 0ll;
ret[0] = 1ll;
for(int i=1; i<=(dgr>>1); i<<=1)
{
polyln((i<<1),ret,tmp9);
for(int j=0; j<(i<<1); j++) tmp9[j] = (-tmp9[j]+poly[j]+P)%P; tmp9[0]++;
polymul((i<<2),ret,tmp9,tmp10);
for(int j=0; j<(i<<1); j++) ret[j] = tmp10[j];
}
}
}
5.2.3. 多项式带余除法
描述: 时间复杂度\(O(n\log n)\).
代码:
void polydiv(int dgr1,int dgr2,llong poly1[],llong poly2[],llong ret1[],llong ret2[])
{
int _dgr1 = getdgr(dgr1),_dgr2 = getdgr(dgr2);
polyrev(dgr2,poly2,tmp5); polyrev(dgr1,poly1,tmp9);
polyinv(_dgr1,tmp5,tmp6);
for(int i=dgr1-dgr2+1; i<(_dgr1<<1); i++) tmp6[i] = 0ll;
ntt(_dgr1<<1,1,tmp9,tmp7); ntt(_dgr1<<1,1,tmp6,tmp8);
for(int i=0; i<(_dgr1<<1); i++) tmp7[i] = tmp7[i]*tmp8[i]%P;
ntt(_dgr1<<1,-1,tmp7,tmp8);
for(int i=dgr1-dgr2+1; i<(_dgr1<<1); i++) tmp8[i] = 0ll;
polyrev(dgr1-dgr2+1,tmp8,ret1);
ntt(_dgr1<<1,1,poly2,tmp7); ntt(_dgr1<<1,1,ret1,tmp8);
for(int i=0; i<(_dgr1<<1); i++) tmp7[i] = tmp7[i]*tmp8[i]%P;
ntt(_dgr1<<1,-1,tmp7,ret2);
for(int i=dgr2; i<(_dgr1<<1); i++) ret2[i] = 0ll;
for(int i=0; i<dgr2-1; i++) ret2[i] = (poly1[i]-ret2[i]+P)%P;
}
5.2.4. 常系数线性递推
描述: 时间复杂度\(O(k\log k\log n)\).
代码:
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#define llong long long
#define modinc(x) {if(x>=P) x-=P;}
using namespace std;
const int N = 1<<21;
const int LGN = 21;
const int G = 3;
const int P = 998244353;
llong a[N+3];
llong b[N+3];
llong c[N+3];
llong f[N+3],g[N+3],f_t[N+3],invg[N+3],invg_t[N+3];
llong tmp1[N+3],tmp2[N+3],tmp3[N+3],tmp4[N+3];
llong tmp5[N+3],tmp6[N+3],tmp7[N+3],tmp8[N+3],tmp9[N+3];
llong tmp10[N+3],tmp11[N+3],tmp12[N+3],tmp13[N+3];
llong sexp[N+3];
int fftid[N+3];
int n,dgr; llong m;
llong quickpow(llong x,llong y)
{
llong cur = x,ret = 1ll;
for(int i=0; y; i++)
{
if(y&(1ll<<i))
{
ret = ret*cur%P;
y-=(1ll<<i);
}
cur = cur*cur%P;
}
return ret;
}
llong mulinv(llong x) {return quickpow(x,P-2);}
void init_fftid(int dgr)
{
int len = 0; for(int i=0; i<=LGN; i++) if(dgr==(1<<i)) {len = i; break;}
fftid[0] = 0ll;
for(int i=1; i<dgr; i++) fftid[i] = ((fftid[i>>1])>>1)|((i&1)<<(len-1));
}
int getdgr(int x)
{
int ret = 1; while(ret<x) ret<<=1;
return ret;
}
void ntt(int dgr,int coe,llong poly[],llong ret[])
{
init_fftid(dgr);
if(poly!=ret) {for(int i=0; i<dgr; i++) ret[i] = poly[fftid[i]];}
else {for(int i=0; i<dgr; i++) if(i<fftid[i]) swap(ret[i],ret[fftid[i]]);}
for(int i=1; i<=(dgr>>1); i<<=1)
{
llong tmp = quickpow(G,(P-1)/(i<<1));
if(coe==-1) tmp = mulinv(tmp);
sexp[0] = 1ll; for(int j=1; j<i; j++) sexp[j] = sexp[j-1]*tmp%P;
for(int j=0; j<dgr; j+=(i<<1))
{
int kk=0;
for(llong *k=ret+j; k<ret+j+i; k++)
{
llong y = k[i]*sexp[kk]%P;
k[i] = (*k)-y<0 ? (*k)-y+P : (*k)-y;
*k = (*k)+y>=P ? (*k)+y-P : (*k)+y; kk++;
}
}
}
if(coe==-1)
{
llong tmp = mulinv(dgr);
for(int j=0; j<dgr; j++) ret[j] = ret[j]*tmp%P;
}
}
void polyinv(int dgr,llong poly[],llong ret[])
{
for(int i=0; i<dgr; i++) ret[i] = 0ll;
ret[0] = mulinv(poly[0]);
for(int i=1; i<=(dgr>>1); i<<=1)
{
for(int j=0; j<(i<<2); j++) tmp1[j] = j<i ? ret[j] : 0ll;
for(int j=0; j<(i<<2); j++) tmp2[j] = j<(i<<1) ? poly[j] : 0ll;
ntt((i<<2),1,tmp1,tmp3); ntt((i<<2),1,tmp2,tmp4);
for(int j=0; j<(i<<2); j++) tmp3[j] = (tmp3[j]*tmp3[j]%P)*tmp4[j]%P;
ntt((i<<2),-1,tmp3,tmp4);
for(int j=0; j<(i<<1); j++) ret[j] = (tmp1[j]+tmp1[j]-tmp4[j]+P)%P;
}
for(int i=dgr; i<(dgr<<1); i++) ret[i] = 0ll;
}
void polyrev(int _dgr,llong poly[],llong ret[])
{
for(int i=0; i<_dgr; i++) ret[i] = poly[_dgr-1-i];
}
void polydiv(int _dgr1,llong poly1[],llong ret1[],llong ret2[]) //_dgr1-1 divides n
{
polyrev(_dgr1,poly1,tmp5);
for(int i=_dgr1-n; i<_dgr1; i++) tmp5[i] = 0ll; //Note here: tmp5 is modulo x^{_dgr1-n}. But why it's still correct without this line when the length is (dgr<<2) instead of (dgr<<1)?
ntt((dgr<<1),1,tmp5,tmp8); ntt((dgr<<1),1,invg,invg_t);
for(int i=0; i<(dgr<<1); i++) tmp8[i] = tmp8[i]*invg_t[i]%P;
ntt((dgr<<1),-1,tmp8,tmp9);
for(int i=_dgr1-n; i<_dgr1; i++) tmp9[i] = ret1[i] = 0ll;
polyrev(_dgr1-n,tmp9,ret1);
for(int i=0; i<(dgr<<1); i++) tmp8[i] = f_t[i];
ntt((dgr<<1),1,ret1,tmp9);
for(int i=0; i<(dgr<<1); i++) tmp8[i] = tmp8[i]*tmp9[i]%P;
ntt((dgr<<1),-1,tmp8,ret2);
for(int i=0; i<(dgr<<1); i++) ret2[i] = (i<n) ? (poly1[i]-ret2[i]+P)%P : 0ll;
}
void polyquickpow(llong expn,llong modp[],llong ret[])
{
if(expn<n) {ret[expn] = 1ll; return;}
polyquickpow(expn>>1,modp,ret);
ntt((dgr<<1),1,ret,tmp10);
for(int i=0; i<(dgr<<1); i++) tmp10[i] = tmp10[i]*tmp10[i]%P;
ntt((dgr<<1),-1,tmp10,tmp11);
if(expn&1) {for(int i=(dgr<<1)-1; i>=1; i--) tmp11[i] = tmp11[i-1]; tmp11[0] = 0ll;}
polydiv((n<<1)-1+(expn&1),tmp11,tmp12,tmp13);
for(int i=0; i<(dgr<<1); i++) ret[i] = i<dgr-1 ? tmp13[i] : 0ll;
}
int main()
{
scanf("%lld%d",&m,&n);
for(int i=1; i<=n; i++) {scanf("%lld",&b[i]); b[i] = (b[i]%P+P)%P;}
for(int i=0; i<n; i++) {scanf("%lld",&a[i]); a[i] = (a[i]%P+P)%P;}
dgr = getdgr(n+1);
for(int i=0; i<n; i++) f[i] = (P-b[n-i])%P; f[n] = 1ll;
ntt((dgr<<1),1,f,f_t);
polyrev(n+1,f,g);
polyinv(dgr,g,invg); //Note here: This should be dgr instead of (dgr<<1)!!!!!!!!!!
ntt((dgr<<1),1,invg,invg_t);
polyquickpow(m,f,c);
llong ans = 0ll;
for(int i=0; i<n; i++)
{
ans += c[i]*a[i]%P;
modinc(ans);
}
printf("%lld\n",ans);
return 0;
}
5.2. 多项式乘法的拓展
5.2.1. 三模数NTT
描述: 假设模数为\(P\)多项式次数为\(n\), 则乘积多项式每一项的系数不超过\(nP^2\le 10^{23}\), 使用三种不同的模数CRT合并即可。
代码:
const llong P[3] = {998244353ll,1004535809ll,167772161ll},G[3] = {3,3,3};
llong a[N+3],b[N+3],c[3][N+3],ans[N+3];
int n,m; llong P0;
int main()
{
scanf("%d%d%lld",&n,&m,&P0);
for(int i=0; i<=n; i++) scanf("%lld",&a[i]);
for(int i=0; i<=m; i++) scanf("%lld",&b[i]);
int dgr = FFT::getdgr(max(n,m));
for(int i=0; i<3; i++)
{
FFT::P = P[i]; FFT::G = G[i];
FFT::polymul(dgr,a,b,c[i]);
}
dgr<<=1;
for(int i=0; i<dgr; i++)
{
llong k = (c[1][i]-c[0][i]+P[1])*669690699ll%P[1];
llong x = c[0][i]+k*P[0];
k = (c[2][i]-x%P[2]+P[2])*107800441ll%P[2];
ans[i] = (x+k*P[0]%P0*P[1]%P0)%P0;
}
for(int i=0; i<=n+m; i++) printf("%lld ",ans[i]);
return 0;
}
5.2.2. 一种将两次DFT合并的优化技巧
暂无
5.2.3. 拆系数FFT
描述: 将系数分解为\(a\sqrt P+b\)的形式，然后分别卷积。
代码: 暂无
5.3. 拉格朗日插值法
描述: 给定\(n\)次多项式的\((n+1)\)个点值，唯一确定出这个多项式，时间复杂度\(O(n^2)\).
代码:
namespace Lagrange
{
llong ax[N+3],ay[N+3],poly[N+3];
llong aux[N+3],aux2[N+3];
void lagrange(int n)
{
aux[0] = 1ll;
for(int i=0; i<=n; i++)
{
for(int j=i+1; j>0; j--)
{
aux[j] = (aux[j-1]-aux[j]*ax[i]%P+P)%P;
}
aux[0] = P-aux[0]*ax[i]%P;
}
for(int i=0; i<=n; i++)
{
llong coe = 1ll;
for(int j=0; j<=n; j++)
{
if(i==j) continue;
coe = coe*(ax[i]-ax[j]+P)%P;
}
coe = mulinv(coe);
for(int j=0; j<=n+1; j++) aux2[j] = aux[j];
for(int j=n; j>=0; j--)
{
poly[j] = (poly[j]+ay[i]*aux2[j+1]%P*coe)%P;
aux2[j] = (aux2[j]+aux2[j+1]*ax[i])%P;
}
}
}
void clear(int n)
{
for(int i=0; i<=n+1; i++) aux[i] = aux2[i] = poly[i] = 0ll;
}
}
6. 计算几何
6.1. 计算几何基本模板
代码:
double 型：
inline int dcmp(double x) {return x<-EPS?-1:(x>EPS?1:0);}
struct Point
{
double x,y;
Point() {}
Point(double _x,double _y) {x = _x,y = _y;}
inline double ang() const {return atan2(y,x);}
bool operator ==(const Point &arg) const {return dcmp(x-arg.x)==0&&dcmp(y-arg.y)==0;}
};
typedef Point Vector;
inline Point operator +(const Point &x,const Point &y) {return Point(x.x+y.x,x.y+y.y);}
inline Point operator -(const Point &x,const Point &y) {return Point(x.x-y.x,x.y-y.y);}
inline double Dot(const Point &x,const Point &y) {return x.x*y.x+x.y*y.y;}
inline double Cross(const Point &x,const Point &y) {return x.x*y.y-x.y*y.x;}
inline double EuclidDist2(const Point &x,const Point &y) {return (x.x-y.x)*(x.x-y.x)+(x.y-y.y)*(x.y-y.y);}
long long 型：
struct Point
{
llong x,y;
Point() {}
Point(llong _x,llong _y):x(_x),y(_y) {}
bool operator ==(const Point &arg) const {return x==arg.x&&y==arg.y;}
int quadrant() {return x>=0?(y>=0?0:1):(y>0?0:1);}
};
typedef Point Vector;
Point operator +(Point x,Point y) {return Point(x.x+y.x,x.y+y.y);}
Point operator -(Point x,Point y) {return Point(x.x-y.x,x.y-y.y);}
Point operator *(Point x,llong y) {return Point(x.x*y,x.y*y);}
llong Dot(Vector x,Vector y) {return x.x*y.x+x.y*y.y;}
llong Cross(Vector x,Vector y) {return x.x*y.y-x.y*y.x;}
bool cmp_ang(Point x,Point y) {return x.quadrant()!=y.quadrant()?x.quadrant()<y.quadrant():Cross(x,y)>0;}
llong EuclidDist2(Point x,Point y) {return (x.x-y.x)*(x.x-y.x)+(x.y-y.y)*(x.y-y.y);}
6.2. 求直线交点
描述: 利用叉积和面积关系实现，时间复杂度\(O(1)\).
代码:
Point LineIntersect(Line x,Line y)
{
double t = Cross(y.dir,x.x-y.x)/Cross(x.dir,y.dir);
return x.x+x.dir*t;
}
6.3. 凸包的Minkowski和
描述: 使用类似旋转卡壳的方式，时间复杂度\(O(n)\).
代码:
void MinkowskiSum()
{
int i = 2,j = 2,k = 2; c[1] = a[cha[1]]+b[chb[1]]; cn++; cm++; cha[cn] = chb[cm] = 1;
while(i<=cn && j<=cm)
{
Vector tmp1 = a[cha[i]]-a[cha[i-1]],tmp2 = b[chb[j]]-b[chb[j-1]];
if(Cross(tmp1,tmp2)>0 || (Cross(tmp1,tmp2)==0 && tmp1.length()>tmp2.length())) {c[k] = (c[k-1]+tmp1); i++;}
else {c[k] = (c[k-1]+tmp2); j++;} k++;
}
while(i<=cn) {c[k] = c[k-1]+(a[cha[i]]-a[cha[i-1]]); k++; i++;}
while(j<=cm) {c[k] = c[k-1]+(b[chb[j]]-b[chb[j-1]]); k++; j++;}
ConvexHull(k-1,c,cc,ch);
for(int i=1; i<=cc; i++) chc[i] = c[ch[i]];
}
