[UVa 12581]Divisibility

Let's try to figure out what's the value of $(x_{1},x_{2},x_{3},...)$.
We can discover that the value is the number of different paths from $(0,0,0,...)$ to $(x_{1},x_{2},x_{3},...)$, one step in the path is defined as increase one of the co-ordinates by one.

So we can know the value of $(x_{1},x_{2},x_{3},...)$ is $\frac{(\sum_{k=1}^{N}  x_{k})!}{\prod_{k=1}^{N}  (x_{k}!)}$
If you don't know why, image that there are $x_{1}$ $vectors$ of $(1,0,0,...)$, $x_{2}$ $vectors$ of $(0,1,0,...)$, and so on. Therefore, no matter how you permutate them, they always form a path from $(0,0,0,...)$ to $(x_{1},x_{2},x_{3},...)$, the number of such permutations is $\frac{(\sum_{k=1}^{N}  x_{k})!}{\prod_{k=1}^{N}  (x_{k}!)}$ (Let's call it $M$ later)

The problem is, in which condition $M$ is not divisible by P ?

We can just take only $P$'s multiples into consideration.
The number of factor $P$ in $n!$ is $\left[\frac{n}{P}\right]+\left[\frac{n}{P^{2}}\right]+\left[\frac{n}{P^{3}}\right]+...$
For simplicity, we call the value $f(n)$ (regard $P$ as constant)

Then M is not divisible by P if and only if $f(\sum_{k=1}^{N}  x_{k})=\sum_{k=1}^{N}  f(x_{k})$
i.e., for every $r\in\mathbb{N}$, $\sum_{k=1}^{N}  (x_{k}\mod P^{r})<P^{r}$ must hold.

If we transform every $n$ into $P-dicimal$ numbers, call it $T(n)$, then the above equation holds if and only if $\sum_{k=1}^{N}  (k_{th}$ digit of $T(x_{k}))<P$.

So be question becomes, how many set, $(x_{1},x_{2},x_{3},...)$, are there, such that for each $i\in\mathbb{N}$, $\sum_{k=1}^{N}  (i_{th}$ digit of $T(x_{k}))<P$ ?

Let's try to design a Dynamic Programming method to get the answer.
Suppose that we can enumerate from high digit to low digit.
A state includes the following information:
1. $i$, represent we are considering the $i_{th}$ digit now.
2. For every $T(x_{k})$, the range of its $(i-1)_{th}$ digit should be.(whether the $(i-1)_{th}$ digit has low bound or up bound)

For example, if the $i_{th}$ digit has a low bound of $a$, a up bound of $b$, the $(i-1)_{th}$ digit has no bounds if the $i_{th}$ digit$\in[a+1,b-1]$, has a low bound if the $i_{th}$ digit is $a$, has a up bound if the $i_{th}$ digit is $b$
Thus, we can update the values by enumerate the subsets of its set of bounds.
It's convenient to use bitmask to represent the sets.
Partial time complexity: $O(3^{2N}\log_{P}\max(x_{k}))$
However, we can optimize it by skipping if the value is 0.

Use another DP to calculate the number of cases, such that $\sum_{k=1}^{N} \{ i_{th}$ digit of $x_{k}$, within the given bound$\}<P$.
Partial time complexity: $O(NP)$.

Total time complexity: $O(3^{2N}NP\log_{P}\max(x_{k}))$.

code:

