PI值計算

HDU 2179 pi值計算 

先發上大數版本的程式（java水的，不想寫高精度了。。）

import java.math.BigDecimal;import java.math.BigInteger;import java.util.Scanner;public class Main {public static void main(String[] args) {BigDecimal TWO=BigDecimal.valueOf(2);BigDecimal ans=BigDecimal.valueOf(2);for(int i=5000;i>=1;i--){ans=TWO.add(ans.multiply(td(i)).divide(td(2*i+1),1800,BigDecimal.ROUND_HALF_UP));}String stra=ans.toString();Scanner sc=new Scanner(System.in);while(sc.hasNext()){int n=sc.nextInt();if(n==0)break;System.out.println("3.");for(int i=0;i<n;i++){System.out.print(" ");System.out.print(stra.substring(i*5+2,i*5+7));if((i+1)%10==0)System.out.println();}if(n%10!=0)System.out.println();}}public static BigDecimal td(int n){return BigDecimal.valueOf(n);}}

公式是 PI=2+（1/3×（2+2/5×（2+...）））大概要迭代到5000層左右才能精確到1500位

不知道為什麼之前用冪級數展開一直不夠精度，但這個公式就可以了。。

下面是從網上摘抄的一段關於計算PI的演算法，同樣的公式，但沒有用高精度，代碼真的很難理解。。

