连分数近似算法

　　(2019年2月19日注：这篇文章原先发在自己github那边的博客，时间是2017年2月5日)

　　这道题源自数学实验上面的一组实验，当时困扰了我特别久，题目的内容是用matlab求出π的连分数展开及每层迭代的值。

　　因为matlab的数值精度的问题，当你运行3+16450/16421时，就算你设置了format rat，也没有办法显示准确的有理分数的值。原书是用mathmatica做的，直接一句命令就出来了。

　　其实matlab里面也是有相应的连分数展开的命令的，详细可以查看以下：

Continued fraction expansion

　　但是显然精度高时，就会遇到那天叙森写的学习笔记的另外一面，精度要求太高，函数无法满足要求，举个例子。

>> rat(pi,1e-31)
ans =
3 + 1/(7 + 1/(16 + 1/(-294 + 1/(3 + 1/(-4 + 1/(5 + 1/(-15)))))))
>> rat(pi,1e-18)
ans =
3 + 1/(7 + 1/(16 + 1/(-294 + 1/(3 + 1/(-4 + 1/(5 + 1/(-15)))))))

　　其实后来我测试自己的算法也是有问题的，但是至少比之前迭代不到10次就直接误差为0了要好得多，究其原因还是因为前面运算的时候用的是double型数，这里我尝试过，如果在一开始就使用符号运算的话，不仅运算时间无法忍受，而且很容易迭代到20次以后就直接程序堆栈溢出或者NaN或者是Inf了。

　　好了，说了这么多，该介绍下连分数是什么了，其实也不难，据说拉马努金对这个非常的了解。$$S_{n} = frac{1}{a_{1}+frac{1}{a_{2}+frac{1}{a_{3}+…}}}$$

　　为了算法的简便，我们设定的是这组数列全部都是正整数。对于有理数，我们知道，这组数列一定是有限数列，对于无理数，这组数列大部分情况下是无限数列。我们可以从某一个数开始，设定它为0，按照这样子的顺序，加出来一个有理分数，这个有理分数就是我们需要的值，举个例子。$$pi = frac{1}{3+frac{1}{7+frac{1}{15+…}}}$$

　　算法的大致步骤如下：

　　　　(1) 根据输入的Number值，取不足近似整数(matlab里面就是fix函数)，比如π就是取3，同时将3存入数组P的第一个值$P[1]=fix(Number)$(matlab的数组从1开始的)

　　　　(2) 做运算$π-3$，作为数组A的第一个值$A[1]$

　　　　(3) $P[2]=frac{1}{A[1]}$

　　　　(4) $A[2]=fix(P[2])$

　　下一个P和A重复上述过程

function result=test(Number,n)
%Number=(sqrt(5)-1)/2;	%测试用
%n=100;					%测试用
result=0;
P=zeros(1,n);A=zeros(1,n);
P(1)=fix(Number);
A(1)=Number-fix(Number);
for i=2:n-1
    P(i)=1/A(i-1);
    A(i)=P(i)-fix(P(i));
end
P=sym(fix(P));    %转成符号计算
for j=n-1:-1:2
    result=1/(P(j)+result);
end
result=result+P(1);

　　测试以后发现，这组算法还是有问题的，前面也说过，最大问题就是数组A是一个double型的数组，把精度截断了。我尝试过用sym型来写，但是会在某个数值的时候因为大数相乘(分子分母比较大的时候，通分造成的原因)直接超过1e308变成inf，然后就NaN了。综合上面，只能放弃掉一些精度。经过测试，直接sym型的时候连分数能写到20层，如果将A进行double化，目前没发现层数的问题，但是有效的层数，100层的话大约在70层左右，vpa以后精度在1e-15左右，还是可以的。计算处理速度上基本没有大的问题，100层的时间大约是2秒钟左右，300层的时间大约是3秒钟左右。