#include<cstdio>
#include<cassert>
#include<algorithm>
using namespace std;
typedef long long LL;
const LL MOD=1000000009LL;
void getmax(int &a,const int b){if(b>a)a=b;}
int N,P;
vector<vector<LL> >A,B;
void ReadInput()
{
    scanf("%d%d",&N,&P);
    LL *a=new LL[N],*b=new LL[N];
    for(int i=0;i<N;i++)scanf("%lld",&a[i]);
    for(int i=0;i<N;i++)scanf("%lld",&b[i]);
    for(int i=0;i<N;i++)if(a[i]>b[i])swap(a[i],b[i]);
    A.clear(),B.clear();
    for(int i=0;i<N;i++)
    {
        vector<LL>va,vb;
        while(a[i])va.push_back(a[i]%P),a[i]/=P;
        while(b[i])vb.push_back(b[i]%P),b[i]/=P;
        A.push_back(va),B.push_back(vb);
    }
    delete[]a;delete[]b;
}
LL Solve(const vector<LL>&l,const vector<LL>&r)//O(7*20)
{
    assert((int)l.size()==N&&l.size()==r.size());
    static LL dp[2][19];
    for(int i=0;i<P;i++)dp[0][i]=0LL;
    dp[0][0]=1LL;
    int d=0;
    for(int i=0;i<N;i++,d^=1)
    {
        for(int j=0;j<P;j++)dp[d^1][j]=0LL;
        for(int j=0;j<P;j++)
        {
            for(int k=l[i];k<=r[i]&&j+k<P;k++)(dp[d^1][j+k]+=dp[d][j])%=MOD;
        }
    }
    LL ans=0LL;
    for(int i=0;i<P;i++)(ans+=dp[d][i])%=MOD;
    return ans;
}
LL Get(const vector<vector<LL> >&s,const int i,const int j){return j<(int)s[i].size()?s[i][j]:0LL;}
LL DP[51][1<<14];
bool Get(const int v,const int i){return (v&(1<<i))>0;}
LL Solve(const vector<vector<LL> >&A,const vector<vector<LL> >&B,const int M)//O(50*(2187*(21+140)+(2187*7)))
{
    assert(M<=50);
    for(int u=0;u<(1<<(N*2));u++)DP[M][u]=0LL;
    DP[M][0]=1LL;//all incomplete
    for(int i=M-1;i>=0;i--)
    {
        for(int u=0;u<(1<<(2*N));u++)
        {
            DP[i][u]=0LL;
            for(int a=u;;a=(a-1)&u)//complete can come from incomplete
            {
                if(DP[i+1][a])
                {
                    vector<LL>ls,rs;
                    for(int j=0;j<N;j++)
                    {
                        LL l,r;
                        if(!Get(a,j*2)&&!Get(a,j*2+1))
                        {
                            const int v1=Get(A,j,i),v2=Get(B,j,i);
                            if((v1==v2)!=(!Get(u,j*2)&&!Get(u,j*2+1)))
                            {
                                l=0,r=-1;
                            }
                            else if(!Get(u,j*2)&&!Get(u,j*2+1))l=r=v1;
                            else if(!Get(u,j*2))l=r=v1;
                            else if(!Get(u,j*2+1))l=r=v2;
                            else l=v1+1,r=v2-1;
                        }
                        else if(!Get(a,j*2))
                        {
                            if(Get(u,j*2))l=Get(A,j,i)+1,r=P-1;
                            else l=r=Get(A,j,i);
                        }
                        else if(!Get(a,j*2+1))
                        {
                            if(Get(u,j*2+1))l=0,r=Get(B,j,i)-1;
                            else l=r=Get(B,j,i);
                        }
                        else l=0,r=P-1;
                        ls.push_back(l),rs.push_back(r);
                    }
                    (DP[i][u]+=DP[i+1][a]*Solve(ls,rs))%=MOD;
                }
                if(a==0)break;
            }
        }
    }
    LL ans=0LL;
    for(int u=0;u<(1<<(N*2));u++)(ans+=DP[0][u])%=MOD;
    return ans;
}
LL Solve()
{
    int m=0;
    for(int i=0;i<N;i++)getmax(m,A[i].size()),getmax(m,B[i].size());
    return Solve(A,B,m);
}
int main()
{
//    freopen("in.txt","r",stdin);
    int testcount;scanf("%d",&testcount);
    while(testcount--)
    {
        ReadInput();
        static int kase=1;
        printf("Case %d: %lld\n",kase++,Solve());
    }
    return 0;
}