CN104636113A - 一种计算机处理大整数的算法 - Google Patents

Abstract

本发明公开了一种计算机处理大整数的算法，包括以下步骤：S1：采用手工或从文件中输入数据到计算机中存储；S2：通过预处理模块对S1中输入的数据进行预处理；S3：对S2中预处理后得到的数据，通过计算模块进行计算；S4：对S3中计算得出的结果通过结果处理模块进行处理。本发明的有益效果：本发明采用动态分配的或者固定长度的字符串作为大整数的存储结构，可以进行4000位以上各类运算，解决了现有计算机进行大整数计算时不精确的问题。

Description

一种计算机处理大整数的算法

技术领域

本发明涉及大整数计算技术，具体涉及一种计算机处理大整数的算法。

背景技术

现有程序设计语言提供的整数数据类型一般最多为8字节(64位二进制整数)，可表示的非负整数范围为0～2⁶⁴-1，而2⁶⁴-1＝18446744073709551615只有20位十进制数字。所以，针对更大整数的处理和计算问题，必须设计专门的软件来实现有关功能。

Windows操作系统提供的计算器可以进行32位以内十进制整数的运算。但是当运算结果超过十进制32位时，只能用科学计数法表示，其小数部分有效数字最多为30位，不能保证得到正确结果。

美国Wolfram Research公司用C语言开发的数学系统软件Mathematica是一个功能强大的计算机数学系统软件。它不但可以解决数学中的常用的加、减、乘、除、乘方、开方、阶乘、Euclid除法等数值计算问题，而且可以自动地完成积分、多项式的因式分解等复杂的计算工作，还可以解决符号演算问题，方便地绘出各种函数图形。

美国PTC公司的MathCAD是类似于Mathematica的工程软件，但不如Mathematica功能强大。它的代数运算模块可以进行整数的基本四则运算和乘方、开方，也可以解四次以下的代数方程，但能够处理的整数范围不大。

另外还有一款Maple软件也类似于MathCAD和Mathematica。

这些软件都没有提供Stirling数、Catalan数等特殊数的计算能力，也不能求解不定方程和同余方程组，不能计算一个正整数关于另外一个正整数的乘法阶和乘法逆。

1994年美国Peter Shor提出了在世界上有重大影响的基于量子计算的大整数因子分解算法。Shor算法成功的关键正是计算一个正整数关于另外一个正整数的乘法阶。可见美国等西方国家也在研究大数计算问题，并且在电子计算条件下也没有大的突破。

国内未见提供不定方程求解、同余方程组求解和特殊数计算等功能且拥有自主知识产权的大整数计算技术，国外相应技术也未完全提供这些功能。

在国家军事、政治、外交、经济、商业等领域的信息安全问题中，在计算机科学、数学等领域的科学研究中，常常涉及或者必须考察“大”到数十位甚至上百位十进制整数的有关计算，有时候甚至要求计算精度达到小数点后数十位甚至上百位(多位小数的有关计算都可以化为整数的有关计算)。比如我国“两弹一星”理论设计和实际试验中有些数据要计算到小数点后数十位。又如在军事信息安全领域中，有的地方要用到RSA、ELGamal或ECDLP密码算法或数字签名算法，这些算法都离不开大整数的运算，甚至涉及有关大素数的整数运算，诸如大整数的加法、乘法、乘法阶、乘法逆、最大公约数、因子分解以及同余等问题，其它密码算法或数字签名算法也会涉及大大整数的运算。再如军事领域很多地方要用到伪随机数，而伪随机数一般是用伪随机序列产生或者用同余式计算得到的，这二者都离不开有关大整数的计算。现有的大整数计算技术不能完全满足这些领域的专门需要或者在一定程度上制约着这些领域有关课题的研究和发展。

故现有技术亟需发展和改进。

发明内容

本发明要解决的技术问题是提供一种能够基于字符串而处理大整数的计算机算法。

本发明的技术解决方案是：

一种计算机处理大整数的算法，包括以下步骤：

S1：采用手工或从文件中输入数据到计算机中存储；

S2：通过预处理模块对S1中输入的数据进行预处理；

S3：对S2中预处理后得到的数据，通过计算模块进行计算；

S4：对S3中计算得出的结果通过结果处理模块进行处理。

S1输入的数据采用动态分配的或者固定长度的字符串作为大整数的存储结构；输入的数据中每一位十进制数字直接用其ASCII码存储，负整数的负号“-”写在字符串首位，非负整数不带负号。

在面向对象技术编程时，可将S1中输入的数据单独设计为数据成员。

所述预处理模块是用来删除S1输入数据中所有的空格和删除输入数据前面多余的数字“0”以及判断输入的数据是不是整数。

所述计算模块包括基本运算模块、乘方与开方模块、阶乘与因子模块、素数模块、最大公约数模块、特殊数模块、不定方程模块、同余式模块。

所述计算模块在计算时不改变S1中输入的数据。

所述结果处理模块是用于对S3的计算结果数据的整数部分从后向前每隔一个固定的位数就插入一个空格，对计算结果数据的小数部分从前向后每隔一个固定的位数就插入一个空格。

所述S4得到的数据长度应大于S1输入数据长度的3倍。

本发明的有益效果：

本发明采用动态分配的或者固定长度的字符串作为大整数的存储结构，可以进行4000位以上十进制整数的加法、减法、乘法、乘方、开平方、Euclid除法、除法和整除性测试，除法运算可精确到小数点后1023位以上，其中开平方运算可精确到小数点后3000位以上，乘方运算可以计算出2的13604次方、3的8583次方、5的5859次方、5的5859次方、7的4846次方等等；可以计算1493以内十进制正整数的阶乘、4000位以上十进制正整数的阶乘因子和最大公约数、2000位以上十进制正整数的最小公倍数；可以计算参数在2506以内的二项式系数、参数在458以内的Stirling数和参数在1254以内的Catalan数；可以计算斐波那契数到19598项以内；可以求解系数在十进制2000位以上的二元一次、三元一次不定方程和系数在十进制16位以内的一元二次方程；可以计算十进制2000位以内正整数的乘法逆，求解系数在十进制2000位以内的一次同余方程和方程组；可以成功地对16位以内的十进制正整数进行素性判定和因子分解。一元二次方程、素性判定和因子分解等问题的计算复杂性高，速度一般较慢，因此这里对其范围限制得很小。

附图说明

图1为本发明步骤方框图；

图2是本发明主界面图；

图3是“基本运算”菜单；

图4是“乘方与开方”菜单；

图5是“阶乘与因子”菜单；

图6是“素数”菜单；

图7是“最大公约数”菜单；

图8是“特殊数”菜单；

图9是“不定方程”菜单；

图10是“同余式”菜单。

图11是一个111位正整数与一个123位正整数乘法的计算实例；

图12是图11中两个整数除法的计算实例，精确到小数点后1023位；

图13是上述两个整数Euclid除法的计算实例；

图14是2的13604次方的计算实例，计算结果有4096位；

图15是1234567890987654321的因子分解的计算实例，分解结果为1234567890987654321＝32×7×19×928163×1111211111其中928163和1111211111都是素数；

图16是30位数123456789098765432101234567890的平方根的计算实例，精确到小数点后3000位；

图17是1493！的计算实例，计算结果有4093位，末尾有370个0；

图18是14位数42949672971917的最小素因子的计算实例，计算结果表明42949672971917是一个素数；

图19是计算图11中两个整数最大公约数的计算实例，计算结果表明那两个数是互素的；

图20是计算图11中两个整数最小公倍数的计算实例，可以用乘法验证这个结果；

图21是二项式系数的计算实例，计算结果有753位；

图22是第19598个斐波那契数的计算实例，计算结果有4096位；

图23是第一类Stirling数s(458，234)的计算实例，计算结果有656位；

