Topcoder SRM 569 1000pt

Problem Statement

The factorial of the k-th order of a number n is denoted n!k and defined by the following recurrences:
1) n!k = n!(k-1) * (n-1)!k for n > 0 and k > 0
2) n!k = 1 for n = 0
3) n!k = n for k = 0

For example, 7!1 = 7! (the traditional factorial), and 5!3 = 5!2 * 4!3 = (5!1 * 4!2) * 4!3.

You are given s N, K and B. Count the number of trailing zeros in the number N!K when it is written in base B and return it modulo 1,000,000,009.

Definition

Class: MegaFactorial
Method: countTrailingZeros
Parameters: int, int, int
Returns: int
Method signature: int countTrailingZeros(int N, int K, int B)
(be sure your method is public)

Limits

Time limit (s): 840.000
Memory limit (MB): 64
Constraints
- N will be between 1 and 1,000,000,000 (10^9), inclusive.
- K will be between 1 and 16, inclusive.
- B will be between 2 and 10, inclusive.

Examples

0) 6 1 4
Returns: 2
6!1 = 6! = 23100 in base 4.
1) 4 2 6
Returns: 2
4!2 = 4!1 * 3!2 = … = 4! * 3! * 2! = 1200 in base 6.
2) 10 3 10
Returns: 22
3) 50 10 8
Returns: 806813906
4) 1000000000 16 2
Returns: 633700413

This problem statement is the exclusive and proprietary property of TopCoder, Inc. Any unauthorized use or reproduction of this information without the prior written consent of TopCoder, Inc. is strictly prohibited. (c)2003, TopCoder, Inc. All rights reserved.

我A的Topcoder第一道题啊，纪念一下。

题目非常简洁啊（似乎某人说过“题目越简洁，就越是暗藏杀机”），题目中定义了一种函数f(i,j)，我们不妨类比组合数的公式，考虑把f(i,j)不断地展开，那么f(i,j)=i∗∏jk=1f(i−1,k)，如果我们把f(n,k)唯一分解，那么它在B进制下末尾0的个数就等于几个相乘=B的素数的最小的指数，比如B=10的时候，答案就取决于f(n,k)中min(2的指数,5的指数)（应该算是比较好理解吧，相当于就是看可以凑出来多少个B就可以了），考虑每个质数出现次数，便易知最小的指数是在B的唯一分解下最大的质数（称其为答案质数P）得到的，这样的话我们就可以只考虑某一个质数的指数了（你可能会怀疑有问题，比如8=23这种P的指数不为1的情况怎么办？我们一会再说），同时考虑到指数的性质（乘法=对应指数相加），乘法成功转变成加法！
解法一：
不妨单独考虑第一行中1、2、3...每一个对某一个f(n,k)的贡献，由于是指数相加，我们可以得到转移方程g(i,j)=g(i−1,j)+g(i,j−1)，这个转移方程的闭合形式显然是一个组合数！而分析一下可以得到i的指数就是C(k−1+n−i,n−i)，也就是(1,i)走到(n,k)的方案数。所以f(n,k)=∏ni=1iC(k−1+n−i,n−i)，这样的话最终答案就是几个k项组合数相加的形式，（经过漫长的划式子过程，基本上就是把累乘多项式展开，然后再统一处理系数的问题）我们发现这居然是一个伯努利数的形式，直接求和。（ORZ TA学长，ORZ reflash）
解法二：
我们来考虑暴力一点的做法，观察到n非常大，考虑把n给矩乘，但是由于转移的时候还需要加上当前行号i中P的指数，似乎没法矩乘啊。我们把1...n中p的指数都写出来，便会发现一个非常优美的性质：1...pk−1和pk+1...2∗pk−1的指数是对应相等的！矩乘的时候我们是把行号中指数对应相加，而这种重复出现无疑为我们提供了便利，相当于是同一个规则重复出现了多次，考虑倍增！设矩阵A[i]表示当行号中P的指数为i时的转移矩阵，差不多也就是

A[i]=⎛⎝⎜⎜⎜⎜⎜⎜100...i110i111i............0001⎞⎠⎟⎟⎟⎟⎟⎟

