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
2015/03/29
Project Euler與Python:Problem 21~30
位於 09:29
標籤: Project Euler, Python
Subscribe to:
Post Comments (Atom)
No comments:
Post a Comment