Project Euler針對數學和程式愛好者，設計了一系列的問題，至今已有五百多道，供大家挑戰，每一題都要求應在1分鐘內計算出答案。
我將以Python語言來解題，原始碼放在GitHub（yehnan/project_euler_python），基本上都是先採用暴力法，若運算量太大，無法在1分鐘內算出答案，才著手研究想出更快的演算法。當題目的數字很小時，隨便寫就能在1分鐘內算出答案，但是，也就只有幾題而已，大部分的題目都需要大量計算，可沒那麼好混。
我使用硬體Raspberry Pi 2與軟體Raspbian，作為執行速度的依據。
注意：不提供詳細的說明與分析。
Problem 41: Pandigital prime，n位數若含有1、2、3、...、到n，稱為pandigital。請問既是pandigital且是質數的數字，最大的是多少。
若直接檢查1到999999999，耗時甚久。經過分析，所有符合pandigital條件的2、3、5、6、8、9位數，都可被3整除，所以不是質數，不必檢查。
首先是檢查質數的函式。
def is_prime(n):
if n < 2:
return False
if n == 2:
return True
if n % 2 == 0:
return False
for x in range(3, int(n ** 0.5) + 1, 2):
if n % x == 0:
return False
return True
然後使用Python內建模組itertools的permutations，幫我們產生出所有符合pandigital的n位數，然後只需檢查是否為質數即可。呼叫底下函式便可算出答案。
from itertools import permutations as perm
def pp():
n_max = 0
dcheck = [4, 7]
for d in dcheck:
ds = '123456789'[:d]
for nst in perm(ds):
n = int(''.join(nst))
if is_prime(n) and n > n_max:
n_max = n
return n_max
Problem 42: Coded triangle numbers，給定檔案中含有許多英文單字，其中某些單字的值等於三角數（triangle number），請問有幾個。
根據定義，底下函式可算出三角數。
def tn(n):
return n * (n+1) // 2
然後是判斷某值是否為三角數。有記憶功能，當欲判斷的值不在裡頭時，就繼續算出更大的三角數放進去。
memo = [0, 1, 3, 6, 10, 15, 21, 28, 36, 45, 55]
def is_tn(n):
mn = len(memo)
mmax = memo[-1]
while n > mmax:
mmax = tn(mn)
memo.append(mmax)
mn += 1
if n in memo:
return True
return False
然後根據定義，算出某單字的值。
def word_value(word):
value = 0
for c in word:
value += ord(c) - ord('A') + 1
return value
最後是問題本身，開檔讀檔，逐一判斷裡頭的單字是否符合條件。呼叫底下函式，傳入檔名，便可得到答案。
from io import open
def ctn(filename):
result = 0
with open(filename, 'r', encoding='ascii') as fin:
words = eval(''.join(('[', fin.read(), ']')))
for word in words:
if is_tn(word_value(word)):
result += 1
return result
Problem 43: Sub-string divisibility，呃，題目很長，請自己看。
若用暴力法，將會超過1分鐘，必須進行分析，縮減需要處理的量。發現只有尾端為....952867或....357289的數字才可能符合條件。
判斷某數（字串形式）是否符合條件。
def is_ssd(ns):
divs = (2, 3, 5, 7, 11, 13, 17)
shift = 1
for i in range(0, len(divs)):
n_sub = int(ns[i+shift:i+shift+3])
if n_sub % divs[i] != 0:
return False
return True
使用內建模組itertools的permutations，產生所有的排列。呼叫底下函式即可求得答案。
from itertools import permutations as perm
def ssd():
result = 0
for np in perm('0134'):
if np[0] == '0':
continue
ns = ''.join(np) + '952867'
if is_ssd(ns):
result += int(ns)
for np in perm('0146'):
if np[0] == '0':
continue
ns = ''.join(np) + '357289'
if is_ssd(ns):
result += int(ns)
return result
Problem 44: Pentagon numbers，找出一對五角數（entagonal number）Pj與Pk，使得其和與差也都是五角數，請問最小的Pj - Pk的絕對值是多少。
假定Pj小於Pk，根據題目要求，Pk - Pj = Pi是五角數，Pk + Pj = Pn是五角數，所以Pn - Pk是五角數，兩式相加得到2Pk = Pi + Pn，所以Pn - 2Pk = Pi是五角數，想要找的就是Pi = Pk - Pj。
def pn():
pn_set = set()
n = 1
while True:
pn = n * (3*n-1) // 2
for pk in pn_set:
if pn-pk in pn_set and pn-2*pk in pn_set:
return pn-2*pk
pn_set.add(pn)
n += 1
Problem 45: Triangular, pentagonal, and hexagonal，請找出第三個三角數、五角數、六角數皆相同的數字。
底下這個產生器函式，參數f是個函式，將會傳入能計算三角數、五角數、六角數的函式，不斷給出下一個值。
def gen_n(f):
i = 1
while True:
n = f(i)
yield n
i += 1
比較兩個函數的值，相同時才給出。
def sub_gf(fa, fb):
a = next(fa)
b = next(fb)
while True:
if a < b:
a = next(fa)
elif a > b:
b = next(fb)
else:
yield a
a = next(fa)
題目本身，暴力法，不斷比較。呼叫底下函式傳入3，可得到答案。
def tph(n):
tn = gen_n(lambda i: i * (i+1) // 2)
pn = gen_n(lambda i: i * (3*i - 1) // 2)
hn = gen_n(lambda i: i * (2*i - 1))
sub1 = sub_gf(tn, pn)
sub2 = sub_gf(sub1, hn)
result = 0
for i in range(n):
result = next(sub2)
return result
Problem 46: Goldbach's other conjecture，Goldbach說每個奇數組合數，都可被寫作質數加上兩倍的某數平方。但這是錯的，請找出無法滿足條件的最小的奇數組合數。
底下函式有一半會不斷找出下一個質數，2的倍數不必看，而且3的倍數也不必看，觀察需要檢查的數5、7、(9)、11、13、(15)、17、19、(21)、23、25、(27)、...，其中括起來的可跳過，發現有規律性，由底下程式的n += 3-flag與flag = -flag實作。
當n不是質數時，就檢查它是否不能被寫作質數加上兩倍的某數平方。呼叫底下函式可求得答案。
def goc():
n = 5
flag = +1
primes = set([2, 3])
while True:
if all(n%p!=0 for p in primes):
primes.add(n)
else:
if not any((n - 2*i*i) in primes for i in range(1, n)):
return n
n += 3-flag
flag = -flag
Problem 47: Distinct primes factors，找出連續四個數字，每個數字都有四個不同的質因數。
底下的trial_division與prime_factors取自他人寫的輔助函式，譬如呼叫prime_factors(14)會得到((2, 1), (7, 1))。
def trial_division(n, bound=None):
if n == 1:
return 1
for p in [2, 3, 5]:
if n % p == 0:
return p
if bound == None:
bound = n
dif = [6, 4, 2, 4, 2, 4, 6, 2]
m = 7
i = 1
while m <= bound and m*m <= n:
if n % m == 0:
return m
m += dif[i%8]
i += 1
return n
def prime_factors(n):
if n == 1:
return []
factors = []
while n != 1:
p = trial_division(n)
e = 1
n /= p
while n%p == 0:
e += 1
n /= p
factors.append((p, e))
factors.sort()
return tuple(factors)
有了能做質因數分解的函式後，題目本身就簡單了，不斷測試連續的四個數字即可。呼叫底下函式便可求得答案。
def dpf():
n = 647
while True:
if (len(prime_factors(n)) == 4 and
len(prime_factors(n+1)) == 4 and
len(prime_factors(n+2)) == 4 and
len(prime_factors(n+3)) == 4):
return n
n += 1
Problem 48: Self powers，1**1 + 2**2 + 3**3 + ... + 1000**1000的最後十個位數是多少。
這題對Python來說很簡單，直接暴力法。底下函式傳入1000與10，可得到答案。
def sp(n, d=None):
result = 0
for i in range(1, n+1):
result += i ** i
result = str(result)
if d is not None:
result = result[-d:]
return result
Problem 49: Prime permutations，4位數的等差數列：1487、4817、8147，等差是3330，每個數都是質數，而且互為重新排列後的數。4位數的等差數列中，還有一組（三個數字）也符合此條件，請找出來。
只要有判斷質數與判斷是否互為重新排列的函式，便可輕鬆解題。
def is_prime(n):
if n < 2:
return False
if n == 2:
return True
if n % 2 == 0:
return False
for x in range(3, int(n ** 0.5) + 1, 2):
if n % x == 0:
return False
return True
def is_perm(a, b):
return sorted(str(a)) == sorted(str(b))
然後是題目本身，底下函式傳入1001與1488，可求得1487、4817、8147。若傳入1489與9999，便可算出答案。
def pp(n_from, n_to):
for n in range(n_from, n_to+1, 2):
n2, n3 = n+3330, n+6660
if (is_prime(n) and is_prime(n2) and is_prime(n3) and
is_perm(n, n2) and is_perm(n2, n3)):
return str(n) + str(n2) + str(n3)
return None
Problem 50: Consecutive prime sum，質數41可以寫成連續6個質數相加：2 + 3 + 5 + 7 + 11 + 13。請問小於一百萬的質數當中，可以寫成連續質數相加而且項目個數最多的是哪一個。
先找出小於一百萬的質數，然後累加，再使用雙層迴圈，讓累加和相減，測試是否符合條件。底下是找出質數的函式。
def prime_sieve(n):
sieve = [True] * (n//2)
for i in range(3, int(n**0.5)+1, 2):
if sieve[i//2]:
sieve[i*i//2::i] = [False] * ((n-i*i-1) // (2*i)+1)
return [2] + [2*i+1 for i in range(1, n//2) if sieve[i]]
然後是累加的函式。
def cumulative_sum(li):
result = [0]
tmp = 0
for e in li:
if tmp > 1000000: # ! 注意！
break # ! 注意！
tmp += e
result.append(tmp)
return result
接下來是問題本身，其中primes_x多了個[0]，因為這樣比較方便。底下函式傳入1000000，可算出答案。
def cps(m):
primes = prime_sieve(m)
primes_set = set(primes)
primes_sum = cumulative_sum(primes)
p_max = 2
p_len_max = 1
for i in range(p_len_max, len(primes_sum)):
for j in range(i-1-p_len_max, -1, -1):
p = primes_sum[i] - primes_sum[j]
if p > m:
break
if p in primes_set and p > p_max:
p_max = p
p_len_max = i - j
return (p_max, p_len_max)
雖然上面程式可以在1分鐘內算出答案，但其實有偷吃步，地方在標示紅色「注意！」之處，把大於1000000的累加和砍掉了，因為這題答案不會落在那一區，所以沒問題。若拿掉「注意！」的程式碼，將需要超過5分鐘的計算時間。我在網路上找了許多解答（主要是以Python寫的），都有類似的部分，才能在1分鐘內完成。有人把上限提高到100,000,000，而結果是99,819,619，也就是說，到了那麼大的地方，答案仍會小於上限，因此不管大於上限的累加和。
我看了你的帖子。 这对我很有帮助，也很有帮助。 我很欣赏您在文章中提供的有价值的信息。 谢谢你再次发帖！
ReplyDelete