2013年辽宁省赛 NEU1311 Rightmost non-zero Digit in N!

1311: Rightmost non-zero Digit in N!

时间限制: 2 Sec   内存限制: 128 MB
提交: 67   解决: 3

题目描述

The expression N!, read as "N factorial," denotes the product of the first N positive integers, where N is nonnegative. So, for example

N N!

0 1

1 1

2 2

3 6

4 24

5 120

10 3628800

For this problem, you are to write a program that can compute the rightmost non-zero digit of the factorial for N in base B.

输入要求

In the first line there is an integer T (T <= 100), indicates the number of test cases.

In each case, there is a single line with integers N and B

(0<=N<=10^7, 2<=B<=10^4).

输出要求

For each case, the output format is “Case #c: ans”. 

c is the case number start from 1.

ans is the answer of this problem. The ans is in base 10.

See the sample for more details.

假如输入

2 13 3 11 11

应当输出

Case #1: 2 Case #2: 10

提示

Case #1:

13*12*11*10*9*8*7*6*5*4*3*2*1=6227020800 base 10, which in base 3 is

121001222020102200000, so the right-most non-zero digit is 2.

 

警告：本人语文是体育老师教的。。水平有限。。可能导致你看了之后对世界失去信心。。故请慎看。。

省赛的数论题。。貌似封榜之前没有人交。。但比赛之后看有六支队过了这题。。。。（貌似我赛前偷窥的时候还看到admin没有过这题。。OTZ。。）

废话少说了。。来来说说这题的思路吧。。

有一个很直接的思路，就是将N的阶乘用唯一分解定理写开，找出决定0的个数的因子（对B进行因式分解，可以依次枚举。。方法是假设这个因子决定零的个数。。算出应该会有的零的个数。。取其中最小值。。）然后按比例除掉多余的因子。。剩下的因子用快速幂计算。。。我估计了一下，素数的个数差不多是N/lnN,故这个算法应该是一个N/lgN*lgB的算法。。。由于有大量乘法所以实际上的效率比想象中的要低。。所以交了之后果断TLE。。。

木有办法。。只好再想一个新的呗。。可是新的也不好想啊。。注意到威尔逊定理和高斯对威尔逊定理的概括。。我们可以设计一种在B是素数情况的lgN+B的算法

具体是这样的

由于B是素数，注意到k!  = (B+1)*(B+2)*(B+3)...*(B+k-1) mod B 且当k<B时我们有k! mod B 不等于 0

故在计算N！的最右非零位时可以这样算: 1*2*3...*(B-1)*(B+1)*...(B+B-1)*.......*((N/B)*B+N%B)*B*2B*...*(N/B)B再将所有的B去掉，利用上面的同余关系在O(B)时间内算出k! mod B的值lgN时间内可以将N！的最右非零位计算完。。。。

然后呢。。我就把抱着试一试的心态。。把素数的情况分开来单独做，再交了一次。。居然过了。。。orz。。。虽然程序跑的有点慢的说。。。。。。。。（貌似是一秒多来着。。。）

恩。。这就只能归结于数据水了。。可是我们不应该满足于此嘛。。。对于一般情况的高效做法肯定是存在的。。。。（为什么？感觉啊！）

嗯嗯。。。可是到底是怎么样的呢？又是一种比较直接的方法是仿造上面的做法，对于一般的B因式分解。。将所有B中的素因子拿出来，剩余部分可以仿造之前的算法。。但是拿出来的部分就不好做了。。。为什么呢？举个例子。。比如说要计算100！在12进制下的最右非零位。。我们将12分解成12 = 2 ^2 * 3这样就会有问题了。。我拿2的倍数的时候我把6拿了出来。。可是我再拿3的倍数的时候怎么办？若是不拿。。拿出的3的部分将3去掉后就不是一个连续的阶乘序列了。。要是拿了什么时候把多拿的数补回去？

所以捏。。我们需要再换个思路想想。。有这么一个结论叫做乘性函数的函数值由素数的函数值决定。。这就启示我们了对于某个进制B =  p1^a1 * p2^a2 * p3^a3... * pk*ak我们是不是可以把它看成是在

B1 = p1^a1 B2 = p2^a2 ...Bk = pk^ak进制下的解合并而成的呢？答案当然是可以啦。。。。

具体好麻烦哦。。不想敲了。。就是先算在Bi进制下的最右非零位将少乘的pi乘回来，没有去掉的其他因子乘上逆（注意到gcd(Bi, B/Bi) = 1）最后用中国剩余定理将结果合并就行啦！

