为了描述其算法原理，请我们首先以二项式定理开始。二项式定理给出两个数之和的整数次幂诸如(x + y)ⁿ 展开为类似ax^by^c 项之和的恒等式，公式如下，

                                                                                                        (1)

其中   称为二项式系数，亦可记作  ,见文献[1]。特别的，当x=1，公式(1)可写成如下形式，

                         (2)

值得注意的是，不但a为正整数时，上述的公式成立，当为实数甚至复数时，上面公式亦同样成立，见文献[2]。当a=1/2，公式2 变为

上面 (2K-3)!!和2k!!中的“!!”符号表示双阶乘，当n为奇数时，n!!=n ×(n-2)×(n-4)×* …1, 当n为偶数时，n!!=n×(n-2)×(n-4)× …2, 当n=0 或者-1时，n!!=1，见文献[3]。按照定义有 ，上面推导出来的公式可最终表示为：                                              (3)

其x的k次项的系数的递推公式为： ，

在文献[4] 给出如下公式，它实际上和(3)是一致的，不过(3)更简单些。

在文献[5] 和 文献[6]， 给出更多的这个级数的系数，他们是1, 1/2, -1/8, 1/16, -5/128, 7/256, -21/1024, 33/2048, -429/32768, 715/65536, -2431/262144, 4199/524288, -29393/4194304, 52003/8388608。

替换公式3中的“+x”为”-x”,可以发现，多项式中的所有系数都是负的，公式如下

                                                                                                   (4)

本文中的算法就是基于(4)，其主要方法如下

1. 找出x 的平方根的一个近似值 base，使得base稍大于x的平方根, base可表示为Q/P的形式。依下面的步骤，可推出N的平方根的一种表达式

                                                                                                        (5)

例如：

2. 若 ，对根式按照(4)展开，逐项求和，可得出一个稍小于1的数，然后乘以base，就得到最终的结果。

下面分别介绍步骤1和步骤2的算法

步骤1的目标是求出N的平方根的近似值的一个分数表示，记作Q/P，为了方便计算级数的和，Q² 不能大于计算机内置的整数类型，对于32位系统，要求Q²不能大于2³²-1。求1个浮点数的分数表示采用如下方法进行，逐步求精的迭代过程。

1.1 在迭代之初，求出N的平方根为s，s用c数学库函数sqrt来计算，可精确到15-16位10进制数字，然后计算s的一个下界low和上界high。1.1到1.6的代码，请参照下文C代码中的函数findFraction。

1.2 计算mid,取low和high的分母之和作为mid的分母，取low和high的分子之和作为mid的分子，容易证明，low<mid<high.

1.3 比较mid的分子，若其平方超过内置整数的最大值，则转步骤1.5.

1.4 比较s和mid的值，若mid>s, 用mid的值替换high原有的值，记作high←mid，否则low←mid

1.5 转1.2，继续下一次迭代。

1.6 base←high

容易看出，步骤1.2到1.5是一个迭代的过程，整个过程类似2分查找。在整个迭代过程中，区间high-low越来越小，high也就越来越近似于N的平方根

得到base=P/Q的值，利用(5)就可求出rem的值。

2. 计算 的值，可使用公式(4)来求，rem具有这样的特点，分子很小而分母较大，容易知道，rem越小，级数收敛的就速度越快。在计算过程中使用递推法。用res表示级数的和，暂时不考虑其正负，用item表示级数第i项的值。重复如下步骤，直到item足够小,其和便达到指定精度。下面的代码中rem.d表示rem的分母，rem.n表示rem的分子.

2.1 item←1/2

2.2 item=item*rem, 即item *=rem.n, item/=rem.d

2.3 res += item

2.4 i←2

2.5 item *= rem.n, item/=rem.d

2.6 item *= (2i-3), item /= 2i;

2.7 res+=item

2.8 检查item的值，如果item足够小，转2.11

2.9 i++

2.10 转2.5，继续求下一项item和item的和res

2.11 将res乘以-1，这里采用求补的方法，类似于对数表示法，将小数部分求补，整数部分减1，如-0.1235 的补数表示为-1+(1-0.1235)=（-1）+0.8765。

2.12 到此为止，res等于级数中除了第一项以外各项的和，即 。

2.13 res+=1, 其整数部分从-1变为0

2.14 res*=base.

至此，N的平方根就计算完毕。

下面的内容叙述任意精度数的表示法和其四则运算

本文中的C程序用动态数组表示一个任意精度的定点小数，采用1000000000进制，数组元素的类型为C语言的内置整数类型unsigned long，数组的每个元素表示9位10进制数。结构体LONG_FLOAT定义了任意精度数的表示，data指针表示数组的地址。任意精度数采用定点格式，整数部分用数组data的第一个元素data[0]表示，小数部分用数组的其它元素表示。这样的表示刚好满足需要，因为小整数的平方根的整数部分恰好可用1个元素来表示。结构体LONG_FLOAT的成员len表示数组的长度，成员first_idx表示数组中第一个非0元素的偏移地址。在计算过程中，item的值会越来越小，其数组前部会出现连续的0，为了减少计算量，仅从第一个非0元素算起。first_idx>=len-1视为任意精度数的值为0。为保证精度，在创建数组时，实际分配的数组的长度要比所需的精度多2个单元。任意精度数的存储顺序和显示的顺序一致。在计算乘法和加法时，采用从右至左的顺序计算，在计算除法时，采用从左至右的顺序计算。

