多项式相关
我学科技?真的假的?
起因是在做斯特林数,发现没事就要佛佛塔或者呢塔塔,遂学之。
FFT
现在要求两个多项式的卷积。
卷积
现在给你两个序列 \(a_i,b_i\) 并定义 \(c_i\) 为
可以认为 \(c_i\) 就是 \(a_i,b_i\) 的卷积。换句话说就是给你两个多项式 \(F(x)=\sum_{i=0}^n a_ix^i,G(x)=\sum_{i=0}^nb_ix^i\),然后 \(H(x)=F(x)*G(x)=\sum_{i=0}^{2*n}c_ix^i\),手玩一下二次式乘二次式,算出来的多项式就是卷积结果,多项式的 \(i\) 次式就是 \(c_i\)。
点值表示法
现在让你算两个多项式的乘法就可以 \(O(n^2)\) 算了,但是不够快,考虑优化这个过程,但是发现难以优化。
前人的智慧指出:对任意 \(n\) 个数点(x不同)都可以唯一地插值出来一个对应的 \(n\) 次多项式,既然 \(H=(F*G)\) 则对于原先的数点必定有 \(H(x_0)=F(x_0)*G(x_0)\)。确立这些 \(H(x)\) 的数点的过程应当是 \(O(n)\) 的,如果现在可以用一个比较高速的方法插值出 \(H(x)\) 就可以实现复杂度优化了。
单位根
如果你在看这段建议先学复数及其复平面上的三角表示和加减乘除于复平面上的表示。
考虑在复平面内给单位圆劈成 \(n\) 份,\(x\) 正方向向上的第一刀对应的向量 \(\omega_n\) 肯定满足 \(\omega_n^n=1\),因为转了一圈,就把这个东西叫单位根。
这个东西肯定就有一些性质比如
砍两倍次数每次拼两个上去和砍一倍每次拼一个一样...
因为转了一半过去。
分治法法塔(FFT)
现在有这样一个多项式
我直接给他奇偶劈两半,不妨认为 \(n\) 是偶数。
再设
然后就有
那很明显这个 \(x\) 太丑了直接给他换成幂次等效于旋转的单位根。
同理
你知道这个玩意有多牛逼吗就是说 \(F(\omega_n^k)\) 和 \(F(\omega_n^{k+n/2})\) 是可以同时算出来的,加号改减号即可而且这个玩意跨度正好是 \(n/2\)。
直接考虑一个递归分治状物然后根据6定理得到时间复杂度为 \(O(nlogn)\)。
蛋是这个玩意常数大到飞起是无法通过板题的,至少复杂度对了常数优化都是后事。
先考虑这样一个问题:这样kuku算完算出来的是一堆 \(F(x)\) 啊插值的事情还没解决。
挨法法塔(IFFT)
现在要解决如何把点值表示法转换为一般的系数表示法,而且要充分利用单位根优势,要不然不如去考虑拉插求系数。
然后前人的智慧再次发力提出一个牛逼的构造,定义
\(y_i\) 就是刚才算出来的那坨 \(F()\)。
然后这个东西可以化。
然后先放在这,设 \(S(x)=\sum_{i=0}^{n-1}x^i\) 并带入 \(\omega_n^k\) 然后构造等比数列求和式
但是这是建立在 \(\omega_n^{kn}=1\) 的基础上的,\(k=0\) 时原式为 \(n\)
所以这个
当且仅当 \(j=k\) 时为 \(n\),否则为 \(0\)。
进而 \(a_k=\frac{c_k}{n}\)。
所以流程就是:FFT,点乘计算,反着FFT(带个-1),然后除以 \(n\),即次数。
放一下板子代码方便以后照搬。
#include<bits/stdc++.h>
#define db double
#define MAXN 10000005
using namespace std;
int n,m,pw,len=1,loc[MAXN];
const db pi=acos(-1.0);
struct comp{
db x,y;
comp(db a=0,db b=0){x=a,y=b;}
}a[MAXN],b[MAXN];
comp operator+(comp a,comp b){return comp(a.x+b.x,a.y+b.y);}
comp operator-(comp a,comp b){return comp(a.x-b.x,a.y-b.y);}
comp operator*(comp a,comp b){return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
inline void FFT(comp *A,int typ){
for(int i=0;i<len;i++)if(i<loc[i])swap(A[i],A[loc[i]]);
for(int mid=1;mid<len;mid<<=1){
comp wn(cos(pi/mid),typ*sin(pi/mid));
for(int r=mid<<1,j=0;j<len;j+=r){
comp w(1,0);
for(int k=0;k<mid;k++,w=w*wn){
comp x=A[j+k],y=w*A[j+k+mid];
A[j+k]=x+y;
A[j+k+mid]=x-y;
}
}
}
}
signed main(){
scanf("%d%d",&n,&m);
for(int i=0;i<=n;i++)scanf("%lf",&a[i].x);
for(int i=0;i<=m;i++)scanf("%lf",&b[i].x);
while(len<=n+m)len<<=1,++pw;
for(int i=0;i<len;i++)loc[i]=(loc[i>>1]>>1)|((i&1)<<(pw-1));
FFT(a,1);
FFT(b,1);
for(int i=0;i<=len;i++)a[i]=a[i]*b[i];
FFT(a,-1);
for(int i=0;i<=n+m;i++)printf("%d ",(int)(a[i].x/len+0.5));
return 0;
}
NTT
就是把单位根换成原根。这个东西是用来解决带取模卷积的,一般是mod998244353。
namespace POLY{
int loc[MAXN];
const int Gx=3,GI=332748118;
struct Polynomial{
int len;
vector<int>a;
int& operator[](int x){return a[x];}
inline void fix(int x){
len=x;
a.resize(x+1);
}
inline void reset(){
a.clear();len=0;
}
inline void NTT(int n,int typ){
fix(n-1);
for(int i=0;i<n;i++)if(i<loc[i])swap(a[i],a[loc[i]]);
for(int len=1;len<n;len<<=1){
int step=qpow(typ==1?Gx:GI,(mod-1)/(len<<1));
for(int l=0;l<n;l+=(len<<1)){
int w=1;
for(int k=0;k<len;k++,w=w*step%mod){
int x=a[l+k],y=w*a[l+len+k]%mod;
a[l+k]=(x+y)%mod;
a[l+len+k]=(x-y)%mod;
}
}
}
if(typ==-1){
int inv=qpow(n,mod-2);
for(int i=0;i<n;i++)a[i]=a[i]*inv%mod;
}
}
};
Polynomial operator*(Polynomial F,Polynomial G){
Polynomial res;
int n=F.len+G.len;
int len=1,pw=0;
while(len<=n)len<<=1;
res.fix(len-1);
for(int i=0;i<len;i++)loc[i]=(loc[i>>1]>>1)|((i&1)*(len>>1));
F.NTT(len,1);G.NTT(len,1);
for(int i=0;i<len;i++)res[i]=F[i]*G[i]%mod;
res.NTT(len,-1);
res.fix(n);
return res;
}
}
题
多项式相关的题感觉主角都不是多项式而是推式子...
求
这个长得就特别卷积啊想想怎么转化成卷积形式。就是说
这个加号就很恶心所以考虑设 \(a'_i=a_{n-1-i}\),这个reverse一下就好了。
进而
然后直接开卷就行了。
乍一看和多项式没有任何关系啊,然后就要开始考虑构造了。
首先原式的上下指标都能稍微改一改,改成到 \(j\) 和从 \(j\) 开始的,不过不重要。
然后考虑构造 \(f(x)=q_x,g(x)=\frac{1}{i^2}\)
则原式变为
前半部分可以直接卷,考虑后半部分的处理。
根据上一题的经验尝试替换指标和反转
然后也可以卷了。
神秘综合题。
首先进行一波推式子。考虑对于一个点 \(x\) 什么时候会贡献答案。相当于对于当前的分治中心 \(y\),当且仅当 \(x\rightarrow y\) 的路径上还没有点成为过分治中心时(即在子树中)会贡献答案,进而要求路径上不能选点,分治中心要选其他随意,这时有 \(\frac{1}{dis(x,y)+1}\) 的概率对答案做出 \(siz_y\) 上的 \(x\) 的 \(1\) 贡献。
所以答案即为
因为距离跑到分母上了所以答案没法简单统计...
考虑点分治,对于当前分治中心 \(u\) 会得到若干距离不同的路径,如果知道不同距离路径的个数就可以简单计算了。
但是点分治不同子树枝杈间可以靠 \(u\) 两两合并路径,发现这一过程恰好是多项式卷积的形式(\(a*dis_x merge b*dis_y\rightarrow a*b*dis_{x+y}\))。但是还没完,如果有多个子树直接就假掉了...前两个子树卷积完再来一个就要卷积一个跨三个子树的路径了,显然是错的。
转而考虑容斥,用在 \(u\) 中任意选两个的方案减去所有单个子树中任选两个的不合法方案就是当前分治中心的贡献,发现正好可以变为 \(u\) 子树自己卷积减去所有它子树内自己卷积的结果,用fft处理。
然后分析复杂度,因为每次卷积的多项式大小和点分树子树大小一致所以复杂度是 \(T(n)=O(nlogn)+2(T(\frac{n}{2})+O(\frac{n}{2}logn))=O(nlog^2n)\) 的。
#include<bits/stdc++.h>
#define db double
#define MAXN 65005
using namespace std;
int n,pw,len=1,loc[MAXN];
db ans;
const db pi=acos(-1.0);
struct comp{
db x,y;
comp(db a=0,db b=0){x=a,y=b;}
}a[MAXN],f[MAXN],b[MAXN];
db P1[MAXN],P2[MAXN];
comp operator+(comp a,comp b){return comp(a.x+b.x,a.y+b.y);}
comp operator-(comp a,comp b){return comp(a.x-b.x,a.y-b.y);}
comp operator*(comp a,comp b){return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
struct Polynomial{
int len;
vector<comp>a;
comp& operator[](int x){return a[x];}
inline void fix(int x){
len=x;
a.resize(x+1);
}
inline void reset(){
a.clear();len=0;
}
inline void FFT(int n,db typ){
fix(n-1);
for(int i=0;i<n;i++)if(i<loc[i])swap(a[i],a[loc[i]]);
for(int mid=1;mid<n;mid<<=1){
comp wn(cos(pi/mid),typ*sin(pi/mid));
for(int r=mid<<1,j=0;j<n;j+=r){
comp w(1,0);
for(int k=0;k<mid;k++,w=w*wn){
comp x=a[j+k],y=w*a[j+k+mid];
a[j+k]=x+y;
a[j+k+mid]=x-y;
}
}
}
if(typ==-1)for(int i=0;i<n;i++)a[i].x/=(db)n;
}
};
Polynomial operator*(Polynomial F,Polynomial G){
Polynomial res;
int n=F.len+G.len;
int len=1,pw=0;
while(len<=n)len<<=1;
res.fix(len-1);
for(int i=0;i<len;i++)loc[i]=(loc[i>>1]>>1)|((i&1)*(len>>1));
F.FFT(len,1);G.FFT(len,1);
for(int i=0;i<len;i++)res[i]=F[i]*G[i];
res.FFT(len,-1);
res.fix(n);
return res;
}
struct node{
int v,nxt;
}edge[MAXN<<1];
int h[MAXN],tmp;
inline void add(int u,int v){
edge[++tmp]=(node){v,h[u]};
h[u]=tmp;
}
int siz[MAXN],dp[MAXN],sum,root;
bool vis[MAXN];
inline void getrt(int u,int fa){
siz[u]=1;
dp[u]=0;
for(int i=h[u];i;i=edge[i].nxt){
int v=edge[i].v;
if(v==fa||vis[v])continue;
getrt(v,u);
siz[u]+=siz[v];
dp[u]=max(dp[u],siz[v]);
}
dp[u]=max(dp[u],sum-dp[u]);
if(dp[u]<dp[root])root=u;
}
inline void getdis(int u,int fa,int d,Polynomial &F){
++F[d].x;
for(int i=h[u];i;i=edge[i].nxt){
int v=edge[i].v;
if(v==fa||vis[v])continue;
getdis(v,u,d+1,F);
}
}
inline void work(int u){
for(int i=h[u];i;i=edge[i].nxt){
int v=edge[i].v;
if(vis[v])continue;
Polynomial son;
son.fix(siz[v]);
getdis(v,u,1,son);
son=son*son;
for(int i=0;i<=son.len;i++)ans-=son[i].x/(i+1);
}
Polynomial son;
son.fix(siz[u]);
getdis(u,0,0,son);
son=son*son;
for(int i=0;i<=son.len;i++)ans+=son[i].x/(i+1);
}
inline void solve(int u){
vis[u]=1;
//printf("solving nowu=%d\n",u);
work(u);
for(int i=h[u];i;i=edge[i].nxt){
int v=edge[i].v;
if(vis[v])continue;
sum=siz[v],root=0;
getrt(v,u);
solve(root);
}
}
signed main(){
scanf("%d",&n);
for(int i=1,u,v;i<n;i++){
scanf("%d%d",&u,&v);
++u,++v;
add(u,v);
add(v,u);
}
root=0,sum=dp[root]=n;
getrt(1,0);
solve(root);
printf("%.4f",ans);
return 0;
}
发现丢失斧头的价值求和可以用多项式卷积的指数表示(恰好斧头的价值都是唯一的),系数则为方案。
考虑去重,分讨丢了一把,两把和三把斧头的情况。
对于丢了一把斧头的情况则为给定的 \(a_i\)
对于两把斧头的情况:因为不能重复地丢一把斧头,设输入给定的数据表示的多项式为 \(F\),令 \(G\) 为仅重复选两次表示出的多项式(即将\(F\)的所有次数乘二)则贡献为 \(F*F-G\)。
对于丢三把斧头的情况,这时候不合法的情况可以为重复选三把斧头或者重复选两把斧头,保留上文 \(F,G\) 的基础上设 \(H\) 是仅重复选三次得到的多项式。总方案为 \(F*F*F\),减去选两个一样且第三个不一样的方案 \(F*G-H\),但是这样的情况在三个f中会重复出现三次所以乘三,最后减去选三个一样的,贡献为 \(F*F*F-3*(F*G-H)-H\)。
这个题自排列是不合法的所以丢多于一个的给方案除个自排列就行。
有点卡精,建议统计答案全用double最后输出的时候再转int。
众所周知第二类斯特林数有通项公式
考虑转化原式为
考虑拆分后面的分式为卷积形式,设
预处理一下开卷就完事了。
然后可以顺手捏掉第二类斯特林数行。
挺有意思的。
答案相当于是回文子序列个数-回文串个数,后者好说主要看回文子序列怎么求。
两个字符回文的前提是相等且关于轴线对称,所有满足对称的位置对 \((x_1,x_2)\) 都有 \(x_1+x_2=2_0\),恰好满足卷积的下标构成。
考虑分别处理出所有 a 和 b 位置的多项式并自己卷积,则卷积后的第 \(i\) 项系数 \(a_i\) 就应该是关于 \(i/2\) 位置对称的字符对个数,对答案的贡献就应是 \(2^{\lceil\frac{a_i}{2}\rceil}\)。
设现在往上又给了 \(x\)。
发现除了 \(\sum_{i=1}^na_ib_i\) 都可以 \(O(m)\) 求出。
根据以往的经验把式子变成 \(a'_i=a_{n-i}\)
这样做的原因是对手环的旋转等效于将任一元素倍长然后旋转 \(k\) 就变为
所以翻转并倍长一个序列然后卷积求 \([n,2n)\) 的系数极大值即可。
对单点的贡献:相当于枚举度数计算,剩下 \(n-1\) 个点随便连边即可,\(n\) 个点就是式子乘 \(n\)。
然后就是处理一行斯特林数
ntt卷一卷就行。

浙公网安备 33010602011771号