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，也就是說，到了那麼大的地方，答案仍會小於上限，因此不管大於上限的累加和。
2015/03/29
Project Euler與Python：Problem 41~50
Project Euler與Python：Problem 31~40
Project Euler針對數學和程式愛好者，設計了一系列的問題，至今已有五百多道，供大家挑戰，每一題都要求應在1分鐘內計算出答案。
我將以Python語言來解題，原始碼放在GitHub（yehnan/project_euler_python），基本上都是先採用暴力法，若運算量太大，無法在1分鐘內算出答案，才著手研究想出更快的演算法。當題目的數字很小時，隨便寫就能在1分鐘內算出答案，但是，也就只有幾題而已，大部分的題目都需要大量計算，可沒那麼好混。
我使用硬體Raspberry Pi 2與軟體Raspbian，作為執行速度的依據。
注意：不提供詳細的說明與分析。
Problem 31: Coin sums，湊硬幣問題，英國硬幣面額有1、2、5、10、20、50、100、200，請問湊出總值200的湊法共有幾種。
若使用遞迴形式，也可以，但較慢；此處使用動態規劃法（dynamic programming）。
底下函式傳入200與coins_england，即可算出答案。
coins_england = (1, 2, 5, 10, 20, 50, 100, 200)
def coin_sum(total, coins):
ways = [1] + ([0] * total)
for coin in coins:
for i in range(coin, total+1):
ways[i] += ways[i - coin]
return ways[total]
Problem 32: Pandigital products，a × b = c，而a、b、c的位數恰恰含有數字1到9，譬如39 × 186 = 7254。請找出所有c，求總和。
使用Python內建程式庫的排列函式。
from itertools import permutations as perm
底下函式判斷某數字（字串形式）是否含有數字1到9，而且各只有一個。
ds = '123456789'
def is_pandigital(num):
ins = sum(1 for d in ds if d in num)
return ins == len(ds)
分析後，a、b、c的位數只有兩種可能，分別是1、4、4，以及2、3、4，採用暴力法處理。底下函式會算出答案。
def pd():
result = set()
answers = set()
for i in ds: # 1 4 4
for jli in perm(ds[:int(i)] + ds[int(i)+1:], 4):
j = ''.join(jli)
k = str(int(i) * int(j))
if len(k) == 4 and is_pandigital(i + j + k):
result.add((i, j, k))
answers.add(int(k))
for ili in perm(ds, 2): # 2 3 4
for jli in perm(''.join(set(ds)-set(ili)), 3):
i = ''.join(ili)
j = ''.join(jli)
k = str(int(i) * int(j))
if len(k) == 4 and is_pandigital(i + j + k):
result.add((i, j, k))
answers.add(int(k))
return sum(answers)
Problem 33: Digit cancelling fractions，找出誤算也正確的約分。
沒什麼好說的，看清楚題目，努力寫出程式解題。
使用Python內建程式庫的gcd，求最大公因數。
from fractions import gcd
判斷分子n與分母d，是否符合題目的條件。
# n: numerator, d: denominator
def is_it(n, d):
ns = str(n)
ds = str(d)
if ns[0] == ds[0]:
x = int(ns[1])
y = int(ds[1])
if (n // gcd(n, d) == x // gcd(x, y) and
d // gcd(n, d) == y // gcd(x, y)):
return True
elif ns[0] == ds[1]:
x = int(ns[1])
y = int(ds[0])
if (n // gcd(n, d) == x // gcd(x, y) and
d // gcd(n, d) == y // gcd(x, y)):
return True
elif ns[1] == ds[0]:
x = int(ns[0])
y = int(ds[1])
if (n // gcd(n, d) == x // gcd(x, y) and
d // gcd(n, d) == y // gcd(x, y)):
return True
return False
判斷某數字是否為能輕易判斷的。
def is_dd(num): # double digit, like 11, 22...
return (num // 10) == (num % 10)
def is_m10(num): # mutiple of 10, like 10, 20...
return (num % 10) == 0
題目本身，用暴力法迭代，找出答案。呼叫底下函式即可算出答案，會回傳乘積約分後的分子與分母，以及共有幾組組合符合條件，其中分母是此題答案。
def dcf():
prod_n = 1
prod_d = 1
count = 0
for d in range(10, 99+1): # two digits
for n in range(10, d): # less than 1
if (not is_dd(n) and not is_dd(d) and
not is_m10(n) and not is_m10(d) and
is_it(n, d)):
prod_n *= n
prod_d *= d
count += 1
gg = gcd(prod_n, prod_d)
return prod_n//gg, prod_d//gg, count
Problem 34: Digit factorials，某些數字，把一個個位數拿出來，作階乘運算後相加，又會等於原本的數字。請找出這些數字、求總和。
使用Python內建程式庫的階乘函式，先算出0到9的階乘，加快速度。
from math import factorial
facts = [factorial(i) for i in range(10)]
找出上限，7位數已經遠遠超過其位數階乘和。
def limit(): # 9999999 is way more than its fact-sum
n = 2
while True:
if 10**(n-1) > facts[9]*n:
return n
n += 1
檢查某數字是否符合條件。
def check(n):
tmp = n
total = 0
while tmp > 0:
total += facts[tmp % 10]
tmp //= 10
return total == n
問題本身，暴力法。呼叫底下函式便可算出答案。
def df():
result = 0
for n in range(10, facts[9] * (limit()-1)):
if check(n):
result += n
return result
Problem 35: Circular primes，197、971、719都是質數，稱為環形質數（circular prime），請問一百萬以下有幾個環形質數。
使用埃氏篩法找出小於n的質數。
def prime_sieve(n):
primes = set(range(2, n+1))
for i in range(2, (n+1) // 2):
if i in primes:
m = 2
while i*m <= n:
primes.discard(i*m)
m += 1
return primes
給定某數n，找出所有的環形排列。
def circulars(n):
result = []
s = str(n)
for i in range(len(s)):
result.append(int(s[i:] + s[0:i]))
return result
判斷是否都是質數。
def all_in(iterable, primes):
for x in iterable:
if x not in primes:
return False
return True
問題本身，底下函式傳入1000000，可算出答案。
def cp(max):
primes = prime_sieve(max)
count = 0
for p in primes:
cs = circulars(p)
if all_in(cs, primes):
count += 1
return count
Problem 36: Double-base palindromes，數字585，不論是十進位表示法還是二進位表示法1001001001，都是迴文數（palindromic ）。小於一百萬的數字中，找出符合此條件的數字，求總和。
判斷十進位、二進位是不是迴文數。
def is_palindrome_10(n):
s = str(n)
return s == s[::-1]
def is_palindrome_2(n):
s = bin(n)[2:]
return s == s[::-1]
def is_palindrome(n):
return is_palindrome_10(n) and is_palindrome_2(n)
問題本身，暴力法，底下函式傳入1000000，可算出答案。
def dbp(max):
return sum(i for i in range(1, max) if is_palindrome(i))
Problem 37: Truncatable primes，3797若從左邊逐一拿掉一個位數，會變成797、97、7，通通都是質數，從右邊逐一拿掉也是。這種數字總共只有11個，不包括2、3、5、7，請找出這些數字求總和。
首先是判斷質數的函式。
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_truncatable_left(n):
s = str(n)
s_len = len(s)
for i in range(1, s_len):
x = s[i:]
if not is_prime(int(x)):
return False
return True
def is_truncatable_right(n):
s = str(n)
s_len = len(s)
for i in range(1, s_len):
x = s[:-i]
if not is_prime(int(x)):
return False
return True
為了簡便，把以上的檢查通通放在一個函式內。
def check_all(n):
return is_prime(n) and is_truncatable_left(n) and is_truncatable_right(n)
最後是問題本身，因已知道共有11個，所以使用內建程式庫itertools的count，從10開始，判斷是否符合條件。底下函式傳入11，可求出答案。
import itertools
def tp(cnt_max):
result = 0
cnt = 0
for n in itertools.count(10):
if check_all(n):
result += n
cnt += 1
if cnt == cnt_max:
break
return result
Problem 38: Pandigital multiples，以192為例，乘上1、2、3後得到192、384、576，此三數連接起來得到192384576，稱為連接乘積，恰好是9位數，而且剛好含有1到9。請問連接乘積中最大的是哪個。
首先撰寫函式判斷某數字是否為9位數、而且剛好含有1到9。
def is_pandigital(num):
ds = '123456789'
s = str(num)
ins = sum(1 for d in ds if d in s)
return ins == len(ds)
然後，底下函式給定某數字參數num，會乘上1、2、3、等等，然後連接，若位數是9個才回傳，若不是則回傳False。回傳值有兩個，前一個是連接乘積，後一個是最大乘上了多少，譬如192最大需乘上3。
def prod(num):
s = str(num)
n = 2
while True:
s += str(num * n)
if len(s) == 9:
return s, n
elif len(s) > 9:
return False, False
n += 1
return False, False
最後是問題本身，迭代1到4位數，因為5位數以上的連接乘積不會是9位數。然後以暴力法迭代每個數字，算出連接乘積，判斷是否符合條件，找出最大的。呼叫底下函式即可算出答案。
def pm():
num_max = 0
for i in range(1, 4+1):
for num in range(10**(i-1), 10**i - 1):
s, n = prod(num)
if s and is_pandigital(s) and int(s) > num_max:
num_max = int(s)
return num_max
Problem 39: Integer right triangles，給定直角三角形周長p為120的話，可算出三組整數邊長的解{20,48,52}、{24,45,51}、{30,40,50}。在周長小於1000（包含）的情況下，請問哪個周長的直角三角形的整數邊長解最多。
用暴力法，會超過1分鐘，必須進行分析，縮減需要處理的計算量。
def irt(p_upperbound):
p_max = 0
t_max = 0
for p in range(p_upperbound//4*2, p_upperbound+1, 2):
t = 0
for a in range(2, int(p / (2**0.5+1)) + 1):
if p * (p-2*a) % (2* (p-a)) == 0:
t += 1
if t > t_max:
t_max = t
p_max = p
return p_max
Problem 40: Champernowne's constant，把從1開始的數字連接起來，變成123456789101112131415161...，請問第1、10、100、1000、10000、100000、1000000個位數的乘積為何。
把從1開始的數字不斷連接起來，直到超過給定長度dn。
def d(dn):
li = []
i = 1
n = 0
while True:
s = str(i)
li.append(s)
n += len(s)
if n > dn:
break
i += 1
return ''.join(li)
然後，拿出想要的位數，相乘。底下函式傳入7，可算出答案。
def cc(pow):
s = d(10**pow)
result = 1
for i in range(pow):
result *= int(s[10**i-1])
return result
Project Euler與Python：Problem 21~30
Project Euler針對數學和程式愛好者，設計了一系列的問題，至今已有五百多道，供大家挑戰，每一題都要求應在1分鐘內計算出答案。
我將以Python語言來解題，原始碼放在GitHub（yehnan/project_euler_python），基本上都是先採用暴力法，若運算量太大，無法在1分鐘內算出答案，才著手研究想出更快的演算法。當題目的數字很小時，隨便寫就能在1分鐘內算出答案，但是，也就只有幾題而已，大部分的題目都需要大量計算，可沒那麼好混。
我使用硬體Raspberry Pi 2與軟體Raspbian，作為執行速度的依據。
注意：不提供詳細的說明與分析。
Problem 21: Amicable numbers，找出小於10000的所有親和數（amicable number），求總和。
首先是找出數字n的所有真因數，不包括n本身。
def divisors(n):
result = [1]
result_2 = []
for i in range(2, int(sqrt(n))+1):
if n % i == 0:
result.append(i)
x = n // i
if x != i:
result_2.insert(0, n//i)
result.extend(result_2)
return result
然後，算出所有因數的和。
def sum_divisors(n):
return sum(divisors(n))
那麼便可判斷某數是否為親和數。
def amicable(a):
b = da = sum_divisors(a)
db = sum_divisors(b)
if a != b and a == db and b == da:
return b
else:
None
最後是題目本身，底下函式傳入10000，即可求出答案。
def an(m):
set_amicable = set()
for a in range(1, m):
if a not in set_amicable:
b = amicable(a)
if b:
set_amicable.add(a)
set_amicable.add(b)
return sum(set_amicable)
Problem 22: Names scores，算出檔案裡英文名字的分數的總和。
首先是算出英文名字的字母次序值。
from io import open
def alphabetical_value(s):
v = 0
for c in s:
v += ord(c) - ord('A') + 1
return v
開啟檔案，使其成為合乎語法的Python串列，運用eval得到串列物件。然後計算每個英文名字的分數，求總和。
def ns(filename):
with open(filename, 'r', encoding='ascii') as fin:
names = sorted(eval(''.join(('[', fin.read(), ']'))))
scores = 0
for i, name in enumerate(names):
scores += alphabetical_value(name) * (i + 1)
return scores
Problem 23: Non-abundant sums，找出所有不會是兩個盈數（abundant number）之和的正整數，求總和。
盈數的定義如下，其中sum_divisors請見Problem 21。
def is_abundant(n):
return sum_divisors(n) > n
然後，從1迭代到28123，測試能否成為兩個盈數之和，不能的話就加起來。
def nas():
abn = set()
result = 0
for n in range(1, 28123):
if is_abundant(n):
abn.add(n)
if not any((n-a in abn) for a in abn):
result += n
return result
Problem 24: Lexicographic permutations，給定數字0到9，按照數值順序來作排列，請問第一百萬個會是多少。
可以撰寫排列函式，然後再找出第幾個。底下的產生器函式，若輸入'0123456789'，便會照順序給出'0123456789'、'0123456798'、'0123456879'、'0123456897'、等等。
def perm_gf(s):
if len(s) == 1:
yield s[0]
else:
for i, c in enumerate(s):
stmp = s[:i] + s[i+1:]
for sub in perm_gf(stmp):
yield c + sub
然後是題目本身，傳入想要第幾個nth（從0起跳）。底下函式傳入'0123456789'與10**6 - 1可得到答案。
def lp_perm(s, nth):
for i, x in enumerate(perm_gf(s)):
if i == nth:
return x
不過，上面這種寫法雖能在1分鐘內求出，但此題只需要找出某一個而已，不需要產生出所有的排列。我們知道，以0開頭的個數有9!個，以1開頭的個數也有9!個，這麼一來便可找出第nth個位在哪一區，然後運用遞迴。
底下使用能計算階乘的函式factorial。算出想要的在哪一區，把該區的開頭數字拿出來，調整剩餘的排列元素與索引，遞迴呼叫。同樣的，底下函式傳入'0123456789'與10**6 - 1可得到答案，快得多了。
def lp(s, nth):
if nth == 0:
return s
else:
seg = factorial(len(s) - 1)
idx = nth // seg
c = s[idx]
stmp = s[:idx] + s[idx+1:]
return c + lp(stmp, nth - idx*seg)
Problem 25: 1000-digit Fibonacci number，找出第一個含有1000個位數的費氏數。
具有記憶功能的費氏數計算函式，遞迴形式。
memo = {0: 0, 1: 1}
def fib_m(n):
if n not in memo:
memo[n] = fib_m(n-1) + fib_m(n-2)
return memo[n]
然後，不斷計算下一個費氏數，直到位數個數達到1000。底下函式傳入1000可算出答案。
def ndfn(nd):
i = 1
while True:
fn = fib_m(i)
if len(str(fn)) >= nd:
return i
i += 1
Problem 26: Reciprocal cycles，在小於1000的正整數d中，哪個1/d擁有最長的循環小數。
運用費馬小定理的話，很快就能求出答案。首先是找出質數的函式。
def prime_sieve(n):
primes = set(range(2, n+1))
for i in range(2, (n+1) // 2):
if i in primes:
m = 2
while i*m <= n:
primes.discard(i*m)
m += 1
return primes
然後是問題本身，如果10**c-1 % p（質數）為0，那麼1/p會有c位數的循環小數，而且分母若有質數p，最多會有p-1個位數的循環小數。
def rc(n):
for p in reversed(sorted(prime_sieve(n))):
c = 1
while pow(10, c, p) - 1 != 0:
c += 1
if p-1 == c:
return p
Problem 27: Quadratic primes，找出a與b，使得二次方程式n**2 + a*n + b能夠給出最多個質數。其中a與b的絕對值小於1000，n從0開始代入。
底下使用埃氏篩法找出質數，但因為並不知道需要判斷多大的數字，所以必須能依情況延伸。
def prime_sieve(n, initials=set(), nmax=2):
primes = set(range(nmax, n+1))
if initials != set():
primes |= initials
for i in range(2, (n+1) // 2):
if i in primes:
m = 2
while i*m <= n:
primes.discard(i*m)
m += 1
return primes
n_max = 100
primes = prime_sieve(n_max)
def is_prime(n):
global n_max
global primes
if n_max < n:
n_max_new = n + (n//2)
primes = prime_sieve(n_max_new, primes, n_max)
n_max = n_max_new
return n in primes
底下函式計算某二次方程式能給出幾個質數。
def howmany(a, b):
n = 0
while True:
if is_prime(n**2 + a*n + b):
n += 1
continue
else:
break
return n
然後是問題本身，暴力法，從-1000到1000迭代a與b。呼叫底下函式可得到答案。
def qp():
prod = None
cnt_max = 0
for a in range(-1000, 1000+1):
for b in range(-1000, 1000+1):
cnt = howmany(a, b)
if cnt > cnt_max:
cnt_max = cnt
prod = a * b
return prod
Problem 28: Number spiral diagonals，從中間開始，按照順時鐘方向擺上數字1、2、3、等等，形成1001×1001的正方形，請問兩條對角線上的數字總和是多少。
沒什麼好說的，根據題目找出規則，然後逐一相加。底下函式傳入1001可得到答案。
def nsd(size):
result = 1 # 中央
for i in range(3, size+2, 2):
result += i**2 # 右上
tmp = (i-2)**2 + 1 + i - 2
result += tmp # 右下
result += tmp + i - 1 # 左下
tmp = tmp + i - 1
result += tmp + i - 1 # 左上
return result
Problem 29: Distinct powers，a與b是從2到100的正整數，請問a**b共會出現幾個不相同的數字。
不多說，直接暴力法。底下函式傳入2、100、2、100，可算出答案。
def dp(al, ah, bl, bh):
ps = set()
for a in range(al, ah+1):
for b in range(bl, bh+1):
ps.add(a**b)
return len(ps)
Problem 30: Digit fifth powers，把數字的每個位數作5次方後相加，會等於原本的數字。請問符合此條件的數字總和是多少。
使用暴力法，但必須知道檢查到哪個上限就可停止。譬如底下函式若傳入代表次方的參數5，若是6位數，每位數作5次方後相加後得到354294（6位數），所以只需檢查到此就夠了。而7位數的次方相加無法得到7位數，也就不必檢查。
使用底下函式得知應檢查到哪裡。
def limit(power):
w = 2
while True:
num_lower = 10 ** (w - 1)
sop_upper = (9**power) * w
if num_lower > sop_upper:
return sop_upper
w += 1
判斷是否符合條件。
def is_sop(n, power):
sop = sum(int(c)**power for c in str(n))
return n == sop
問題本身，逐一檢查。底下函式傳入5，可算出答案。
def dfp(power):
result = 0
for n in range(2, limit(power)+1):
if is_sop(n, power):
result += n
return result
Project Euler與Python：Problem 11~20與67
Project Euler針對數學和程式愛好者，設計了一系列的問題，至今已有五百多道，供大家挑戰，每一題都要求應在1分鐘內計算出答案。
我將以Python語言來解題，原始碼放在GitHub（yehnan/project_euler_python），基本上都是先採用暴力法，若運算量太大，無法在1分鐘內算出答案，才著手研究想出更快的演算法。當題目的數字很小時，隨便寫就能在1分鐘內算出答案，但是，也就只有幾題而已，大部分的題目都需要大量計算，可沒那麼好混。
我使用硬體Raspberry Pi 2與軟體Raspbian，作為執行速度的依據。
注意：不提供詳細的說明與分析。
Problem 11: Largest product in a grid，20×20格子，每格有個數字，把相鄰的四個數相乘，找出乘積最大的。所謂相鄰是指同一方向的連續四個數，方向有上下、左右、左上右下的對角線、右上左下的對角線。
把題目給的20×20格子，放進檔案裡，然後讀取放進二維串列。
from io import open
grid = []
with open('p011_data.txt', 'r', encoding='ascii') as fin:
for line in fin:
grid.append(list(map(int, line.split())))
給定含有數字的串列，算出乘積。
def product(li):
result = 1
for n in li:
result *= n
return result
給定串列、個數count，找出相鄰count個數字相乘後乘積最大的。
def product_max(li, count):
max = 0
for i in range(0, len(li)-count + 1):
m = product(li[i:i+count])
if m > max:
max = m
return max
底下函式可拿出左上右下或右上左下的對角線，由參數slant決定。
def diag(grid, row, col, slant):
d = []
while True:
try:
d.append(grid[row][col])
row += 1
col += slant
except IndexError:
break
return d
終於到了問題本身，底下函式可找出相鄰四個數乘積最大的，但只會尋找一半：每一列、右上半部對角線為左上右下、左上半部對角線為右上左下。
def find_max(grid, count):
max = 0
for row in grid:
m = product_max(row, count)
if m > max:
max = m
for col in range(len(grid[0])):
m = product_max(diag(grid, 0, col, 1), count)
if m > max:
max = m
for col in range(len(grid[0])):
m = product_max(diag(grid, 0, col, -1), count)
if m > max:
max = m
return max
不過，只要把格子旋轉90度，再丟進上面的函式，便可處理另外一半。
def lpiad(grid, count):
m = find_max(grid, count)
grid_rotated = list(zip(*grid))
mr = find_max(grid_rotated, count)
return max(m, mr)
Problem 12: Highly divisible triangular number，找出第一個擁有超過500個因數的三角數（triangle number）。
暴力法太慢，應該先作質因數分解，根據質因數的次方數，便可算出因數的個數。不過此處並未作質因數分解，而是作2與奇數的分解。底下函式，給定參數n，求出n的因數個數。
def divisors_count(n):
dc = 1
cnt = 0
while n % 2 == 0:
n //= 2
cnt += 1
dc *= (cnt + 1)
p = 3
while n != 1:
cnt = 0
while n % p == 0:
n //= p
cnt += 1
dc *= (cnt + 1)
p += 2
return dc
下面這個產生器函式，會不斷給出下一個三角數。
def triangle_number_gf():
tn = 0
i = 1
while True:
tn += i
i += 1
yield tn
最後是問題本身，傳入參數500，便可算出答案。
def hdtn(dc):
for tn in triangle_number_gf():
if divisors_count(tn) > dc:
return tn
Problem 13: Large sum，給定100個數字，皆為50位數，算出總和後，取出前10個位數。
因為只要記憶體夠大，Python的整數便能不斷變大，所以這題對Python來說很簡單。把給定的東西放在檔案裡，呼叫底下函式傳入檔名，即可算出答案。
from io import open
def ls(filename):
with open(filename, 'r', encoding='ascii') as fin:
total = sum(int(line) for line in fin)
return str(total)[:10]
Problem 14: Longest Collatz sequence，根據奇偶歸一猜想，每個正整數經過數次Collatz轉換後，最終會得到1。在小於一百萬的數字中，哪個數需要最多次轉換才能到1。
首先，底下函式給定參數n（起始數），回傳Collatz轉換的次數。
memo = {1: 1}
def collatz(n):
if n not in memo:
if n % 2 == 0:
memo[n] = 1 + collatz(n // 2)
else:
memo[n] = 1 + collatz(3 * n + 1)
return memo[n]
然後，以迴圈迭代所有起始數，找出轉換次數最多的那一個。底下函式傳入1000000即可求出答案。
def lcs(m):
n = 1
n_len = 1
for i in range(2, m+1):
i_len = collatz(i)
if n_len < i_len:
n = i
n_len = i_len
return n
Problem 15: Lattice paths，20×20的格子，從最左上走到最右下，只能向左或向下移動，請問有幾種路徑。
此題等同於：有40個位置，手上有20個紅球和20個綠球，請問共有幾種擺法。也就是40! / 20! / 20!。所以先寫階乘函式。
def factorial(n):
result = 1
for i in range(2, n+1):
result *= i
return result
然後是問題本身，底下函式傳入20與20，即可求出答案。
def lp(row, col):
return factorial(row+col) // factorial(row) // factorial(col)
Problem 16: Power digit sum，求2**1000的所有位數的總和。
因Python支援無限大的整數（只要記憶體足夠的話），所以這題對Python來說很簡單。底下函式傳入1000即可求出答案。
def pds(n):
num = 2 ** n
result = 0
for d in str(num):
result += int(d)
return result
Problem 17: Number letter counts，把1到1000寫成英文，總共有幾個字母。
嗯，只能根據規則把數字轉成英文字串，底下函式只能處理1到1000的情況。
d3 = {1: 'one thousand'}
d2a = { 10: 'ten', 11: 'eleven', 12: 'twelve', 13: 'thirteen',
14: 'fourteen', 15: 'fifteen', 16: 'sixteen',
17: 'seventeen', 18: 'eighteen', 19: 'nineteen' }
d2b = { 2: 'twenty', 3: 'thirty', 4: 'forty', 5: 'fifty',
6: 'sixty', 7: 'seventy', 8: 'eighty', 9: 'ninety' }
d1 = { 1: 'one', 2: 'two', 3: 'three', 4: 'four', 5: 'five',
6: 'six', 7: 'seven', 8: 'eight', 9: 'nine' }
def i2s(i): # from 1 to 1000
a4 = i // 1000
a3 = (i % 1000) // 100
a2 = (i % 100) // 10
a1 = i % 10
result = []
if a4 == 1:
result.append('one thousand')
if a3 != 0 or a2 != 0 or a1 != 0:
result.append(' and ')
if a3 != 0:
result.append(d1[a3])
result.append(' hundred')
if a2 != 0 or a1 != 0:
result.append(' and ')
if a2 == 1:
result.append(d2a[a2*10 + a1])
elif a2 >= 2:
result.append(d2b[a2])
if a1 != 0:
result.append('-')
if a2 != 1 and a1 != 0:
result.append(d1[a1])
return ''.join(result)
然後是問題本身，底下函式傳入1與1000，即可算出答案。
def nlc(n, m=None):
if m is None:
m = n
result = 0
for i in range(n, m+1):
s = i2s(i)
s = s.replace(' ', '').replace('-', '')
result += len(s)
return result
Problem 18: Maximum path sum I，Problem 67: Maximum path sum II，給定呈現三角形的數字層，從頂端開始走，往下一層走時可往左或往右，走到底部時形成一條路徑，把走過的數字加總，請問路徑加總和最大的是多少。
首先寫函式把代表三角形的資料從檔案載入，變成多維串列。
from io import open
def get_triangle(filename):
triangle = []
with open(filename, 'r', encoding='ascii') as fin:
for line in fin:
triangle.append([int(n) for n in line.split()])
return triangle
然後是尋找最大路徑加總和的函式，屬於從上往下的遞迴形式。當你站在某一個位置時，只要算出左下三角形的最大路徑加總和、以及右下三角形的最大路徑加總和，挑出比較大的那一條，便可決定往哪走。
def max_path_sum(tri, row, col):
if row == len(tri)-1:
return tri[row][col]
else:
sub = max(max_path_sum(tri, row+1, col),
max_path_sum(tri, row+1, col+1))
return tri[row][col] + sub
最後是題目本身，傳入資料檔名，即可求出答案。
def mpsi(filename):
tri = get_triangle(filename)
return max_path_sum(tri, 0, 0)
不過，當三角形層數一多（Problem 67），上述寫法就太慢了，必須由下往上，從倒數第二層開始，位置(r, c)的最大路徑加總和，會是(r, c) + (r+1, c)與(r, c) + (r+1, c+1)兩個之中較大的那一個。一層一層往上，最後便可算出頂端(0, 0)的最大路徑加總和。
def max_path_sum(tri, row, col):
for r in range(len(tri)-1 -1, -1, -1):
for c in range(len(tri[r])):
tri[r][c] += max(tri[r+1][c], tri[r+1][c+1])
return tri[row][col]
Problem 19: Counting Sundays，二十世紀每個月第一天是星期日的日子，有幾天。
我們知道1900年1月1日是星期一，定其天數為1；而1900年1月7日是星期日，天數為7，所以只要算出某日期（日月年）的天數，能被7整除的話即是星期日。
首先是判斷潤年的函式。
def is_leapyear(year):
if year%4 == 0 and year%100 != 0 or year%400 == 0:
return 1
else:
return 0
給定月份與年份，回傳該月份有幾天。
month = [31, 28, 31, 30, 31, 30,
31, 31, 30, 31, 30, 31]
def days_of_month(m, y):
return month[m-1] + (is_leapyear(y) if m == 2 else 0)
給定年份，回傳該年份有幾天。
def days_of_year(y):
return sum(month) + is_leapyear(y)
給定日期（日月年），算出其天數。
def date_to_days(date):
dy, mn, yr = date
days = dy
for y in range(1900, yr):
days += days_of_year(y)
for m in range(1, mn):
days += days_of_month(m, yr)
return days
能被7整除，則為星期日。
def is_sunday(days):
return days % 7 == 0
迭代所有年份的每月第一天。
def cs():
count = 0
for y in range(1901, 2000+1):
for m in range(1, 12+1):
days = date_to_days((1, m, y))
if is_sunday(days):
count += 1
return count
Problem 20: Factorial digit sum，加總100階乘的所有位數。
沒什麼好說的，算出階乘，然後把每個位數加起來。
def fact(n):
result = 1
for i in range(1, n):
result *= i
return result
底下函式傳入100，便可求得答案。
def fds(n):
s = str(fact(n))
result = 0
for c in s:
result += int(c)
return result
Project Euler與Python：Problem 1~10
Project Euler針對數學和程式愛好者，設計了一系列的問題，至今已有五百多道，供大家挑戰，每一題都要求應在1分鐘內計算出答案。
我使用硬體Raspberry Pi 2與軟體Raspbian，作為執行速度的依據。
十八世紀數學家Leonhard Euler。
注意：並不提供詳細的說明與分析。
Problem 1: Multiples of 3 and 5，小於1000的正整數中，找出3或5的倍數，相加求總和。
底下函式直接使用暴力法，傳入參數1000可算出此題的答案。
def m35(n):
result = 0
for i in range(1, n):
if i % 3 == 0 or i % 5 == 0:
result += i
return result
也可縮成一行，其中運用了產生器運算式（generator expression）。
def m35(n):
return sum(i for i in range(1, n) if i%3==0 or i%5==0)
Problem 2: Even Fibonacci numbers，找出小於四百萬而且是偶數的費波那契數，相加求總和。
首先是計算費波那契數（fibonacci number）的函式，採取遞迴形式，但會記憶算過的值，免得重複計算。
memo = {0: 0, 1: 1}
def fib(n):
if n not in memo:
memo[n] = fib(n-1) + fib(n-2)
return memo[n]
然後是題目本身，使用暴力法，從小到大算出費波那契數，小於四百萬而且是偶數的話就加起來。底下函式傳入參數4000000，即可得到答案。
def efn(max):
result = 0
n = 1
fn = fib(n)
while fn <= max:
if fn % 2 == 0:
result += fn
n += 1
fn = fib(n)
return result
Problem 3: Largest prime factor，找出600851475143最大的質因數。
首先是能從小到大、不斷給出質數的函式，底下prime_gf是產生器函式，prime_f則是產生器。
def prime_gf():
primes = [2]
def is_prime(n):
for i in primes:
if n % i == 0:
return False
return True
while True:
p = primes[-1]
yield p
while True:
p += 1
if is_prime(p):
primes.append(p)
break
prime_f = prime_gf()
然後是題目本身，底下函式傳入參數600851475143，即可得到答案。
def lpf(n):
max = 1
for i in prime_f:
if i <= n:
while n % i == 0:
n /= i
max = i
else:
break
return max
Problem 4: Largest palindrome product，迴文數，譬如9009，從頭端看過來、與從尾端看過去都一樣，而且9009是兩個2位數相乘後的迴文數中最大的（91 × 99），請問兩個3位數相乘後最大的迴文數是多少。
首先，底下函式可判斷是否為迴文數。
def is_palindrome(n):
ns = str(n)
a, b = 0, len(ns)-1
while a < b:
if ns[a] != ns[b]:
return False
else:
a += 1
b -= 1
return True
然後，底下函式的參數d代表幾位數，以暴力法列出所有的兩個d位數組合，相乘後判斷是不是迴文數，通通記錄起來。
def palindromes(d):
result = []
for x in range(10**d-1, 10**(d-1)-1, -1):
for y in range(x, 10**(d-1)-1, -1):
if is_palindrome(x * y):
result.append((x*y, x, y))
return result
最後進行排序，便可得到最大的迴文數。底下函式傳入參數3可得到答案。
def palindromes_max(d):
return sorted(palindromes(d), key=lambda x: x[0])[-1][0]
Problem 5: Smallest multiple，能夠被1到20所有數字整除的正整數裡，最小是多少。
首先是計算最大公因數與最小公倍數的函式。
def gcd(a, b):
while b:
a, b = b, a%b
return a
def lcm(a, b):
return a * b // gcd(a, b)
然後是題目本身，一次處理兩個數，不斷算出兩個數的最小公倍數，1跟2的最小公倍數是2，2跟3的最小公倍數是6，6跟4的最小公倍數是12，依此類推。底下函式傳入參數1與20，即可算出答案。
def sm(x, y):
m = 1
for n in range(x, y+1):
m = lcm(m, n)
return m
Problem 6: Sum square difference，前100個自然數（正整數）的平方和與和平方的差值。
呃，沒什麼好說的，數字小、計算量少，直接暴力法即可。底下函式傳入參數100可得到答案。
def ssd(n):
sum_sq = sum(i**2 for i in range(1, n+1))
sq_sum = sum(range(1, n+1)) ** 2
return sq_sum - sum_sq
Problem 7: 10001st prime，求第10001個質數。
底下函式可判斷參數n是否為質數，檢查小於n的數之中，是否有可以整除n的數存在，但不需要檢查「小於n」的數，只需檢查「小於n開根號+1」的數即可。
def is_prime(n):
if n % 2 == 0:
return False
p = 3
n_sqrt = n**0.5 + 1
while p < n_sqrt:
if n % p == 0:
return False
p += 2
return True
然後是問題本身，從3開始，逐一判斷下一個（+2）是否為質數，直到目標。底下函式傳入參數10001便可算出答案。
def nth_prime(n):
prime = 2
count = 1
p_next = 3
while count < n:
if is_prime(p_next):
prime = p_next
count += 1
p_next += 2
return prime
Problem 8: Largest product in a series，給定一個有100位數的數字，找出其中相鄰的13個位數、位數的乘積最大。
暴力法，逐一算出乘積，比較哪個大。
def lpias(num, n):
num_str = str(num)
num_len = len(num_str)
def prod(ns):
p = 1
for i in ns:
p *= int(i)
return p
max = 1
for i in range(0, num_len - n + 1):
m = prod(num_str[i:i+n])
if m > max:
max = m
return max
Problem 9: Special Pythagorean triplet，只有一組畢氏三元數組a、b、c能使得a+b+c等於1000，請算出a*b*c。
首先，判斷a、b、c是否為畢氏三元數組。
def is_pythagorean_triplet(a, b, c):
return a**2 + b**2 == c**2
然後，底下產生器函式可產生出a + b + c = n而且a < b < c的三數組。注意，其中range的參數並非最佳化。
def abc_gf(n):
for c in range(n//3, n):
for b in range((n-c)//2, n-c):
a = n - c - b
if a < b < c:
yield a, b, c
最後是問題本身，底下函式傳入參數1000，即可找出答案。
def spr(n):
for a, b, c in abc_gf(n):
if is_pythagorean_triplet(a, b, c):
return a*b*c
Problem 10: Summation of primes，找出小於兩百萬的所有質數，加總求和。
使用埃拉托斯特尼篩法找出小於n的質數。
def prime_sieve(n):
primes = set(range(2, n+1))
for i in range(2, (n+1+1) // 2):
if i in primes:
m = 2
while i*m <= n:
primes.discard(i*m)
m += 1
return primes
然後加總即可。底下函式傳入參數2000000即可算出答案。
def sop(n):
return sum(prime_sieve(n))
參考資料：
- Project Euler Solutions - Dreamshire，相當不錯，說明簡短精要。
- Project Euler -- Python，只有程式碼。
- Project Euler Solutions - Zach Denton，提供多種語言的解答程式碼，部分提供分析與說明。
- Project Euler: Python solutions，目標改成在10秒內算出答案。
- GitHub tokland/pyeuler，以純函數式程式設計的方式來解題。