图24是第二类Stirling数S(458，234)的计算实例，计算结果有612位；

图25是Catalan数h(1254)的计算实例，计算结果有749位；

图26是求解以图11中两个整数为系数、常数项为1的二元一次不定方程的计算实例；

图27是求解二元一次不定方程

42949672971917x+1234567890987654321y＝123456789的计算实例；

图28是求解三元一次不定方程

42949672971917x+1234567890987654321y+1111211111z＝123456789的计算实例；

图29是求解一元二次方程

1234567890x2+5432112345x+9876543210＝0的计算实例；

图30是乘法逆的计算实例；

图31是一元一次同余方程的计算实例；

图32是一元一次同余方程组的计算实例。

具体实施方式

实施例：

一种计算机处理大整数的算法，其特征在于，包括以下步骤：

S1：采用手工或从文件中输入数据到计算机中存储；

S2：通过预处理模块对S1中输入的数据进行预处理；

S3：对S2中预处理后得到的数据，通过计算模块进行计算；

S4：对S3中计算得出的结果通过结果处理模块进行处理。

S1输入的数据采用动态分配的或者固定长度的字符串作为大整数的存储结构；输入的数据中每一位十进制数字直接用其ASCII码存储，负整数的负号“-”写在字符串首位，非负整数不带负号。按照以上技术方案如手写整数一样自然，从而回避对各种编程语言现有整型数据类型进行扩展的困难。

在面向对象技术编程时，可将S1中输入的数据单独设计为数据成员。

所述预处理模块是用来删除S1输入数据中所有的空格和删除输入数据前面多余的数字“0”以及判断输入的数据是不是整数。

1.设输入字符串是STR。下面算法删除STR中的所有空格：

(1)令字符型指针变量p指向STR[0]，字符型指针变量q指向R[0]。

(2)判断：若p[0]＝'\0'，则令q[0]＝'\0'，结束。

(3)判断：若p[0]≠″(空格)，则转向第(5)步。

(4)令p＝p+1，转向第(2)步。

(5)令q[0]＝p[0]，p＝p+1，q＝q+1。

(6)转向第(2)步。

算法结束后，结果是字符串R。

2.设输入字符串为STR。下面的算法判定STR是否整数：

(1)删除STR中的所有空格，得到字符串S。

(2)令字符型指针变量p指向S[0]。

(3)判断：若p[0]≠'-'而p[0]>'9'或p[0]<'0'，或者p[0]＝'-'而p[1]>'9'或p[1]<'0'，则令flag＝0，结束。

(4)令p＝p+1。

(5)判断：若p[0]＝'\0'，则令falg＝1，并结束。

(6)判断：若p[0]<'0'或者p[0]>'9'，则令flag＝0，并结束。

(7)令p＝p+1，并转向第(5)步。

算法结束后，若flag的值为1，则STR是整数；否则，若flag的值为0，则STR不是整数。

这个算法允许将整数0表示成负号后面若干个0的形式。

3.设有整数字符串STR。下面的算法删除STR前面所有多余的字符'0'，结果存放在另一个字符串R中：

(1)令字符型指针变量p指向STR[0]，字符型指针变量q指向R[0]。

(2)判断：若p[0]＝'-'，则令q[0]＝'-'，p＝p+1，q＝q+1。

(3)判断：若p[0]≠'0'或p[1]≠'0'，则转向第(5)步。

(4)令p＝p+1，转向第(3)步。

(5)判断：若p[0]＝'0'且p[1]＝'\0'，则令q[0]＝'0'且q[1]＝'\0'，并转向第(10)步。

(6)判断：若p[0]＝'0'，则令p＝p+1。

(7)令q[0]＝p[0]，p＝p+1，q＝q+1。

(8)判断：若p[0]≠'\0'，转向第(7)步。

(9)令q[0]＝'\0'。

(10)判断：若R[0]＝'-'而R[1]＝'0'，则令R[0]＝'0'，R[1]＝'\0'。

(11)结束。

算法结束后，结果是字符串R。

计算模块包括基本运算模块、乘方与开方模块、阶乘与因子模块、素数模块、最大公约数模块、特殊数模块、不定方程模块、同余式模块。计算模块在计算时不改变S1中输入的数据。

基本运算模块：圆周率的近似值，加法，减法，乘法，除法、Euclid除法，整除性测试。

乘方与开方模块：方根化简，乘方，一个正整数在另一正整数中的方次数，平方根。

阶乘与因子模块：阶乘，阶乘因子，正整数的因子分解。

素数模块：最小素因子，素性测试，Euler函数。

最大公约数模块：最大公约数，最小公倍数。

特殊数模块：二项式系数，斐波那契数，第一模块Stirling数，第二模块Stirling数，Catalan数。

不定方程模块：二元一次不定方程ax-by＝1，二元一次不定方程ax+by＝c，三元一次不定方程ax+by+cz＝n，勾股数，整系数二元一次方程。

同余式模块：一个正整数关于另一个正整数的乘法阶，一个正整数关于另一个正整数的乘法逆，一元一次同余方程，一元一次同余方程组。

以下的大小比较和取负算法是为其它算法服务的，不属于计算模块。

1.整数大小比较

设两个整数字符串STRA和STRB的长度分别为lena和lenb。比较整数STRA和STRB大小的算法如下：

(1)令sgna＝0，sgnb＝0，indexa＝0，indexb＝0。

(2)判断：若STRA[0]>'0'，则令sgna＝1；否则，若STRA[0]＝'-'，则令sgna＝-1，indexa＝1。同理，若STRB[0]>'0'，则令sgnb＝1；否则，若STRB[0]＝'-'，则令sgnb＝-1，indexb＝1。

(3)判断：若sgna＝0且sgnb＝0，则令sgn＝0，结束；若sgna>sgnb，则令sgn＝1，结束；若sgna<sgnb，则令sgn＝-1，结束。

(4)当lena≠lenb时，进行如下操作：

若lena<lenb，则令sgn＝-1；若lena>lenb，则令sgn＝1；最后，令sgn＝sgn×sgna，并结束。

(5)令sgn＝0，对于i依次从indexa到lena-1的每一个值，比较STRA[i]和STRB[i]，若首次遇到一个i值使STRA[i]<STRB[i]或STRA[i]>STRB[i]，则对应地令sgn＝-1或sgn＝1，结束比较；最后，令sgn＝sgn×sgna，并结束。

算法结束后，sgn的值等于-1、0或1，相应地，表示STRA<STRB、STRA＝STRB或STRA>STRB。

2.取负运算

设整数字符串STR的长度为len。对STR取负的算法如下：

(1)判断：若STR[0]＝'0'，则令R＝"0"，并结束。

(2)判断：若STR>"0"，则令R[len+1]＝'\0'，并且对于i依次从len到1的每一个值，令R[i]＝STR[i-1]；最后令R[0]＝'-'并结束。

(3)判断：若STR<"0"，则对于i依次从0到len-2的每一个值，令R[i]＝STR[i+1]。最后令R[len-1]＝'\0'，并结束。

算法结束后，R＝-STR。

以下为计算模块：

(一)基本运算类

1.简单加法

设置一个普通整型变量CF作进位标志，模仿CPU的二进制计算过程，也模拟人的手工计算过程。

设两个非负整数字符串STRA和STRB的长度分别为lena和lenb。实现整数STRA、STRB简单加法运算STRA+STRB的算法如下：

(1)判断：若STRA＝"0"，则令R＝STRB，结束；若STRB＝"0"，则令R＝STRA，结束。

(2)判断：若lena≥lenb，则令X＝STRA，Y＝STRB，m＝lena，n＝lenb；否则，令X＝STRB，Y＝STRA，m＝lenb，n＝lena。

(3)令CF＝0，R[m+1]＝'\0'；

