首页 > 其他分享 >多叉分治半在线卷积

多叉分治半在线卷积

时间:2023-02-24 22:46:23浏览次数:49  
标签:typedef 多叉 卷积 分治 poly POLY uint log Mod

唔,把以前一直口胡着的多叉分治半在线卷积实现了一下啊!

这个算法的思想是,我们分治计算贡献时,设目前区间长度为 \(n\),分成 \(B\) 个部分,并且用 cdq 的方式计算贡献,其中 \(B\) 是关于 \(n\) 的一个函数。

考虑平常 \(B=2\) 时最优,是因为计算两两支之间的贡献均为 \(O(n/B\log (n/B))\) 的,也即要进行 \(B^2\) 轮长度为 \(n/B\) 的 FFT,多叉没有优化效果。

此时,我们注意到转移时的乘法中本质不同的 poly 只有 \(O(B)\) 个,我们可以提前 DFT 算出各个多项式的点值,每次把点值对应做贡献,最后分别一轮 IDFT 插值计算即可。这样只用做大约 \(3B-3\) 轮长度为 \(2n/B\) 的 FFT(使用循环卷积!),以及计算点值间乘法贡献的 \(\Theta(B^2)\) 对。

从而单层的复杂度为 \(\Theta(nB+n\log(n/B))\),取 \(B=\log n\),即得单层复杂度为 \(\Theta(n\log n)\) 的。

可以证明递归层数是 \(\Theta(\log n/\log\log n)\) 级别的,总复杂度即为 \(\Theta(n\log^2n/\log\log n)\)。

按照 EI 论文,目前理论界存在一个神秘的 \(O(n\log^{1+\epsilon}n)\) 的半在线卷积算法(\(O(n\log ne^{2\sqrt{\log2\log\log n}})\)),感觉估计在 OI 圈内没啥用,除了指出在线算法往往和迭代法有相近的复杂度以及更小的常数外。

总不会这个东西是对那个 \(O(B^2)\) 对点值转移的过程继续调用当前算法吧?反正我是不会的。。。

放一下不等关系的代码。

封装版:

const ullt Mod=998244353,g=3;
typedef ConstMod::mod_ullt<Mod>modint;
typedef std::vector<modint>modvec;
typedef NTT_POLY::poly_NTT<Mod,g>poly;
typedef NTT_POLY::poly_eval<Mod,g>eval;
typedef NTT_POLY::poly_inter<Mod,g>inter;
typedef NTT_POLY::poly_cpd<Mod,g>cpd;
typedef NTT_POLY::poly_nums<Mod,g>nums;
typedef FWT_MODINT::FWT_Mod<Mod>FWT;
modint P[200005],Q[200005],Dp[100005];
bol A[100005];
voi dfs(uint l,uint r){
    if(r-l<=128){
        for(uint i=l;i<r;i++){
            Dp[i]+=Q[i+1];if(A[i])for(uint j=1;i+j<r;j++)Dp[i+j]-=Dp[i]*Q[j];
        }
        return;
    }
    uint B=1;
    if(r-l<=256)B=32;
    else if(r-l<=512)B=64;
    else if(r-l<=1024)B=128;
    else while(B*8<r-l)B<<=1;
    uint Cnt=(r-l+B-1)/B;
    poly::DIFDIT s(B<<1);
    modvec A[15],C[15];
    for(uint i=0;i<Cnt;i++){
        uint L=l+i*B,R=std::min(r,l+(i+1)*B);
        modvec User(B<<1);
        for(uint j=0;j<i;j++)for(uint k=0;k<(B<<1);k++)User[k]+=A[j][k]*C[i-j-1][k];
        s.dit(User);
        for(uint j=0;j<R-L;j++)Dp[L+j]-=User[j+B];
        dfs(L,R);
        if(i!=Cnt-1){
            A[i].resize(B<<1),C[i].resize(B<<1);
            for(uint j=0;j<(B<<1);j++)A[i][j]=Q[i*B+j];
            for(uint j=L;j<R;j++)if(::A[j])C[i][j-L]=Dp[j];
            s.dif(A[i]),s.dif(C[i]);
        }
    }
}
chr C[100005];
int main()
{
#ifdef MYEE
    freopen("QAQ.in","r",stdin);
    // freopen("QAQ.out","w",stdout);
#endif
    P[0]=1;for(uint i=1;i<=100001;i++)P[i]=P[i-1]*i;
    Q[100001]=P[100001].inv();for(uint i=100001;i;i--)Q[i-1]=Q[i]*i;
    uint n=0;scanf("%s",C);while(C[n])n++;
    uint cnt=0;
    for(uint i=0;i<n;i++)Dp[i]=0,A[i]=C[i]=='>',cnt+=A[i];
    dfs(0,n+1);
    (cnt&1?-Dp[n]*P[n+1]:Dp[n]*P[n+1]).println();
    return 0;
}

