高位数的阶乘计算(简单的matlab，Python以及C语言实现)

今天有同学问我1000万的阶乘用计算机怎么算，所以我花了大概20分钟想了一个思路，分享给大家，主要面向新手，大佬轻喷。
大家都知道，计算机计算的数字精度是有限制的，就整数而言对于不同的CPU，不同的操作系统，不同的语言，不同的编译器等都可能导致计算位数精度的不同。对于一般的编程语言，在常用的平台上处理精度至多不会超过2的63。话不多说，我们来看算法本身：
先说思路，逻辑就是小学教的乘法计算方法。具体操作方法是，第一步把要计算的数的每一位都提出来然后依次放入数组，然后把该数组直接乘以下一个数，对得到的数组进行进位运算，再乘下一个。

function b=jiechen(a)
%将要计算的数的每一位数字依次放入一个一维数组
a_1=a;
len=floor(log10(a))+1;
for i=1:len
    b(i)=mod(a_1,10);
    a_1=floor(a_1/10);
end
%开始计算
for j=a-1:-1:1
    b=b*j;
    len=size(b,2);
    %进位运算
    for i=1:len-1
        if b(i)>=10
            b(i+1)=b(i+1)+floor(b(i)/10);
            b(i)=mod(b(i),10);
        end
    end
    while b(end)>=10
        c=b(end);
        b(end)=mod(b(end),10);
        b=[b,floor(c/10)];
    end
end
%逆序输出
b=b(end:-1:1);

100的阶乘计算结果如下：

使用matlab自带的factorial，prod，gamma函数计算的高位数阶乘并不准确。如果为了简单验证可以使用Windows计算器

按照这个思路，你也可以轻松将matlab代码转换成任何你想要的编程语言代码。

用Python也非常简单，下面给出代码：

import math
import numpy as np

a=10000
a_1=a
b=np.array([a_1%10])
len=math.floor(np.log10(a_1))+1
for i in np.arange(1,len):
    a_1=math.floor(a_1/10)
    b=np.append(b,a_1%10)
print(b)
for j in range(a-1,1,-1):
    b=b*j
    len=np.size(b)
    #进位运算
    for i in range(0,len-1):
        if b[i]>=10:
            b[i+1]=b[i+1]+int(b[i]/10)
            b[i]=b[i]%10

    while b[-1]>=10:
        c=b[-1]
        b[-1]=b[-1]%10
        b=np.append(b,int(c/10))
    
#逆序输出
print(list(reversed(b)),np.size(b))

当然实际上用Python不必这么麻烦，Python可以精确处理长整型的各种运算，而且计算速度非常快。
例如：

>>> math.factorial(100)
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

使用C语言也很简单

#include<stdio.h>
#include<math.h>
#define N 10000

int main()
{
    float a, a_1,len,c ;
    int i, j;
    int b[N]= { 0 };
    printf("请输入要计算阶乘的数字：");
    scanf_s("%f", &a);
    a_1 = a;
    len = floor(log10(a)) + 1;
    for (i = 0;i<len;i++)
    {
        b[i] = (int)a_1%10;
        printf("%d\n", b[i]);
        a_1 = floor(a_1 / 10);
    }
    /*开始计算*/
    for (j = a - 1; j > 0; j--)
    {
        i = N;
        while (b[i - 1] == 0)
        {
            i--;
            if (b[i] != 0)break;
        }
        len = i;
        for (i = 0; i < len; i++)
        {
            b[i] = b[i] * j;
        }
       /*进位运算*/
        for (i = 0; i < N-1 ; i++)
        {
            if (b[i] >= 10)
            {
                b[i + 1] = b[i + 1] + floor(b[i] / 10);
                b[i] = b[i]%10;
            }
        }
    }
    i = N;
    while (b[i - 1] == 0)
    {
        i--;
        if (b[i] != 0)break;
    }
    len = i;
    for (i=len-1;i>=0; i--)
    {
        printf("%d", b[i]);
    }    
}