C语言：阶乘的高精度计算
C语言：阶乘的高精度计算
在编程中，计算一个正整数的阶乘是一个常见的问题。然而，当这个正整数足够大时，普通的整型数据类型就无法储存如此大的结果。本文将介绍如何使用高精度计算算法来解决这一问题，重点讲解了高精度乘法的原理和实现方法，并给出了基于高精度乘法实现阶乘计算的具体代码。
前言
计算一个正整数的阶乘在思路上并不难，难的是当这个正整数足够大时，编程语言的整型数据类型便无法储存如此大的结果。这篇文章将介绍如何用高精度计算算法解决这一问题。
高精度阶乘问题
题目描述
参考题目连接：PTA | 程序设计类实验辅助教学平台 (pintia.cn)
题目要求很简单：输入一个整数，输出它的阶乘。
低精度算法及局限性
按比较直接的思路，用一个for循环即可解决，代码如下
void Print_Factorial ( const int N ){
int s = 1;
for(int i=1;i<=N;i++){
s *= i;
}
printf("%d",s);
}
或者是用一个简单的递归函数：
int Factorial(int N){
if(N == 1) return 1;
return N*Factorial(N-1);
}
void Print_Factorial ( const int N ){
printf("%d",Factorial(N));
}
以上都是一些有一定可行性的做法。但是题目对N的要求为“其值不超过1000”，也就是说，N的值能够做到足够大使得结果无法被C中所有整型数据结构储存。如
100！=93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
这个时候，以上方法便不适用了。
高精度算法
原理
高精度算法，顾名思义，便是比简单数字运算精度更高的一种算法。这种精度就体现在可以进行大数值的运算。
首先要解决的问题便是：怎么才能够储存一个巨大的数值呢？答案是用一个数组储存每一位上的数字。这样一来，我们定义一个长度100一维数组，就可以存放一个位数小于等于100的整数，这已经远远超出了long类型所能存放的大小。
解决了数据存储的问题，另一个问题便随之而来：把数据存储在数组里，该怎么实现两个变量之间的运算呢？其实，不论是高精度加法，减法，还是乘法，其计算模式都与一种小学生都会的方法十分相似——列竖式。
下面以高精度乘法为例。
高精度乘法
对两个数进行高精度乘法的过程就跟对他们进行列竖式计算差不多：
列竖式时，我们将其中一个数A作为“被乘数”，另一个数B作为“乘数”，先用B的第1位（个位）与A的各个位相乘，得到一个结果n1；再用B的第2位（十位）与A的各个位相乘，得到第二个结果n2......重复这个过程直到B的最后一位，将各次结果乘上相印系数后相加即可。
而在C中，用数组实现这一过程更简单，省去了一些不必要的步骤。
下面讨论如何用程序实现。
首先便是要输入两个待计算的数字。由于待乘数可能会很大，故在输入时，直接先将它们存入字符形式数组，再转化位int类型数组（同时进行倒置）。过程函数如下：
int PutInNum(char a[]){
int i=0;
char c;
while(scanf("%c", &c) && c != '\n'){ //遇到换行（回车）时停止输入
a[i++] = c;
}
return i; //返回数字的长度
}
/*转换数组类型并倒置*/
void Transform(char n1[], int n2[], int l){ //此处int类型的数组n2应初始化为{0}
int p = l-1;
for(int i=0;i<l;i++){
n2[p--] = (int)n1[i] - 48; //‘1’的ascii为49，故减去48便得到数字1
}
}
此处进行倒置是为了后面计算方便。此处不展开，后文再细讲。
将数字存入数组后，便是计算部分。我们创建一个数组res来储存结果。res同样应初始化为{0}，为计算答案做准备。
计算思路就是以n1的每一位与n2的每一位做乘法后逐次将结果存入res数组。在数十种，我们知道，十位的数乘以百位的数，结果的最低位数便是10*100 = 1000即千位（2+3-1=4即第四位）。那么在进行高精度乘法时也是一样，n1的第2位n1[1]，与n2的第3位n2[2]相乘，所得位数就应该为2+3-1=4，即res[3]。我们惊奇的发现，其实最后存入res的位置，就是n1,n2对应位数的索引值之和（1+2=3）。
由此可以得出，n1[i]与n2[j]之积的初始位数在res[i+j].（这就是倒置的优势，而且还可以从0开始遍历）
而n1[i]与n2[j]之积可能是一个两位数，我们需要将它们的积的十位部分加入res[i+j+1]，即**res[i+j+1] = (n1[i]n2[j])/10，res[i+j] = (n1[i]n2[j])%10
这样计算后得到的数组，再从后往前遍历输出，便可得到最终答案。计算函数如下
int Print_MULT(int n1[], int n2[], int l1, int l2){
int res[MAXN] = {0}; //储存结果的数组
int len = 0;
/*n2做初始数，n1做乘上的数
即在乘法竖式中，n2在上,n1在下*/
for(int i=0;i<l1;i++){ //i指向n1当前运算位数
for(int j=0;j<l2;j++){ //j指向n2当前运算位数
int tmp = n1[i] * n2[j];
res[i+j] += tmp;
if(res[i+j] >=10){
res[i+j+1] += (int)res[i+j]/10;
res[i+j] %= 10;
len = (i+j+1 > len)?i+j+1:len;
}
else len = i+j;
}
}
for(int i=len; i>=0;i--) printf("%d",res[i]);
}
完整代码如下：
#include<stdio.h>
#define MAXN 10000
/*输入数字*/
int PutInNum(char a[]){
int i=0;
char c;
while(scanf("%c", &c) && c != '\n'){ //遇到换行（回车）时停止输入
a[i++] = c;
}
return i; //返回数字的长度
}
/*转换数组类型并倒置*/
void Transform(char n1[], int n2[], int l){ //此处int类型的数组n2应初始化为{0}
int p = l-1;
for(int i=0;i<l;i++){
n2[p--] = (int)n1[i] - 48; //‘1’的ascii为49，故减去48便得到数字1
}
}
int Print_MULT(int n1[], int n2[], int l1, int l2){
int res[MAXN] = {0}; //储存结果的数组
int len = 0;
/*n2做初始数，n1做乘上的数
即在乘法竖式中，n2在上,n1在下*/
for(int i=0;i<l1;i++){ //i指向n1当前运算位数
for(int j=0;j<l2;j++){ //j指向n2当前运算位数
int tmp = n1[i] * n2[j];
res[i+j] += tmp;
if(res[i+j] >=10){
res[i+j+1] += (int)res[i+j]/10;
res[i+j] %= 10;
len = (i+j+1 > len)?i+j+1:len;
}
else len = i+j;
}
}
for(int i=len; i>=0;i--) printf("%d",res[i]);
}
int main(){
char a[MAXN], b[MAXN];
int la = PutInNum(a); int lb = PutInNum(b);
int n1[MAXN] = {0}; int n2[MAXN] = {0};
Transform(a ,n1,la); Transform(b, n2, lb);
Print_MULT(n1,n2,la,lb);
return 0;
}
高精度算法实现阶乘计算
接下来进入正题：如何基于高精度乘法进行阶乘计算？
PTA | 程序设计类实验辅助教学平台 (pintia.cn)
这其实就是一个进行若干次高精度乘法的过程。求n的阶乘，只需先将n与n-1相乘，其结果再与n-2相乘...如此往复，便可得到结果。
计算n!的过程，仍然是用到n1,n2,res这三个数组。第一次相乘，n1为n，n2为n-1，结果存入res中。再将n1初始化为res，res重置为{0}，n2 = n2 - 1，再进行一次乘法，结果存入res中......直到n2=1时，res中的数便是结果。
实现代码如下：
#include <stdio.h>
void Print_Factorial(const int N);
int main()
{
int N;
scanf("%d", &N);
Print_Factorial(N);
return 0;
}
#define MAXN 10000
int Print_MULT(int n1[], int n2[], int l1, int l2, int res[], int len)
{
/*n2做初始数，n1做乘上的数
即在乘法竖式中，n2在上,n1在下*/
for (int k = 0; k <= l1; k++)
{ // i指向n1当前运算位数
for (int j = 0; j < l2; j++)
{ // j指向n2当前运算位数
int tmp = n1[k] * n2[j];
res[k + j] += tmp;
if (res[k + j] >= 10)
{
res[k + j + 1] += (int)res[k + j] / 10;
res[k + j] %= 10;
len = (k + j + 1 > len) ? k + j + 1 : len;
}
else
len = k + j;
}
}
return len;
// for(int i=len; i>=0;i--) printf("%d",res[i]);
}
void Print_Factorial(const int N)
{
if (N == 0)
printf("1");
else if (N < 0)
printf("Invalid input");
else
{
int n1[MAXN] = {0};
int time = 0;
int last_r = 0;
int p = 0;
for (int i = 10; time == 0; i *= 10)
{
if (N % i == N)
time++;
int r1 = N % i;
int r2 = r1 - last_r;
last_r = r1;
n1[p++] = r2 / (i / 10);
}
int flag = 0;
int n2[MAXN] = {0};
int len = p - 1;
for (int i = 0; i < p; i++)
{
n2[i] = n1[i];
}
for (int i = 0; i < N - 1; i++)
{
int res[MAXN] = {0};
n2[0] -= 1;
for (int j = 0; j < p; j++)
{
if (n2[j] < 0)
{
n2[j] = 9;
n2[j + 1] -= 1;
}
if (n2[p - 1] == 0)
p--;
}
len = Print_MULT(n1, n2, len, p, res, len);
for (int j = 0; j <= len; j++)
{
n1[j] = res[j];
}
}
for (int i = len; i >= 0; i--)
printf("%d", n1[i]);
}
}