【XSY3344】连续段 DP 牛顿迭代 NTT

最新推荐文章于 2021-10-03 11:01:32 发布
最新推荐文章于 2021-10-03 11:01:32 发布
阅读量565 收藏 1
点赞数
文章标签： 数据结构与算法
原文链接：http://www.cnblogs.com/ywwyww/p/10193748.html
版权

题目大意

　　对于一个长度为 \(n\) 的排列 \(p\)，我们称一个区间 \([l,r]\) 是连续的当且仅当 \((\max_{l\leq i\leq r}a_i)-(\min_{l\leq i\leq r}a_i)=r-l\)。

　　对于两个排列 \(p_1,p_2\)，我们称这两个排列是等价的，当且仅当他们的长度相同且连续区间的集合相同。

　　给你 \(n\)，对于 \(1\leq i\leq n\)，所有求 \(i\) 阶排列形成的等价类个数对 \(p\) 取模的值。

　　\(n\leq 100000,p=k\times2^{18}+1,k\in \mathbb {N},p\) 是质数。

题解

　　对于一个区间 \([l,r]\)，如果这个区间的所有子区间都是连续的，就称这个区间为强连续区间。

　　每次我们找一个最大的强连续区间，把这个区间内所有点缩到一起。

　　如果找不到，就找一个最小的连续区间，把这个区间内所有点缩到一起。

　　这个过程就形成了一个树结构。

  • 总共有 \(n\) 个叶子。
  • 对于一个代表强连续区间的点，这个点的儿子个数 \(\geq 2\)。
  • 对于一个代表连续区间的点，这个点的儿子个数 \(\geq 4\)。（可以发现，长度为 \(\leq 3\) 的序列是一定有强连续区间的。）

　　每一个排列 \(p\) 对应了一棵树。

　　对于一棵树，可以构造出一个排列 \(p\)。

　　那么就只用对这棵树计数了。

　　显然可以 \(O(n^2)\) DP。