(4)对于i依次从m-1到m-n以及j依次从n-1到0(i、j同步递减)的每一个值，计算t＝(X[i]-'0')+(Y[j]-'0')+CF，然后判断：若t≤9，则令R[i+1]＝'0'+t，CF＝0；否则令R[i+1]＝'0'+t-10，CF＝1。

(5)判断：若m>n，则对于i依次从m-n-1到0的每一个值，计算t＝(X[i]-'0')+CF，然后判断：若t≤9，则令R[i+1]＝'0'+t，CF＝0；否则，令R[i+1]＝'0'+t-10，CF＝1。

(6)令R[0]＝'0'+CF。

(7)调用预处理算法，删除字符串R前面可能多余的数字'0'，结果仍存放在R中。

(8)结束。

算法结束后，R＝STRA+STRB。

2.简单减法

设置一个普通整型变量CF作借位标志，模仿CPU的二进制计算过程，也模拟人的手工计算过程。

设两个非负整数字符串STRA和STRB的长度分别为lena和lenb。实现整数STRA、STRB简单减法运算STRA-STRB的算法如下：

(1)判断：若STRA＝"0"，则令R＝-STRB，结束；若STRB＝"0"，则令R＝STRA，结束。

(2)判断：若lena≥lenb，则令X＝STRA，Y＝STRB，m＝lena，n＝lenb，sgn＝1；否则，令X＝STRB，Y＝STRA，m＝lenb，n＝lena，sgn＝-1。

(3)令CF＝0，R[m]＝'\0'。

(4)对于i依次从m-1到m-n以及j依次从n-1到0(i、j同步递减)的每一个值，计算t＝X[i]-Y[j]-CF，然后判断：若t<0，则令R[i]＝'0'+t+10，CF＝1；否则，令R[i]＝'0'+t，CF＝0。

(5)判断：若m>n，则对于i依次从m-n-1到0以及k依次从m-n到1(j、k同步递减)的每一个值，计算t＝(X[i]-'0')-CF，然后判断：若t<0，则令R[i]＝'0'+t+10，CF＝1；否则，令R[i]＝'0'+t，CF＝0。

(6)调用预处理算法，删除字符串R前面可能多余的数字'0'，结果仍存放在R中。

(7)判断：若sgn×(1-2×CF)<0，则令R＝-R。

(8)结束。

算法结束后，R＝STRA-STRB。

3.加法

设有两个整数字符串STRA和STRB。实现整数STRA、STRB加法运算STRA+STRB的算法如下：

(1)判断：若STRA＝"0"，则令R＝STRB，结束；若STRB＝"0"，则令R＝STRA，结束。

(2)判断：若STRA>"0"，则令sgna＝1，X＝STRA；否则，若STRA<"0"，则令sgna＝-1，X＝-STRA。同理，若STRB>"0"，则令sgnb＝1，Y＝STRB；否则，若STRB<"0"，则令sgnb＝-1，Y＝-STRB。

(3)令sgn＝sgna×sgnb。

(5)判断：若sgn>0，则根据简单加法计算R＝X+Y。

(6)判断：若sgn<0，则根据简单减法计算R＝X-Y；

(7)判断：若sgna<0，则令R＝-R。

(8)结束。

算法结束后，R＝STRA+STRB。

4.减法

设两个整数字符串STRA和STRB。实现整数STRA、STRB减法运算STRA-STRB的算法如下：

(1)T＝-STRB。

(2)用加法算法计算R＝STRA+T。

(3)结束。

算法结束后，R＝STRA-STRB。

5.简单乘法

设置一个普通整型变量CF作进位标志，模拟手工计算过程。

设非负整数字符串STR的长度为len，数字字符ch满足'0'≤ch≤'9'。实现整数STR乘ch运算的算法如下：

(1)判断：若STR＝"0"或ch＝'0'，则令R＝"0"，结束。

(2)判断：若ch＝'1'，则令R＝STR，结束。

(3)令CF＝0，R[len+1]＝"\0"。

(4)对于i依次从len-1到0的每一个值，计算t＝(STR[i]-'0')×(ch-'0')+CF，然后判断：若t≤9，则令R[i+1]＝'0'+t，CF＝0；否则令R[i+1]＝'0'+t％10，CF＝t/10。

(5)令R[0]＝'0'+CF。

(6)调用预处理算法，删除字符串R前面可能多余的数字'0'。

(7)结束。

算法结束后，R＝STR×ch。

6.乘法

设两个整数字符串STRA和STRB的长度分别为lena和lenb。实现整数STRA、STRB乘法运算STRA×STRB的算法如下：

(1)判断：若STRA＝"0"或STRB＝"0"，则令R＝"0"，结束。

(2)判断：若STRA＝"1"，则令R＝STRB，结束；若STRB＝"1"，则令R＝STRA，结束。

(3)判断：若STRA>"0"，则令sgna＝1，X＝STRA；否则，若STRA<"0"，则令sgna＝-1，X＝-STRA，lena＝lena-1。同理，若STRB>"0"，则令sgnb＝1，Y＝STRB；否则，若STRB<"0"，则令sgnb＝-1，Y＝-STRB，lenb＝lenb-1。

(4)令sgn＝sgna×sgnb，i＝0，R＝""。

(5)判断：若i＝lenb，则转向第(10)步。

(6)计算字符串R的长度len，并令R[len]＝'0'，R[len+1]＝'\0'。

(7)根据简单乘法计算C＝X×Y[i]。

(8)根据加法算法计算R＝R+C；

(9)令i＝i+1，并转向第(5)步。

(10)调用预处理算法，删除字符串R前面可能多余的数字'0'。

(11)判断：若sgn<0，则令R＝-R。

(12)结束。

算法结束后，R＝STRA×STRB。

7.Euclid除法

整数a关于正整数b的Euclid除法指的是求满足a＝bq+r(0≤r＜b)的唯一对整数q和r。

设整数字符串STRA和正整数字符串STRB的长度分别为lena和lenb。实现整数STRA关于整数STRB的Euclid除法运算的算法如下：

(1)判断：若STRA＝"0"，则令Q＝"0"，R＝"0"，结束。

(2)判断：若STRB＝"1"，则令Q＝STRA，R＝"0"，结束。

(3)令Y＝STRB。若STRA[0]>'0'，则令sgna＝1，X＝STRA；否则，令sgna＝-1，X＝-STRA，lena＝lena-1。

(4)当X＝Y时，令R＝"0"，并判断：若sgna>0，则令Q＝"1"，结束；若sgna<0，则令Q＝"-1"，结束。

(5)当X<Y时，判断：若sgna>0，则令Q＝"0"，R＝X，结束；若sgna<0，则令Q＝"-1"，R＝Y-X，结束。

(6)令Q＝"0"。

(7)对于i依次从0到lenb-1的每一个值，令R[i]＝X[i]。然后令R[lenb]＝'\0'。

(8)判断：若R≥Y，则对R反复减Y(结果仍然存放在R中)，对Q反复加"1"(结果仍然存放在Q中)，直到R<Y为止。

(9)令i＝lenb。

(10)判断：若i≥lena，则转向第(14)步。

(11)令R＝R×"10"，R＝R+X[i]；令Q＝Q×"10"。

(12)判断：若R≥Y，则对R反复减Y(结果仍然存放在R中)，对Q反复加"1"(结果仍然存放在Q中)，直到R<Y为止。

(13)令i＝i+1，转向第(10)步。

(14)判断：若sgna≥0，则结束。

(15)判断：若R＝"0"，则令Q＝-Q；否则，令Q＝Q+"1"，Q＝-Q，R＝Y-R。

(16)结束。

算法结束后，Q、R的值满足STRA＝Q×STRB+R，0≤R<STRB。

8.除法

设非负整数字符串STRA和正整数字符串STRB的长度分别为lena和lenb。实现整数STRA÷STRB的除法运算的算法如下：