散装版:

const uint Mod=998244353,g=3;
inline uint chg(uint v){return v<Mod?v:v-Mod;}
class DIFDIT
{
    private:
        uint n;uint*G;
    public:
        DIFDIT():n(0),G(NULL){}
        DIFDIT(uint len)
        {
            n=1;while(n<len)n<<=1;
            G=new uint[n],G[0]=1;
            uint w=power<ullt>(g,(Mod-1)/n,Mod),*End=G+n;
            for(uint*_G=G+1;_G<End;_G++)*_G=(ullt)_G[-1]*w%Mod;
        }
        ~DIFDIT(){if(G!=NULL)delete[]G,G=NULL;}
        inline uint size(){return n;}
        voi dif(std::vector<uint>&x)
        {
            if(x.size()<n)x.resize(n);
            for(uint i=n>>1;i;i>>=1)for(uint j=0;j<n;j+=i<<1)
            {
                uint*w=G;
                for(uint k=0;k<i;k++,w+=n/(2*i))
                {
                    uint u=x[j|k],t=x[i|j|k];
                    x[j|k]=chg(u+t),x[i|j|k]=(ullt)(u+Mod-t)*(*w)%Mod;
                }
            }
        }
        voi dit(std::vector<uint>&x)
        {
            if(x.size()<n)x.resize(n);
            for(uint i=1;i<n;i<<=1)for(uint j=0;j<n;j+=i<<1)
            {
                uint*w=G;
                for(uint k=0;k<i;k++,w+=n/(2*i))
                {
                    uint t=(ullt)x[i|j|k]*(*w)%Mod;
                    x[i|j|k]=chg(x[j|k]+Mod-t),x[j|k]=chg(x[j|k]+t);
                }
            }
            for(uint i=1;i*2<n;i++)std::swap(x[i],x[n-i]);
            uint t=power<ullt>(n,Mod-2,Mod);for(uint i=0;i<n;i++)x[i]=(ullt)x[i]*t%Mod;
        }
};
uint Q[200005],Dp[100005];
bol A[100005];
voi dfs(uint l,uint r){
    if(r-l<=64){
        for(uint i=l;i<r;i++){
            Dp[i]=chg(Dp[i]+Q[i+1]);
            if(A[i])for(uint j=1;i+j<r;j++)
                Dp[i+j]=chg(Dp[i+j]+Mod-(ullt)Dp[i]*Q[j]%Mod);
        }
        return;
    }
    uint B=1;
    if(r-l<=128)B=16;
    else if(r-l<=256)B=32;
    else if(r-l<=512)B=64;
    else if(r-l<=1024)B=128;
    else while(B*8<r-l)B<<=1;
    uint Cnt=(r-l+B-1)/B;
    DIFDIT s(B<<1);
    std::vector<uint>A[7],C[7];
    for(uint i=0;i<Cnt;i++){
        uint L=l+i*B,R=std::min(r,l+(i+1)*B);
        std::vector<uint>User(B<<1);
        for(uint j=0;j<i;j++)for(uint k=0;k<(B<<1);k++)
            User[k]=(User[k]+(ullt)A[j][k]*C[i-j-1][k])%Mod;
        s.dit(User);
        for(uint j=0;j<R-L;j++)Dp[L+j]=chg(Dp[L+j]+Mod-User[j+B]);
        dfs(L,R);
        if(i!=Cnt-1){
            A[i].resize(B<<1),C[i].resize(B<<1);
            for(uint j=0;j<(B<<1);j++)A[i][j]=Q[i*B+j];
            for(uint j=0;j<B;j++)if(::A[j+L])C[i][j]=Dp[j+L];
            s.dif(A[i]),s.dif(C[i]);
        }
    }
}
chr C[100005];
int main()
{
#ifdef MYEE
    freopen("QAQ.in","r",stdin);
    // freopen("QAQ.out","w",stdout);
#endif
    uint n=0,v=1,cnt=0;scanf("%s",C);while(C[n])n++;
    for(uint i=2;i<=n+1;i++)v=(ullt)v*i%Mod;
    Q[n+1]=power<ullt>(v,Mod-2,Mod);
    for(uint i=n+1;i;i--)Q[i-1]=(ullt)Q[i]*i%Mod;
    for(uint i=0;i<n;i++)cnt+=A[i]=C[i]=='>';
    dfs(0,n+1);
    printf("%llu\n",cnt&1?chg(Mod-(ullt)Dp[n]*v%Mod):(ullt)Dp[n]*v%Mod);
    return 0;
}