　　记 \(f_i\) 为长度为 \(i\) 的排列的方案数，\(F(x)=\sum_{i\geq 0}f_ix^i\)，那么就有：
\[ F(x)=\frac{(F(x)^2+F(x)^4)}{1-F(x)}+x \]
　　用牛顿迭代法解这个方程即可。
\[ (F(x)^2+F(x)^4)\frac{1}{1-F(x)}-F(x)+x=0\\ G(y)=(y^2+y^4)\frac{1}{1-y}-y+x\\ G'(y)=(2y+4y^3)\frac{1}{1-y}+(y^2+y^4)(-\frac{1}{(1-y)^2})(-1)-1\\ G(F(x))\equiv 0\pmod {x^n}\\ G(F(x))\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x)) \pmod {x^n}\\ F(x)\equiv F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}\pmod {x^n}\\ \]
　　时间复杂度：\(O(n\log n)\)

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#include<functional>
#include<cmath>
#include<vector>
#include<assert.h>
//using namespace std;
using std::min;
using std::max;
using std::swap;
using std::sort;
using std::reverse;
using std::random_shuffle;
using std::lower_bound;
using std::upper_bound;
using std::unique;
using std::vector;
typedef long long ll;
typedef unsigned long long ull;
typedef double db;
typedef std::pair<int,int> pii;
typedef std::pair<ll,ll> pll;
void open(const char *s){
#ifndef ONLINE_JUDGE
    char str[100];sprintf(str,"%s.in",s);freopen(str,"r",stdin);sprintf(str,"%s.out",s);freopen(str,"w",stdout);
#endif
}
void open2(const char *s){
#ifdef DEBUG
    char str[100];sprintf(str,"%s.in",s);freopen(str,"r",stdin);sprintf(str,"%s.out",s);freopen(str,"w",stdout);
#endif
}
int rd(){int s=0,c,b=0;while(((c=getchar())<'0'||c>'9')&&c!='-');if(c=='-'){c=getchar();b=1;}do{s=s*10+c-'0';}while((c=getchar())>='0'&&c<='9');return b?-s:s;}
void put(int x){if(!x){putchar('0');return;}static int c[20];int t=0;while(x){c[++t]=x%10;x/=10;}while(t)putchar(c[t--]+'0');}
int upmin(int &a,int b){if(b<a){a=b;return 1;}return 0;}
int upmax(int &a,int b){if(b>a){a=b;return 1;}return 0;}
const int N=270000;
ll p,g;
ll fp(ll a,ll b)
{
    ll s=1;
    for(;b;b>>=1,a=a*a%p)
        if(b&1)
            s=s*a%p;
    return s;
}
namespace ntt
{
    const int W=262144;
    int rev[N];
    int *w[20];
    void init()
    {
        ll s=fp(g,(p-1)/W);
        w[18]=new int[1<<17];
        w[18][0]=1;
        for(int i=1;i<W/2;i++)
            w[18][i]=w[18][i-1]*s%p;
        for(int i=17;i>=1;i--)
        {
            w[i]=new int[1<<(i-1)];
            for(int j=0;j<1<<(i-1);j++)
                w[i][j]=w[i+1][j<<1];
        }
    }
    void ntt(ll *a,int n,int t)
    {
        for(int i=1;i<n;i++)
        {
            rev[i]=(rev[i>>1]>>1)|(i&1?n>>1:0);
            if(rev[i]>i)
                swap(a[i],a[rev[i]]);
        }
        for(int i=2,l=1;i<=n;i<<=1,l++)
            for(int j=0;j<n;j+=i)
                for(int k=0;k<i/2;k++)
                {
                    ll u=a[j+k];
                    ll v=a[j+k+i/2]*w[l][k];
                    a[j+k]=(u+v)%p;
                    a[j+k+i/2]=(u-v)%p;
                }
        if(t==-1)
        {
            reverse(a+1,a+n);
            ll inv=fp(n,p-2);
            for(int i=0;i<n;i++)
                a[i]=a[i]*inv%p;
        }
    }
    void add(ll *a,ll *b,ll *c,int n,int m,int l)
    {
        static ll a1[N];
        int k=max(n,m);
        for(int i=0;i<=k;i++)
            a1[i]=0;
        for(int i=0;i<=n;i++)
            a1[i]=(a1[i]+a[i])%p;
        for(int i=0;i<=m;i++)
            a1[i]=(a1[i]+b[i])%p;
        for(int i=0;i<=l;i++)
            c[i]=a1[i];
    }
    void mul(ll *a,ll *b,ll *c,int n,int m,int l)
    {
        static ll a1[N],a2[N];
        int k=1;
        while(k<=n+m)
            k<<=1;
        for(int i=0;i<k;i++)
            a1[i]=a2[i]=0;
        for(int i=0;i<=n;i++)
            a1[i]=a[i];
        for(int i=0;i<=m;i++)
            a2[i]=b[i];
        ntt(a1,k,1);
        ntt(a2,k,1);
        for(int i=0;i<k;i++)
            a1[i]=a1[i]*a2[i]%p;
        ntt(a1,k,-1);
        for(int i=0;i<=l;i++)
            c[i]=(a1[i]+p)%p;
    }
    void getinv(ll *a,ll *b,int n)
    {
        if(n==1)
        {
            b[0]=fp(a[0],p-2);
            return;
        }
        getinv(a,b,n>>1);
        static ll a1[N],a2[N];
        for(int i=0;i<n<<1;i++)
            a1[i]=a2[i]=0;
        for(int i=0;i<n;i++)
            a1[i]=a[i];
        for(int i=0;i<n>>1;i++)
            a2[i]=b[i];
        ntt(a1,n<<1,1);
        ntt(a2,n<<1,1);
        for(int i=0;i<n<<1;i++)
            a1[i]=a2[i]*(2-a1[i]*a2[i]%p)%p;
        ntt(a1,n<<1,-1);
        for(int i=0;i<n;i++)
            b[i]=(a1[i]+p)%p;
    }
}
int e[50],t;
int check(int x)
{
    for(int i=1;i<=t;i++)
        if(fp(x,(p-1)/e[i])==1)
            return 0;
    return 1;
}
void getg()
{
    int _=p-1;
    for(int i=2;i*i<=_;i++)
        if(_%i==0)
        {
            e[++t]=i;
            while(_%i==0)
                _/=i;
        }
    if(_!=1)
        e[++t]=_;
    for(int i=1;;i++)
        if(check(i))
        {
            g=i;
            return;
        }
}
void solve(ll *a,int n)
{
    if(n==1)
        return;
    solve(a,n>>1);
    static ll f1[N],f2[N],f3[N],f4[N],f5[N],a1[N],a2[N],a3[N],a4[N];
    for(int i=0;i<n;i++)
        a1[i]=(p-a[i])%p;
    a1[0]++;
    for(int i=0;i<n;i++)
        f1[i]=a[i];
    ntt::ntt(f1,n<<1,1);
    for(int i=0;i<n<<1;i++)
        f2[i]=f1[i]*f1[i]%p;
    ntt::ntt(f2,n<<1,-1);
    for(int i=0;i<n;i++)
        f4[i]=f2[i];
    ntt::ntt(f4,n<<1,1);
    for(int i=0;i<n<<1;i++)
        f3[i]=f4[i]*f1[i]%p;
    ntt::ntt(f3,n<<1,-1);
    for(int i=0;i<n<<1;i++)
        f4[i]=f4[i]*f4[i]%p;
    ntt::ntt(f4,n<<1,-1);
    for(int i=0;i<n;i++)
        f5[i]=f4[i];
    ntt::ntt(f5,n<<1,1);
    for(int i=0;i<n<<1;i++)
        f5[i]=f5[i]*f1[i]%p;
    ntt::ntt(f5,n<<1,-1);
    for(int i=0;i<n;i++)
        f1[i]=a[i];
    
    for(int i=0;i<n;i++)
        a1[i]=(-f1[i]+3*f2[i]-2*f3[i]+f4[i]-f5[i])%p;
    for(int i=1;i<n;i++)
        a1[i]=(a1[i]-2*f1[i-1]+f2[i-1])%p;
    a1[1]++;
    
    for(int i=0;i<n;i++)
        a2[i]=(4*f1[i]-2*f2[i]+4*f3[i]-3*f4[i])%p;
    a2[0]--;

    ntt::getinv(a2,a3,n);
    ntt::mul(a1,a3,a4,n-1,n-1,n-1);
        
    for(int i=0;i<n;i++)
        a[i]=(a[i]-a4[i]+p)%p;
}
int n;
ll s[N];
int main()
{
    open("b");
    scanf("%d%lld",&n,&p);
    getg();
    ntt::init();
    int k=1;
    while(k<=n)
        k<<=1;
    solve(s,k);
    for(int i=1;i<=n;i++)
    {
        s[i]=(s[i]+p)%p;
        printf("%lld\n",s[i]);
    }
    return 0;
}

转载于:https://www.cnblogs.com/ywwyww/p/10193748.html

weixin_30571465
关注
  • 0
    点赞
  • 踩
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
[XSY] 计数（DP，NTT，分治）
Edward The Bunny 的博客
03-23 294
计数 考虑转化题目，变为网格上有若干个点，要从(0,0)(0,0)(0,0)走到(n,an+1)(n,a_{n+1})(n,an+1​) ，每一步只能往右走一步或往上走一步，且若当前在(i,j)(i,j)(i,j) ，必须满足0≤j≤ai+10\leq j\leq a_{i+1}0≤j≤ai+1​，其中an+1=ana_{n+1}=a_nan+1​=an​ 。要对每个kkk求出有多少条恰走过kkk次形如(i,ai+1)→(i+1,ai+1)(i,a_{i+1})\to(i+1,a_{i+1})(i,ai+
[XSY] 树与图（树形DP、生成函数、分治NTT、重链剖分）
Edward The Bunny 的博客
03-17 230
树与图 如果真的在图上跑算法，那么光建图复杂度就O(n2logn)O(n^2logn)O(n2logn)了，这显然不可行。所以一定要把 在图上的操作 转换成 在树上的操作 在图上删去点u，相当于在树上删去u到根节点的链，并把u的整棵子树删掉 如果我们枚举算法可能产生的所有过程，然后再去求每个过程对应的删点次数，那么光枚举过程就已经可以T爆了 所以我们只能枚举删点次数，然后求出该删掉次数对应多少种过程 这可以用树形DP实现，时间复杂度O(n3)O(n^3)O(n3) 可以想到用生成函数进一步优化
连续段 - dp - 多项式理论 - 牛顿迭代
Mys_C_K的博客
12-30 558
题目大意： 对每个n∈[1,N]n\in[1,N]n∈[1,N]求有多少长为nnn的本质不同的排列，其中两个排列p1,p2p_1,p_2p1​,p2​本质相同，当且仅当S(p1)=S(p2)S(p_1)=S(p_2)S(p1​)=S(p2​)，其中： S(p)={[l,r]∣1≤l≤r≤∣p∣,(min⁡i=lrpi)−(max⁡i=lrpi)=r−l}S(p)=\{[l,r]|1\le l\le...
【清华夏令营模拟2019.5.22】连续段（析合树+多项式牛顿迭代）
Cold_Chair的博客
05-22 678
Description： n<=1e5，P是NTT模数 题解： 析合树见WC2019LCA营员交流讲稿 我们考虑把一个序列划分成本源连续段。 怎么划分呢？就是极大划分，假设划分成了x段。 这x段需要满足任取一个区间的段，要么可以全部可以拼起来(1)，要么除了最大的那一个其它都不行(2)。 仔细思考这样就能表示所有的段了。 (1)就是析点，(2)就是合点。 考虑x=2、3时只能是合点 设F表...
[DP 倍增]区间的连续段
鱼竿钓鱼干
04-23 122
[DP 倍增]区间的连续段 题目 题目链接 思路 代码 在这里插入代码片
zoj 3344 //组合数学那种的DP
hqd_acm的专栏
08-26 1813
import java.io.*; import java.math.*; import java.util.*; /** * * @author acm */ public class Main { public static void main(String[]
UOJ#428. 【集训队作业2018】普通的计数题（牛顿迭代）
DZYO的博客
11-21 3072
传送门 题解： 把0操作看做是叶子，1操作看做非叶节点，一个操作在另一个操作删除，则另一个操作为这个操作的父亲，于是转化成了满足以下条件的nnn个点的树的计数： 1.父亲标号&amp;amp;amp;amp;gt;儿子。 2.若一个点为非叶节点，记其儿子中叶子节点的数量为TTT，则若其儿子中有非叶节点，T∈AT\in AT∈A，否则T∈BT \in BT∈B。 首先可以发现的是BBB集合中有没有0都无所谓（因为必须选非空序列），...
[XSY] 宝藏（LCS，DP）
Edward The Bunny 的博客
03-10 158
宝藏 首先，这个问题等价于给定两个字符串S,T ，每次询问LCS(S[1,n],T[x,y])LCS(S[1,n],T[x,y])LCS(S[1,n],T[x,y])。 对每个询问重新求一遍LCSLCSLCS显然不现实，又因为x,yx,yx,y都是连续变化的，我们考虑探讨 LCS(S[1,n],T[x,y])LCS(S[1,n],T[x,y])LCS(S[1,n],T[x,y])与LCS(S[1,n],T[x,y−1])LCS(S[1,n],T[x,y-1])LCS(S[1,n],T[x,y−1])的关系
[CSA35G][XSY3318]Counting Quests （DP）
Edward The Bunny 的博客
09-28 212
XSY3318 CSA35G 对于一个询问区间的集合SSS，求出每一个数被哪些区间覆盖了，记为SiS_iSi​。 要能保证猜出选中数，当且仅当每个数的SiS_iSi​互不相同。 考虑求出不满足要求的集合SSS的个数。 首先可以观察得到SiS_iSi​的一个性质：若a<b<c<da<b<c<da<b<c<d且Sa=Sc,Sb=SdS_a=S_c,S_b=S_dSa​=Sc​,Sb​=Sd​，则必有Sa=Sb=Sc=SdS_a=S_b=S_c=S_dSa​=
连续段问题小结
Edward The Bunny 的博客
10-03 288
oi-wiki P4747 P6795 CF526F blog
ios app的ipv4连接到hb3344.com
12-02
ios app的ipv4连接到hb3344.com，
数位DP+二分搜索——HDU3943 K-th Nya Number
sanmingzi3344的专栏
07-19 960
题目链接： http://acm.hdu.edu.cn/showproblem.php?pid=3943 题目大意： 给定一个左开右闭区间，(p，q]，求在该区间上第k个Nya Number，如果找不到则输出Nya。 Nya Number的定义：含有x个4和y个7的数字。 比如：x=1，y=1，那么47，74，147。。。。都是Nya Number。 解题思路：
TCP Ports list (3498 ports in list)
热门推荐
zyb378747350的专栏
01-11 43万+
参考地址:https://www.gasmi.net/tcp.php TCP Ports list (3498 ports in list) 1 tcpmux TCP Port Service Multiplexer 2 compressnet Management Utility 3 compressnet Compression Process
python爬虫爬取百度首页
qq_45767476的博客
01-31 6265
import requests #导入requests模块 #利用爬虫代码爬去百度首页 #如果当前python文件作为入口程序执行时，则执行if语句下的代码 if __name__=='__main__': # 指定URL url = 'https://www.baidu.com' # 进行UA伪装，模拟浏览器,注意要将相应的User-Agent封装在一个字典中 headers = { 'User-Agent': 'Mozilla/5.0 (Windows
松下PLC   FP-X和PC通讯
weixin_33736832的博客
10-25 3036
今天刚测试了PLC 和PC 通讯，使用485协议。 前天公司将moxa cp-132串口通讯卡买了回来，这两天都在测试PLC和PC的串口通讯。在进行程序开发前，首先要安装好硬件和测试好各个用到的指令。 1、 moxa串口通讯卡，买回来的是cp-132,通讯协议422/485，实际我想用232会更好，因为用485的话，PLC不得不买一个485的接口模块。串口通讯卡的...
Linux 的 复制命令 【 cp 】 （copy）及其 (常用参数 -fp)
sun0322
10-23 6962
前言: 前一段时间执行一个batch的shell文件， 这个batch会读取一个文件， 但是发现，每次修改文件，再执行shell后 文件都会恢复修改之前的状态， 并且修改时间也是当前时间之前的时间。   于是看了一下shell文件，发现其中有cp命令   -------------------------------------------------- 简单介绍cp命令:
JS获取汉字首字母
weixin_30420305的博客
03-07 5501
//获取 汉字首字母 function makePy(str) { if (typeof (str) != "string") throw new Error(-1, "函数makePy需要字符串类型参数!"); var arrResult = new Array(); for (var i = 0, len = str.length; i ...
FP-growth
温染的笔记
04-27 1788
FP-growth FP-growth 算法能够更有效地挖掘数据，但不能用于发现关联规则。 FP-growth 基于 Apriori 算法构建，但在完成相同任务时采用了一些不同的技术。 Apriori：在每次循环的连接步中都要扫描数据集，来计算当前组合而成的项集的支持度。 FP-growth：只需要对数据库进行两次扫描，并将数据集存储在一个特定的称作 FP 树的数据结构。这种做法能够使得算法的执...
帮我写一段取消抖音全部喜欢的Python代码，我的抖音账号是19905027625，密码是xsy666...
最新发布
02-19
password = "xsy666"url = "https://www.douyin.com/aweme/v1/discover/"data = { 'cursor': 0, 'aweme_id': '', 'count': 20 }s = requests.Session()# 登录 login_url = '...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值