從turtle海龜動畫 學習 Python - 高中彈性課程系列 11.3 連分數演算法與轉轉相除法- 用 turtle 呈現演算法之執行動作

將與數學相關的演算法, 例如轉轉相除法與連分數演算法, 以海龜移動呈現, 讓學生可以了解數學演算法的執行過程, 算是用海龜繪圖, 實現演算法視覺化.

gcd(a,b) = gcd(b,a % b)
將 a 的大小遞降, 最後達到 最小的非零的 r, 就是 a, b 的最大公因數
此種解題思路, 有的叫 無窮遞降法, 用在很多有名的數學大定理的證明中, 例如 John Nash 的一個證明, 也是用遞歸迭代的方式去遞降一個參數, 叫做 Nash iteration method,
看起來很高深, 其實都是類似的思維.

---

“Talk is cheap. Show me the code.”
― Linus Torvalds

老子第41章
上德若谷
大白若辱
大方無隅
大器晚成
大音希聲
大象無形
道隱無名

拳打千遍, 身法自然

“There’s no shortage of remarkable ideas, what’s missing is the will to execute them.” – Seth Godin
「很棒的點子永遠不會匱乏，然而缺少的是執行點子的意志力。」—賽斯．高汀

---

致謝: 非常感謝製作 python turtle demo 的義工群!

本系列文章之連結

• 從turtle海龜動畫學習Python-高中彈性課程1 link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 2 安裝 Python, 線上執行 Python link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 3 烏龜繪圖 所需之Python基礎 link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 4 烏龜開始畫圖 link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 5 用函數封裝重複性指令-呼叫函數令烏龜畫正 n 邊形 link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 6 畫多重旋轉圓,螺旋正方形 link

• 從turtle海龜動畫 學習 Python - 7 遞歸 recursive - 高中彈性課程系列 link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 8 碎形 (分形 fractal) link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 8.1 碎形 L-system link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 9 Python 物件導向介紹 link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 9.1 Python 物件導向的練習 link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 10 藝術畫 自定義海龜形狀 link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 10.1 藝術畫 python繪製天然雪花結晶 https://blog.csdn.net/m0_47985483/article/details/122262036 link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 10.2 藝術畫 Python 製作生成式藝術 link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 11.1 氣泡排序 - 用 turtle 呈現演算法之執行動作 link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 11.2 maze 迷宮 - 用 turtle 呈現演算法之執行動作 link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 11.3 連分數演算法與轉轉相除法- 用 turtle 呈現演算法之執行動作 link

• 從turtle海龜動畫 學習 Python - 高中彈性課程系列 11.4 最短路徑 Dijkstra- 用 turtle 呈現演算法之執行動作 link

---

連分數演算法其實與轉轉相除法是同一個東西,
轉轉相除法一般是用來求最大公因數, gcd(a, b), 他是提取試除的最終步驟所獲得的非零餘數, 這就是 a,b 的最大公因數
連分數演算法則是提取每次試除的步驟所獲得的商數,

以下我們參考 "理解演算法Python初學者的深度歷險 " 這本書的說明及codes.

轉轉相除法

歐基理德輾轉相除法求 gcd(a,b), 主要基於以下這個性質或觀察:

Lemma If a= qb + r, then gcd(a,b) = gcd(b,r).

註: a: dividend 被除數, b: divisior 除數, q: quotient 商數, r: remainder 餘數.
Ref: Burton, Elementary Number Thoery, 7th, Sec 2.4.

基於此性質, 可以透過不斷的 執行

$$a=qb+r$$
$$a\leftarrow b$$
$$b\leftarrow r$$

或是濃縮寫成:
$$gcd(a,b)=gcd(b,a\%b)=gcd(a\%b,b\%(a\%b))=\cdots$$

將 a 的大小遞降, 最後達到 最小的非零的 r, 就是 a, b 的最大公因數

此種解題思路, 有的叫 無窮遞降法, 用在很多有名的數學大定理的證明中, 例如 John Nash 的一個證明, 也是用遞歸迭代的方式去遞降一個參數, 叫做 Nash iteration method,
看起來很高深, 其實都是類似的思維.

另外也可以把 Lemma 改成是相減的, 可以叫做 輾轉相減法求 gcd(m,n), 但是效能會較慢.

以下我們先寫一個用輾轉相除的程式求 gcd,
後面我們再寫一個用輾轉相減的程式求 gcd,

輾轉相除法求 gcd(a,b)

如果用Python 內建的 取餘數的指令 r = a % b 或是 r = mod(a, b)
則可以把 Lemma 改寫成:
Lemma (Python Version) If a>=b and r = a % b, then gcd(a,b) = gcd(b,r).
或是
Lemma (Python Version) If a>=b and r = mod(a, b), then gcd(a,b) = gcd(b,r).

或是更直接
Lemma (Python Version) If a>=b , then gcd(a,b) = gcd(b,a % b).
Lemma (Python Version) If a>=b , then gcd(a,b) = gcd(b, mod(a, b)).

# By Prof. P-J Lai MATH NKNU 2020/2/13
def gcd_Euclid_abr(m,n):
    ''' Find gcd(m.n) with division 用"歐基理德輾轉相除法"求 gcd(m,n)'''

    a = max(m, n)
    b = min(m, n)
    q = a // b
    r = a % b
    
    while r >0:
        a, b = b, r
        r = a % b

    return b