(1)输入除法结果小数部分的最大位数n。

(2)对STRA和STRB施行Euclid除法，得到STRA＝STRB×Q+R。

(3)令S＝""，T＝""，i＝0。

(4)判断：若R＝"0"或i≥n，则结束。

(5)令X＝R×"10"。

(6)对X和STRB施行Euclid除法，得到X＝STRB×T+R。

(7)把字符串T连接在字符串S后面。

(8)令i＝i+1，转向第(4)步。

算法中，T是临时字符串变量。

算法结束后，STRA÷STRB的整数部分为Q，小数部分为S，小数位数至多为n。小数部分不进行四舍五入，各位都是精确数字。若小数位数小于n，说明最后一次除法除尽了，商是精确值。

9.整除性测试

设有正整数字符串STRA和非负整数字符串STRB。判定STRA是否整除STRB的算法如下：

(1)对STRB和STRA施行Euclid除法，得到STRB＝STRA×Q+R。

(2)判断：若R＝"0"，则STRA整除STRB；否则，STRA不整除STRB。

10.圆周率的近似值

计算圆周率π的方法很多。比如可以用下面的数学公式

$\mathrm{\&pi;}=4(1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\frac{1}{9}-\frac{1}{11}+L+\frac{{(-1)}^{n-1}}{2n-1}+L)$

来计算π。但用这样的公式不能得到π的很精确的近似分数。另外，这个级数收敛很慢，取前面2n项，才能保证精度达到若要精确到小数点后m位，就要取前面2×10^m项，时间和内存开销都以指数方式递增，是一般计算机难以承受的。所以本发明采用近似分数来计算π值。

π的第n个近似分数的分子和分母分别按下面的递归公式计算：

P₀＝1，P₁＝q₁，P_n＝q_nP_n-1+P_n-2；

Q₀＝0，Q₁＝1，Q_n＝q_nQ_n-1+Q_n-2。

其中q_i是连分数中的第i个不完全商，可以在公开的文献中查到前31个q_i。第i个近似分数的误差为1/(Q_i-1Q_i)。

设π的连分数的第n(1≤n≤31)个不完全商存放在字符型指针数组q的数组元素q[n-1]中。例如前31个不完全商为

char*q[31]＝{"3","7","15","1","292","1","1","1","21","31","14","2","1","2","2","2","2","1","84","2","1","1","15","3","13","1","4","2","6","6","1"}；

求π的第n个近似值的算法如下：

(1)令PP＝"3"，QQ＝"1"，P0＝"0"，P1＝"1"，Q0＝"0"，Q1＝"0"，TT＝""。

(2)判断：若n＝1，则结束。

(3)令i＝2，并判断：若i>n，则转向第(7)步。

(4)令P0＝P1，P1＝PP，TT＝q[i-1]×P1，PP＝TT+P0，TT＝""。

(5)令Q0＝Q1，Q1＝QQ，TT＝q[i-1]×Q1，QQ＝TT+Q0，TT＝""。

(6)令i＝i+1，转向第(4)步。

(7)令TT＝QQ×Q1，并令i等于TT的长度。

(8)指定小数位数i+1，对PP和QQ施行除法，商的整数部分为"3"，小数部分存放在S中。

(9)结束。

算法结束后，圆周率的第n(1≤n≤31)个近似分数为PP/QQ，其整数部分为I＝"3"，小数部分为S，理论精度小于1/TT。

简单加法、简单减法和简单乘法都是为其它算法服务的，它们都不属于计算模块。

(二)乘方与开方类

1.一个正整数在另一正整数中的方次数。

正整数k＞1在正整数n中的方次数deg(k,n)是使k^t整除n的最大幂指数t。

设有两个正整数字符串K和N。求DEG(K，N)的算法如下：

(1)令X＝N，Q＝""，R＝""，DEG＝"0"。

(2)对X和K施行Euclid除法，得到商Q和余数R。

(3)判断：若R≠"0"，则结束。

(4)令DEG＝DEG+"1"，X＝Q。

(5)对X和K施行Euclid除法，得到商Q和余数R。

(6)转向第(3)步。

算法结束后，DEG＝DEG(K，N)。

2.乘方

设有两个非负整数字符串N和K。求N的K次方的算法如下：

(1)判断：若K＝"0"，则令P＝"1"，结束。

(2)判断：若N<"2"，则令P＝N，结束。

(3)令P＝N，T＝K-"1"。

(4)判断：若T＝"0"，则结束。

(5)令P＝P×N，T＝T-"1"，转向第(4)步。

算法默认0⁰＝1。算法结束后，P就是N的K次方值。

3.方根化简

方根化简就是将根号下能开方的因子开方成为根号前面的系数。例如 $\sqrt[3]{1080}=6\sqrt[3]{5}.$

设有正整数字符串N。化简N的K次方根的算法要调用后面求最小素因子的算法。化简算法如下：

(1)令D＝""，X＝N，P＝""，Q＝""，R＝""，T＝""，A＝"1"，B＝"1"。

(2)判断：若K＝"1"，则令A＝X，结束。

(3)判断：若X＝"1"，则结束。

(4)求X的最小素因子P，并求D＝DEG(P，X)。

(5)对D和K施行Euclid除法，得商Q和余数R。

(6)求P的R次方T，并令B＝B×T，R＝Q-R。

(7)求P的R次方Q，并令Q＝T×Q，A＝A×Q。

(8)求Q的K次方R，并令T＝R×T。

(9)对X和T施行Euclid除法，得商Q和余数R＝"0"。

(10)令X＝Q，转向第(3)步。

算法结束后，N的K次方根被化简为A的形式。

4.平方根

设a是一个正整数，b是十进制数字0～9之一。显然

(10a+b)²-100a²＝(20a+b)b。

这是手工开平方的原理。据此可设计算法模拟手工开方过程。

设非负整数字符串N的长度为lenn。求N的平方根直到小数点后num位的算法如下：

(1)P＝""，Q＝""，R＝""，S＝"0"，T＝""，A＝""，B＝"0"，Root＝""，root＝""，i＝0，len＝0。

(2)判断：若N＝"0"，则令Root＝"0"，结束。

(3)令len＝lenn。

(4)判断：若len是奇数，则令P[0]＝N[0]，P[1]＝'\0'，i＝1；否则，令P[0]＝N[0]，P[1]＝N[1]，P[2]＝'\0'，i＝2。

(5)判断：若S>P，则转向第(8)步。

(6)令A＝B，B＝B+"1"，S＝B×B。

(7)转向第(5)步。

(8)令Root＝A，S＝A×A，R＝P-S。

(9)判断：若i≥len，则转向第(17)步。

(10)令P＝R×"10"，P＝P+N[i]，P＝P×"10"，P＝P+N[i+1]。

(11)令i＝i+2，T＝Root×"20"，B＝"0"，S＝"0"。

(12)判断：若S>P，则转向第(15)步。

(13)令A＝B，R＝P-S，B＝B+"1"，Q＝T+B，S＝Q×B。

(14)转向第(12)步。

(15)令Root＝Root×"10"，Root＝Root+A。

(16)转向第(9)步。

(17)令i＝0，len等于整数字符串Root的当前长度。

(18)判断：若R＝"0"或者i≥num，则转向第(25)步。

(19)令i＝i+1，P＝R×"100"，T＝Root×"20"，B＝"0"，S＝"0"。

(20)判断：若S>P，则转向第(23)步。

(21)令A＝B，R＝P-S，B＝B+"1"，Q＝T+B，S＝Q×B。

(22)转向第(20)步。

(23)令Root＝Root×"10"，Root＝Root+A。

(24)转向第(18)步。

(25)令root等于Root中从Root[len]开始的子字符串，Root[len]＝'\0'。

(26)结束。