的形式，然后矩阵B[i]等于∏Pi−1j=1A[j]，重复出现那我们就要统一处理，最终B[i]=(B[i−1]∗A[i−1])P−1∗B[i−1]（请读者自己画一个图或者打一张表感受一下），最后再处理一下B[i]=B[i]∗A[i]，这样我们把n按照P进制展开，把B矩阵不断累乘转移就好了（转移的时候一定要注意矩阵满足结合律但不满足交换律）！
但是注意还没有结束！我们得到了P的指数（设为X）是没错，但这就是答案吗？不！比如8=23，我们得到了f(n,k)中2的指数，但是三个2次才能凑出一个8，所以我们要把答案进行整除3，可这是在模意义下啊！考虑转化为取模，即（假设P在B中的指数为s）：ans=Xs(mod 109+9)->X−s∗(Xs)=(X mod s) (mod 109+9)->ans=X−(X mod s)s(mod 109+9)。所以计算两遍分别得出X mod 109+9和X mod s，然后计算就好了。
（TC为了时间效率，非常人性化地用了我根本不懂的class和public，贴一下代码防止以后忘了格式）

#include<cstdio>#include<cstring>#include<cmath>#include<algorithm>using namespace std;class MegaFactorial{    typedef long long ll;    int n,m,K,base,aa[50],BB,top;ll mod;    struct data{        int x,y,a[20][20];        data(){x=y=0;memset(a,0,sizeof(a));}        void data1(){x=y=18;memset(a,0,sizeof(a));for (int i=1;i<=18;++i)a[i][i]=1;}    }A[50],B[50],ans;    public:int cheng(int a,int b)    {        ll x=a,y=b,z=(x*y)%mod;        return (int)z;    }    public:ll power(ll a,ll b)    {        ll ans=1;        for (;b;b>>=1,a=(a*a)%mod)        if(b&1)ans=(ans*a)%mod;        return ans;    }    public:data ch(data a,data b)    {        data c;c.x=a.x,c.y=b.y;        for (int i=1;i<=a.x;++i)            for (int j=1;j<=b.y;++j)                for (int k=1;k<=a.y;++k)                c.a[i][j]=(c.a[i][j]+cheng(a.a[i][k],b.a[k][j]))%mod;        return c;    }    public:data calc(data a,int b,int K)    {        data ans,aa=a;        ans.data1();ans.x=ans.y=K+1;        for (;b;b>>=1,aa=ch(aa,aa))        if(b&1)ans=ch(ans,aa);        return ans;    }    public:ll getans(int K)    {        for (int now=0;now<=m;++now)        {            A[now].data1();A[now].x=K+1,A[now].y=K+1;            for (int i=1;i<=K;++i)A[now].a[K+1][i]=now%mod;            for (int i=1;i<=K;++i)                for (int j=1;j<=i;++j)                A[now].a[j][i]=1;        }        B[0].data1();B[0].x=B[0].y=K+1;        B[1]=calc(A[0],BB-1,K);        for (int i=2;i<=m;++i)        {            B[i].data1();B[i].x=B[i].y=K+1;            for (int j=1;j<=BB-1;++j)            {                B[i]=ch(B[i],B[i-1]);                B[i]=ch(B[i],A[i-1]);            }            B[i]=ch(B[i],B[i-1]);        }        for (int i=0;i<=m;++i)B[i]=ch(B[i],A[i]);        ans.x=1,ans.y=K+1;memset(ans.a,0,sizeof(ans.a));ans.a[1][K+1]=1;        for (int i=top;i>=1;--i)ans=ch(ans,calc(B[i-1],aa[i],K));        return ans.a[1][K];    }    public:int countTrailingZeros(int N,int K,int B)    {        int hh[]={0,0,2,3,2,5,3,7,2,3,5};        int hhh[]={0,0,1,1,2,1,1,1,3,2,1};        BB=hh[B];top=0;while(N){aa[++top]=N%BB;N/=BB;}m=top;ll fans;        if(hhh[B]==1){mod=1000000009LL;fans=getans(K);return (int)fans;}        else        {            mod=hhh[B];fans=getans(K);            mod=1000000009LL;fans=(getans(K)-fans+mod)%mod;            fans=fans*power(hhh[B],mod-2)%mod;            return (int)fans;        }    }};

总结：
1、矩阵满足结合律不满足交换律
2、多次处理不要忘记数组的清空
3、观察性质的能力非常重要，而这只能通过不断的锻炼才能逐渐加强
4、当某一规律、某一规则重复出现的时候要考虑倍增