Description
組合數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的倍數。

題解：

相信大家都知道lucas定理的形式：$C_n^m=C_{n/p}^{m/p}*C_{n\%p}^{m\%p}\%p$，實際上，我們還可以這樣理解：把$n,m$看做一個$p$進制的數，那么結果等於$C_{a_0}^{b_0}\times C_{a_1}^{b_1}\times ……C_{a_k}^{b_k}$，$a_i$為$n$分解成$p$進制后的第$i$位，$b$也類似。然后若結果不是$p$的倍數，當且僅當所有$a_i>=b_i$，因為這樣的話所有$a_i、b_i$都是小於$p$的，不會出現$p$的因子；否則就會有幾項為$0$，然后根據這個數位DP即可，轉移有點惡心，要細心。

代碼：

#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define pa pair<int,int>
const int Maxn=110;
const LL mod=1000000007LL,inv=500000004LL;
LL read()
{
    LL x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
    return x*f;
}
LL n,m,a[70],b[70],la,lb;
LL f[70][2][2];
//位數 是否小於n 是否小於m
LL sum(LL l,LL r){l%=mod;r%=mod;return (l+r)%mod*((r-l+1LL)%mod)%mod*inv%mod;}   
int main()
{
    int T=read();
    LL k=read();
    while(T--)
    {
        memset(a,0,sizeof(a));
        memset(b,0,sizeof(b));
        memset(f,0,sizeof(f));
        n=read(),m=read();la=lb=0;
        if(m>n)m=n;
        LL t,ans;
        if(n==m)ans=sum(1,n+1);
        else ans=(sum(1,m+1)+(n-m)%mod*((m+1)%mod)%mod)%mod;
        t=n;
        while(t)a[++la]=t%k,t/=k;
        for(int i=1;i<=(la>>1);i++)swap(a[i],a[la-i+1]);
        t=m;
        while(t)b[++lb]=t%k,t/=k;
        for(int i=1;i<=(lb>>1);i++)swap(b[i],b[lb-i+1]);
        int len=max(la,lb);
        if(la!=len)for(int i=la;i;i--)a[i+len-la]=a[i],a[i]=0;
        if(lb!=len)for(int i=lb;i;i--)b[i+len-lb]=b[i],b[i]=0;
        f[0][0][0]=1;
        for(int i=0;i<len;i++)
        for(int pn=0;pn<2;pn++)
        for(int pm=0;pm<2;pm++)
        if(f[i][pn][pm])
        {
            LL t=f[i][pn][pm];
            if(pn&&pm)f[i+1][1][1]=(f[i+1][1][1]+t*sum(1,k)%mod)%mod;
            if(pn&&!pm)
            {
                f[i+1][1][1]=(f[i+1][1][1]+t*(sum(1,b[i+1])+(k-b[i+1])*b[i+1]%mod)%mod)%mod;
                f[i+1][1][0]=(f[i+1][1][0]+t*(k-b[i+1])%mod)%mod;
            }
            if(!pn&&pm)
            {
                f[i+1][1][1]=(f[i+1][1][1]+t*sum(1,a[i+1])%mod)%mod;
                f[i+1][0][1]=(f[i+1][0][1]+t*(a[i+1]+1)%mod)%mod;
            }
            if(!pn&&!pm)
            {
                if(a[i+1]&&b[i+1])
                {
                    if(a[i+1]<=b[i+1])f[i+1][1][1]=(f[i+1][1][1]+t*sum(1,a[i+1])%mod)%mod;
                    else f[i+1][1][1]=(f[i+1][1][1]+t*(sum(1,b[i+1])+(a[i+1]-b[i+1])*b[i+1]%mod)%mod)%mod;
                }
                if(b[i+1])f[i+1][0][1]=(f[i+1][0][1]+t*min(a[i+1]+1,b[i+1])%mod)%mod;
                if(a[i+1]&&a[i+1]-1>=b[i+1])f[i+1][1][0]=(f[i+1][1][0]+t*(a[i+1]-b[i+1])%mod)%mod;
                if(a[i+1]>=b[i+1])f[i+1][0][0]=(f[i+1][0][0]+t)%mod;
            }
        }
        for(int i=0;i<2;i++)
        for(int j=0;j<2;j++)
        ans=((ans-f[len][i][j])%mod+mod)%mod;
        printf("%lld\n",ans);
    }
}