算法中整数字符串A和B专门用来存放十进制数字，长度声明为2就行了。但整数字符串Root要长到足以存放N平方根的整数部分和小数部分。算法结束后，N的平方根的整数部分是整数字符串Root，小数部分是整数字符串root。若root是空字符串，则说明N本身是平方数。

(三)阶乘与因子类

1.阶乘

设有非负整数字符串N。求N阶乘的算法如下：

(1)令X＝N，F＝"1"。

(2)判断：若X<"2"，则结束。

(3)判断：若X＝"1"，则结束。

(4)令F＝F×X，X＝X-"1"，转向第(3)步。

算法结束后，N阶乘在F中。

2.阶乘因子

阶乘因子问题就是对于一个素数p和一个正整数n求deg(p,n！)。根据数学上的结论，有

$\mathrm{deg}(p,n!)=\underset{k=1}{\overset{\left[{\log }_{p}n\right]}{\mathrm{\&Sigma;}}}{m}_k,m_0=n,m_k=\left[\frac{m_{k-1}}{p}\right].$

设有素数字符串P和正整数字符串N。求DEG(P，N！)的算法如下：

(1)令M＝N，D＝"0"，Q＝""，R＝""。

(2)判断：若M<P，则结束。

(3)对M和P施行Euclid除法，得到商Q和余数R。

(4)令D＝D+Q，M＝Q，转向第(2)步。

算法结束后，D＝DEG(P，N！)。

3.正整数的因子分解

正整数的因子分解问题就是求正整数n＞1的标准分解式其中诸e_k都是正整数，诸p_k都是素数，且p₁＜p₂＜L＜p_K。

设有正整数字符串N。分解N的算法要调用后面求最小素因子的算法。算法如下：

(1)令X＝N，T＝""，P＝""，D＝""，Q＝""，R＝""，len＝0。

(2)判断：若X<"4"，则令A[0]＝"1"，A[1]＝X，A[2]＝"1"，结束。

(3)令A[0]＝"0"，len＝0。

(4)判断：若X＝"1"，则结束。

(5)令A[0]＝A[0]+"1"，len＝len+1。

(6)求X的最小素因子P，并求D＝DEG(P，X)。

(7)令A[2×len-1]＝P，A[2×len]＝D。

(8)求P的D次方T。

(9)对X和T施行Euclid除法，得到商Q和余数R＝"0"。

(10)令X＝Q，转向第(4)步。

算法中A是一个字符串数组。

算法结束后，A[0]中存放的是N的不同素因子个数len，A[2×k-1]和A[2×k]中分别存放的是N的第k(1≤k≤len)个素因子及其在N中的方次数。

方根化简问题和素性判定问题本质上是正整数的因子分解问题。正整数的因子分解问题是很难的数学问题。由于对于较大的正整数，因子分解问题的计算时间长得不能忍受，因此人们转而考虑概率算法。而概率算法不能确保成功，因此这里未采用概率算法。

(四)素数类

1.最小素因子

设正整数字符串N>1的长度为len。求N的最小素因子的算法如下：

(1)令D＝""，SQRT＝""，Q＝""，R＝""。

(2)判断：若N[len-1]∈{'0'，'2'，'4'，'6'，'8'}，则令D＝"2"，结束。

(3)判断：若N[len-1]＝'5'，则令D＝"5"，结束。

(4)求N的平方根的整数部分SQRT，并令I＝"3"。

(5)判断：若I>SQRT，则令D＝N，结束。

(6)对N和I施行Euclid除法，得到商Q和余数R。

(7)判断：若R＝"0"，则令D＝I，结束。

(8)令I＝I+"2"，转向第(5)步。

算法结束后，D是N的最小素因子。

2.素性测试

设有整数字符串N。判定N是否素数的算法如下：

(1)判断：若N<"2"，则N不是素数，结束。

(2)求N的最小素因子D。

(3)判断：若D<N，则N不是素数，结束；否则N是素数，结束。

3.Euler函数

Euler函数φ(n)的值等于集合{1,2,L,n}中与n互素的整数的个数。计算方法是先求得n的标准分解式然后按照下面公式计算：

$\mathrm{\&phi;}\left(n\right)=\left((p_1-1)p_1^{e_1-1}\right)\left((p_2-1)p_2^{e_2-1}\right)L\left((p_K-1)p_K^{e_K-1}\right)$

设有正整数字符串N。求Euler函数值φ(N)的算法如下：

(1)令E＝""，PHI＝"1"，DEG＝"0"。

(2)判断：若N<"3"，则结束。

(3)调用因子分解算法分解N，结果存放在字符串数组A中。

(4)令I＝"1"，i＝1。

(5)判断：若I>A[0]，则结束。

(6)令E＝A[2×i-1]，E＝E-"1"，DEG＝A[2×i]-"1"。

(7)令J＝"1"，j＝1。

(8)判断：若J>DEG，则转向第(11)步。

(9)令E＝E×A[2×i-1]。

(10)令J＝J+"1"，j＝j+1，转向第(8)步。

(11)令PHI＝PHI×E。

(12)令I＝I+"1"，i＝i+1，转向第(5)步。

算法中A是一个字符串数组，A的每一个元素都初始化为空字符串。

算法结束后，PHI＝φ(N)。

(五)最大公约数类

1.最大公约数

设有正整数字符串M和N。求M、N的最大公约数GCD(M，N)的辗转相除算法如下：

(1)令A＝""，B＝M，Q＝""，R＝N。

(2)令A＝B，B＝R。

(3)对A和B施行Euclid算法，得到商Q和余数R。

(4)判断：若R＝"0"，则转向第(2)步。

算法结束后，GCD(M，N)＝B。

两个非零整数X、Y的指的就是它们绝对值的最大公约数。

2.最小公倍数

设有正整数字符串M和N。求M、N的最小公倍数LCM(M，N)的算法如下：

(1)令LCM＝""，Q＝""，R＝""，。

(2)令LCM＝GCD(M，N)。

(3)对M和LCM施行Euclid算法，得到商Q和余数R＝"0"。

(4)令LCM＝Q×N，结束。

算法结束后，LCM就是M、N的最小公倍数。

(六)特殊数类

1.二项式系数

根据下面的一组公式进行计算：

$C_n^k=\frac{n(n-1)L(n-k+1)}{k!},C_n^k=C_n^{n-k},C_n^0=1,C_n^1=n$

设有非负整数字符串N和K。求二项式系数C(N，K)的算法如下：

(1)令F＝"1"，I＝"0"，M＝""，P＝"1"，Q＝""，R＝""。

(2)判断：若N<K，则令Q＝"0"，结束。

(3)对N和"2"施行Euclid除法，得到商Q和余数R。

(4)判断：若K>Q，则令M＝N-K。

(5)判断：若M＝"0"，则令Q＝"1"，结束。

(6)判断：若M＝"1"，则令Q＝N，结束。

(7)令F＝N-I，P＝P×F，I＝I+"1"。

(8)判断：若I<M，则转向第(7)步。

(9)计算M的阶乘，结果为F。

(10)对P和F施行Euclid除法，得到商Q和余数R＝"0"。

算法结束后，Q＝C(N，K)。

2.斐波那契数

根据递归公式

F₁＝1，F₂＝1，F_n＝F_n-1+F_n-2，n≥3

进行计算。

设有正整数字符串N。求斐波那契数F(N)的算法如下：

(1)令K＝N，F＝"0"，F1＝"1"，F2＝"0"，T＝"0"。

(2)判断：若K＝"0"，则结束。

(3)令F＝F1+F2，F2＝F1，F1＝F。

(4)令K＝K-"1"，转向第(2)步。

算法结束后，F＝F(N)。

3.第一类Stirling数

主要根据递归公式

s(n,k)＝0，n＜k；

s(n,n)＝1，n≥0；

