[JZOJ5330]密码

题目大意

给定n,p,k,求 i<=j<=n[pk|Cji]
其中, n<=101000p<=109pk<=109

分析

首先我们看一个组合数,怎么样才能被 pk 整除呢?我们有 Cnm=m!n!(mn)! ,我们可以分别算出上下阶乘因子含有多少个p,记为 cntp() 再相减。对于一个阶乘 n!cntp(n)=i=1[npi] 。那么如果cnt(m)-cnt(n)-cnt(m-n)>=k,该组合数能被 pk 整除。
这样不能看出什么,我们把cnt都合起来,原式变成 i=1[mpi][npi][mnpi] ,可以发现,对于一个i,多项式的值只有可能是0或者1,1就意味着有因子p了嘛。有lucas定理,我们容易想到把n,m,m-n看成p进制数,那么多项式的值为1的时候,就意味着二进制下 n,mni1m=n+(mn)im1
这个其实是库默尔定理。
那么我们尝试构造n和m-n,问题就变成:有多少对(x,y),满足十进制时 x,y<=n ,p进制时 x+y 有K次或更多次进位。
然后就可以用数位DP做了,P进制下的数位DP。考虑设计状态。
状态[i][j][0/1][0/1]分别表示从高位往低位考虑到第i位,一共有j位进位了,x+y是否贴着上界N,i+1位是否进位。f[i][j][0,1][0,1]就表示此时的方案数。然后算的时候,对于每个状态我们计算x和y在i+1位中的满足条件方案数,乘上当前的f,转移到新的f。转移十分恶心,我推了好久,系数是一堆等差数列。

代码

#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#include<bitset> 
using namespace std;
typedef long long ll;
typedef double db;
#define fo(i,j,k) for(i=j;i<=k;i++)
#define fd(i,j,k) for(i=j;i>=k;i--)
const int N=1005,mo=1e9+7;
ll f[N][N][2][2],a[N],b[N],a00,a01,a10,a11,ax,ay,bx,by,siz,o00,o10,i,p,K,n,j,ans;
char s[N];
void predo()
{
    while (a[0])
    {
        fd(i,a[0],2)
        {
            a[i-1]+=a[i]%p*10;
            a[i]/=p;
        }
        b[++b[0]]=a[1]%p;
        a[1]/=p;
        while (a[0]&&a[a[0]]==0) a[0]--;
    }
}
int main()
{
    freopen("t3.in","r",stdin);
    //freopen("password.out","w",stdout);
    scanf("%s %d %d",s+1,&p,&K);
    a[0]=strlen(s+1);
    fo(i,1,a[0]) a[a[0]-i+1]=s[i]-'0';
    predo();
    n=b[0];
    fo(i,1,n/2) swap(b[i],b[n-i+1]);
    f[1][0][0][0]=(b[1]+1)*b[1]/2%mo;
    f[1][0][1][0]=b[1]+1;
    f[1][1][0][1]=(b[1])*(b[1]-1)/2%mo;//等下重新推 
    f[1][1][1][1]=b[1];
    o00=a00=(1+p)*p/2;
    a01=p*p-a00;
    o10=a10=p*(p-1)/2;
    a11=p*p-a10;
    a00%=mo;a01%=mo;a10%=mo;a11%=mo;
    // 第i+1位还没有考虑方案,只考虑第i位。 
    fo(i,1,n)
        fo(j,0,i)
        {
                //f[i+1][j+1][0][1], 0~p-1,
                //考虑i+1位放置情况 
                ax=(1+b[i+1])*b[i+1]/2;
                siz=(p-b[i+1]-1)+1;
                ay=(b[i+1]+1+p)*siz/2+(p-siz)*p-o00;     
                ax%=mo;ay%=mo;
                bx=b[i+1]*(b[i+1]-1)/2;
                siz=(p-b[i+1])+1;
                by=(b[i+1]+p)*siz/2+(p-siz)*p-o10;
                bx%=mo;by%=mo;
                (f[i+1][j][0][0]+=f[i][j][0][0]*a00+f[i][j][0][1]*a01+f[i][j][1][0]*ax+f[i][j][1][1]*ay)%=mo;
                (f[i+1][j+1][0][1]+=f[i][j][0][0]*a10+f[i][j][0][1]*a11+f[i][j][1][0]*bx+f[i][j][1][1]*by)%=mo;
                (f[i+1][j][1][0]+=f[i][j][1][0]*(b[i+1]+1)+f[i][j][1][1]*(p-1-(b[i+1]+1)+1))%=mo;
                (f[i+1][j+1][1][1]+=f[i][j][1][0]*(b[i+1]-1+1)+f[i][j][1][1]*(p-1-b[i+1]+1))%=mo;
        }
    fo(j,K,n)
    {
        ans=(ans+f[n][j][0][0]+f[n][j][1][0])%mo;
    }
    printf("%lld\n",ans);
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值