UOJ#50 【UR#3C】链式反应 FFT求解多项式线性常微分方程

最新推荐文章于 2021-03-17 09:13:37 发布
最新推荐文章于 2021-03-17 09:13:37 发布
阅读量2.6k 收藏 1
点赞数 2
分类专栏： 其它题库 快速傅里叶变换 文章标签： UOJ FFT 多项式 线性常微分方程
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/popoqqq/article/details/53172703
版权

题目大意：给定 n 和集合C，对于 i=1..n 求多少 i 个节点有标号的多叉树满足：
1.父亲节点的标号大于子节点
2.一个点如果有儿子，则有两个无序的α型儿子，有 c 个无序的β型儿子，其中 c∈C
3.如果一个点是根节点或 α 型儿子，那么它可以有儿子或者是一个叶节点；如果一个点是 β 型儿子，那么它只能是一个叶节点

由于有标号，所以这里显然要使用指数级生成函数
设 fi 表示 i 个节点的多叉数数量，其指数级生成函数为：
F(x)=∑fixii!
那么考虑：
根节点是固定的，不参与标号的排列，首先把根节点刨掉，即：
F′(x)=∑fixi−1(i−1)!
然后一棵树可能只有一个叶节点，刨掉之后就什么都不剩了；
也可能有儿子，这时候儿子分为三部分：
α1,α2,β ，三者的指数级生成函数分别为 F(x),F(x),C(x)
故直观上来看应该是 F2(x)C(x) ，但是这样不对，因为 α1 和 α2 没有顺序，所以应该是 12F2(x)C(x)
可以列出方程：
F′(x)=12F2(x)C(x)+1

这个方程怎么解？

牛顿迭代。

