LISP Programming 1062_ex5

http://erdos.csie.ncnu.edu.tw/~klim/scheme/lisp-1062.html

exercise 5: calculating continued-fraction
•Refer to exercise 1.37.(參考這篇：黃金分割)
•You have to write two procedures. One is (cont-frac n d k) and another is (pi epsilon).
•The "cont-frac" should be implemented in linear iterative way.
•To calculate pi, the sequence N_i is: 4, 1^2, 2^2, 3^2, ...
 and the sequence D_i is: 1, 3, 5, 7, ...
•To approximate pi, calculate (cont-frac n d 1), (cont-fract n d 2), .. until two consecutive values
 are  close enough, that is, the distance between them is less than epsilon.

$ scheme48
> ,load ex5.scm
> (cont-frac (lambda (i) 1.0) (lambda (i) 1.0) 100)
0.6180339887498948
> (pi 0.001)
3.141463414634146
> (pi 0.0001)
3.1415888250921244
> (pi 0.0000000000001)
3.141592653589791
> ,exit
$
解答說明請參考：【Scheme】連分法求pi

解答:

;continued fraction calc pi
(define (cont-frac n d k)
   (define (frac-iter i result)
       (if (= i 0)
           result
           (frac-iter (- i 1.0) (/ (n i) (+ (d i) result)))))
   (frac-iter (- k 1.0) (/ (n k) (d k))))
(define (calc-pi k)
    (define (d i)
        (- (* 2 i) 1))
    (define (n i)
        (if (= i 1)
            4
            (* (- i 1) (- i 1))))
    (cont-frac n d k))
(define (calc-enough index accuracy)
    (define next_i index) ;next_i=index
    (if (< (abs (- (calc-pi index) (calc-pi (+ next_i 1)))) accuracy)
        (calc-pi (+ next_i 1))               ;true
        (calc-enough (+ index 1) accuracy))) ;false
(define (pi accuracy)
    (calc-enough 1 accuracy))

執行結果：

【Scheme】連分法求Pi

題目要求:
•Refer to exercise 1.37.
•You have to write two procedures. One is (cont-frac n d k) and another is (pi epsilon).
•The "cont-frac" should be implemented in linear iterative way.
•To calculate pi, the sequence N_i is: 4, 1^2, 2^2, 3^2, ...
 and the sequence D_i is: 1, 3, 5, 7, ...
•To approximate pi, calculate (cont-frac n d 1), (cont-fract n d 2), .. until two consecutive values
 are  close enough, that is, the distance between them is less than epsilon.

參考SICP ex 1.37，寫兩個程序
1.  (cont-frac n d k) and (pi epsilon)
2. "cont-frac" 使用迭代法
3. 計算 Pi，序列N_i 是4，1^2，2^2，3^2 . . .D_i是 1，3，5，7 . . .
4. 近似求Pi，直到(cont-frac n d 1), (cont-fract n d 2)值小於精確值

首先 π 的連分公式如下:

其中4,1,4,9,16就是N_i 項; 而 1,3,5,7,9就是Ｄ_i項,
參考上一篇的(cont-frac)，迭代法程式:

(define (cont-frac-iter n d k)
   (define (frac-iter i result)
       (if (= i 0)
           result
           (frac-iter (- i 1) (/ (n i) (+ (d i) result)))))
   (frac-iter (- k 1) (/ (n k) (d k))))

接著解釋 n項 與 d項:
這裡要注意的是 
d項就是迭代 i值的(2*i-1)，例如第一項是(2*1-1)=1，第二項是(2*2-1)=3 . . .依此類推
而n項就是 迭代 i值的平方(i^2)，但要注意第一項是4，第二項是((2-1)^2)，第三項是((3-1)^2) . . .

(define (calc-pi k)
    (define (d i)
        (- (* 2 i) 1)) ;((2*i)-1)
    (define (n i)
        (if (= i 1)
            4
            (* (- i 1) (- i 1))))
    (cont-frac n d k))

這樣就可以求Pi，例如(calc-pi 10)，求出聯分法10項的Pi值: 
(calc-pi 10)
3.14159254044654 
(calc-pi 20)
3.14159265358979
用這種方式求Pi因為收斂很快所以不會花費太多時間。

接著要給精確值求出Pi，就是用迭代法每一次求出的Pi值與下一個項次比較，直到小於精確值為止
如:
(calc-pi 5)-(calc 6) < accuracy ?
如果是則停止
如果不是則(calc-pi 6)-(calc-pi 7)

(define (inc x)
    (+ x 1))
(define (enough-pi index accuracy)
    (define next_i index) ;next_i=index
    (if (< (abs (- (calc-pi index) (calc-pi (inc next_i)))) accuracy)
        (calc-pi (inc next_i))            ;true
        (enough-pi (inc index) accuracy)) ;false
)
(define (pi accuracy)
    (enough-pi 1 accuracy))

輸入與結果:
(pi 0.001)
3.14146341463415
(pi 0.0001)
3.14158882509212
(pi 0.0000000000001)
3.14159265358979
如果將各項的輸出印出

index:5 next_i:5 err:0.001
3.14146341463415
在第5項的時候收斂到0.001

index:7 next_i:7 err:0.0001
3.14158882509212
在第7項的時候收斂到0.0001

index:19 next_i:19 err:1e-13
3.14159265358979
在第19項的時候收斂到0.0000000000001


參考:
1.維基百科圓周率
2. SICP Ex 1.37

【Scheme】SICP Ex.1.37

Exercise 1.37.  a. An infinite continued fraction is an expression of the form
在SICP練習1.37 中要求寫出連分法的表示，公式如下:

而這個除法式可以用一個遞迴表達式來表示如下：

其中函數 cf(i) 表示連分式的第 i 個項。

因此Scheme程式遞迴寫法如下:

(define (cont-frac n d k)
  (define (frac i)
     (if (< i k)
         (/ (n i) (+ (d i) (frac (+ i 1))))
         (/ (n i) (d i))))
  (frac 1))

迭代寫法如下:

(define (cont-frac-iter n d k)
   (define (frac-iter i result)
       (if (= i 0)
           result
           (frac-iter (- i 1) (/ (n i) (+ (d i) result)))))
   (frac-iter (- k 1) (/ (n k) (d k))))

可以用下面的近似驗證:

> (cont-frac (lambda (i) 1.0) (lambda (i) 1.0) 100)0.618033988749895

這就是黃金分割比例值。


參考:
1. SCIP 教材