s(n,0)＝0，n≥1；

s(n,k)＝s(n-1,k-1)-(n-1)s(n-1,k)，n≥k≥1

进行计算。计算过程中，可以利用以下性质：

$s(n,1)={(-1)}^{n-1}!,s(n,n-1)=-C_n^2.$

设有非负整数字符串N和K，N的长度为len。第一类Stirling数S(N，K)的算法如下：

(1)令T＝""，T0＝""，T1＝""，I＝""，STR＝""，i＝0，j＝0。

(2)判断：若N<K，或者N>0而K＝0，则令STR＝"0"，结束。

(3)判断：若N＝K，则令STR＝"1"，结束。

(4)令I＝N-"1"。

(5)判断：若I≠K，则转向第(9)步。

(6)令STR＝N×I。

(7)对STR和"2"施行Euclid除法，得到商T0和余数T1＝"0"。

(8)令STR＝-T0，并结束。

(9)判断：若K≠"1"，则转向第(13)步。

(10)计算I的阶乘T0。

(11)判断：若(N[len-1]-'0')％2＝1，则令STR＝T0；否则，令STR＝-T0。

(12)结束。

(13)令S[0]＝"0"，S[1]＝"1"，T＝"0"，I＝"2"，i＝2。

(14)判断：若I>N，则结束。

(15)令S[i]＝"1"，T＝I-"1"，J＝T，j＝i-1。

(16)判断：若j≤0，则转向第(20)步。

(17)令T0＝S[j-1]，T1＝S[j]，T1＝T1×T，S[j]＝T0-T1。

(18)判断：若I＝N或者J＝K，则令STR＝S[j]，并转向第(20)步。

(19)令J＝J-"1"，j＝j-1，转向第(16)步。

(20)令I＝I+"1"，i＝i+1，转向第(14)步。

算法中S是一个字符串数组，S的每一个元素都初始化为空字符串。

算法结束后，STR＝S(N，K)。

4.第二类Stirling数

主要根据递归公式

S(n,k)＝0，n＜k；

S(n,n)＝1，n≥0；

S(n,0)＝0，n≥1；

S(n,k)＝S(n-1,k-1)+kS(n-1,k)，n≥k≥1

进行计算。计算过程中，可以利用以下性质：

$s(n,1)=1,s(n,2)=2^{n-1}-1,s(n,n-1)=C_n^2.$

设有非负整数字符串N和K，N的长度为len。第二类Stirling数S(N，K)的算法如下：

(1)令T0＝""，T1＝""，N＝""，K＝""，I＝""，J＝""，i＝0，j＝0。

(2)判断：若N<K，或者N>0而K＝0，则令STR＝"0"，结束。

(3)判断：若N＝K或者K＝"1"，则令STR＝"1"，结束。

(4)令J＝N-"1"。

(5)判断：若J≠K，则转向第(9)步。

(6)令STR＝N×J。

(7)对STR和"2"施行Euclid除法，得到商T0和余数T1＝"0"。

(8)令STR＝T0，并结束。

(9)判断：若K≠"2"，则转向第(11)步。

(10)计算2的J次方STR，并令STR＝STR-1，结束。

(11)令S[0]＝"0"，S[1]＝"1"，T＝"0"，I＝"2"，i＝2。

(12)判断：若I>N，则结束。

(13)令S[i]＝"1"，T＝I-"1"，J＝T，j＝i-1。

(14)判断：若j≤0，则转向第(18)步。

(15)令T0＝S[j-1]，T1＝S[j]，T1＝T1×J，S[j]＝T0+T1。

(16)判断：若I＝N或者J＝K，则令STR＝S[j]，并转向第(18)步。

(17)令J＝J-"1"，j＝j-1，转向第(14)步。

(18)令I＝I+"1"，i＝i+1，转向第(12)步。

算法中S是一个字符串数组，S的每一个元素都初始化为空字符串。

算法结束后，STR＝S(N，K)。

5.Catalan数

Catalan数就是满足递归关系

$h_1=1,h_n=\underset{k=1}{\overset{n-1}{\mathrm{\&Sigma;}}}{h}_k{h}_{n-k}$

的数列。其公式为

$h_n=\frac{1}{n}\left(\begin{array}{c}2n-2\\ n-1\end{array}\right)$

设有正整数字符串N。求Catalan数H(N)的算法如下：

(1)令C＝""，M＝""，K＝""，H＝""。

(2)判断：若N≤"2"，则令C＝"1"，结束。

(3)令K＝N-1，M＝2×K。

(4)计算二项式系数C＝C(M，K)。

(5)对C和N施行Euclid除法，得到商H和余数M＝"0"。

算法结束后，H＝H(N)。

(七)不定方程类

要求整数解的方程称为不定方程。

1.二元一次不定方程ax-by＝1

设a＞0，b＞0。二元一次不定方程ax-by＝1有解的充分必要条件是GCD(a,b)＝1。当这一条件满足时，若(x₀,y₀)是此方程的一个特解，则此方程的所有解为x＝x₀+bt，y＝y₀+at。其中t为任意整数。因此，只须设计算法解决求一个特解的问题。

求特解的方法是用辗转相除法把分数表示成连分数然后按照下面的递归公式计算P_k和Q_k：

P₀＝1，P₁＝q₁，P_k＝q_kP_k-1+P_k-2；

Q₀＝0，Q₁＝1，Q_k＝q_kQ_k-1+Q_k-2。

则得到一个特解x₀＝(-1)ⁿQ_n-1，y₀＝(-1)ⁿP_n-1。实践中一般还要求(x₀,y₀)满足0≤x₀＜b，0≤y₀＜a，这时候可以对求得的x₀和y₀分别加上或减去b和a的相同倍数。

设正整数字符串A和B满足GCD(A，B)＝"1"。求二元一次不定方程A×X-B×Y＝"1"的一个特解的算法如下：

(1)PP＝"0"，P1＝"0"，P2＝"1"，QQ＝"0"，Q1＝"1"，Q2＝"0"，T＝"0"，Q＝"0"，R＝"0"，TA＝"0"，TB＝"0"，X＝"0"，Y＝"0"，sgn＝0。

(2)对A和B施行Euclid除法，得到商Q和余数R。

(3)令P1＝Q，TA＝B，TB＝R，sgn＝-1。

(4)判断：若R＝"0"，则转向第(9)步。

(5)对TA和TB施行Euclid除法，得到商Q和余数R。

(6)令TA＝TB，TB＝R，T＝Q×P1，PP＝T+P2，T＝Q×Q1，QQ＝T+Q2，sgn＝-sgn，

(7)令P2＝P1，Q2＝Q1，P1＝PP，Q1＝QQ。

(8)转向第(4)步。

(9)判断：若sgn>0，则令X＝Q2，Y＝P2；否则，令X＝-Q2，Y＝-P2。

(10)对X和B施行Euclid除法，得到商Q和余数R。

(11)令T＝A×Q，X＝R，Y＝Y-T，结束。

算法结束后，(X，Y)是所求的一个特解。

2.二元一次不定方程ax+by＝c

设ab≠0。二元一次不定方程ax+by＝c有解的充分必要条件是GCD(|a|,|b|)整除c，因此可设GCD(a,b)＝1。当这一条件满足时，若(x₀,y₀)是此方程的一个特解，则此方程的所有解为x＝x₀-bt，y＝y₀+at。其中t为任意整数。因此，只须设计算法解决求一个特解的问题。

求特解的方法是先求二元一次不定方程|a|x-|b|y＝1的一个特解(x₀,y₀)，得到二元一次不定方程ax+by＝1的一个特解(x₀sgna,-y₀sgnb)，然后将这个特解乘以c，得到原二元一次不定方程ax+by＝c的一个特解，再对这个特解进行化简。