假设现在我们有待定多项式 x(t) 、常多项式 a(t) 、已知函数 f(x)=12ax2+1 ，求解方程：
ddtx=f(x)
老办法，倍增处理。
假设我们已经知道了 x(t) 的前 n 项xn，求 x(t) 的前 2n 项 x2n
泰勒展开：
ddtx2n=f(x2n)=f(xn)+f′(xn)(x2n−xn)+f′′(xn)(x2x−xn)2+...
ddtx2n=f(xn)+f′(xn)(x2n−xn) (mod t2n)
ddtx2n−x2nf′(xn)=f(xn)−xnf′(xn) (mod t2n)
看到这里学过微积分的人已经会构造了吧。
令 r=e−∫f′(xn)dt ，则 ddtr=−f′(xn)r
等式两边同乘 r 得到：
rddtx2n+x2nddtr=(f(xn)−xnf′(xn))r (mod t2n)
ddt(x2nr)=(f(xn)−xnf′(xn))r (mod t2n)
x2n=∫(f(xn)−xnf′(xn))r dtr (mod t2n)
x2n=∫(1−12ax2n)r dtr (mod t2n)
里面的一切计算都是 O(nlogn) 的，总时间复杂度 T(n)=O(nlogn)+T(n2)=O(nlogn)
注意常数……

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define M 530000
#define MOD 998244353
#define G 3
using namespace std;
int n,m,d;
long long inv[M];
int A[M],B[M];
void Linear_Shaker()
{
    int i;
    for(inv[1]=1,i=2;i<=d<<1;i++)
        inv[i]=(MOD-MOD/i)*inv[MOD%i]%MOD;
}
long long Quick_Power(long long x,int y)
{
    long long re=1;
    while(y)
    {
        if(y&1) (re*=x)%=MOD;
        (x*=x)%=MOD; y>>=1;
    }
    return re;
}
void FFT(int a[],int n,int type){
    static int rev_bit[M];
    int i,j,k,w,wn,t,bit;
    for(bit=0;1<<bit<n;++bit);
    static int b[M];
    for(i=0;i<n;i++)
        rev_bit[i]=(rev_bit[i>>1]>>1)|((i&1)<<bit-1);
    for(i=0;i<n;i++)
        if(i<rev_bit[i])
            swap(a[i],a[rev_bit[i]]);
    for(k=2;k<=n;k<<=1)
        for(wn=Quick_Power(G,(long long)(MOD-1)/k*type%(MOD-1)),i=0;i<n;i+=k)
            for(w=1,j=0;j<k>>1;++j,w=(long long)w*wn%MOD)
            {
                t=(long long)w*a[i+j+(k>>1)]%MOD;
                a[i+j+(k>>1)]=a[i+j]-t<0?a[i+j]-t+MOD:a[i+j]-t;
                a[i+j]=a[i+j]+t>=MOD?a[i+j]+t-MOD:a[i+j]+t;
            }
    if(type!=1)
    {
        for(i=0;i<n;i++)
            a[i]=(long long)a[i]*inv[n]%MOD;
    }
}
void Get_Inv(int a[],int b[],int n) 
//求a的逆，长度为n，结果储存在b中。
//要求传入时b[0...n<<1]为空。
{
    static int temp[M];
    int i;
    if(n==1)
    {
        b[0]=Quick_Power(a[0],MOD-2);
        return ;
    }
    Get_Inv(a,b,n>>1);
    memcpy(temp,a,sizeof(a[0])*n);
    memset(temp+n,0,sizeof(a[0])*n);
    FFT(temp,n<<1,1);
    FFT(b,n<<1,1);
    for(i=0;i<n<<1;i++)
        temp[i]=(long long)b[i]*(2-(long long)temp[i]*b[i]%MOD+MOD)%MOD;
    FFT(temp,n<<1,MOD-2);
    memcpy(b,temp,sizeof(a[0])*n);
    memset(b+n,0,sizeof(a[0])*n);
}
void Get_Ln(int a[],int b[],int n)
//求a的ln，长度为n，结果储存在b中。
//要求a的常数项为1。
{
    static int a_[M],a_inv[M];
    int i;
    Get_Inv(a,a_inv,n);
    for(i=0;i<n-1;i++)
        a_[i]=(long long)a[i+1]*(i+1)%MOD;
    FFT(a_,n<<1,1);
    FFT(a_inv,n<<1,1);
    for(i=0;i<n<<1;i++)
        b[i]=(long long)a_[i]*a_inv[i]%MOD;
    FFT(b,n<<1,MOD-2);
    for(i=n-1;i;i--)
        b[i]=b[i-1]*inv[i]%MOD;
    b[0]=0;
    memset(b+n,0,sizeof(b[0])*n);
    memset(a_,0,sizeof(a_[0])*n<<1);
    memset(a_inv,0,sizeof(a_inv[0])*n<<1);
}
void Get_Exp(int a[],int b[],int n)
//求a的exp，长度为n，结果储存在b中。
//要求传入时b[0...n<<1]为空。
//要求a的常数项为0。
{
    static int temp[M];
    int i;
    if(n==1)
    {
        b[0]=1;
        return ;
    }
    Get_Exp(a,b,n>>1);
    Get_Ln(b,temp,n);
    for(i=0;i<n;i++)
        temp[i]=((i==0)+MOD-temp[i]+a[i])%MOD;
    FFT(temp,n<<1,1);
    FFT(b,n<<1,1);
    for(i=0;i<n<<1;i++)
        b[i]=(long long)b[i]*temp[i]%MOD;
    FFT(b,n<<1,MOD-2);
    memset(b+n,0,sizeof(b[0])*n);
}
void Newton_Method(int a[],int b[],int n)
//求解常微分方程dx/dt=ax^2/2+1，长度为n，结果储存在b中。
//要求传入时b[0...n<<1]为空。
//x[2n]=int( (1-ax^2/2)*r )/r
//r=exp(-int(ax))
{
    static int temp[M],temp_temp[M],r[M],r_temp[M];
    int i;
    if(n==1)
    {
        b[0]=0;
        return ;
    }
    Newton_Method(a,b,n>>1);
    memcpy(temp,a,sizeof(a[0])*n);
    memset(temp+n,0,sizeof(a[0])*n);
    FFT(temp,n<<1,1);
    FFT(b,n<<1,1);
    for(i=0;i<n<<1;i++)
        temp_temp[i]=(long long)temp[i]*b[i]%MOD;
    FFT(temp_temp,n<<1,MOD-2);
    for(i=n-1;i;i--)
        temp_temp[i]=temp_temp[i-1]*inv[i]%MOD*(MOD-1)%MOD;
    temp_temp[0]=0;
    memset(r,0,sizeof(a[0])*n<<1);
    Get_Exp(temp_temp,r,n);

    for(int i=0;i<n<<1;i++)
        temp[i]=(1+inv[2]*(MOD-1)%MOD*temp[i]%MOD*b[i]%MOD*b[i]%MOD)%MOD;
    FFT(temp,n<<1,MOD-2);
    memset(temp+n,0,sizeof(a[0])*n);

    memcpy(r_temp,r,sizeof(a[0])*n<<1);
    FFT(temp,n<<1,1);
    FFT(r_temp,n<<1,1);
    for(int i=0;i<n<<1;i++)
        temp[i]=(long long)temp[i]*r_temp[i]%MOD;
    FFT(temp,n<<1,MOD-2);

    for(i=n-1;i;i--)
        temp[i]=temp[i-1]*inv[i]%MOD;
    temp[0]=0;
    memset(temp+n,0,sizeof(a[0])*n);

    memset(r_temp,0,sizeof(a[0])*n<<1);
    Get_Inv(r,r_temp,n);

    FFT(temp,n<<1,1);
    FFT(r_temp,n<<1,1);
    for(int i=0;i<n<<1;i++)
        b[i]=(long long)temp[i]*r_temp[i]%MOD;
    FFT(b,n<<1,MOD-2);

    memset(b+n,0,sizeof(a[0])*n);
    memset(temp+n,0,sizeof(a[0])*n);
}
int main()
{
    cin>>n;
    for(d=1;d<=n+1;d<<=1);
    Linear_Shaker();

    long long temp=1;
    for(int i=0;i<n;i++)
    {
        scanf("%1d",&A[i]);
        A[i]=A[i]*temp;
        (temp*=inv[i+1])%=MOD;
    }

    Newton_Method(A,B,d);

    temp=1;
    for(int i=1;i<=n;i++)
    {
        (temp*=i)%=MOD;
        printf("%d\n",int(B[i]*temp%MOD));
    }

    return 0;
}
PoPoQQQ
关注
专栏目录
UOJ#188. 【UR #13】Sanrd（min25筛）
m0_52639539的博客
05-02 363
UOJ#188. 【UR #13】Sanrd（min25筛） 题目链接：https://uoj.ac/problem/188 题目大意：定义数论函数 f(x)f(x)f(x) 的值为 xxx 的第二大质因子，其中质因子可以重复贡献，例如 f(36)=2f(36)=2f(36)=2 ，
UOJ #33. 【UR #2】树上GCD
Freopen的博客
04-28 423
题目 先求得GCD是iii的倍数的点对数再反演。 用长链剖分求出子树内各个深度的节点个数，然后在合并子树的时候统计答案。 先求出深度小的子树vvv的每个iii的倍数的点数，再O(maxdepv)O(maxdep_v)O(maxdepv​)枚举iii，在之前的子树中求出iii的倍数的点数，然后乘起来即可。 在之前的子树中可以如果i>Si>Si>S就暴力，i<Si<Si&...
[UOJ50]链式反应
weixin_34081595的博客
11-15 276
这个题意说人话就是：一棵带标号的有根树，编号满足堆性质，根节点有$x$个儿子是叶子（$x\in A$），另外的$2$个儿子也是这样的一棵树，求不同的树的个数 设$f_n$为答案，枚举那两棵子树的大小$i,j$，那么$f_i=\frac12\sum\limits_{1\leq j,k\lt i}[i-1-j-k\in A]\binom{i-1}j\binom{i-1-j}kf_jf_k$，$f_1...
uoj50【UR#3】链式反应
weixin_34148340的博客
03-07 173
题解： 令$a(x)$为破坏死光的$EFG$,$f(x)$为方案的$EGF$：$f(x) = x + \int \ \frac{1}{2} f^2(x) a(x) \ dt$; 注意到$f(0)=0$，所以考虑如何解：$f'(x) = \frac{1}{2} a(x) f(x)^2 +1$ 设$g(f) = 1 + \frac{1}{2}af^2$，即求解$f' = ...
【UR #3】链式反应
weixin_34062469的博客
04-18 322
http://uoj.ac/problem/50 %炮姐 好博客 树形结构 枚举根节点的儿子是哪两个 然后列出方程： 然后有EGF的影子！ 倍增？ 泰勒展开可以把未知数从函数里拿出来！并且变成1次项， 只要计算h(F0(x))以及h'(F0(x)) 考虑把F(x)有关项移到左边 想办法把导数或者积分放到等号右边 乘上一个关键的v(x) 因为这...
【UOJ#50】【UR #3】链式反应（常微分方程）（牛顿迭代）
zxyoi_dreamer的博客（不定期诈尸）
09-20 581
传送门 一般当牛顿迭代里面套了Exp的时候，能写分治FFT就写分治FFT，多半跑得比牛顿迭代快。 ——zxyoi 分治FFT题解：https://blog.csdn.net/zxyoi_dreamer/article/details/101037455 题解： 首先将问题转化为如下形式： 对满足如下条件的nnn个点的树计数，其父亲儿子的标号满足堆性质，非叶节点有c+2c+2c+2个儿...
uoj #181. 【UR #12】密码锁 搜索
beginend
04-11 348
题意 给出一张nnn个点的完全图，现在要给这个完全图的每一条边随机定向成一个有向图。对于一条边 (i,j)(i&lt;j)(i,j)(i&lt;j)(i,j)(iiii到jjj的概率是numi,j/10000numi,j/10000num_{i,j}/10000，numi,jnumi,jnum_{i,j} 指这条边旁边的数字，否则就是jjj到iii在随机定向后，设这张有向图的强连通分量数目为xx...
【UOJ#33】【UR #2】树上GCD（长链剖分/根号分类讨论）
xajd1906的博客
01-18 171
【UOJ#33】【UR #2】树上GCD 求解树上两个点到lca的距离的最大公约数是k的对数 首先我们很容易就想到莫比乌斯反演，那么利用倍数形式，我们只需要求解是i的倍数的对数。 考虑枚举lca，这个问题就和深度有关，那么可以长链剖分，然后每次枚举倍数是O(dlogd)的，但是这里的深度不止是轻链还有重链，那么复杂度就不对了。 但是我们发现当这个倍数很大时复杂度是正确的，每次维护轻链的倍数，然后对于每个倍数跳，总复杂度就是O(nn)O(n\sqrt{n})O(nn​)。但是很小的时候就不对了，但是我们利用根
python中finfo函数_python – 使用scipy fft和ifft以数字方式求解常微分方程
weixin_29131533的博客
01-12 483
正弦曲线的导数或复指数与其频率成正比,相移π/ 2.对于复指数,相移相当于乘以j.例如,d / dt exp(j *Ω* t)== j *Ω* exp(j *Ω* t)==Ω* exp(j *π/ 2)* exp(j *Ω* t)= =Ω* exp(j *(Ω*tπ/ 2)).因此,如果你有傅里叶变换对u(t)< => U(Ω),然后du / dt< => jΩ* U(Ω)...
【UOJ #50】【UR #3】—链式反应（生成函数+分治NTT/多项式Exp+常微分方程）
Stargazer的博客
11-05 439
传送门 首先这个题面就很胃疼 而且感觉讲的不是很清楚 实际上是要求满足如下条件的树的个数： 对于每个非叶节点，有c+2c+2c+2个儿子，其中有c∈Ac\in Ac∈A个叶子节点和2个非叶节点 且点的编号满足父亲小于儿子 按照套路设f[i]f[i]f[i]为iii个点的答案可以列出dpdpdp式 f[i]=12∑j∑k[i−j−k−1∈A](i−1j)(i−j−1k)fjfkf[i]=\frac...
【UOJ#50】【UR #3】链式反应（CDQ分治）（NTT）
zxyoi_dreamer的博客（不定期诈尸）
09-19 310
传送门 题解： 首先把这个长且扯淡的题面先转成人话。 请你计算有多少棵有根树，父亲儿子标号满足堆性质，并且非叶节点有c+2c+2c+2个儿子，其中c∈Ac\in Ac∈A，ccc个儿子全部都是叶子，剩下两个儿子不是叶子，它们的子树符合上述性质（请递归理解）。 一个很显然的DP就出来了，设fif_ifi​表示iii个点的方案数，则 fi=12∑j=1∑k=1[i−1−j−k∈A](i−1j)(i−...
【UOJ#50】【UR #3】链式反应（分治FFT，动态规划）
weixin_30827565的博客
06-19 223
【UOJ#50】【UR #3】链式反应（分治FFT，动态规划） 题面 UOJ 题解 首先把题目意思捋一捋，大概就是有\(n\)个节点的一棵树，父亲的编号大于儿子。 满足一个点的儿子有\(2+c\)个，其中\(c\in A\)，且\(c\)个儿子是叶子，另外\(2\)个存在子树，且两种点的链接的边是不同的，求方案数。 那么就考虑一个暴力\(dp\)，设\(f[i]\)表示有\(i\)个节点的树的个数...
fft的时移特性MATLAB,求解Matlab微分方程组中的时移问题！！！
weixin_39572794的博客
03-17 316
//https://www.ilovematlab.cn/thread-61233-1-2.html!using["XSLSF"]; //使用命名空间XSLSF//数组xArray存放x的值；ti为当前有效值的个数；tmax为ti对应的时间；tmin为起始时间。xt(t:k:xArray,ti,tmax,tmin)={k=(t-tmin)/(tmax-tmin)*ti...
洛谷P5205 【模板】多项式开根（FFT）
weixin_34248023的博客
03-20 103
题面 传送门 题解 考虑分治 假设我们已经求出\(A'^2\equiv B\pmod{x^n}\)，考虑如何计算出\(A^2\equiv B\pmod{x^{2n}}\) 首先肯定存在\(A^2\equiv B\pmod{x^n}\) 然后两式相减\[A'^2-A^2\equiv 0\pmod{x^n}\] \[(A'-A)(A'+A)\equiv 0\pmod{x^n}\] 我们假设\(A'-A...
FFT算法学习笔记
dehukun0483的博客
04-17 329
写在前边 　　1.辣鸡RRRR_wys之前csdn的博客，千年不更。。。还很水。。。于是开了这个Blog。。。妄图拯救一下自己 　　2.最近接触了一些多项式理论。于是翘掉了愉快的高频自控，通过《算导》稍稍学习了一下 　　3.算法竞赛中，FFT主要解决多项式的乘法等问题 FFT基本概念 　　1.FFT即快速傅里叶变换, 是离散傅里叶变换的加速算法。可以在o(nlogn)的时间内，...
FFT求解卷积
ftx456789的博客
01-26 1487
卷积 已知两个多项式 A(x)=a0+a1x1+a2x2+a3x3+⋯+anxn" role="presentation" style="position: relative;">A(x)=a0+a1x1+a2x2+a3x3+⋯+anxnA(x)=a0+a1x1+a2x2+a3x3+⋯+anxnA(x)=a_{0}+a_{1}x^1+a_{2}x^2+a_{3}x^3+\cdot
两个矩阵是否相交的算法_数学、运筹学、计算机等领域的36个重要算法
weixin_39957461的博客
11-26 341
哈哈哈教你一句话激怒各个专业！！！中国大学的现实是：大学层次越低，上课越多，学生读书越少30省份开学时间一览！仅剩1地未公布开学时间重要算法要学会1. A *搜索算法图搜索算法，用于查找从给定初始节点到给定目标节点的路径。它采用启发式估计，通过估计通过该节点的最佳路径对每个节点进行排名。它按照此启发式估计的顺序访问节点。因此，A *算法是最佳优先搜索的示例。2. 波束搜索...
FFT（快速傅里叶变换）详解
cqbzcsq的博客
12-26 1975
FFT简介 一种在O(nlogn)时间复杂度下把多项式（给定系数表达式）在系数表达式和点值表达式中互换的算法。 FFT的相关知识 一、多项式 1、定义：由若干个单项式相加组成的代数式叫做多项式（如：） 2、多项式的次数：组成多项式的单项式中，次数最高的单项式的次数，为该多项式的次数 3、多项式的次数界：多项式的次数加1 4、多项式的表示方法： （1）、系数表示法： 如：、 （2...
UOJ#435题目“Simple Tree”
最新发布
05-30
UOJ#435题目“Simple Tree”要求实现一种数据结构，支持以下两种操作： 1. 修改某个节点的权值； 2. 查询某个节点到根节点路径上权值的最大值。 因为需要查询某个节点到根节点的路径，所以可以考虑使用树链剖分算法来解决这个问题。具体而言，可以对树进行一遍深度优先遍历（DFS），在遍历到每个节点的时候，记录它所在的链的编号、它在链中的位置以及链的信息，包括链的最大值和次大值等。这样，在查询某个节点到根节点路径上权值的最大值时，就可以分解成若干个链的查询，并且对于每个链的查询，可以直接使用预处理好的信息进行处理，而不需要暴力地查询每个节点的权值。 对于修改某个节点的权值，可以使用线段树来维护每个节点的信息，包括节点的权值以及节点到根节点路径上的最大值和次大值等。由于每个节点只会被访问一次，所以时间复杂度为 O(nlogn)。 总的时间复杂度为 O(nlogn)，其中 n 表示树的节点数。
评论 2
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值