操作 - 分治NTT

最新推荐文章于 2024-01-26 02:49:59 发布

原创最新推荐文章于 2024-01-26 02:49:59 发布 · 707 阅读

1 ·

CC 4.0 BY-SA版权

快速傅里叶变换同时被 2 个专栏收录

23 篇文章

订阅专栏

分治

10 篇文章

订阅专栏

题目大意：
有n个操作和一个数字初始为0，每个操作形如有p的概率加a，1-p的概率乘b。现在随机打乱这个操作序列，问操作完数值的期望是多少。 $n≤105n\le10^5$
题解：
每次就是让x变成kx+b。手玩一下发现答案就是：
$1n!∑i=1nbi∑j=0n−1j!(n−j−1)![xj]∏t=1n(x+kt)x+ki\frac1{n!}\sum_{i=1}^nb_i\sum_{j=0}^{n-1}j!(n-j-1)![x^j]\frac{\prod_{t=1}^n(x+k_t)}{x+k_i}$
后面那一坨交换求和顺序可以直接分治NTT。

 
#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define Rep(i,v) rep(i,0,(int)v.size()-1)
#define lint long long
#define p 998244353
#define ull unsigned lint
#define db long double
#define pb push_back
#define mp make_pair
#define fir first
#define sec second
#define gc getchar()
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
using namespace std;
typedef pair<int,int> pii;
typedef set<int>::iterator sit;
inline int inn()
{
    int x,ch;while((ch=gc)<'0'||ch>'9');
    x=ch^'0';while((ch=gc)>='0'&&ch<='9')
        x=(x<<1)+(x<<3)+(ch^'0');return x;
}
const int N=100010;int k[N],b[N];
inline int show(int *a,int n) { debug(n)sp<<", ";rep(i,0,n-1) cerr<<a[i]sp;cerr ln;return 0; }
inline int fast_pow(int x,int k,int ans=1) { for(;k;k>>=1,x=(lint)x*x%p) (k&1)?ans=(lint)ans*x%p:0;return ans; }
inline int *newInt(int n) { int *t=new int[n];memset(t,0,sizeof(int)*n);return t; }
inline int clr(int *a,int n) { return memset(a,0,sizeof(int)*n),0; }
inline int cpy(int *a,int *b,int n) { return memcpy(a,b,sizeof(int)*n),0; }
namespace NTT_space{
    const int N=131072+100,LOG=20;
    int *dwg[LOG],*dwgi[LOG];
    struct dwg_preluder{
        dwg_preluder()
        {
            for(int i=2,c=1;i+5<N;i<<=1,c++)
            {
                dwg[c]=newInt(i>>1),dwgi[c]=newInt(i>>1);
                int *d=dwg[c],*di=dwgi[c];d[0]=di[0]=1;
                int w=fast_pow(3,(p-1)/i),wn=fast_pow(w,p-2);
                rep(j,1,(i>>1)-1) d[j]=(lint)d[j-1]*w%p,di[j]=(lint)di[j-1]*wn%p;
            }
        }
    //  ~dwg_preluder() { for(int i=2,c=1;i+5<N;i<<=1,c++) delete dwg[c],delete dwgi[c]; }
    }__dwg_preluder;
    namespace NTT_SPACE{
        int r[N];
        inline int pre(int m)
        {
            int L=0,n=1;while(n<m) n<<=1,L++;
            rep(i,1,n-1) r[i]=(r[i>>1]>>1)|((i&1)<<(L-1));
            return n;
        }
        inline int NTT(int *a,int n,int s)
        {
            rep(i,1,n-1) if(i<r[i]) swap(a[i],a[r[i]]);
            for(int i=2,c=1;i<=n;i<<=1,c++)
            {
                int *d=(s>0?dwg[c]:dwgi[c]);
                for(int j=0,t=i>>1;j<n;j+=i) for(int k=0,x,y;k<t;k++)
                    x=a[j+k],y=(lint)d[k]*a[j+k+t]%p,
                    a[j+k]=(x+y>=p?x+y-p:x+y),a[j+k+t]=(x-y<0?x-y+p:x-y);
            }
            if(s<0) for(int i=0,v=fast_pow(n,p-2);i<n;i++) a[i]=(lint)a[i]*v%p;return 0;
        }
    }using NTT_SPACE::pre;using NTT_SPACE::NTT;
    namespace TMS_space{
        int A[N],B[N],as[N],bs[N];
        inline int trans(int *a,int *A,int m,int n,int s) { return cpy(A,a,m),clr(A+m,n-m),NTT(A,n,s); }
        inline int tms(int *a,int m1,int *b,int m2,int* &c,int m3=-1)
        {
            if(m3<0) m3=m1+m2-1;int n=pre(m1+m2-1);
            trans(a,A,m1,n,1),trans(b,B,m2,n,1);
            rep(i,0,n-1) A[i]=(lint)A[i]*B[i]%p;
            if(c==NULL) c=newInt(m3);
            return NTT(A,n,-1),cpy(c,A,m3),0;
        }
    }using TMS_space::tms;
    inline int pls(int *a,int n,int *b,int m,int* &c)
    {
        if(c==NULL) c=newInt(max(n,m));
        rep(i,0,max(n,m)-1) c[i]=(i<n?a[i]:0)+(i<m?b[i]:0),(c[i]>=p?c[i]-=p:0);
        return 0;
    }
}using NTT_space::tms;using NTT_space::pls;
int *ans[N],*pd[N],ansl[N],pdl[N],fac[N];
inline int solve(int l,int r)
{
    if(l==r) return ans[l]=newInt(ansl[l]=1),ans[l][0]=b[l],pd[l]=newInt(pdl[l]=2),pd[l][0]=k[l],pd[l][1]=1;
    int mid=(l+r)>>1,lc=l,rc=mid+1,*t1,*t2,t1l,t2l;solve(l,mid),solve(mid+1,r);
     
    t1=NULL,tms(ans[lc],ansl[lc],pd[rc],pdl[rc],t1,t1l=ansl[lc]+pdl[rc]-1),//delete ans[lc],
    t2=NULL,tms(pd[lc],pdl[lc],ans[rc],ansl[rc],t2,t2l=pdl[lc]+ansl[rc]-1),//delete ans[rc],
    ans[lc]=NULL,pls(t1,t1l,t2,t2l,ans[lc]),ansl[lc]=max(t1l,t2l);
     
    t1=NULL,tms(pd[lc],pdl[lc],pd[rc],pdl[rc],t1),//delete pd[lc],delete pd[rc],
    pd[lc]=t1,pdl[lc]+=pdl[rc]-1;return 0;
}
int main()
{
    int n=inn(),x,y,z,res=0;
    rep(i,1,n) x=inn(),y=inn(),z=inn(),k[i]=(x+(1ll-x)*z)%p,(k[i]<0?k[i]+=p:0),b[i]=(lint)x*y%p;
    solve(1,n);rep(i,fac[0]=1,n) fac[i]=(lint)i*fac[i-1]%p;
    rep(i,0,n-1) res=(res+(lint)fac[i]*fac[n-i-1]%p*ans[1][i])%p;
    res=(lint)res*fast_pow(fac[n],p-2)%p;return !printf("%d\n",res);
}