BZOJ4737 组合数问题 Lucas定理 + 数位dp

Posted 2020-11-01 Mychael

题目

组合数C(n,m)表示的是从n个物品中选出m个物品的方案数。举个例子，从(1,2,3)三个物品中选择两个物品可以有(

1,2),(1,3),(2,3)这三种选择方法。根据组合数的定义，我们可以给出计算组合数C(n,m)的一般公式：
C(n,m)=n!/m!*(n?m)!
其中n!=1×2×?×n。（额外的，当n=0时，n!=1）
小葱想知道如果给定n,m和k，对于所有的0≤i≤n，0≤j≤min(i,m)有多少对(i,j)满足C(i,j)是k的倍数。

输入格式

第一行有两个整数t,k，其中t代表该测试点总共有多少组测试数据，k的意义见。
接下来t行每行两个整数n,m，其中n,m的意义见。

输出格式

t行，每行一个整数代表所有的0≤i≤n，0≤j≤min(i,m)中有多少对(i,j))满足C(i,j)是k的倍数

答案对10^9+7取模。

输入样例

3 23

23333333 23333333

233333333 233333333

2333333333 2333333333

输出样例

851883128

959557926

680723120

提示

1≤n,m≤10^18，1≤t,k≤100，且 k 是一个质数

题解

根据\(Lucas\)定理我们知道，在模质数\(k\)下
\[{n \choose m} \equiv \prod_{i = 1} {n \mod k^i \choose m \mod k^i} \pmod k\]
由此，结果为\(0\)，当且仅当存在一个\(i\)，使得\({n \mod k^i \choose m \mod k^i} \equiv 0 \pmod k\)
一个组合数在模质数意义下为0，当且仅当\(n < m\)
一个组合数在模质数意义下不为0，那么就是\(n >= m\)

那么我们将\(n\)和\(m\)拆分为\(k\)进制数，就可以设一个dp:
\(f[i][0/1][0/1]\)表示第\(i\)位之前(\(n\)是否达到上界) (\(m\)是否达到上界) 的\(n >= m\)方案数，即不为\(0\)的方案数
最后用总方案减去就可以了
转移就自己推推

要注意乘法可能溢出，要用快速乘

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define LL long long int
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<‘ ‘; puts("");
#define cls(s) memset(s,0,sizeof(s))
using namespace std;
const int maxn = 105,maxm = 100005,INF = 1000000000,P = 1e9 + 7;
LL n,m,f[maxn][2][2],a[maxn],b[maxn],ai,bi,N,K,v2;
LL qpow(LL a,LL b){
    LL ans = 1;
    for (; b; b >>= 1,a = a * a % P)
        if (b & 1) ans = ans * a % P;
    return ans;
}
void add(LL& a,LL b){
    a += b;
    if (a >= P) a -= P;
}
LL Mul(LL a,LL b){
    LL re = 0;
    for (; b; b >>= 1,a = (a + a) % P) if (b & 1) re = (re + a) % P;
    return re;
}
LL S(LL x){return Mul(x,x + 1) * v2 % P;}
int main(){
    int T; scanf("%d%lld",&T,&K); v2 = qpow(2,P - 2);
    while (T--){
        cls(a); cls(b); cls(f);
        scanf("%lld%lld",&n,&m); ai = bi = 0;
        LL ans = S(min(n,m) + 1);
        ans = (ans + Mul(m + 1,max(n - m,0ll))) % P;
        while (n) a[++ai] = n % K,n /= K;
        while (m) b[++bi] = m % K,m /= K;
        N = max(ai,bi);
        f[N][1][1] = 1;
        for (int i = N; i; i--){
            //0 0
            add(f[i - 1][0][0],S(K) * f[i][0][0] % P);
            //0 1
            add(f[i - 1][0][0] ,(S(b[i]) + b[i] * (K - b[i]) % P) % P * f[i][0][1] % P);
            add(f[i - 1][0][1],(K - b[i]) * f[i][0][1] % P);
            //1 0
            add(f[i - 1][0][0],S(a[i]) * f[i][1][0] % P);
            add(f[i - 1][1][0],(a[i] + 1) * f[i][1][0] % P);
            //1 1
            if (a[i] >= b[i]){
                add(f[i - 1][0][0],(S(b[i]) + b[i] * ((a[i] - b[i] + P) % P) % P) % P * f[i][1][1] % P);
                add(f[i - 1][1][0],b[i] * f[i][1][1] % P);
                add(f[i - 1][0][1],(a[i] - b[i] + P) % P * f[i][1][1] % P);
                add(f[i - 1][1][1],f[i][1][1]);
            }
            else {
                add(f[i - 1][0][0],S(a[i]) * f[i][1][1] % P);
                add(f[i - 1][1][0],(a[i] + 1) % P * f[i][1][1] % P);
            }
        }
        for (int i = 0; i < 2; i++)
            for (int j = 0; j < 2; j++)
                ans = (ans - f[0][i][j] + P) % P;
        printf("%lld\n",(ans % P + P) % P);
    }
    return 0;
}