2015/03/29

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

No comments:

Post a Comment