目前 loj 最优解(301ms):链接

Update:被 Alpha 重测无了。/kel

标签:typedef,多叉,卷积,分治,poly,POLY,uint,log,Mod
From: https://www.cnblogs.com/myee/p/multi-div-semi-relaxed-mul.html

相关文章

  • 题解 北大集训2018 点分治
    题意给定\(n\)个点的树,求点分治方案数,对\(10^9+7\)取模。两种方案不同当且仅当点分树不同。\(1\leqn\leq5000\)题解考虑怎样合并两个点分治方案。如果我们有......
  • 【YBT2023寒假Day15 C】缺口一样(数论)(莫队)(根号分治)
    缺口一样题目链接:YBT2023寒假Day15C题目大意给你一个序列,多次询问,每次问你一个区间这里面所有非空点集的最大公约数之积,对质数取模。思路首先质因子之间是独立的,考虑......
  • 点分治练习题单(动态更新)
    传送门有点难,慢慢做。1.P2634[国家集训队]聪聪可可比板子要简单一点,分治时寻找路径时用桶记录模数为\(0,1,2\)的个数,再进行\(01\)背包即可。统计答案时由于两点可......
  • 浅谈点分治
    非常有趣的一个知识点。所谓点分治,也就是在树上进行分治的操作。分治可以处理许多问题,各种区间类问题很多都可以利用分治思想,例如线段树就是利用分治处理许多区间性问题......
  • 微服务拆分治理最佳实践
    作者:京东零售徐强黄威张均杰背景部门中维护了一个老系统,功能都耦合在一个单体应用中(300+接口),表也放在同一个库中(200+表),导致系统存在很多风险和缺陷。经常出现问题:如......
  • 基于python 的 分治算法 解决 大数乘法问题 -- 浅浅记录
    分治方法总的来说还是递归的调用,将大的问题分解为小的问题1#通过分治计算n个0的字符串2defzero(n):3ifn==0:4return""5ifn==1:......
  • 微服务拆分治理最佳实践
    作者:京东零售徐强黄威张均杰背景部门中维护了一个老系统,功能都耦合在一个单体应用中(300+接口),表也放在同一个库中(200+表),导致系统存在很多风险和缺陷。经常出现问题......
  • 多叉分治半在线卷积
    内容咕了。放一下不等关系的代码。封装版:constulltMod=998244353,g=3;typedefConstMod::mod_ullt<Mod>modint;typedefstd::vector<modint>modvec;typedefNTT_PO......
  • Solution Luogu6097 子集卷积
    其实是暴力。因为这是模板题,所以模板的前置知识也要讲。前置知识:FWT计算或卷积。这里只需要掌握快速计算或卷积的方法,所以内容较少。如果向了解更多(比如异或卷积)的话......
  • Tensorflow中常用的卷积函数
    卷积函数(1)计算N维卷积的和tf.nn.convolution(input,filter,padding,strides=None,dilation_rate=None,name=None,data_format=None)(2)对一个四维的输入数据input和卷积核......