设非零整数字符串A、B和C满足GCD(A，B)＝"1"的条件，并设A、B的符号分别为sgna和sgnb。求二元一次不定方程A×X+B×Y＝C的一个特解的算法如下：

(1)令TA＝A，TB＝B，X＝"0"，Y＝"0"。

(2)判断：若sgna<0，则令TA＝-TA；若sgnb<0，则令TB＝-TB。

(3)求TA×X-TB×Y＝1的一个特解(X，Y)。

(4)判断：若sgna<0，则令X＝-X；若sgnb>0，则令Y＝-Y。

(5)令X＝X×C，Y＝Y×C。

(6)对X和TB施行Euclid除法，得到商Q和余数R。

(7)令T＝A×Q，并判断：若sgnb<0，则令T＝-T。

(8)令X＝R，Y＝Y+T，结束。

算法结束后，(X，Y)是所求的一个特解。

3.三元一次不定方程ax+by+cz＝n

设abc≠0。三元一次不定方程ax+by+cz＝n有解的充分必要条件是GCD(a,b,c)整除n，因此可设GCD(a,b,c)＝1。当这一条件满足时，若(x₀,y₀,z₀)是此方程的一个特解，则此方程的所有解为x＝x₀-b's-u₀ct，y＝y₀+a's-v₀ct，z＝z₀+dt。其中d＝GCD(a,b)，a'＝a/d，b'＝b/d，(u₀,v₀)是二元一次不定方程a'u+b'v＝1的一个特解，t为任意整数。因此，只须设计算法，求得所有这些系数，并求得一个特解。

设非零整数字符串A、B、C和整数字符串N满足GCD(A，B，C)＝"1"的条件。求三元一次不定方程A×X+B×Y+C×Z＝N的一个特解的算法如下：

(1)令AP＝""，BP＝""，U＝""，V＝""，W＝""，X＝""，Y＝""，Z＝""。

(2)求D＝GCD(A，B)。

(3)对A和D施行Euclid算法，得到商U和余数V＝"0"，并令AP＝U。

(4)对B和D施行Euclid算法，得到商U和余数V＝"0"，并令BP＝U。

(5)求二元一次不定方程AP×U+BP×V＝"1"的一个特解(U，V)。

(6)求二元一次不定方程D×W+C×Z＝N的一个特解(W，Z)。

(7)令X＝U×W，Y＝V×W。

(8)令U＝U×C，V＝V×C，结束。

算法结束后，(X，Y，Z)就是所求的一个特解(x₀,y₀,z₀)，D的值就是d＝gcd(a,b)，AP、BP的值分别就是a'、b'，U、V的值分别就是u₀c、v₀c。

4.勾股数

满足方程x²+y²＝z²的正整数组(x,y,z)称为勾股数。若(x,y,z)是勾股数，则(y,x,z)显然也是勾股数。下面是一个勾股数公式

x＝2ab，y＝a²-b²，y＝a²+b²，a＞b＞0。

现在要求对于输入的正整数a＞b＞0求相应的勾股数。

设整数字符串A、B满足A>B的条件。下面的算法按上面公式求相应的勾股数：

(1)令X＝""，Y＝""，Z＝""，S＝""，T＝""。

(2)令S＝2×A，X＝S×B。

(3)令S＝A×A，T＝B×B。

(4)令Y＝S-T，Z＝S+T，结束。

算法结束后，(X，Y，Z)就是所求的勾股数。

5.整系数一元二次方程

设a＞0，GCD(a,b,c)＝1。对于整系数二元一次方程ax²+bx+c＝0，要求按求根公式

$x=\frac{-b\&PlusMinus;\sqrt{b^2-4\mathrm{ac}}}{2a}$

求得其根，并对根号和分式进行化简。

(1)令Delta＝""，T＝""，P＝""，Q＝""，R＝""，I＝""，i＝0，sgn＝0。

(2)令Delta＝B×B，T＝4×A，T＝T×C，Delta＝Delta-T。

(3)令S[0]＝-B，S[2]＝2×A，T＝GCD(S[0]，S[2])。

(4)若Delta≠"0"，则转向第(9)步。

(5)对S[0]和T施行Euclid除法，得到商Q和余数R＝"0"。

(6)令S[0]＝Q，S[1]＝"0"，S[3]＝"0"。

(7)对S[2]和T施行Euclid除法，得到商Q和余数R＝"0"。

(8)令S[2]＝Q，结束。

(9)令S[1]＝"1"，S[3]＝"1"，sgn＝1，并判断：若Delta<"0"，则令Delta＝-Delta，sgn＝-1。

(10)用因子分解算法分解Delta，结果存放在字符串数组D中。

(11)令I＝"1"，i＝1。

(12)判断：若I>D[0]，则转向第(18)步。

(13)对D[2×i]和"2"施行Euclid除法，得到商Q和余数R。

(14)求D[2×i-1]的Q次方P。

(15)令S[1]＝S[1]×P。

(16)判断：若R>0，则令S[3]＝S[3]×D[2×i-1]。

(17)令I＝I+"1"，i＝i+1，转向第(12)步。

(18)令T＝GCD(T，S[1])。

(19)对S[0]和T施行Euclid除法，得到商Q和余数R。

(20)令S[0]＝Q。

(21)对S[1]和T施行Euclid除法，得到商Q和余数R。

(22)令S[1]＝Q。

(23)对S[2]和T施行Euclid除法，得到商Q和余数R。

(24)令S[2]＝Q，结束。

算法中S、D都是字符串数组，它们的每一个元素都初始化为空字符串。S有4个数组元素。算法结束后，原方程的根的形式是

