高精度计算
多精度计算
许剑伟 2006-10-31
一、多（高）精度数据表示法：
用字符型数组表示一个高精度的数，以下示范数据结构，左边为数组底端（或说内存底端），下表以底端高位（或说高端低位）权方式表示多精度数，反之，也可用底端低位（或说高端高位）方式表示多精度数。
|
数组脚标 |
0 |
1 |
2 |
3 |
4 |
5 |
|
对应位权 |
个位 |
10分位 |
100分位 |
1000分位 |
10的-4 |
10的-5 |
|
数组值 |
3 |
1 |
4 |
5 |
7 |
9 |
如上表所示，该数组表示的是3.14579。10分位的位权是0.1，100分位的位权是0.01。
二、使用中国小学算术教育指明的算法进行多精度四则运算。
（一）做乘法和加减法时，从建议最底位权处开始计算，这点很重要，不然做进位时会遇到一些麻烦。做除法时，从最高位权处开始计算。
（二）关于截断误差：用数组表示一个高精度数，仍然可能存在误差，由于数组位数有限，只好对小数5位以后的数做截断处理，从而产生误差，如数1/3=0.3333…. ，用上图数组表示为0.33333，总共6位精度（实有效数只有5个）。
（三）多精度算法。
1、多精度数加法（多对多加法）：逐位加法并做进位处理(取和与进位可同步处理，以提高速度)。
自从有了10进制等有效的数制问世之后，加法问题也随之得到有效的解决。人们通过九九加法表，很快的知道两个1位数的加法结果，如果是多位数，只需逐位相加并进位即可。
|
位权 |
个位 |
10分位 |
100分位 |
1000分位 |
10的-4 |
10的-5 |
|
数组1 |
3 |
1 |
4 |
5 |
7 |
9 |
|
数组2 |
6 |
6 |
5 |
8 |
2 |
2 |
|
数组结果 |
3+6=9 |
1+6=7 |
4+5=9 |
5+8=13 |
7+2=9 |
9+2=11 |
|
进位处理 |
9 |
8 |
0 |
4 |
0 |
1 |
2、多精度数与整数的乘法（1对多乘法）：逐位相乘并做进位处理。
这种法与加法差不多，但需要乘法口决表，对电脑而言，当然不需要口决表，它使用二进制计算(0和1)，所谓的口决只是有一个“乘数是1则结果是被乘数否则就是0”，这种乘法简单得不能再简单，如果是多位数乘法，则在相加之前要考虑移位。
|
位权 |
个位 |
10分位 |
100分位 |
1000分位 |
10的-4 |
10的-5 |
|
数组 |
3 |
1 |
4 |
5 |
7 |
9 |
|
整数 |
6 |
|||||
|
数组结果 |
3*6=18 |
1*6=6 |
4*6=24 |
5*6=30 |
7*6=42 |
9*6=54 |
|
进位处理 |
8(溢出1) |
8 |
7 |
4 |
7 |
4 |
3、多精度数与整数的除法（多对1除法）：逐位除法并做余数处理(逐位求精)。
逐位求精算法可以比较完美的解决除数为整数(int类型)的问题，如果除数是一个高精度数，逐位求精算法不一定可行。
4、多精度数与多精度数的乘法：逐位相乘并做进位处理——硬乘法。
多位数与多位数的乘法要复杂很多。如果要考虑快速相乘，则会更加复杂。
这种算法是小学老师告诉我们的，那时教材里也给出这样的算法，该算法的明显好处是算法比较单，并且用笔和纸两个工具便可解决问题。然而离开了笔和纸，这种笔算与口算混合的算法的优势几乎完全失去，几乎成了最臭的乘法，远不如古代使用歌诀或珠算等方法，所以要在电脑中实现乘，需对算法稍做修改。小学算法主要问题在：每行都要作进位处理，最后取总和时还要做进位处理，再次说，该算法占用至少了占用了N*N的空间（设两个乘数的位数都是N），做个10万位的二数相乘，可能把你的内存耗尽，再三，该算法总是从最最低位开始，乘法结果的最低位往往是要舍去的，我们却花了大量时间计算(如两个有效数字为5的数相乘，结果的精度也应该是5位，不是10位)，所以强硬将算法原原本本搬给电脑计算是很不科学的。
计算时请注意，按表格中示意纵向优先做乘法计算并取和进位，计算时从高端（底位权）开始，即第1轮在右边，这样你将不需要额外的数组空间，结果可直接填回原数组的相应位置中，另外，这么做还可以很方便的做一次性纵列加法及进位，运算过程大大简化。原因很简单，第n位的值总是由乘数的第n1位和被乘数的第n2位相乘得到，且n1+n2=n+1，易得n1<=n且n2<=n，这就意未着要计算第n位的结果，需要的是n左边的数，不需要n右边的数，所以计算左边各位结果将不再需要当前位及其右边位的乘数与被乘数，当前位的空间可被利用。
这种算法的计算量为N*N/2，如果你进行全宽度完全乘法计算，其计算量为N*N。
乘法卷积
卷积过程就是交两两相乘后，同位权的取和，这个乘法过程与上述的纵列乘法计算相同
以下示例卷积过程
循环卷积(圆卷积)，循环卷积要做循环移位操作：
卷积的另一种直观方法是使用行列式，比上面的图示法更简洁一些。
上面二种方法卷积都取得相同的结果，必须在做圆卷积时a、b后半段全部置0。有了这种乘法规律，应用数论的一些方法做乘法，可成千上成倍提高乘法速度。
5、多精度与多精度除法(多对多除法)
除法工作则比乘法更加复杂。数学家们想方设法简化除法，并取得了较好的效果，使用现行大众除法得到了普及，这里对这种除法做进一步剖析，并把它应用在大数相除，虽然这种普及的算法效率低下，但在处理位数较少的除法仍然是可取的。
例：设有大数a和b，a=9.1457978734，b=1.98669，c=a/b，求c的值。
注意，在做除法时，如果a、b前面前有0，则应移位补齐，保证a、b的第1位大于0。
我们已假设a和b为大数，这时不能一次性准确逐步求精，原因是我们在计算法第1位除法结果时就已经遇到困难（以后各位遇到的困难与第1位的相同），a除b的第1数位应如何计算呢？先考虑[9/1]=9作为第1位除法结果，显然不对。接下来取前两位做除法，即[91/19]=4，这个结果也不一定是准确值，但肯定更接进准确值。是[914/198]=4吗？也不一定准确，但此结果会更可靠一些！然而为了准确得到第1位，总不能把b的所以位数都拿来除，b是个大数，我们除不了，这正是要解决的多对多的除法问题！现在用以下方法解决：
设a的前1位为a1，b的前1位为b1，a的前2位为a2，b的前2位为b2，b的总位数为m，a的前m位为am，c1=am/b，第1位正确结果则为[am/b]。
很明显，用a1和b1整除作为结果误差太大，先做试探cx=a2/b2，cx与准确值c1=am/b会差多少呢？
三、分治法快速大数相乘：
1、分治法乘法减少乘法次数的原理
硬乘法算法简单，计算量超大，如果每个被乘数有N，，要做N*N次乘法，N次加法。下文价绍的分治乘法（也称二分乘法）可大大提高乘法速度。
设有大数A和B，长度为N。把大数分成二段，前后各段长均为N/2，A前段为a后段为b，B前段为c后段为d。为了表述方便，设前段的数权为k。
如此算法，只进行了次3次加法和3次乘法，而硬法则是4次乘法和2次加法。既然对分后可减少乘法次数，我们可对这3个乘法再对半分，计算9个乘法，如此层层对分计算，以获得更好的速度。这种分段处理的思想属于分治算法一族（古人云分而治之）。
2、分治法乘法的计算量
四、快速傅里叶变换实现大数相乘
阅读本节有一定的难度，因为这种变换本是一种神秘而深奥的数学。我尽力将这种方法讲清楚，但我知道可能难达到这个目标，即便把它写成一本书与未必能够写完整，所以能够让想使用该算法的人们掌握本算法实现过程，就算完成任务吧。
快速傅里叶变换简称FFT。
（一）离散傅里叶变换的形式：
按传统的方法是从离散转向连续，即使用极限的一些方法把离散傅里叶转为连续傅里变换。考虑到很多人可
对于时间连续信号，可利用傅里叶变换获得其频谱函数；或由其频谱函数通过反变换得到原时间函数。用公式表示为：
在离散信号处理中，应将傅里叶变换的积分形式改变为离散傅里叶变换的求和形式，把连续傅里叶变换的积分区间化成离散傅里叶变换的求和区间。
对时域的有限区间0～L内的信号x(t)按等时间间隔ΔT 抽样，t=nΔT，得到N=L/ΔT个抽样值x(nΔT)，n=0,1,2,…,N-1。
在频域的有限带宽0～F内的频谱函数x(jω)按等频域间隔Δf＝1/L抽样，即Δf=1/NΔT=F/N，或Δω=2π/NΔT，ω= k2π/NΔT，得到N个频率抽样值X(jkΔω)=X(jk2πΔf)，其中k=0，1，2，…，N-1。
通过以上变换后得ωt=(k2π/NΔT)(nΔT)= 2πnk/N。
令ΔT为1，得Δω=2π/N，ω=2π/NΔT =k2π/N，t=nΔT=n
时间序列；并将X(jk2πΔf)记作X(k)，称为频谱序列。
做离散变换：积分变换为取和,dt变为ΔT,dω变为Δω，代入后得到
由上式可知,当给定时域信号序列x(n)的长度为N时,计算频域上的一个抽样值X(k)，就需要进行N次复数乘法运算和N-1次复数加法运算;要得到X(k)的N个抽样值,就需要进行N2次复数乘法运算和N(N-1)次复数加法运算。
从物理意义上看，非周期时间信号的频谱是连续的和非周期的。时间信号抽样之后成为离散时间信号，它的频谱就变为周期性的连续谱。对频谱函数进行抽样，则对应的时间函数就变为周期性的连续信号。同时对时间信号和相应的频谱函数进行抽样，则得到离散的和周期的时间信号函数和频谱函数，这样就构成了上述离散傅里叶变换对。
（二）基2的快速傅里叶变换：
计算离散傅里叶变换的一种快速算法,简称FFT。快速傅里叶变换是1965年由美国人J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少,特别是被变换的抽样点数N越多,FFT算法计算量的节省就越显著。
计算离散傅里叶变换的快速方法，有按时间抽取的FFT算法和按频率抽取的FFT算法。前者是将时域信号序列按偶奇分排,后者是将频域信号序列按偶奇分排。利用Wn的以下特性可大大减少计算量。
1、按时间抽取算法
令信号序列的长度为N＝2的M次方,其中M是正整数，可以将时域信号序列x(n)分解成两部分，一是偶数部分x（2n），另一是奇数部分x（2n+1），其中n=0,1,2…N/2-1。于是信号序列x(n)的离散傅里叶变换可以用两个 N/2抽样点的离散傅里叶变换来表示和计算。考虑到和离散傅里叶变换的周期性，离散傅里叶变换可写成:
这种方法本质上也上分治算法，事实上算到这一步，已经将一个序列对分为2个小序列计算，计算量减少了，可是G（k）和H(k)的取和范围变窄了一半，直接通过G(k)和H(k)代入上式计算，只能得到一半的X(k)序列。为了得到更一半序列，还得作进一步变换。
到此我们看到了什么？熟悉电脑算法的人们很快就会意识到，我们将大幅度减少计算！求X(k)前半段和后半段本来是长度为N的序列变换，已转为长度2个长度为N/2的变换，N/2的变换的计算量为(N/2)的平方，那么X(N)的计算理为N平方/2，计算量减少了一半！如果对G(k)、H(k)也对分计算，总计算量将再次变小，按上一节关于分治算法计算量的计算，层层对分后，计算量将减少为NlogN(2为底数)
这个计算过程是递归过程序。可得X(k)，需求得偶序列变换H(k)各奇G(k)，H(k)、G(k)再分奇偶…，最后应如何x(n)序列应如何放置才便于求解？把X(k)序列写成二进制形式，
x(000),x(001),…,x(111),以N=8为例
一层变换要求奇偶分开
偶x(000),x(010),x(100),x(110)…
奇x(001),x(011),x(101),x(111)…
即把k末位是0的放一组,是1的放一组
如做二层变换
对以上2序列整除2后得以为1步进单位的序列,对同样按末位0、1分组，在2进制中无需除2，除2就是右移1位，所以直接看倒数第2位即可。
第三层及以后同上
最后得到的二进序列顺序的规规律为末位按0和1分2组，末第2位按0和1分4组，末第3位按0和1分8组…正好与原序列反过来,原序列是第1位分2组,第2位分4组…显然重置序列顺序为2进制的倒序
(0)000→000(0)
(1)001→100(4)
(2)010→010(2)
(3)011→110(6)
(4)100→001(1)
(5)101→101(5)
(6)110→011(3)
(7)111→111(7)
以下N=8时的信息号流程简图
可见，这一变换要求x(k)倒序
第二种时域抽取
要求X(k)的奇部分和偶部分，只需求前半段的变换X1和后半段X2，同样要求X1、X2，也只需再往左求前后半段。层层递归，最后需取x(k)的前半段与后半段。
这一变换不要求x(k)倒序，但得到的X(k)是倒序的
2、按频率抽取(逆变换)
3、某些性质
离散傅里叶变换具有连续傅里叶变换的性质：
离散傅里叶变换除有周期性之外，还具有一般线性变换的性质。
①线性：若组合信号为几个时域信号之和，其离散傅里叶变换等于各个信号的离散傅里叶变换之和。
②选择性：离散傅里叶变换的算法可以等效为一个线性系统的作用。频域变换值X(k)代表不同频率的谱线输出，这意味着离散傅里叶变换算法对频率具有选择性。
③循环移位性：有限长度的序列x(n)可以扩展为周期序列B(n),而x(n)可以看作是周期序列中主值区间内的主值序列,它的各个抽样序列好像放在一个N等分的圆周上，序列的移位就相当于它在圆周上旋转，由此可依次重复地看到周期序列x(n)。这种序列的移位称为循环移位，或圆周移位。这种性质对计算循环卷积和循环相关很有用。
④其他：如序列的离散傅里叶变换对称性和循环卷积性（即圆周卷积性）等。
离散傅里叶变换
离散傅里叶变换有与傅里叶变换相类似的作用和性质，在离散信号分析和数字系统综合中占有极其重要的地位。它不仅建立了离散时域与离散频域之间的联系，而且由于它存在周期性，还兼有连续时域中傅里叶级数的作用，与离散傅里叶级数有着密切联系。在计算速度方面，FFT快速计算的算法，使离散傅里叶变换的应用更为普遍，在实现各种数字信号处理系统中起着核心的作用。例如，通过计算信号序列的离散傅里叶变换可以直接分析它的数字频谱等。
关于卷积
电压、电流时域采样序列u(n)和i(n)，功率p(n)则为u(n)点积i(n)
则P(k)=U(k)卷积I(k)，因为
关于对称
关于乘法
五、牛顿迭代法(切线法)
当重复以上迭代，则每迭代一次将得到双倍精度
减少求开方程中不断做多对多大数相除的方法：用例3的方法求开方的倒数，最后再将结果用例1方法做一次倒数计算即可。
用级数计算圆周率
一、用最简单的公式：
以上列出了几个PI的公式，公式一比较简单，收敛速度也比较快。
//------------------------------
//----使用公式一计算
//------------------------------
#include<stdio.h>
#include<conio.h>
void add(char *a,char *b,int wn){//算a+b,结果放在a
for(int i=wn-1,y=0;i>=0;i--){
a[i]+=b[i]+y;
if(a[i]>=10) a[i]-=10,y=1; else y=0;
}
}
void mul_div(char *a,int b,int b2,int wn){ //乘除,算a*(b/b2),结果放在a
for(int i=wn-1,y=0,c;i>=0;i--){//乘b
c=a[i]*b+y;
y=c/10;
a[i]=c%10;
}
for(int i=0,y=0,c;i<wn;i++){//除b2
c=a[i]+y*10;
a[i]=c/b2;
y=c%b2;
}
}
void main(int argc, char* argv[])
{ int i,N=8*1024,N2=3.4*N;
char *a=new char[N*3],*b=a+N;
for(int i=0;i<N;i++) a[i]=b[i]=0;
for(i=1,a[0]=b[0]=2;i<N2;i++){
mul_div(b,i,2*i+1,N);
add(a,b,N);
}
for(i=0;i<N;i++){
printf("%d",a[i]);
if(i%50==0) printf(" : %d\r\n",i);
}
delete[]a;
getch();
}
