iT邦幫忙

第 11 屆 iT 邦幫忙鐵人賽

DAY 16
3
Software Development

活用python- 路遙知碼力,日久練成精系列 第 16

Day16- Project1- 數字世界的朋友與戀人 (趣味數論問題: 友好數、婚約數)

路遙知碼力,日久練成精- 只要在程式之路鑽研的夠深,便能夠充分發揮程式碼的力量; 練習的日子夠久,便能夠練成寫出精簡代碼的能力。

大家好,我是心原一馬,內心原來一心喜歡打程式碼。
經歷了十五天python精華語法介紹的前奏後,
終於要來做第一個程式小專題囉。

古希臘人相信每個數字都有自己代表的意義,
例如1表示創造,2表示女性,3代表男性,
2+3變成的5代表結婚,
接下來則是特別的6,
6的因數有1, 2, 3, 6,
6除了自己之外的因數和恰好是1+2+3=6,
因此6被視為是一個「完美數」。

今天的主題主要是探討兩種特別的數對:「友好數」、「婚約數」。
在以前計算機、電腦還不發達的年代,
要找到諸如此類特別的數字,
只能靠數學家用大量的數學技巧去慢慢發現。

今日,藉著電腦高速的運算速度,
即便只有中學的數學知識,
你我都有能力找出這些特別的數,
以滿足人類自古以來的好奇心

完美數在數學界的名聲比較大,
數學上已經有大量研究了,
今天我們主要來找友好數(amicable pair)、婚約數(betrothed numbers)

友好數、婚約數

首先給個定義:
友好數就是兩個正整數,
彼此除了自身之外的正因數之和與對方相等。
例如 (220, 284) 是一對最小的友好數。
220 自身外的正因數和為1+2+4+5+10+11+20+22+44+55+110=284
284 自身外的正因數和為1+2+4+71+142=220
所以友好數就像一對好朋友一樣,
彷彿一種「你中有我,我中有你」的感情。

婚約數的定義跟友好數非常相似,
婚約數是彼此除了1與自身之外的正因數之和與對方相等。
例如 (48, 75) 是一對最小的婚約數。
48除了1和本身的正因數相加是:2+3+4+6+8+12+16+24=75,
75除了1和本身的正因數相加是:3+5+15+25=48。

相信這時候可能部分讀者看到有點頭痛了,
小馬彷彿感受到讀者們的心聲:
蛤,寫程式還要會算數學哦,
我以前最不擅長的科目就是數學了

由於人們希望藉助電腦程式高速的運算能力來解決問題,
數學能力確實有助於寫出更有效率的演算法。
但是不要緊,本系列文所需要的數學能力大概就到中學數學而已,
小馬會一步步帶大家理解的。
大家堅持下去。

第一步: 求出整數的正因數和

希望大家還記得什麼是因數,
n的因數就是可以整除n的數字。
我們運用之前學過的列表生成式,
達到方便求和的效果。
以下程式可以求出除了自身外的正因數和。

def divisorsum(n):
    return sum([i for i in range(1,n) if n%i==0])

相信這是最直覺的方法,
我們從1開始搜索,
利用 n%i==0 這個判斷式判斷是否為因數,
如果是因數就添加到列表中,
最後用sum()函數求和。

第二步: 雖然非主軸但順便教大家求完美數

有了求正因數和的函數,
我們可以很容易判斷一個數是否為完美數。

def isPerfect(n):
    return divisorsum(n) == n

我們定義一個函數isPerfect(n)
判斷一個數是不是完美數,
即判斷一個數自身外的正因數和是否等於自己。
加入以下程式找出小於10000的完美數

for i in range(1,10000):
    if isPerfect(i):
        print(i)

結果:
6
28
496
8128
你可以參考維基百科
確認我們找到的完全數正確無誤。

第三步: 更有效率的求出整數的正因數和

求正因數和的程式碼很精簡,
但目前我們仍然不太滿意,
當我們將上述程式改為尋找小於十萬的完美數:

for i in range(1,100000):
    if isPerfect(i):
        print(i)

你可以試著在你的電腦執行看看,
執行時間很可能會跑個幾十秒。
因為我們的計算因數和的函數是沒效率的。
我們如何找到n的所有因數呢?
原本我們讓程式去檢查所有1~n-1的數字是否可以整除它
但實際上,我們只要搜索到根號n就可以找到所有因數。

對n來說,
我們把小於根號n 的數字稱為「小數字」,
大於根號n的數字稱為「大數字」。
例如根號30大約是5點多,

30 = 1*30
   = 2*15
   = 3*10
   = 5*6

每個分解都必然是一個「小數字」乘上「大數字」。
你硬是讓兩個「大數字」相乘,
乘出來的數字就太大了。
(兩個大於根號n的數相乘一定大於n)

所以說,我們只要能夠找到小於根號30的因數「1, 2, 3, 5」,
便能透過除法去算出大於根號30的因數「30, 15, 10, 6」。

