zoukankan      html  css  js  c++  java
  • HDU3439 Sequence

    今天下午学习了二项式反演,做了一道错排的题,开始了苦逼的经历。

    显然答案是C(︀n,k)︀*H(n − k).
    其中H(i)为长度为i的错排序列

    然后经过课件上一番二项式反演的推导

    我就写了个扩展卢卡斯然后交上去了。

    一直t啊.....

    我算了算复杂度差不多是O(T*P*log^3P)

    后来剪了剪枝,应该低了点。

    还是t啊.....

    我搜了搜题解发现没有我这么写的。

    看了一下错排是有规律的,果然还是打表大法吼啊。

    发个正解

    #include<bits/stdc++.h>
    using namespace std;
    typedef long long ll;
    ll qmod(ll n,ll m,ll p)
    {
        ll ans=1;
        while(m)
        {
            if(m&1)ans=ans*n%p;
            n=n*n%p;m>>=1;
        }
        return ans;
    }
    void exgcd(ll a,ll b,ll &x,ll &y)
    {
        if(!b){
            x=1;y=0;return;
        }
        exgcd(b,a%b,y,x);y-=a/b*x;
    }
    ll inv(ll n,ll p)
    {
        if(!n)return 0;
        ll a=n,b=p,x=0,y=0;
        exgcd(a,b,x,y);
        x=(x%b+b)%b;
        if(!x)x+=b;
        return x;
    }
    ll mul(ll n,ll pi,ll pk)
    {
        if(!n)return 1ll;
        ll ans=1;
        if(n/pk)
        {
            for(ll i=2;i<=pk;++i)
                if(i%pi)ans=ans*i%pk;
            ans=qmod(ans,n/pk,pk)%pk;
        }
        for(ll i=2;i<=n%pk;++i)
            if(i%pi)ans=ans*i%pk;
        return ans*mul(n/pi,pi,pk)%pk;
    }
    int C(ll n,ll m,ll p,ll pi,ll pk)
    {
        if(m>n)return 0;
        ll a=mul(n,pi,pk),b=mul(m,pi,pk),c=mul(n-m,pi,pk);
        ll k=0,ans;
        for(ll i=n;i;i/=pi)k+=i/pi;
        for(ll i=m;i;i/=pi)k-=i/pi;
        for(ll i=n-m;i;i/=pi)k-=i/pi;
        ans=a*inv(b,pk)%pk*inv(c,pk)%pk*qmod(pi,k,pk)%pk;
        return ans*(p/pk)%p*inv(p/pk,pk)%p;//CRT
    }
    bool v[100005];
    int pp[100005],cnt;
    void pri()
    {
        for(int i=2;i<=100000;++i)
        {
            if(!v[i])
            {
                pp[++cnt]=i;
            }
            for(int j=1;j<=cnt&&i*pp[j]<=100000;++j)
            {
                v[i*pp[j]]=1;
                if(i%pp[j]==0)break;
            }
        }
    }
    ll calc(ll n,ll m,ll p)
    {
        ll ans=0;
        for(ll x=p,i=1;i<=cnt&&x;++i)
        {
            if(x==1)break;
            if(x%pp[i]==0)
            {
                ll num=1;
                while(x%pp[i]==0)x/=pp[i],num*=pp[i];
                ans=(ans+C(n,m,p,pp[i],num))%p;
            }  
        }
        return ans;
    }
    ll F(ll x,ll p)
    {
           ll ans=0;
        if(x==0)return 1;
        x=x%(2*p);
        if(x==0)x=2*p;
        for(int i=2;i<=x;++i)
        ans=(ans*i+(i%2==0?1:-1))%p; 
        return (ans+p)%p;
    }
    int main()
    {
        ll n,m,p,t,ans=0;
        scanf("%I64d",&t);pri();v[1]=1;
        for(int ii=1;ii<=t;++ii)
        {
            scanf("%I64d%I64d%I64d",&n,&m,&p);
            ans=calc(n,m,p)%p;
            ans=ans*F(n-m,p)%p;
            printf("Case %d: %I64d
    ",ii,ans);
        }
        return 0;
    }

    再补个我的辣鸡程序,路过的dalao帮忙看看也中啊、

    #include<bits/stdc++.h>
    using namespace std;
    typedef long long ll;
    ll qmod(ll n,ll m,ll p)
    {
        ll ans=1;
        while(m)
        {
            if(m&1)ans=ans*n%p;
            n=n*n%p;m>>=1;
        }
        return ans;
    }
    void exgcd(ll a,ll b,ll &x,ll &y)
    {
        if(!b){
            x=1;y=0;return;
        }
        exgcd(b,a%b,y,x);y-=a/b*x;
    }
    ll inv(ll n,ll p)
    {
        if(!n)return 0;
        ll a=n,b=p,x=0,y=0;
        exgcd(a,b,x,y);
        x=(x%b+b)%b;
        if(!x)x+=b;
        return x;
    }
    ll mul(ll n,ll pi,ll pk)
    {
        if(!n)return 1ll;
        ll ans=1;
        if(n/pk)
        {
            for(ll i=2;i<=pk;++i)
                if(i%pi)ans=ans*i%pk;
            ans=qmod(ans,n/pk,pk)%pk;
        }
        for(ll i=2;i<=n%pk;++i)
            if(i%pi)ans=ans*i%pk;
        return ans*mul(n/pi,pi,pk)%pk;
    }
    int C(ll n,ll m,ll p,ll pi,ll pk)
    {
        if(m>n)return 0;
        ll a=mul(n,pi,pk),b=mul(m,pi,pk),c=mul(n-m,pi,pk);
        ll k=0,ans;
        for(ll i=n;i;i/=pi)k+=i/pi;
        for(ll i=m;i;i/=pi)k-=i/pi;
        for(ll i=n-m;i;i/=pi)k-=i/pi;
        ans=a*inv(b,pk)%pk*inv(c,pk)%pk*qmod(pi,k,pk)%pk;
        return ans*(p/pk)%p*inv(p/pk,pk)%p;//CRT
    }
    bool v[100005];
    int pp[100005],cnt;
    void pri()
    {
        for(int i=2;i<=100000;++i)
        {
            if(!v[i])
            {
                pp[++cnt]=i;
            }
            for(int j=1;j<=cnt&&i*pp[j]<=100000;++j)
            {
                v[i*pp[j]]=1;
                if(i%pp[j]==0)break;
            }
        }
    }
    ll calc(ll n,ll m,ll p)
    {
        ll ans=0;
        for(ll x=p,i=1;i<=cnt&&x;++i)
        {
            if(x==1)break;
            if(x%pp[i]==0)
            {
                ll num=1;
                while(x%pp[i]==0)x/=pp[i],num*=pp[i];
                ans=(ans+C(n,m,p,pp[i],num))%p;
            }  
        }
        return ans;
    }
    int main()
    {
        ll n,m,p,t,ans=0;
        scanf("%I64d",&t);pri();v[1]=1;
        for(int ii=1;ii<=t;++ii)
        {
            scanf("%I64d%I64d%I64d",&n,&m,&p);
            ans=calc(n,m,p)%p;n-=m;ll pre=0;ll num=1,pos=max(0ll,n-p);
            if(!ans)
            {
                 printf("Case %d: %I64d
    ",ii,pre*ans%p);continue;
            }
            for(ll k=n;k>=pos;--k)
            {
                if(n!=k)num=num*(n-k)%p;
                if(k&1ll)pre=(pre-calc(n,k,p)%p*num%p+p)%p;
                else pre=(pre+calc(n,k,p)%p*num%p)%p;
                if(!num)break;
            }
            printf("Case %d: %I64d
    ",ii,pre*ans%p);
        }
        return 0;
    }

    好吧,蒟蒻苦逼的一下午。

    同时纪念衡水人民224起义。

  • 相关阅读:
    python操作elasticsearch
    php源码的编译
    linux 访问windows 共享文件
    list容器排除重复单词的程序
    求组合数m_n
    简单堆排序
    快速排序
    判断点在直线左侧或者右侧
    求取点到直线的距离
    求给定三个点的夹角
  • 原文地址:https://www.cnblogs.com/nbwzyzngyl/p/8474596.html
Copyright © 2011-2022 走看看