## 2016年1月25日 星期一

### [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;
{
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--)
{
}