$X=\left\{\begin{array}{cc}\frac{S\left[0\right]}{S\left[2\right]},& \mathrm{sgn}=0;\\ \frac{S\left[0\right]\&PlusMinus;S\left[1\right]\sqrt{S\left[3\right]}}{S\left[2\right]}& \mathrm{sgn}=1;\\ \frac{S\left[0\right]\&PlusMinus;\mathrm{iS}\left[1\right]\sqrt{S\left[3\right]}}{S\left[2\right]}& \mathrm{sgn}=-1.\end{array}\right.$

其中sgn是Delta的符号值，i是虚数单位。

(八)同余式类

1.一个正整数关于另一个正整数的乘法阶

一个正整数a关于另一个正整数m＞1的乘法阶ORD(a,m)等于使m整除a^k-1的最小正整数。这里a与m必须满足GCD(a,m)＝1。一般地，ORD(a,m)是Euler函数值φ(m)的一个因子。但是计算φ(m)要对m进行因子分解，徒增计算量。因此，这里采用直接搜索法。

设两个正整数字符串A和M>"1"满足GCD(A，M)＝"1"。计算ORD(A，M)的算法如下：

(1)令P＝A，Q＝""，R＝""，D＝"1"。

(2)对A和M施行Euclid除法，得到商Q和余数R。

(3)令AR＝R。

(4)判断：若R＝"1"，则结束。

(5)令D＝D+"1"，P＝R×AR。

(6)对P和M施行Euclid除法，得到商Q和余数R。

(7)转向第(4)步。

算法结束后D＝ORD(A，M)。

2.一个正整数关于另一个正整数的乘法逆

一个正整数a关于另一个正整数m＞1的乘法逆等于使m整除ab-1的最小正整数b。这里a与m必须满足GCD(a,m)＝1。因此，可以采用解二元一次不定方程ax-my＝1的方法解决。

设两个正整数字符串A和M>"1"满足GCD(A，M)＝"1"。计算A关于M乘法逆的算法如下：

(1)令X＝""，Y＝""。

(2)用求解二元一次不定方程的算法求二元一次不定方程A×X-M×Y＝"1"的一个特解(X，Y)。

(3)结束。

算法结束后，X即为所求。

3.一元一次同余方程

一元一次同余方程是指形如ax≡b(modm)的方程。其中x是未知数，m＞1称为模。它表示m整除ax-b。满足这个条件的任意一个整数都称为它的解。使m整除x₁-x₂的任何两个解x₁和x₂被视为本质上相同的。一元一次同余方程ax≡b(modm)有解的充分必要条件是d＝GCD(a,m)整除b。当这一条件满足时，所有d个本质上不同的解是0≤k＜d。其中x₀是简化的一元一次同余方程的任意一个解。由此可见，解一元一次同余方程的关键是求满足GCD(a,m)＝1的一元一次同余方程的ax≡b(modm)的一个特解。而这又等价于求解二元一次不定方程ax-my＝b。

设两个正整数字符串A和M>"1"满足GCD(A，M)＝"1"。B是一个整数字符串。求解一元一次同余方程AX≡B(mod M)的算法如下：

(1)令X＝""，Y＝""。

(2)用求解二元一次不定方程的算法求二元一次不定方程A×X-M×Y＝B的一个特解(X，Y)。

(3)结束。

算法结束后，X即为所求。

4.一元一次同余方程组

若干个一元一次同余方程联立组成的方程组

$\left\{\begin{array}{c}x\&equiv;b_1\left({\mathrm{mod}m}_1\right)\\ x\&equiv;b_2\left({\mathrm{mod}m}_2\right)\\ \mathrm{L\; L\; L\; L\; L\; L\; L}\\ x\&equiv;b_n\left(\mathrm{mod}{m}_n\right)\end{array}\right.$

称为一元一次同余方程组，其中诸m_k大于1且两两互素。若整数x₀满足上面同余方程组中的每一个同余方程，则称x≡x₀(modm)为该同余方程组的解，m＝m₁m₂L m_n。

根据中国剩余定理，上面一元一次同余方程组的解为

$x\&equiv;\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{b}_i{M}_i{N}_i\left(\mathrm{mod}m\right).$

其中M_i＝m/m_i，N_i是M_i关于m_i的一个乘法逆。

因此，求解一元一次同余方程组的问题化为计算M_i、N_i和的问题。

设有正整数字符串N表示正整数n>1，B是整数字符串数组，MK是正整数字符串数组，B和MK共有N个数组元素，MK中的所有N个正整数字符串都大于1且两两互素。求解形如

$\left\{\begin{array}{c}X\&equiv;B\left[0\right]\left(\mathrm{mod}M\right[0\left]\right)\\ X\&equiv;B\left[1\right]\left(\mathrm{mod}M\right[1\left]\right)\\ \mathrm{L\; L\; L\; L\; L\; L\; L\; L}\\ X\&equiv;B[N-1]\left(\mathrm{mod}M\right[N-1\left]\right)\end{array}\right.$

的一元一次同余方程组的算法如下：

(1)令M＝"1"，Q＝""，R＝""，T＝""，X＝""，I＝"0"，i＝0，J＝"0"，j＝0。

(2)判断：若I＝N，则转向第(10)步。

(3)令MP[i]＝"1"，NP[i]＝"1"，J＝"0"，j＝0。

(4)判断：若J＝N，则转向第(9)步。

(5)判断：若J≠I，则令MP[i]＝MP[i]×MK[j]。

(6)令J＝J+"1"，j＝j+1，转向第(4)步。

(7)令M＝M×MK[i]。

(8)求MP[i]关于MK[i]的乘法逆NP[i]。

(9)令I＝I+"1"，i＝i+1，转向第(2)步。

(10)令Q＝""，R＝""，T＝""，X＝"0"，I＝"0"，i＝0。

(11)判断：若I＝N，则转向第(16)步。

(12)令T＝B[i]×MP[i]，T＝T×NP[i]。

(13)对T和M施行Euclid除法，得到商Q和余数R。

(14)令T＝R，X＝X+R。

(15)令I＝I+"1"，i＝i+1，转向第(11)步。

(16)对X和M施行Euclid除法，得到商Q和余数R。

(17)令X＝R，结束。

算法中，MP、NP是有N个数组元素的字符串数组。

算法结束后，X即为原一元一次同余方程组的解，M的值即为m。

结果处理模块是用于对S3的计算结果数据的整数部分从后向前每隔一个固定的位数优选为10就插入一个空格，对计算结果数据的小数部分从前向后每隔一个固定的位数优选为10就插入一个空格。具体算法如下：

1.设输入计算结果的整数部分为字符串STR。下面的算法对STR从后向前每隔一个固定的位数n就插入一个空格：

(1)判断：若n<1，则结束。

(2)令R＝""，并令字符型指针变量p指向STR[0]，字符型指针变量q指向R[0]。

(3)判断：若p[0]＝'-'，则令p＝p+1，q＝q+2，q[0]＝p[0]，q[1]＝″(空格)。

(4)令len等于p所指字符串的长度，并令r＝len％n。

(5)令i＝0，j＝0。

(6)判断：若p[i]＝'\0'，则令q[j]＝'\0'，结束。

(7)令q[j]＝p[i]，i＝i+1，j＝j+1。

(8)判断：若i＝r或者i％n＝0，则令q[j]＝″(空格)，j＝j+1。

(9)转向第(6)步。

算法结束后，结果是字符串R。

2.设计算结果的小数部分为字符串STR。下面的算法对STR从前向后每隔一个固定的位数n就插入一个空格：

(1)判断：若n<1，则结束。

(2)令R＝""，i＝0，j＝0。

(3)判断：若STR[i]＝'\0'，则令R[j]＝'\0'，结束。

(4)令R[j]＝STR[i]，i＝i+1，j＝j+1。

(5)判断：若i％n＝0，则令R[j]＝″(空格)，j＝j+1。

(6)转向第(3)步。

算法结束后，结果是字符串R。

S4得到的数据长度应大于S1输入数据长度的3倍。

上列详细说明是针对本发明可行实施例的具体说明，该实施例并非用以限制本发明的专利范围，凡未脱离本发明所为的等效实施或变更，均应包含于本案的专利范围中。

Claims (8)

1.一种计算机处理大整数的算法，其特征在于，包括以下步骤：

S1：手工或从文件中输入数据到计算机中存储；

S2：通过预处理模块对S1中输入的数据进行预处理；

S3：对S2中预处理后得到的数据，通过计算模块进行计算；

S4：对S3中计算得出的结果通过结果处理模块进行处理。

2.根据权利要求1所述的计算机处理大整数的算法,其特征在于：S1输入的数据采用动态分配的或者固定长度的字符串作为大整数的存储结构；输入的数据中每一位十进制数字直接用其ASCII码存储，负整数的负号“-”写在字符串首位，非负整数不带负号。

3.根据权利要求2所述的计算机处理大整数的算法,其特征在于：在面向对象技术编程时，可将S1中输入的数据单独设计为数据成员。

4.根据权利要求1所述的计算机处理大整数的算法,其特征在于：所述预处理模块是用来删除S1输入数据中所有的空格和删除输入数据前面多余的数字“0”以及判断输入的数据是不是整数。

5.根据权利要求1所述的计算机处理大整数的算法,其特征在于：所述计算模块包括基本运算模块、乘方与开方模块、阶乘与因子模块、素数模块、最大公约数模块、特殊数模块、不定方程模块、同余式模块。

6.根据权利要求5所述的计算机处理大整数的算法,其特征在于：所述计算模块在计算时不改变S1中输入的数据。

7.根据权利要求1所述的计算机处理大整数的算法,其特征在于：所述结果处理模块是用于对S3的计算结果数据的整数部分从后向前每隔一个固定的位数就插入一个空格，对计算结果数据的小数部分从前向后每隔一个固定的位数就插入一个空格。

8.根据权利要求1-7任意一项中所述的计算机处理大整数的算法,其特征在于：所述S4得到的数据长度应大于S1输入数据长度的3倍。