print("By gcd_Euclid_abr")
print(gcd_Euclid_abr(4,2))
print(gcd_Euclid_abr(7,3))
print(gcd_Euclid_abr(5,20))
print(gcd_Euclid_abr(72511,1229))
# 楊克昌趣味數學編程及拓展, P.59: 72511 = 59 * 1229
#>>> #gcd_Euclid(7787340436465476554757657804244,7341600567657341#2345766867345250282)
#34
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

輾轉相減法求 gcd(a,b)

用"相減法"時,
Be careful!
假設 a >= b 且 d = a - b,
It maybe b > d or b < d. (相減之後的 d 不一定比 b 小!)
所以每次都要用 a, b = max(b, d), min(b, d) 重新取大與小,

如果用"歐基理德輾轉相除法", 則 b > r (商數一定比餘數大), 則不須max, min 的動作

# gcd_substract
# By Prof. P-J Lai MATH NKNU 2020/2/13
# 用"相減法"求 gcd(m,n)
# Use max(), min() no need to import math
# 用"相減法"時, # Be careful! It maybe b > d or b < d
# 所以要用 a, b = max(b, d), min(b, d)
# 如果用"歐基理德輾轉相除法", 則 b > r (商數一定比餘數大), 
# 則不須max, min 的動作

def gcd_substract_abd(m,n):
    ''' Find gcd(m.n) with division 用"相減法"求 gcd(m,n)'''

    a = max(m, n)
    b = min(m, n)
    d = a - b
      
    while d >0:
        a, b = max(b, d), min(b, d)
        # Be careful! It may b > d or b < d
        d = a - b

    return b

print("gcd_substract_abd")
print(gcd_substract_abd(4,2))
print(gcd_substract_abd(7,3))
print(gcd_substract_abd(5,20))
print(gcd_substract_abd(72511,1229))
# 楊克昌趣味數學編程及拓展, P.59: 72511 = 59 * 1229
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

連分數演算法

$$\frac{33}{105}=3+\frac{1}{5+\frac{1}{2}}$$
一些特殊的數字, 會有很漂亮的連分數表法,
例如:
黃金分割比例 $\varphi =\frac{1+\sqrt{5}}{2}$,
$\varphi$ 有很漂亮的連分數表法:
$$\varphi =1+\frac{1}{1+\frac{1}{1+\frac{1}{1+\frac{1}{1+\cdots }}}}-------(0)$$

$$\displaystyle
 \phi= 1+
 \frac{1} {1+
 \frac{1} {  1
 +\frac{1}  {1
 +\frac{1}{1+\cdots}
  } 
  } 
  }$$
1
2
3
4
5
6
7
8
9

因為 $\varphi$ 會滿足方程式 $x^2-x-1=0$, 所以有
$\varphi =1+\frac{1}{\varphi }$ -------(1)
繼續以(1) 帶入(1)中的 $\varphi$, 得到
$\varphi =1+\frac{1}{1+\frac{1}{\varphi }}$ -------(2)
如此持續代入, 就得到 (0)式.

---

以下是原書的codes, 我們加上說明,
後面會再用 divmod(), 寫一個更精簡易讀的版本

# Noted by P-J Lai MATH NKNU 20221110
# Bradford Tuckfield, 陳仁和譯, 
# 理解演算法Python初學者的深度歷險, 基峰出版, 2021.
# chapter5, P.96
# 連分數演算法其實就是輾轉相除法,
# gcd(a,b) 是取最後非零的餘數為最大公因數
# 連分數則是取每次得出之商數
##>>> help(divmod)
##Help on built-in function divmod in module builtins:
##
##divmod(x, y, /)
##    Return the tuple (x//y, x%y).  Invariant: div*y + mod == x.
    
import math
def continued_fraction(x,y,length_tolerance):
    output = []
    big = max(x,y)
    small = min(x,y)
    # 取大小就是做倒數的動作, 發現不是,
    # 其實連分數演算法與輾轉相除法是完全一樣的, 每次都是前一次之除數除以前一次之餘數.
    # 所以之後不需要再取大小, 因為除數一定比餘數大.
    
    while small > 0 and len(output) < length_tolerance:
        quotient = math.floor(big/small)
        # 以上可以用 big//small 效果一樣
        
        output.append(quotient)
        new_small = big % small
        big = small
        small = new_small

        # 以上可以用 (q, r) = divmod(big, small) 一次得到商數跟餘數 20221111
    return(output)


print(continued_fraction(105,33,10))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

按 F5 執行編譯:

最後一行為
print(continued_fraction(105,33,10))
得

>>> 
=== RESTART: C:\Users\user\Desktop\202302上課\Listing05-4.0 連分數演算法其實就是輾轉相除法.py ===
[3, 5, 2]
1
2
3

附錄

(q, r) = divmod(big, small)

如果用Python 內建的 同時取商與餘數的指令 (q, r) = divmod(big, small)

References

• Bradford Tuckfield, 陳仁和譯, 理解演算法Python初學者的深度歷險, 基峰出版, 2021.