求解n的阶乘的十进制表示的最右非零位上的数字[转]

求解n的阶乘的十进制表示的最右非零位上的数字[转]

 
下面是算法说明

如果用 C 表示 [n/5] + [n/25] + [n/125] + ..., 那么需要求的是下面的同余方程

         n ! ≡ x * 10^c (mod 10^(c+1))

的解 x , 0< x ≤ 9.
上面这个同余方程等价于下面的方程组

        n! ≡ x * 5^c * 2^c (mod 2^(c+1)), n! ≡ x * 2^c * 5^c (mod 5^(c+1))

当 n > 1 时所有的偶数都是上面的方程组中第一个方程的解,而且, n > 1 时该方程组的
第一个方程没有奇数解,因此 n > 1 时只需要考虑下面这个方程(即方程组的第二个方程)
的符合 0 < x <9 的偶数解:

         n! ≡ x * 2^c * 5^c (mod 5^(c+1)).                         (a)

用 h(n) 表示 所有与 5 互素且不大于 n 的正整数的连乘积, 则 n! 可以表为

h(n) * 5^[n/5] * h([n/5]) * 5^[n/25] * h([n/25]) * 5^[n/125] * h([n/125]) *... ,

代入 (a) 式,消去 5 的乘方后得到下面的同余方程

         h(n) * h([n/5]) * h([n/25]) * ... ≡ x * 2^c (mod 5), (b)

由于 3 * 2 ≡ 1 (mod5), 因此 (b) 式变为
  
         3^c * h(n) * h([n/5]) * h([n/25]) * ...≡ x (mod5),         (c)

由 Euler-Fermat 公式知( % 表示求模运算 )

         3^c ≡ 3 ^ (c % 4) (mod 5),

由 Wilson 定理有

         h(n) ≡ (-1)^([n/5]) * (( n % 5 )! ) (mod 5),

把上面两式代入 (c) 就得到了

         3^(c%4) * (-1)^c * ( n % 5 )! * ([n/5] % 5)! * ([n/25] % 5)!*...

≡ x (mod 5)                                                  (d)

由 c 的表达式可知, 我们可以在求 [n/5],[n/25],[n/125],... 的过程中
求出 c 和 (n % 5)!,([n/5] % 5)!,([n/25] % 5)!,..., 这样就可以求得
(d) 式的左边.把最后的结果 模 5 后即可以求得 x.

这个过程实际上是用连除法求 n 的 5 进制表示. 由于 5 进制表示的 各个 位上的数字是任意的,因此 (d) 式的左边已不能再简化.由于任意 求进制表示的 方法本质上都是连除法,因此这个算法 本质上 已是最优算法..

这是一个时间复杂度 O(log n) 的算法

用 C 语言写的该算法的实现

为了确保跨平台性,我使用了
C 99 规定的类型 uint_fast32_t, 为使用这个类型,请 包含 c 语言的 标准头文件 stdint.h. 下面是我的函数

uint_fast32_t        GetRightestNoZeroDigitA(uint_fast32_t        n)
{
          uint_fast32_t        m = 1,    c = 0,   i;
          uint_fast32_t        cr[] = { 1, 3, 4, 2 };
          uint_fast32_t        rr[] = { 1, 6, 2, 8, 4 };
          while( n > 1){
              if ( 2 == (i = n % 5) )
                  m <<= 1;
              else if ( 4 == i )
                  m <<= 2;
              c += n /= 5;
          }
          if ( c   & 1 )
              return        rr[ ( cr[ c % 4 ] * 4 * m ) % 5 ];
          return        rr[ ( cr[ c % 4 ] * m ) % 5 ];   
}