其实呢这题肯定不用那么麻烦。。可是简单的办法由于我智商有限没有想到。。orz。。。要是有更好的办法欢迎交流！。。最后附上13msAC代码。。

#include <cstdio>
using namespace std;
typedef long long LL;
int T, N, B;
bool bp[10001];
int p[10001][20], num[10001];
LL fb[20], fn[20], nnum[20], cop[20], coprnz[20], mini;
void getp()
{
    for (int i=2; i<10001; i++)
    {
        if (!bp[i])for (int j=i; j<10001; j+=i)
        {
            p[j][num[j]++] = i;
            if(i!=j)bp[j] = 1;
        }
    }
}
LL mod_mult ( LL a, LL b, LL n )
{
    LL ret = 0;
    a = a % n;b = b % n;
    while ( b >= 1 )
    {
        if ( b & 1 )
        {
            ret += a;
            if ( ret >= n ) ret -= n;
        }
        a = a << 1;
        if ( a >= n ) a -= n;
        b = b >> 1;
    }
    return ret;
}
 
LL mod_exp ( LL a, LL b, LL n )
{
    LL ret = 1;
    a = a % n;
    if (b==1) return a;
    while ( b >= 1 )
    {
        if ( b & 1 )
           ret =  mod_mult(ret,a,n);
        a = mod_mult(a,a,n);
        b = b >> 1;
    }
    return ret;
}
void exgcd(LL a, LL b, LL &d, LL &x, LL &y)
{
    if (b==0)
    {
        d = a;x = 1;y = 0;
        return;
    }
    exgcd(b, a%b, d, y, x); y-=x*(a/b);
}
LL inv(LL a, LL n)
{
    LL x, d, y;
    exgcd(a, n, d, x, y);
    if (d!=1) return -1;
    return (x+n)%n;
}
void gcoprnz(int cur)
{
    int tB = cop[cur];
    int base[tB+1];base[0] = 1;
    for (int i=1; i<=tB; i++)
    {
        if (i%fb[cur])
        {
            base[i] = i * base[i-1] % tB;
        }
        else
        {
            base[i] = base[i-1];
        }
    }
    LL ans = 1;
    for (int t=N; t>0; t/=p[B][cur])
    {
        ans = ans * base[t%tB] % tB;
        ans = ans * mod_exp(base[tB], t/tB, tB)%tB;
    }
    ans = ans * mod_exp(fb[cur], nnum[cur]%fn[cur], tB) % tB;
    coprnz[cur] = ans;
    //printf ("%I64d\n", ans);
}
LL china()
{
    LL c[20];
    for (int i=0; i<num[B]; i++)
    {
        c[i] = inv(B/cop[i], cop[i]) * B/cop[i] % B;
        c[i] = c[i] * coprnz[i] % B;
    }
    LL res = 0;
    for (int i=0; i<num[B]; i++)
    {
        res += c[i];
    }
    return res%B;
}
int main()
{
    getp();
    scanf ("%d", &T);
    for (int ci=1; ci<=T; ci++)
    {
        scanf ("%d%d", &N, &B);
        int b = B;
        for (int i=0; i<num[B]; i++)
        {
            fb[i] = p[B][i];fn[i] = 0;cop[i] = 1;
            while (b%p[B][i]==0)
            {
                b /= p[B][i];
                fn[i]++;
                cop[i] *= p[B][i];
            }
            //printf ("p[B][i] = %d, time = %d\n", fb[i], fn[i]);
        }
        mini = 0;
        for (int i=0; i<num[B]; i++)
        {
            nnum[i] = 0;
            for (int j=N/p[B][i]; j>0; j/=p[B][i])
            {
                nnum[i] += j;
            }
            if (nnum[i]/fn[i]<nnum[mini]/fn[mini]) mini = i;
        }
        for (int i=0; i<num[B]; i++)
        {
            gcoprnz(i);
        }
        int u = nnum[mini]/fn[mini];
        for (int i=0; i<num[B]; i++)
        {
            //补数
            int t = (nnum[i]/fn[i]-u)*fn[i];
            coprnz[i] = coprnz[i] * mod_exp(fb[i], t, cop[i]) % cop[i];
            //乘逆
            int x = B/cop[i], v = inv(x, cop[i]);
            coprnz[i] = coprnz[i] * mod_exp(v, u, cop[i]) % cop[i];
        }
        LL res = china();
        printf ("Case #%d: ", ci);
        printf ("%d\n", (int)res);
    }
    return 0;
}
 
/**************************************************************
    Problem: 1311
    User: 20124867
    Language: C++
    Result: 正确
    Time:13 ms
    Memory:1616 kb
****************************************************************/