(2019年2月19日注：這篇文章原先發在自己github那邊的博客，時間是2017年2月5日)
這道題源自數學實驗上面的一組實驗，當時困擾了我特別久，題目的內容是用matlab求出π的連分數展開及每層迭代的值。
因為matlab的數值精度的問題，當你運行3+16450/16421時，就算你設置了format rat，也沒有辦法顯示准確的有理分數的值。原書是用mathmatica做的，直接一句命令就出來了。
其實matlab里面也是有相應的連分數展開的命令的，詳細可以查看以下：
但是顯然精度高時，就會遇到那天敘森寫的學習筆記的另外一面，精度要求太高，函數無法滿足要求，舉個例子。
>> 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秒鍾左右。