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