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

最新推荐文章于 2019-11-05 19:48:54 发布
最新推荐文章于 2019-11-05 19:48:54 发布
阅读量2.5k 收藏 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
关注
  • 2
    点赞
  • 踩
  • 1
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
专栏目录
一些常见数列的生成函数推导
世界
11-14 1万+
曾经有人问过我：“斐波那契数列的生成函数长啥样？” 。。。所以这东西我还是写一发吧 它有什么用？它没啥用。。。1.齐次线性递推数列定义：给定常数k,a1,a2,...,ak,h0,h1,...,hk−1k,a_1,a_2,...,a_k,h_0,h_1,...,h_{k-1}，构造如下数列： hn={hna1hn−1+a2hn−2+...+akhn−kn<kn≥kh_n= \begin{c
python中finfo函数_python – 使用scipy fft和ifft以数字方式求解常微分方程
weixin_29131533的博客
01-12 386
正弦曲线的导数或复指数与其频率成正比,相移π/ 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】链式反应（分治FFT，动态规划）
weixin_30827565的博客
06-19 205
【UOJ#50】【UR #3】链式反应（分治FFT，动态规划） 题面 UOJ 题解 首先把题目意思捋一捋，大概就是有\(n\)个节点的一棵树，父亲的编号大于儿子。 满足一个点的儿子有\(2+c\)个，其中\(c\in A\)，且\(c\)个儿子是叶子，另外\(2\)个存在子树，且两种点的链接的边是不同的，求方案数。 那么就考虑一个暴力\(dp\)，设\(f[i]\)表示有\(i\)个节点的树的个数...
uoj50【UR#3】链式反应
weixin_34148340的博客
03-07 153
题解： 令$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' = ...
[UOJ50]链式反应
weixin_34081595的博客
11-15 251
这个题意说人话就是：一棵带标号的有根树，编号满足堆性质，根节点有$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...
【UOJ#50】【UR #3】链式反应（常微分方程）（牛顿迭代）
zxyoi_dreamer的博客（不定期诈尸）
09-20 498
传送门 一般当牛顿迭代里面套了Exp的时候，能写分治FFT就写分治FFT，多半跑得比牛顿迭代快。 ——zxyoi 分治FFT题解：https://blog.csdn.net/zxyoi_dreamer/article/details/101037455 题解： 首先将问题转化为如下形式： 对满足如下条件的nnn个点的树计数，其父亲儿子的标号满足堆性质，非叶节点有c+2c+2c+2个儿...
uoj-data-converter.zip
09-21
UOJ Data Converter
南京大学服务器上搭建UOJ系统的保姆级步骤
最新发布
09-30
从申请南京大学服务器开始，然后搭建docker系统，之后创建UOJ系统用于发布课堂作业，用于帮助助教进行UOJ系统搭建，帮助助教更好的完成老师布置的任务
UOJ-System:通用在线法官（社区版）
03-11
环球在线法官 一款通用的在线评论系统。 特性 ... 多种部署方式，各取所需，省心省力，方便快捷。 各组成部分可单点部署，也可分离部署...将UOJ代码 向原项目或本项目贡献代码的人 给我们启发与灵感以及提供意见和建议的人
UOJ228基础数据结构练习题.cpp
01-11
UOJ228基础数据结构练习题.cpp
unh6uoj.github.io
03-07
欢迎来到GitHub Pages 您可以使用的来维护和预览Markdown文件中网站的内容。 每当您提交到该存储库时，GitHub Pages都将运行从Markdown文件中的内容重建站点中的页面。 降价促销 Markdown是一种轻巧且易于使用的...
用傅里叶变换解偏微分方程.ppt
12-02
该课件详细讲解了用傅里叶变换求解偏微分方程的过程，以及热传导等具体问题的方程求解，希望有所帮助...
【UOJ #50】【UR #3】—链式反应（生成函数+分治NTT/多项式Exp+常微分方程）
Stargazer的博客
11-05 387
传送门 首先这个题面就很胃疼 而且感觉讲的不是很清楚 实际上是要求满足如下条件的树的个数： 对于每个非叶节点，有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...
【UR #3】链式反应
weixin_34062469的博客
04-18 303
http://uoj.ac/problem/50 %炮姐 好博客 树形结构 枚举根节点的儿子是哪两个 然后列出方程： 然后有EGF的影子！ 倍增？ 泰勒展开可以把未知数从函数里拿出来！并且变成1次项， 只要计算h(F0(x))以及h'(F0(x)) 考虑把F(x)有关项移到左边 想办法把导数或者积分放到等号右边 乘上一个关键的v(x) 因为这...
【UOJ#50】【UR #3】链式反应（CDQ分治）（NTT）
zxyoi_dreamer的博客（不定期诈尸）
09-19 261
传送门 题解： 首先把这个长且扯淡的题面先转成人话。 请你计算有多少棵有根树，父亲儿子标号满足堆性质，并且非叶节点有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−...
FFT求解卷积
ftx456789的博客
01-26 1407
卷积 已知两个多项式 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
深入浅出的讲解傅里叶变换
热门推荐
LYJ个人博客，记录技术、八卦、生活，分享各种心得，期待与朋友们一同进步
10-21 1万+
文章来自：http://www.elecfans.com/engineer/blog/20140527344277.html
微积分拾遗——链式法则
weixin_30321449的博客
07-22 1104
链式法则（chain rule）微积分中求导法则，用于求复合函数的导数； 　　链式法则应用广泛，比如神经网络中的反向传播算法就是已链式法则为基础演变的；接下来先说说链式法则的概念然后通过链式法则的两种形式学习链式法则； 　　链式法则：两个函数组合起来的复合函数，导数等于里面函数代入外函数值的导乘以里面函数之导数； 　　链式法则有两种形式： 　　　　 　　应用链式法则求复合函数导数最基础的还是函数的...
BZOJ 3160 万径人踪灭 Manacher算法+快速傅里叶变换
世界
12-27 5841
题目大意：给定一个由'a'和'b'构成的字符串，求不连续回文子序列的个数 首先回文一定是将字符串倍增 由于求的是不连续回文子序列的个数 因此我们可以求出总回文子序列的个数，然后减掉连续的 连续的就是回文子串 用Manacher算法可以O(n)求解 不连续的就有些难搞了 首先我们令f[i]表示以i为中心的对称字符对个数 比如s[]=$#a#b#a 那么s[4]='b' f[4]=2 那么
UOJ#435题目“Simple Tree”
05-30
UOJ#435题目“Simple Tree”要求实现一种数据结构，支持以下两种操作： 1. 修改某个节点的权值； 2. 查询某个节点到根节点路径上权值的最大值。 因为需要查询某个节点到根节点的路径，所以可以考虑使用树链剖分...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值