#include <iostream.h>long a=10000,b,c=2800,d,e,f[2801],g;void main(){for(;b-c;)f[b++]=a/5;for(;d=0,g=c*2;c-=14,cout<<e+d/a,e=d%a)for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}一、來源程式本文分析下面這個很流行的計算PI的小程式。下面這個程式初看起來似乎摸不到頭腦，不過不用擔心，當你讀完本文的時候就能夠基本讀懂它了。程式一：很牛的計算Pi的程式int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5;for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}二、數學公式數學家們研究了數不清的方法來計算PI,這個程式所用的公式如下：        1         2         3                  kpi=2+(----- *(2+----- *(2+----- *(2+ ... *(2+----- *(2+...))...)))      2*1+1     2*2+1     2*3+1              2*k+1至於這個公式為什麼能夠計算出PI，已經超出了本文的能力範圍。下面要做的事情就是要分析清楚程式是如何?這個公式的。我們先來驗證一下這個公式：程式二：Pi公式驗證程式#include "stdio.h"void main(){        float pi=2;        int i;        for(i=100;i>=1;i--)                pi=pi*(float)i/(2*i+1)+2;        printf("%f\n",pi);        getchar();}上面這個程式的結果是3.141593。三、程式展開在正式剖析器之前，我們需要對程式一進行一下展開。我們可以看出程式一都是使用for迴圈來完成計算的，這樣做雖然可以使得程式短小，但是卻很難讀懂。根據for迴圈的運行順序，我們可以把它展開為如下while迴圈的程式：程式三：for轉換為while之後的程式int a=10000,b,c=2800,d,e,f[2801],g;main() {        int i;        for(i=0;i<c;i++)                f[i]=a/5;        while(c!=0)        {                d=0;                g=c*2;                b=c;                while(1)                {                        d=d+f[b]*a;                        g--;                        f[b]=d%g;                        d=d/g;                        g--;                        b--;                        if(b==0) break;                        d=d*b;                }                c=c-14;                printf("%.4d",e+d/a);                e=d%a;        }}註：for([1];[2];[3]) {[4];}的運行順序是[1],[2],[4],[3]。如果有逗號操作符，例如：d=0,g=c*2，則先運行d=0,然後運行g=c*2,並且最終的結果是最後一個運算式的值，也就是這裡的c*2。下面我們就針對展開後的程式來分析。四、程式分析要想計算出無限精度的PI，我們需要上述的迭代公式運行無數次，並且其中每個分數也是完全精確的，這在電腦中自然是無法實現的。那麼基本實現思想就是迭代足夠多次，並且每個分數也足夠精確，這樣就能夠計算出PI的前n位來。上面這個程式計算800位，迭代公式一共迭代2800次。int a=10000,b,c=2800,d,e,f[2801],g;這句話中的2800就是迭代次數。由於float或者double的精度遠遠不夠，因此程式中使用整數類型（實際是長整型），分段運算（每次計算4位）。我們可以看到輸出語句 printf("%.4d",e+d/a); 其中%.4就是把計算出來的4位輸出，我們看到c每次減少14（ c=c-14;），而c的初始大小為2800，因此一共就分了200段運算，並且每次輸出4位，所以一共輸出了800位。由於使用整型數運算，因此有必要乘上一個係數，在這個程式中係數為1000，也就是說，公式如下：               1          2          3                    k1000*pi = 2k+ --- * (2k+ --- * (2k+ --- * (2k+ ... (2k+ ---- * (2k+ ... ))...)))               3          5          7                  2k+1這裡的2k表示2000，也就是f[2801]數組初始化以後的資料,a=10000,a/5=2000,所以下面的程式把f中的每個元素都賦值為2000：for(i=0;i<c;i++)        f[i]=a/5;你可能會覺得奇怪，為什麼這裡要把一個常數儲存到數組中去，請繼續往下看。我們先來跟蹤一下程式的運行：while(c!=0)             //假設這是第一次運行，c=2800，為迭代次數{        d=0;        g=c*2;          //這裡的g是用來做k/(2k+1)中的分子                b=c;    //這裡的b是用來做k/(2k+1)中的分子                while(1)                {                        d=d+f[b]*a; //f中的所有的值都為2000，這裡在計算時又把係數擴大了a=10000倍。                                                //這樣做的目的稍候介紹，你可以看到輸出的時候是d/a,所以這不影                                                //計算                                g--;                        f[b]=d%g;       //先不管這一行                                d=d/g;  //第一次啟動並執行g為2*2799+1，你可以看到g做了分母                                g--;                        b--;                        if(b==0) break;                        d=d*b;          //這裡的b為2799，可以看到d做了分子。                }                c=c-14;                printf("%.4d",e+d/a);                e=d%a;}只需要粗略的看看上面的程式，我們就大概知道它的確是使用的那個迭代公式來計算Pi的了，不過不知道到現在為止你是否明白了f數組的用處。如果沒有明白，請繼續閱讀。d=d/g,這一行的目的是除以2k+1，我們知道之所以程式無法精確計算的原因就是這個除法。即使用浮點數，答案也是不夠精確的，因此直接用來計算800位的Pi是不可能的。那麼不精確的成分在哪裡？很明顯：就是那個餘數d%g。程式用f數組把這個誤差儲存起來，在下次計算的時候使用。現在你也應該知道為什麼d=d+f[b]*a;中間需要乘上a了吧。把分子擴大之後，才好把誤差精確的算出來。d如果不乘10000這個係數，則其值為2000，那麼運行d=d/g；則是2000/(2*2799+1)，這種整數的除法答案為0，根本無法迭代下去了。現在我們知道程式就是把餘數儲存起來，作為下次迭代的時候的參數，那麼為什麼這麼做就可以使得下次迭代出來的結果為接下來的數字呢？這實際上和我們在紙上作除法很類似：0142/——------7 / 1107---------------3028---------------2014---------------60.....我們可以發現，在做除法的時候，我們通常把餘數擴大之後再來計算，f中既然儲存的是餘數，而f[b]*a;則正好把這個餘數擴大了a倍，然後如此迴圈下去，可以計算到任意精度。這裡要說明的是，事實上每次計算出來的d並不一定只有4位元，例如第一次計算的時候，d的值為31415926，輸出4位時候，把低四位的值儲存在e中間，e=d%a，也就是5926。最後，這個c=c-14不太好理解。事實上沒有這條語句，程式計算出來的仍然正確。只是因為如果迭代2800次，無論分數如何精確，最後Pi的精度只能夠達到800。你可以把程式改為如下形式嘗試一下：for(i=0;i<800;i++){        d=0;        g=c*2;        b=c;        while(1)        {                d=d+f[b]*a;                g--;                f[b]=d%g;                d=d/g;                g--;                b--;                if(b==0) break;                d=d*b;        }        // c=c-14; //不要這句話。        printf("%.4d",e+d/a);        e=d%a;}最後的答案仍然正確。不過我們可以看到內迴圈的次數是c次，也就是說每次迭代計算c次。而每次計算後續位元的時候，迭代次數減少14，而不影響精度。為什麼會這樣，我沒有研究。另外最後的e+d/a,和e=d/a的作用就由讀者自己考慮吧。