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

题目大意

　　对于一个长度为 \(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

时间： 2024-10-22 19:06:01

【XSY3344】连续段 DP 牛顿迭代 NTT的相关文章

【XSY2680】玩具谜题 NTT 牛顿迭代

题目描述 小南一共有\(n\)种不同的玩具小人,每种玩具小人的数量都可以被认为是无限大.每种玩具小人都有特定的血量,第\(i\)种玩具小人的血量就是整数\(i\).此外,每种玩具小人还有自己的攻击力,攻击力可以是任意非负整数,且两种不同的玩具小人的攻击力可以相同.我们把第\(i\)种玩具小人的血量和攻击力表示成\(a_i\)和\(b_i\). 为了让玩具小人们进行战斗,小南打算把一些小人选出来,编成队伍.一个队伍可以表示成一个由玩具小人组成的序列:\((p_1,p_2,\ldots,p_l)\)

【loj6538】烷基计数 加强版 加强版 Burnside引理+多项式牛顿迭代

别问我为啥突然刷了道OI题,也别问我为啥花括号不换行了... 题目描述 求含 $n$ 个碳原子的本质不同的烷基数目模 $998244353$ 的结果.$1\le n\le 10^5$ . 题解 Burnside引理+多项式牛顿迭代 不考虑同构的话,很容易想到dp方程 $\begin{cases}f_0=1\\f_i=\sum\limits_{j+k+l+1=i}f_jf_kf_l\end{cases}$ . 考虑同构,可以通过容斥原理,大力讨论一下容斥系数.一个更简单的方法是考虑Burnside

HDU 1231 最大连续子序列 DP题解

典型的DP题目,增加一个额外要求,输出子序列的开始和结尾的数值. 增加一个记录方法,nothing special. 记录最终ans的时候,同时记录开始和结尾下标: 更新当前最大值sum的时候,更新开始节点. const int MAX_N = 10001; long long arr[MAX_N]; int N, sta, end; long long getMaxSubs() { long long sum = 0, ans = LLONG_MIN; int ts = 0; for (int

求方程解，牛顿迭代和二分

牛顿迭代 #include<iostream> #include<string.h> #include<math.h> using namespace std; float f(float x){ return (pow(x,3)-5*pow(x,2)+16*x+80); } float f1(float x){ return (3*pow(x,2)-5*x+16); } int main(){ // x*x*x-5*x*x+16*x+80; float x=1,x1,

HLG 2116 Maximum continuous product （最大连续积 DP）

链接:  http://acm.hrbust.edu.cn/index.php?m=ProblemSet&a=showProblem&problem_id=2116 Description Wind and his GF(game friend) are playing a small game. They use the computer to randomly generated a number sequence which only include number 2,0 and -

POJ 3764 The xor-longest Path ( 字典树应用—— 求连续段相异或最大最小的线性算法)（好题）

题意:已知:给出n个结点的树,定义:两结点间的权值为两点之间所有边相异或的值.求:树中的某两点间的最大权值. 思路:先说简单一点的题:有道CowXor,是一串线性序列,求某连续段异或的最大值,这题的思路是先求前i项序列相异或的值Si,所以x到y的连续异或就是Sx^Sy ,因为a^b = (a ^ c) ^ (b ^ c). 这题同样是这个思路把线性拓展到树上,先求任何点到某一定点的连续异或值,比如选根结点0,所以这时候有两种情况,1.x,y的路径通过了根结点,显然正确.2.x,y的路径不通过根结

poj 3111 K Best ，二分，牛顿迭代

poj 3111  K Best 有n个物品的重量和价值分别是wi和vi.从中选出k个物品使得单位重量的价值最大. 题解: 1.二分做法 2.牛顿迭代 效率比较: 二分做法: 转换成判断是否存在选取K个物品的集合S满足下面的条件: sigma(vi) / sigma(wi) >= x   {vi∈S, wi∈S} -->   simga(vi - x*wi) >= 0 这样我们对  yi= vi - x*wi {1<=i<=n}从大到小排序,计算sum(yi) {1<=

牛顿迭代（开方，板子）

就是正常的牛顿迭代求开方,牛顿迭代日后更新 开方: #include <bits/stdc++.h> using namespace std; const double eps=1e-9; double sqr(double x) { double k=x; while(k*k-x>eps) k=0.5*(k+x/k); return k; } int main() { double x; while(cin>>x) printf("%.9f\n",sqr

HDU 2899 Strange fuction(牛顿迭代)

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)Total Submission(s): 9333    Accepted Submission(s): 6352 Problem Description Now, here is a fuction:  F(x) = 6 * x^7+8*x^6+7*x^3+5*x^2-y*x (0 <= x <=100)Can you find