zoukankan      html  css  js  c++  java
  • loj6077

    题解:

    网上的做法好像都是容斥

    那就先说一下容斥

    首先问题等价于求下面这个式子的方案数

    $$sum_{i=1}^{n} ai (0<ai<i) =k$$

    直接$dp$复杂度是$nk$的,无法接受

    我们考虑容斥

    答案=所有-至少1个不满足+至少2个不满足-....

    而不满足等价于$ai>=i$ 那么我们只需要令$b[i]=a[i]-i$就可以与其他等价了

    实际上这样子看上去这个容斥并没有用处

    因为要$dp$出容斥系数看上去至少是$O(n^2)$

    现在我们需要算出,用i个不相同的数相加=k的方案数

    我们注意一件事情,i的取值只有$sqrt{k}$,所以如果状态设计得好,是可以的

    要求数不同有一个经典套路就是 把操作转化成

    1.给当前数组内全体数+1

    2.给当前数组内全体数+1再在末尾放一个1

    于是$f[i][j]$表示数组内有i个数,和为j的方案数 $$f[i][j]=f[i-1][j-i]+f[i][j-i]$$复杂度$nsqrt{n}$

    然后你写一下,肯定会发现出问题了(我没写这种做法但下面这种的时候遇到了同样的问题)

    因为这么做可能有数>n了

    怎么保证没有数>n呢,我们只需要令每个时刻的$f[i][j]$都是<=n的

    那么出现大于n的只可能是n+1,$f[i][j]-=f[i-1][j-n-1]$就可以了

    还是比较难想的。。

    这题还有个做法是生成函数,在思维上就简单很多

    $$f(x)=frac{prod_{i=1}^{n} (1-x^i)}{(1-x)^n}$$ 

    然后要求这个东西,暴力是n^2logn的

    然后乘法比较显然的是两边取ln,得到

    $$ln(f(x))=sum_{i=1}^{n} {ln(1-x^i)} - n*ln(1-x)$$

    这样子我们发现要解决的就是$ln(1-x^i)$

    这个东西就两项但我们要用$nlogn$的多项式求逆是不是太浪费了

    于是打表找规律,可以发现(不过要是不知道结论考场谁有时间去打个多项式求ln找规律啊。。)

    $ln(1-x)=frac{1}{1}x+frac{1}{2}x^2+frac{1}{3}x^3+...$

    $ln(1-x^2)=frac{1}{1}x^2+frac{1}{2}x^4+frac{1}{3}x^6+...$

    下面的规律同理

    于是我们可以利用筛法在$nlogn$时间内求出

    当然由于$$f(x)=sum {d|x} {}{ g1(d)g2(frac{x}{d}) }$$ 其中$g1(i)=frac{1}{i},g2(i)=1$

    因为$g1(i),g2(i)$是完全积性函数(只要是积性就可以了),所以他们的卷积也是积性函数

    所以可以用线性筛做到$O(n)$ 我们类似于做约数和再维护一个h(i)表示不包括i的最小素因子的g1的和

    于是就是多项式exp模板题了

    #pragma GCC optimize(2)
    #include <bits/stdc++.h>
    using namespace std;
    #define rint register int
    #define IL inline
    #define rep(i,h,t) for(int i=h;i<=t;i++)
    #define dep(i,t,h) for(int i=t;i>=h;i--)
    #define ll long long
    #define me(x) memset(x,0,sizeof(x))
    #define mep(x,y) memcpy(x,y,sizeof(y))
    #define mid (t<=0?(h+t-1)/2:(h+t)/2)
    namespace IO{
        char ss[1<<24],*A=ss,*B=ss;
        IL char gc()
        {
            return A==B&&(B=(A=ss)+fread(ss,1,1<<24,stdin),A==B)?EOF:*A++;
        }
        template<class T> void read(T &x)
        {
            rint f=1,c; while (c=gc(),c<48||c>57) if (c=='-') f=-1; x=(c^48);
            while (c=gc(),c>47&&c<58) x=(x<<3)+(x<<1)+(c^48); x*=f; 
        }
        char sr[1<<24],z[20]; ll Z,C1=-1;
        template<class T>void wer(T x)
        {
            if (x<0) sr[++C1]='-',x=-x;
            while (z[++Z]=x%10+48,x/=10);
            while (sr[++C1]=z[Z],--Z);
        }
        IL void wer1()
        {
            sr[++C1]=' ';
        }
        IL void wer2()
        {
            sr[++C1]='
    ';
        }
        template<class T>IL void maxa(T &x,T y) {if (x<y) x=y;}
        template<class T>IL void mina(T &x,T y) {if (x>y) x=y;} 
        template<class T>IL T MAX(T x,T y){return x>y?x:y;}
        template<class T>IL T MIN(T x,T y){return x<y?x:y;}
    };
    using namespace IO;
    const double ee=1.00000000000000;
    const double pi=acos(-1.0);
    const int N=4e5+10;
    const int mo=1e9+7;
    int r[N],n,m,l,inv[N],n1[N],n2[N];
    struct cp{
        double a,b;
        IL cp operator +(const cp &o) const
        {
            return (cp){a+o.a,b+o.b};
        }
        IL cp operator -(const cp &o) const
        {
            return (cp){a-o.a,b-o.b};
        }
        IL cp operator *(register const cp &o) const
        {
            return (cp){a*o.a-b*o.b,o.a*b+o.b*a};
        }
    }a[N],b[N],c[N],d[N];
    IL int fsp(int x,int y)
    {
        ll now=1;
        while (y)
        {
            if (y&1) now=now*x%mo;
            x=1ll*x*x%mo;
            y>>=1;
        }
        return now;
    }
    IL void clear()
    {
        for (int i=0;i<=n;i++) a[i].a=a[i].b=b[i].a=b[i].b=c[i].a=c[i].b=d[i].a=d[i].b=0;
    }
    cp *w[N],tmp[N*2];
    int p;
    IL void init()
    {
        cp *now=tmp;
        for (int i=1;i<=p;i<<=1)
        {
            w[i]=now;
            for (int j=0;j<i;j++) w[i][j]=(cp){cos(pi*j/i),sin(pi*j/i)};
            now+=i;
        }
    }
    IL void fft_init()
    {
        l=0; for (n=1;n<=m;n<<=1) l++;
        for (int i=0;i<n;i++) r[i]=(r[i/2]/2)|((i&1)<<(l-1));
    }
    void fft(cp *a,int o)
    {
        for (int i=0;i<n;i++) if (i>r[i]) swap(a[i],a[r[i]]);
        for (int i=1;i<n;i<<=1)
            for (int j=0;j<n;j+=(i*2))
            {
              cp *x1=a+j,*x2=a+i+j,*W=w[i];
              for (int k=0;k<i;k++,x1++,x2++,W++)
              {
                  cp x=*x1,y=(cp){(*W).a,(*W).b*o}*(*x2); 
                  *x1=x+y,*x2=x-y;
              }
            }
        if (o==-1) for(int i=0;i<n;i++) a[i].a/=n;
    }
    IL void getcj(int *A,int *B,int len)
    {
        rep(i,0,len)
        {
            A[i]=(A[i]+mo)%mo,B[i]=(B[i]+mo)%mo;
        }
        for (int i=0;i<len;i++)
        {
           a[i]=(cp){A[i]&32767,A[i]>>15};
           b[i]=(cp){B[i]&32767,B[i]>>15};
        }
        m=len*2; fft_init();
        fft(a,1); fft(b,1);
        for (int i=0;i<n;i++)
        {
            int j=(n-1)&(n-i);
            c[j]=(cp){0.5*(a[i].a+a[j].a),0.5*(a[i].b-a[j].b)}*b[i];
            d[j]=(cp){0.5*(a[i].b+a[j].b),0.5*(a[j].a-a[i].a)}*b[i];
        }
        fft(c,1); fft(d,1);
        double inv=ee/n;
        rep(i,0,n) c[i].a*=inv,c[i].b*=inv;
        rep(i,0,n) d[i].a*=inv,d[i].b*=inv;
        rep(i,0,len)
        {
            ll a1=c[i].a+0.5,a2=c[i].b+0.5;
            ll a3=d[i].a+0.5,a4=d[i].b+0.5;
            B[i]=(a1+((a2+a3)%mo<<15)+((a4%mo)<<30))%mo;
        }
        clear();
    }
    void getinv(int *A,int *B,int len)
    {
        if (len==1) { B[0]=fsp(A[0],mo-2); return;};
        getinv(A,B,(len+1)/2);
        int C[N]={};
        rep(i,0,len-1) C[i]=A[i];
        getcj(B,C,len);
        getcj(B,C,len);
        for (int i=0;i<len;i++) B[i]=((2ll*B[i]-C[i])%mo+mo)%mo;
    }
    IL void getDao(int *a,int *b,int len)
    {
        for (int i=1;i<len;i++) b[i-1]=1ll*i*a[i]%mo;
        b[len-1]=0;
    }
    IL void getjf(int *a,int *b,int len)
    {
        for (int i=0;i<len;i++) b[i+1]=1ll*a[i]*inv[i+1]%mo;
        b[0]=0;
    }
    IL void getln(int *A,int *B,int len)
    {
        int C[N]={},D[N]={};
        getDao(A,C,len);
        getinv(A,D,len);
        getcj(C,D,len);
        getjf(D,B,len);
    }
    IL void getexp(int *A,int *B,int len)
    {
        if (len==1) {B[0]=1; return;}
        getexp(A,B,(len+1)>>1);
        int C[N]={};
        getln(B,C,len);
        for(int i=0;i<len;i++) C[i]=(-C[i]+A[i])%mo;
        C[0]=(C[0]+1)%mo;
        getcj(C,B,len);
    }
    int main()
    {
        freopen("1.in","r",stdin);
        freopen("1.out","w",stdout);
        inv[1]=1;
        rep(i,2,1e5+20) inv[i]=(1ll*inv[mo%i]*(mo-(mo/i)))%mo;
        int n,k;
        read(n); read(k); k++;
        n1[0]=1; n1[1]=-1;
        p=k<<1;
        init();
        getln(n1,n2,k);
        rep(i,0,k) n1[i]=-1ll*n*n2[i]%mo;
        rep(i,1,n)  
        {
          for(int j=1;j*i<=k;j++)
            (n1[i*j]-=inv[j])%=mo;
        }
        me(n2); 
        getexp(n1,n2,k);
        cout<<(n2[k-1]+mo)%mo<<endl;
        return 0;
    }
  • 相关阅读:
    win10 uwp 弹起键盘不隐藏界面元素
    win10 uwp 存放网络图片到本地
    win10 uwp 存放网络图片到本地
    sublime Text 正则替换
    sublime Text 正则替换
    win10 uwp 绘图 Line 控件使用
    win10 uwp 绘图 Line 控件使用
    AJAX 是什么?
    什么是 PHP SimpleXML?
    PHP XML DOM:DOM 是什么?
  • 原文地址:https://www.cnblogs.com/yinwuxiao/p/9741568.html
Copyright © 2011-2022 走看看