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
No comments:
Post a Comment