本文中所需的任意精度计算仅仅涉及以下几种种运算。

1. 两个任意精度定点数相加，函数LF_Add实现了两个同样长度的任意精度定点数的加法运算。

2.任意精度定点数乘以一个单精度整数，单精度整数的类型定义为unsigned long，函数LF_Mul用于计算任意精度定点数以一个单精度整数。

3.任意精度定点数除以一个单精度整数，函数LF_Div用于实现这一功能。

4.求补运算，即计算-1乘以任意精度定点数的积，结果采用补码表示，函数LF_NEG用于实现这一功能。

5.任意精度定点数的创建和初始化，函数LF_Init用于创建动态数组，并对各个元素清0. 其他函数包括GCD32等，函数GCD32用于计算2个单精度整数的最大公约数。

下面为完整的源代码。

#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#include <math.h>
#define RADIX   1000000000
#define LOG10_RAD   9   //log10(RADIX)
#define MAX_WORD 0xffff#if _MSC_VER > 1000       //The compiler is VC
typedef unsigned __int64 UINT64;
#else                   //The compiler is GCC
typedef unsigned long long  UINT64;
#endif typedef unsigned long UINT32;
typedef struct _frac
{UINT32 n;  //numeratorUINT32 d;  //denominator
}FRAC;/*定义一个定点小数类型，数的有效数字存储在数组data[]中，每个数组元素data[i]可表示9位10进制数字。数组的长度为len,故可表示len*9位10进制数字,其中data[0]表示整数部分，data[1]到data[len-1]表示小数部分数字的存储顺序和显示顺序相同，如data[0]=2,data[1]=222212345,data[2]=111198765,则表示2.222212345111198765数组的前几个存储单元可能为0，为了节省计算量，可从数组的第一个非0元素开始算起。第一个非0数组单元的顺序号用first_idx表示，如first_idx=3,则表示data[0],data[1],data[2]全部为0
*/
typedef struct _long_float
{UINT32 *data;int len;int first_idx;
}LONG_FLOAT;// 功能: 求n1,n2的最大公约数
int GCD32(UINT32 n1, UINT32 n2) //Greatest Common Divisor
{UINT32 modNum;UINT32 bigNum=  (n1>n2)?n1:n2;UINT32 smallNum=(n1<n2)?n1:n2;if (smallNum==0)return bigNum;while (1){modNum=bigNum % smallNum;if (modNum==0)break;bigNum=smallNum;smallNum=modNum;}return smallNum;
}void FRAC_set(FRAC *frac,UINT32 n,UINT32 d)
{frac->n=n; frac->d=d;
}void LF_Init(LONG_FLOAT *pNum,int len)
{pNum->data=(UINT32 *)malloc(len*sizeof(UINT32));memset(pNum->data,0,sizeof(UINT32)*len);pNum->first_idx=0;pNum->len=len;
}void LF_Free(LONG_FLOAT *pNum)
{free(pNum->data);pNum->len=0;
}void LF_Add(LONG_FLOAT *pTarget,const LONG_FLOAT *pSrc)
{int i;UINT32 t,c;for (c=0,i=pSrc->len-1;i>=pSrc->first_idx;i--){t=pTarget->data[i] + pSrc->data[i] +c;c=0;if (t>=RADIX){ t-=RADIX; c=1;}pTarget->data[i]=t;}while (c>0 && i>=0){t=pTarget->data[i] + c;c=0;if (t>=RADIX){ t-=RADIX; c=1;}pTarget->data[i--]=t;}pTarget->first_idx=i+1;if (i<0){printf("Error in %d line\n",__LINE__);return ;}
}void LF_Mul(LONG_FLOAT *pTarget,UINT32 k)
{int i;UINT64 t,c;for (c=0,i=pTarget->len-1;i>=pTarget->first_idx;i--){t=(UINT64)(pTarget->data[i]) * (UINT64)k +c;pTarget->data[i]=(UINT32)(t % RADIX);c=(UINT32)(t / RADIX);}if (c>0){if (i<0){printf("Error in %d line\n",__LINE__);return;}pTarget->data[i]=(UINT32)c;pTarget->first_idx=i;}
}//return remainder
UINT32 LF_Div(LONG_FLOAT *pTarget,UINT32 k)
{int i;UINT64 t;UINT32 m=0;for (i=pTarget->first_idx;i<pTarget->len;i++){t=(UINT64)m * (UINT64)RADIX + pTarget->data[i];pTarget->data[i]= (UINT32)(t / k);m=(UINT32)(t % k);}if (pTarget->data[pTarget->first_idx]==0)pTarget->first_idx++;return m;
}void LF_NEG(LONG_FLOAT *pTarget)
{int i=pTarget->len-1;while (pTarget->data[i]==0)  i--;if (i>0){pTarget->data[i]=RADIX-pTarget->data[i];i--;while (i>=1){    pTarget->data[i]=(RADIX-1)-pTarget->data[i]; i--;}}pTarget->data[0]--;pTarget->first_idx=0;
}/*Find a fraction x, such that x equal to f approximatively and x>f
*/
FRAC findFraction(double f)
{FRAC low,mid,high;UINT32 c,root_n;double f1=f- floor(f);root_n=(UINT32)(floor(f));if ( f1 <= 0.5)  {  c=(UINT32)(1.0/f1);  FRAC_set(&low, 1,c+1);FRAC_set(&high,1,c);}  else  {  c=(UINT32)(1.0/(1.0-f1));  FRAC_set(&low, c-1,c);FRAC_set(&high,c,c+1);}low.n  +=  root_n * low.d;high.n +=  root_n * high.d;if ( high.n> MAX_WORD){FRAC_set(&low,  root_n,1);FRAC_set(&high, root_n+1,1);}while (1){FRAC_set(&mid,low.n+high.n,low.d+high.d);if ( mid.n  >= MAX_WORD)break;if ( (double)(mid.n)/(double)(mid.d) > f)high=mid;elselow=mid;}return high;
}void calc_sqrt(UINT32 N, LONG_FLOAT *pItem, LONG_FLOAT *pResult )
{int i;UINT32 max_UINT32,tGCD;FRAC base,rem;double f;f=sqrt((double)(N));/* assume base approximatively equal to sqrt(N) and base>sqrt(N),sqrt(N)= base * sqrt(1- rem) = base * sqrt(1-(N * base.d * base.d)/( base.n * base.n))*/base= findFraction(f);FRAC_set(&rem,(base.n * base.n)-N*base.d*base.d,base.n * base.n);tGCD=GCD32(rem.n,rem.d);if (tGCD>1) {  rem.n /= tGCD;  rem.d /= tGCD; }printf("sqrt(%u)=(%u/%u)/sqrt(1-%u/%u)\n",N, base.n,base.d,rem.n,rem.d);// item[1] = rem * 1/2// item[i] = item[i-1] * rem * (2*i-1)/(2*i+1))    // result  = 1-sum(item[i]) where i=2 up to x and item[x]<10^-d)pItem->data[1]=RADIX/2;    // calculating item[1], now item=1/2LF_Mul(pItem,rem.n);   // item[1]= item[1] * rem, now the item=1/2*remLF_Div(pItem,rem.d);   LF_Add(pResult,pItem);  //sum= sum + item[i]i=2;max_UINT32=~0;while (pItem->first_idx < pItem->len){if (  rem.n * (2*i-3) < max_UINT32){LF_Mul(pItem,rem.n*(2*i-3));    // item *= rem;LF_Div(pItem,rem.d);    // }else{LF_Mul(pItem,rem.n);   // item *=remLF_Div(pItem,rem.d);  // LF_Mul(pItem,2*i-3); // item *= (2i-3)/2i}LF_Div(pItem,2*i);        LF_Add(pResult,pItem);      //sum= sum + item[i]i++;}LF_NEG(pResult);       //pResult= - pResult;pResult->data[0] +=1;    //pResult=pResult+1;LF_Mul(pResult,base.n);       //result=result * baseLF_Div(pResult,base.d);
}//calculating sqrt(n) to d decimal digtal and print the result
void calc_print(UINT32 n,int d)
{int i,buffLen;LONG_FLOAT item,result;UINT32 r;r=(UINT32)(sqrt(n));if (r * r ==n){ printf("sqrt(%u)=%u\n",n,r);return ; }buffLen=(d+LOG10_RAD-1)/LOG10_RAD+3;LF_Init(&item,buffLen);LF_Init(&result,buffLen);calc_sqrt(n,&item,&result);printf("\t=%u.",result.data[0]);for (i=1;i<buffLen-2;i++)printf("%09u",result.data[i]);printf("\n");LF_Free(&item);LF_Free(&result);}int main(int argc, char* argv[])
{UINT32 i;UINT32 base;base=2;for (i=base;i<base+100;i++){calc_print(i,100);}return 0;
}

版权所有，转载请注明出处。

参考文献：

[1] “Binomial theorem” from Wikipedia,http://en.wikipedia.org/wiki/Binomial_theorem.

[2] “Binomial series” from Wikipedia,http://en.wikipedia.org/wiki/Binomial_series

[3] Weisstein, Eric W., "Double Factorial " from MathWorld,http://mathworld.wolfram.com/DoubleFactorial.html

[4] “Square root” From Wikipedia,http://en.wikipedia.org/wiki/Square_root