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;
    }
  • 相关阅读:
    适用于实数范围的中缀表达式的 +
    Django官方文档学习2——数据库及模板
    github命令
    千行代码入门Python
    Notepad++配置Python运行环境
    Python常用网页字符串处理技巧
    requests设置headers,proxies,cookies
    Django官方文档学习1——第一个helloworld页面
    笔记本键盘上没有break键的解决方案
    Python beautifulsoup模块
  • 原文地址:https://www.cnblogs.com/yinwuxiao/p/9741568.html
Copyright © 2011-2022 走看看