2015/03/29

Project Euler與Python:Problem 41~50

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,也就是說,到了那麼大的地方,答案仍會小於上限,因此不管大於上限的累加和。

1 comment:

  1. 我看了你的帖子。 这对我很有帮助,也很有帮助。 我很欣赏您在文章中提供的有价值的信息。 谢谢你再次发帖!

    ReplyDelete