先教大家如何在python開根號,
開根號在數學上相當於「0.5次方」,
例如:

>>> 2 ** 0.5
1.4142135623730951
>>> 4 ** 0.5
2.0
>>> 9 ** 0.5
3.0

但是注意用次方運算出來的數是帶有小數點的浮點數,
我們可以取int()把小數點去掉,例如:

>>> 30 ** 0.5
5.477225575051661
>>> int(30 ** 0.5)
5

根據上述想法,
我們可以把求正因數和的計算式改寫如下:

def divisorsum(n):
    return sum([sum([i,n//i]) for i in range(1,int(n**0.5)+1) if n%i==0])-n

看這段代碼可能不容易理解,
我們試著印出列表中的內容看看:
(假設n是30)

print([[i,30//i] for i in range(1,int(30**0.5)+1) if 30%i==0])

結果為 [[1, 30], [2, 15], [3, 10], [5, 6]]
這不就是上例分解30為兩個整數乘積的狀況嗎?
我們把列表中的數字全部相加就是因數和了。
但是注意若我們把30改為9 (9本身是完全平方數),
結果為 [[1, 9], [3, 3]]
「根號9」這個因數會被重複計算,
簡潔起見,我們把列表內的列表改為集合以去除重複

print([{i,9//i} for i in range(1,int(9**0.5)+1) if 9%i==0])

結果變為[{1, 9}, {3}]

因此原程式改為:

def divisorsum(n):
    return sum([sum({i,n//i}) for i in range(1,int(n**0.5)+1) if n%i==0])-n

注意因為我們這邊要計算除了「自身」外的正因數和,
計算的結果還要再減去n。

第四步: 求出小於十萬的友好數數對

先做個小測試,
我們用簡單的列表生成式算出1~6除了自身外的正因數之和:

print([(i,divisorsum(i)) for i in range(1,7)])

結果為:
[(1, 0), (2, 1), (3, 1), (4, 3), (5, 1), (6, 6)]
意思是
1(除了自身)的正因數和=0,
2(除了自身)的正因數和=1,

6(除了自身)的正因數和=6。
大家可以自行驗算。

還記得什麼是友好數嗎?
友好數的條件是自身外的正因數和等於彼此,
例如(220, 284) 是一對最小的友好數,
divisorsum(220)=284divisorsum(284)=220
也就是說,若n滿足divisorsum(divisorsum(n))==n的條件,
則(n, divisorsum(n))便是一組友好數。

噢,對了,divisorsum這個函數的名字有點長,
這邊故且以DS這個名字代替。
改寫如下:

def DS(n):
    return sum([sum({i,n//i}) for i in range(1,int(n**0.5)+1) if n%i==0])-n

另外我們定義amicablePair(k)函數,
計算小於k的友好數,
判斷條件就用 DS(DS(i))==i 來判斷,如下:

def amicablePair(k):
    return [(i,DS(i)) for i in range(1,k+1) if DS(DS(i))==i]

我們先試試看效果如何,
測試找300以下的友好數:

print(amicablePair(300))

結果為[(6, 6), (28, 28), (220, 284), (284, 220)]
疑疑? 好像有不速之客哦,
6和28本身是完全數也被收納進來了,
看到這邊,
小馬覺得6這個數自己和自己當朋友,
乾脆6不要叫「完全數」,
改名叫「自戀數」算了。
(6, 6)因為是同一個數,
並不是我們想要的答案。
另外(220, 284), (284, 220)其實算同一組數對,
被重複計算了。
我們還要再加上一條判斷條件i< DS(i)才完整,
調整如下:

def amicablePair(k):
    return [(i,DS(i)) for i in range(1,k+1) if i< DS(i) and DS(DS(i))==i]

重新測試結果:

print(amicablePair(300))

結果為[(220, 284)]
這次已經能正確找到一組數對了。
試試找出小於十萬的友好數吧:

print(amicablePair(100000))

結果為[(220, 284), (1184, 1210), (2620, 2924), (5020, 5564), (6232, 6368), (10744, 10856), (12285, 14595), (17296, 18416), (63020, 76084), (66928, 66992), (67095, 71145), (69615, 87633), (79750, 88730)]

第五步: 舉一反三- 求出小於十萬的婚約數數對

婚約數與友好數定義很類似,
婚約數是彼此除了1與自身之外的正因數之和與對方相等。
我們定義函數DS2(n) 計算n除了1與自身之外的正因數之和:

def DS2(n):
    return sum([sum({i,n//i}) for i in range(1,int(n**0.5)+1) if n%i==0])-n-1

只要剛剛的講解你都有跟上的話,
計算婚約數的函數也只要微調一下即可:

def betrothedPair(k):
    return [(i,DS2(i)) for i in range(1,k+1) if i< DS2(i) and DS2(DS2(i))==i]

再附上小於十萬的婚約數結果:

print(betrothedPair(100000))

結果: [(48, 75), (140, 195), (1050, 1925), (1575, 1648), (2024, 2295), (5775, 6128), (8892, 16587), (9504, 20735), (62744, 75495)]

第六步: 大功告成,總結所有程式碼

我們將今天講解的所有程式碼統整在一起,
並順便計算出小於十萬的完美數。
完整的程式碼如下:

def DS(n):
    return sum([sum({i,n//i}) for i in range(1,int(n**0.5)+1) if n%i==0])-n

def DS2(n):
    return sum([sum({i,n//i}) for i in range(1,int(n**0.5)+1) if n%i==0])-n-1

def findPerfect(k):
    return [i for i in range(1,k+1) if DS(i) == i]

def amicablePair(k): #求小於k的友好數對
    return [(i,DS(i)) for i in range(1,k+1) if i< DS(i) and DS(DS(i))==i]

def betrothedPair(k): #求小於k的婚約數對
    return [(i,DS2(i)) for i in range(1,k+1) if i< DS2(i) and DS2(DS2(i))==i]
    
print(findPerfect(100000))
print(amicablePair(100000))
print(betrothedPair(100000))

結果:
[6, 28, 496, 8128]
[(220, 284), (1184, 1210), (2620, 2924), (5020, 5564), (6232, 6368), (10744, 10856), (12285, 14595), (17296, 18416), (63020, 76084), (66928, 66992), (67095, 71145), (69615, 87633), (79750, 88730)]
[(48, 75), (140, 195), (1050, 1925), (1575, 1648), (2024, 2295), (5775, 6128), (8892, 16587), (9504, 20735), (62744, 75495)]

名字彩蛋

結果解完了,不知讀者們有發現友好數婚約數的命名由來嗎?
目前已知的友好數對全部都是「兩個奇數」或同為「兩個偶數」的組合,
而目前已知的婚約數則全部都是「一奇一偶」的組合。
或許因為這樣的特性(相同奇偶性 => 同性; 相異奇偶性 => 異性),
使得這樣的數對各自稱為友好數婚約數呢。

以上就是今天的內容了,
原本看似不好解的數字問題,
經過python巧手一變之後,
可以求出小於十萬的完美數、友好數及婚約數
完整的程式附帶印出結果總共就十三行而已,
相信如果你本來使用其它程式語言,
光寫判斷質數這件事可能就不只十三行了吧?
Python真是太有魅力了。
疑? 說到判斷質數,
Day13的範例13-3好像留有判斷質數如何更有效率的問題未解。
機會難得,給個課後練習好了。

課後練習

在範例13-3中,
給出了兩種判斷質數的解法:

解法一

def isPrime(n):
    return len([i for i in range(1,n+1) if n%i==0])==2

解法二

def isPrime(n):
    return sum([1 if n%i==0 else 0 for i in range(1,n+1)])==2

目前這兩個解法都有兩個問題可以優化:

  1. 希望檢測到1以外的因數提早結束程式 (例如100可以被2整除,便不必再檢查3,4,5...可否整除100)
  2. 檢查到根號n以下的數能否整除n即可 (原因與今日教的優化找正因數和的效率相同)

我們的課後練習要做的是優化第2點,
第1點於明日介紹新函數時一併解答。

習題: 優化找質數效率

請修改上述「解法一」或「解法二」的isPrime函數,
使程式僅檢查到根號n以下的數能否整除n即完成判斷n是否為質數。

小馬貼心給你一個檢驗程式,
讓你在修改完函數後,
可以自我檢查結果是否正確。

for i in range(1,20):
    print(f'{i}是質數?: {isPrime(i)}')

預期結果:

1是質數?: False
2是質數?: True
3是質數?: True
4是質數?: False
5是質數?: True
6是質數?: False
7是質數?: True
8是質數?: False
9是質數?: False
10是質數?: False
11是質數?: True
12是質數?: False
13是質數?: True
14是質數?: False
15是質數?: False
16是質數?: False
17是質數?: True
18是質數?: False
19是質數?: True

上一篇
Day15- 想要「程」長,不只語法
下一篇
Day17- 存在所有人喜歡的文章?認識any(), all()函數
系列文
活用python- 路遙知碼力,日久練成精30

1 則留言

0
sowhat1124
iT邦新手 5 級 ‧ 2020-09-28 18:17:43

感謝您的教學,有個地方想請問一下
我把 amicablePair(k)、betrothedPair(k) 中 if 後的兩個條件順序對調,想說應該不會影響結果

def amicablePair(k): #求小於k的友好數對
    return [(i,DS(i)) for i in range(1,k+1) if DS(DS(i))==i and i< DS(i)]

def betrothedPair(k): #求小於k的婚約數對
    return [(i,DS2(i)) for i in range(1,k+1) if DS2(DS2(i))==i and i< DS2(i)]

結果 amicablePair(k) 可以正常運作,但 betrothedPair(k) 則會報錯,錯誤訊息為 can't convert complex to int ,請問這問題出在哪裡呢?

我要留言

立即登入留言