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

Project Euler與Python:Problem 31~40

Project Euler針對數學和程式愛好者,設計了一系列的問題,至今已有五百多道,供大家挑戰,每一題都要求應在1分鐘內計算出答案。

我將以Python語言來解題,原始碼放在GitHub(yehnan/project_euler_python),基本上都是先採用暴力法,若運算量太大,無法在1分鐘內算出答案,才著手研究想出更快的演算法。當題目的數字很小時,隨便寫就能在1分鐘內算出答案,但是,也就只有幾題而已,大部分的題目都需要大量計算,可沒那麼好混。

我使用硬體Raspberry Pi 2與軟體Raspbian,作為執行速度的依據。

注意:提供詳細的說明與分析。

Problem 31: Coin sums,湊硬幣問題,英國硬幣面額有1、2、5、10、20、50、100、200,請問湊出總值200的湊法共有幾種。

若使用遞迴形式,也可以,但較慢;此處使用動態規劃法(dynamic programming)。

底下函式傳入200與coins_england,即可算出答案。

coins_england = (1, 2, 5, 10, 20, 50, 100, 200)

def coin_sum(total, coins):
    ways = [1] + ([0] * total)
    for coin in coins:
        for i in range(coin, total+1):
            ways[i] += ways[i - coin]
    return ways[total]


Problem 32: Pandigital products,a × b = c,而a、b、c的位數恰恰含有數字1到9,譬如39 × 186 = 7254。請找出所有c,求總和。

使用Python內建程式庫的排列函式。

from itertools import permutations as perm

底下函式判斷某數字(字串形式)是否含有數字1到9,而且各只有一個。

ds = '123456789'

def is_pandigital(num):
    ins = sum(1 for d in ds if d in num)
    return ins == len(ds)


分析後,a、b、c的位數只有兩種可能,分別是1、4、4,以及2、3、4,採用暴力法處理。底下函式會算出答案。
           
def pd():
    result = set()
    answers = set()
    for i in ds: # 1 4 4
        for jli in perm(ds[:int(i)] + ds[int(i)+1:], 4):
            j = ''.join(jli)
            k = str(int(i) * int(j))
            if len(k) == 4 and is_pandigital(i + j + k):
                result.add((i, j, k))
                answers.add(int(k))
    for ili in perm(ds, 2): # 2 3 4
        for jli in perm(''.join(set(ds)-set(ili)), 3):
            i = ''.join(ili)
            j = ''.join(jli)
            k = str(int(i) * int(j))
            if len(k) == 4 and is_pandigital(i + j + k):
                result.add((i, j, k))
                answers.add(int(k))
    return sum(answers)


Problem 33: Digit cancelling fractions,找出誤算也正確的約分。

沒什麼好說的,看清楚題目,努力寫出程式解題。

使用Python內建程式庫的gcd,求最大公因數。

from fractions import gcd

判斷分子n與分母d,是否符合題目的條件。

# n: numerator, d: denominator
def is_it(n, d):
    ns = str(n)
    ds = str(d)
    if ns[0] == ds[0]:
        x = int(ns[1])
        y = int(ds[1])
        if (n // gcd(n, d) == x // gcd(x, y) and
           d // gcd(n, d) == y // gcd(x, y)):
           return True
    elif ns[0] == ds[1]:
        x = int(ns[1])
        y = int(ds[0])
        if (n // gcd(n, d) == x // gcd(x, y) and
           d // gcd(n, d) == y // gcd(x, y)):
           return True
    elif ns[1] == ds[0]:
        x = int(ns[0])
        y = int(ds[1])
        if (n // gcd(n, d) == x // gcd(x, y) and
           d // gcd(n, d) == y // gcd(x, y)):
           return True
    return False


判斷某數字是否為能輕易判斷的。
   
def is_dd(num): # double digit, like 11, 22...
    return (num // 10) == (num % 10)
   
def is_m10(num): # mutiple of 10, like 10, 20...
    return (num % 10) == 0


題目本身,用暴力法迭代,找出答案。呼叫底下函式即可算出答案,會回傳乘積約分後的分子與分母,以及共有幾組組合符合條件,其中分母是此題答案。

def dcf():
    prod_n = 1
    prod_d = 1
    count = 0
    for d in range(10, 99+1):  # two digits
        for n in range(10, d): # less than 1
            if (not is_dd(n) and not is_dd(d) and
                not is_m10(n) and not is_m10(d) and
                is_it(n, d)):
                prod_n *= n
                prod_d *= d
                count += 1
    gg = gcd(prod_n, prod_d)
    return prod_n//gg, prod_d//gg, count 


Problem 34: Digit factorials,某些數字,把一個個位數拿出來,作階乘運算後相加,又會等於原本的數字。請找出這些數字、求總和。

使用Python內建程式庫的階乘函式,先算出0到9的階乘,加快速度。

from math import factorial
facts = [factorial(i) for i in range(10)]


找出上限,7位數已經遠遠超過其位數階乘和。

def limit(): # 9999999 is way more than its fact-sum
    n = 2
    while True:
        if 10**(n-1) > facts[9]*n:
            return n
        n += 1


檢查某數字是否符合條件。

def check(n):
    tmp = n
    total = 0
    while tmp > 0:
        total += facts[tmp % 10]
        tmp //= 10
    return total == n


問題本身,暴力法。呼叫底下函式便可算出答案。 
 
def df():
    result = 0
    for n in range(10, facts[9] * (limit()-1)):
        if check(n):
            result += n
    return result


Problem 35: Circular primes,197、971、719都是質數,稱為環形質數(circular prime),請問一百萬以下有幾個環形質數。

使用埃氏篩法找出小於n的質數。

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


給定某數n,找出所有的環形排列。
  
def circulars(n):
    result = []
    s = str(n)
    for i in range(len(s)):
        result.append(int(s[i:] + s[0:i]))
    return result


判斷是否都是質數。 

def all_in(iterable, primes):
    for x in iterable:
        if x not in primes:
            return False
    return True


問題本身,底下函式傳入1000000,可算出答案。
   
def cp(max):
    primes = prime_sieve(max)
    count = 0
    for p in primes:
        cs = circulars(p)
        if all_in(cs, primes):
            count += 1
    return count


Problem 36: Double-base palindromes,數字585,不論是十進位表示法還是二進位表示法1001001001,都是迴文數(palindromic )。小於一百萬的數字中,找出符合此條件的數字,求總和。

判斷十進位、二進位是不是迴文數。

def is_palindrome_10(n):
    s = str(n)
    return s == s[::-1]
   
def is_palindrome_2(n):
    s = bin(n)[2:]
    return s == s[::-1]

   
def is_palindrome(n):
    return is_palindrome_10(n) and is_palindrome_2(n)


問題本身,暴力法,底下函式傳入1000000,可算出答案。
   
def dbp(max):
    return sum(i for i in range(1, max) if is_palindrome(i))


Problem 37: Truncatable primes,3797若從左邊逐一拿掉一個位數,會變成797、97、7,通通都是質數,從右邊逐一拿掉也是。這種數字總共只有11個,不包括2、3、5、7,請找出這些數字求總和。

首先是判斷質數的函式。

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_truncatable_left(n):
    s = str(n)
    s_len = len(s)
    for i in range(1, s_len):
        x = s[i:]
        if not is_prime(int(x)):
            return False
    return True
       
def is_truncatable_right(n):
    s = str(n)
    s_len = len(s)
    for i in range(1, s_len):
        x = s[:-i]
        if not is_prime(int(x)):
            return False
    return True


為了簡便,把以上的檢查通通放在一個函式內。

def check_all(n):
    return is_prime(n) and is_truncatable_left(n) and is_truncatable_right(n)


最後是問題本身,因已知道共有11個,所以使用內建程式庫itertools的count,從10開始,判斷是否符合條件。底下函式傳入11,可求出答案。
 
import itertools

def tp(cnt_max):
    result = 0
    cnt = 0
    for n in itertools.count(10):
        if check_all(n):
            result += n
            cnt += 1
            if cnt == cnt_max:
                break
    return result

   

Problem 38: Pandigital multiples,以192為例,乘上1、2、3後得到192、384、576,此三數連接起來得到192384576,稱為連接乘積,恰好是9位數,而且剛好含有1到9。請問連接乘積中最大的是哪個。

首先撰寫函式判斷某數字是否為9位數、而且剛好含有1到9。

def is_pandigital(num):
    ds = '123456789'
    s = str(num)
    ins = sum(1 for d in ds if d in s)
    return ins == len(ds)


然後,底下函式給定某數字參數num,會乘上1、2、3、等等,然後連接,若位數是9個才回傳,若不是則回傳False。回傳值有兩個,前一個是連接乘積,後一個是最大乘上了多少,譬如192最大需乘上3。
           
def prod(num):
    s = str(num)
    n = 2
    while True:
        s += str(num * n)
        if len(s) == 9:
            return s, n
        elif len(s) > 9:
            return False, False
        n += 1
    return False, False


最後是問題本身,迭代1到4位數,因為5位數以上的連接乘積不會是9位數。然後以暴力法迭代每個數字,算出連接乘積,判斷是否符合條件,找出最大的。呼叫底下函式即可算出答案。

def pm():
    num_max = 0
    for i in range(1, 4+1):
        for num in range(10**(i-1), 10**i - 1):
            s, n = prod(num)
            if s and is_pandigital(s) and int(s) > num_max:
                num_max = int(s)
    return num_max


Problem 39: Integer right triangles,給定直角三角形周長p為120的話,可算出三組整數邊長的解{20,48,52}、{24,45,51}、{30,40,50}。在周長小於1000(包含)的情況下,請問哪個周長的直角三角形的整數邊長解最多。

用暴力法,會超過1分鐘,必須進行分析,縮減需要處理的計算量。

def irt(p_upperbound):
    p_max = 0
    t_max = 0
    for p in range(p_upperbound//4*2, p_upperbound+1, 2):
        t = 0
        for a in range(2, int(p / (2**0.5+1)) + 1):
            if p * (p-2*a) % (2* (p-a)) == 0:

                t += 1
        if t > t_max:
            t_max = t
            p_max = p
   
    return p_max


Problem 40: Champernowne's constant,把從1開始的數字連接起來,變成123456789101112131415161...,請問第1、10、100、1000、10000、100000、1000000個位數的乘積為何。

把從1開始的數字不斷連接起來,直到超過給定長度dn。

def d(dn):
    li = []
    i = 1
    n = 0
    while True:
        s = str(i)
        li.append(s)
        n += len(s)
        if n > dn:
            break
        i += 1
    return ''.join(li)


然後,拿出想要的位數,相乘。底下函式傳入7,可算出答案。

def cc(pow):
    s = d(10**pow)
    result = 1
    for i in range(pow):
        result *= int(s[10**i-1])
    return result

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

Project Euler與Python:Problem 11~20與67

Project Euler針對數學和程式愛好者,設計了一系列的問題,至今已有五百多道,供大家挑戰,每一題都要求應在1分鐘內計算出答案。

我將以Python語言來解題,原始碼放在GitHub(yehnan/project_euler_python),基本上都是先採用暴力法,若運算量太大,無法在1分鐘內算出答案,才著手研究想出更快的演算法。當題目的數字很小時,隨便寫就能在1分鐘內算出答案,但是,也就只有幾題而已,大部分的題目都需要大量計算,可沒那麼好混。

我使用硬體Raspberry Pi 2與軟體Raspbian,作為執行速度的依據。

注意:提供詳細的說明與分析。

Problem 11: Largest product in a grid,20×20格子,每格有個數字,把相鄰的四個數相乘,找出乘積最大的。所謂相鄰是指同一方向的連續四個數,方向有上下、左右、左上右下的對角線、右上左下的對角線。

把題目給的20×20格子,放進檔案裡,然後讀取放進二維串列。

from io import open
grid = []
with open('p011_data.txt', 'r', encoding='ascii') as fin:
    for line in fin:
        grid.append(list(map(int, line.split())))

給定含有數字的串列,算出乘積。

def product(li):
    result = 1
    for n in li:
        result *= n
    return result

給定串列、個數count,找出相鄰count個數字相乘後乘積最大的。

def product_max(li, count):
    max = 0
    for i in range(0, len(li)-count + 1):
        m = product(li[i:i+count])
        if m > max:
            max = m
    return max

底下函式可拿出左上右下或右上左下的對角線,由參數slant決定。

def diag(grid, row, col, slant):
    d = []
    while True:
        try:
            d.append(grid[row][col])
            row += 1
            col += slant
        except IndexError:
            break
    return d

終於到了問題本身,底下函式可找出相鄰四個數乘積最大的,但只會尋找一半:每一列、右上半部對角線為左上右下、左上半部對角線為右上左下。

def find_max(grid, count):
    max = 0
    for row in grid:
        m = product_max(row, count)
        if m > max:
            max = m
    for col in range(len(grid[0])):
        m = product_max(diag(grid, 0, col, 1), count)
        if m > max:
            max = m
    for col in range(len(grid[0])):
        m = product_max(diag(grid, 0, col, -1), count)
        if m > max:
            max = m
    return max

不過,只要把格子旋轉90度,再丟進上面的函式,便可處理另外一半。

def lpiad(grid, count):
    m = find_max(grid, count)
    grid_rotated = list(zip(*grid))
    mr = find_max(grid_rotated, count)
    return max(m, mr)

Problem 12: Highly divisible triangular number,找出第一個擁有超過500個因數的三角數(triangle number)。

暴力法太慢,應該先作質因數分解,根據質因數的次方數,便可算出因數的個數。不過此處並未作質因數分解,而是作2與奇數的分解。底下函式,給定參數n,求出n的因數個數。

def divisors_count(n):
    dc = 1
    cnt = 0
    while n % 2 == 0:
        n //= 2
        cnt += 1
    dc *= (cnt + 1)
   
    p = 3
    while n != 1:
        cnt = 0
        while n % p == 0:
            n //= p
            cnt += 1
        dc *= (cnt + 1)
        p += 2
    return dc

下面這個產生器函式,會不斷給出下一個三角數。
   
def triangle_number_gf():
    tn = 0
    i = 1
    while True:
        tn += i
        i += 1
        yield tn

最後是問題本身,傳入參數500,便可算出答案。

def hdtn(dc):
    for tn in triangle_number_gf():
        if divisors_count(tn) > dc:
            return tn

Problem 13: Large sum,給定100個數字,皆為50位數,算出總和後,取出前10個位數。

因為只要記憶體夠大,Python的整數便能不斷變大,所以這題對Python來說很簡單。把給定的東西放在檔案裡,呼叫底下函式傳入檔名,即可算出答案。

from io import open

def ls(filename):
    with open(filename, 'r', encoding='ascii') as fin:
        total = sum(int(line) for line in fin)
        return str(total)[:10]


Problem 14: Longest Collatz sequence,根據奇偶歸一猜想,每個正整數經過數次Collatz轉換後,最終會得到1。在小於一百萬的數字中,哪個數需要最多次轉換才能到1。

首先,底下函式給定參數n(起始數),回傳Collatz轉換的次數。

memo = {1: 1}
def collatz(n):
    if n not in memo:
        if n % 2 == 0:
            memo[n] = 1 + collatz(n // 2)
        else:
            memo[n] = 1 + collatz(3 * n + 1)
    return memo[n]


然後,以迴圈迭代所有起始數,找出轉換次數最多的那一個。底下函式傳入1000000即可求出答案。

def lcs(m):
    n = 1
    n_len = 1
    for i in range(2, m+1):
        i_len = collatz(i)
        if n_len < i_len:
            n = i
            n_len = i_len
    return n


Problem 15: Lattice paths,20×20的格子,從最左上走到最右下,只能向左或向下移動,請問有幾種路徑。

此題等同於:有40個位置,手上有20個紅球和20個綠球,請問共有幾種擺法。也就是40! / 20! / 20!。所以先寫階乘函式。

def factorial(n):
    result = 1
    for i in range(2, n+1):
        result *= i
    return result
 

然後是問題本身,底下函式傳入20與20,即可求出答案。

def lp(row, col):
    return factorial(row+col) // factorial(row) // factorial(col)


Problem 16: Power digit sum,求2**1000的所有位數的總和。

因Python支援無限大的整數(只要記憶體足夠的話),所以這題對Python來說很簡單。底下函式傳入1000即可求出答案。

def pds(n):
    num = 2 ** n
    result = 0
    for d in str(num):
        result += int(d)
    return result


Problem 17: Number letter counts,把1到1000寫成英文,總共有幾個字母。

嗯,只能根據規則把數字轉成英文字串,底下函式只能處理1到1000的情況。

d3 = {1: 'one thousand'}
d2a = { 10: 'ten', 11: 'eleven', 12: 'twelve', 13: 'thirteen',
        14: 'fourteen', 15: 'fifteen', 16: 'sixteen',
        17: 'seventeen', 18: 'eighteen', 19: 'nineteen' }
d2b = { 2: 'twenty', 3: 'thirty', 4: 'forty', 5: 'fifty',
        6: 'sixty', 7: 'seventy', 8: 'eighty', 9: 'ninety' }
d1 = { 1: 'one', 2: 'two', 3: 'three', 4: 'four', 5: 'five',
       6: 'six', 7: 'seven', 8: 'eight', 9: 'nine' }

def i2s(i): # from 1 to 1000
    a4 = i // 1000
    a3 = (i % 1000) // 100
    a2 = (i % 100) // 10
    a1 = i % 10
    result = []
    if a4 == 1:
        result.append('one thousand')
        if a3 != 0 or a2 != 0 or a1 != 0:
            result.append(' and ')
    if a3 != 0:
        result.append(d1[a3])
        result.append(' hundred')
        if a2 != 0 or a1 != 0:
            result.append(' and ')
    if a2 == 1:
        result.append(d2a[a2*10 + a1])
    elif a2 >= 2:
        result.append(d2b[a2])
        if a1 != 0:
            result.append('-')
    if a2 != 1 and a1 != 0:
        result.append(d1[a1])
    return ''.join(result)


然後是問題本身,底下函式傳入1與1000,即可算出答案。
   
def nlc(n, m=None):
    if m is None:
        m = n
    result = 0
    for i in range(n, m+1):
        s = i2s(i)
        s = s.replace(' ', '').replace('-', '')
        result += len(s)
    return result


Problem 18: Maximum path sum IProblem 67: Maximum path sum II,給定呈現三角形的數字層,從頂端開始走,往下一層走時可往左或往右,走到底部時形成一條路徑,把走過的數字加總,請問路徑加總和最大的是多少。

首先寫函式把代表三角形的資料從檔案載入,變成多維串列。

from io import open

def get_triangle(filename):
    triangle = []
    with open(filename, 'r', encoding='ascii') as fin:
        for line in fin:
            triangle.append([int(n) for n in line.split()])
    return triangle


然後是尋找最大路徑加總和的函式,屬於從上往下的遞迴形式。當你站在某一個位置時,只要算出左下三角形的最大路徑加總和、以及右下三角形的最大路徑加總和,挑出比較大的那一條,便可決定往哪走。

def max_path_sum(tri, row, col):
    if row == len(tri)-1:
        return tri[row][col]
    else:
        sub = max(max_path_sum(tri, row+1, col),
                  max_path_sum(tri, row+1, col+1))

        return tri[row][col] + sub


最後是題目本身,傳入資料檔名,即可求出答案。
       
def mpsi(filename):
    tri = get_triangle(filename)
    return max_path_sum(tri, 0, 0)


不過,當三角形層數一多(Problem 67),上述寫法就太慢了,必須由下往上,從倒數第二層開始,位置(r, c)的最大路徑加總和,會是(r, c) + (r+1, c)與(r, c) + (r+1, c+1)兩個之中較大的那一個。一層一層往上,最後便可算出頂端(0, 0)的最大路徑加總和。

def max_path_sum(tri, row, col):
    for r in range(len(tri)-1 -1, -1, -1):
        for c in range(len(tri[r])):
            tri[r][c] += max(tri[r+1][c], tri[r+1][c+1])
    return tri[row][col]


Problem 19: Counting Sundays,二十世紀每個月第一天是星期日的日子,有幾天。

我們知道1900年1月1日是星期一,定其天數為1;而1900年1月7日是星期日,天數為7,所以只要算出某日期(日月年)的天數,能被7整除的話即是星期日。

首先是判斷潤年的函式。

def is_leapyear(year):
    if year%4 == 0 and year%100 != 0 or year%400 == 0:
        return 1
    else:
        return 0


給定月份與年份,回傳該月份有幾天。
       
month = [31, 28, 31, 30, 31, 30,
         31, 31, 30, 31, 30, 31]
        
def days_of_month(m, y):
    return month[m-1] + (is_leapyear(y) if m == 2 else 0)


給定年份,回傳該年份有幾天。

def days_of_year(y):
    return sum(month) + is_leapyear(y)


給定日期(日月年),算出其天數。

def date_to_days(date):
    dy, mn, yr = date
    days = dy
    for y in range(1900, yr):
        days += days_of_year(y)
    for m in range(1, mn):
        days += days_of_month(m, yr)
    return days


能被7整除,則為星期日。

def is_sunday(days):
    return days % 7 == 0


迭代所有年份的每月第一天。
 
def cs():
    count = 0
    for y in range(1901, 2000+1):
        for m in range(1, 12+1):
            days = date_to_days((1, m, y))
            if is_sunday(days):
                count += 1
    return count


Problem 20: Factorial digit sum,加總100階乘的所有位數。

沒什麼好說的,算出階乘,然後把每個位數加起來。

def fact(n):
    result = 1
    for i in range(1, n):
        result *= i
    return result


底下函式傳入100,便可求得答案。

def fds(n):
    s = str(fact(n))
    result = 0
    for c in s:
        result += int(c)
    return result


Project Euler與Python:Problem 1~10

Project Euler針對數學和程式愛好者,設計了一系列的問題,至今已有五百多道,供大家挑戰,每一題都要求應在1分鐘內計算出答案。

我使用硬體Raspberry Pi 2與軟體Raspbian,作為執行速度的依據。

十八世紀數學家Leonhard Euler。

根據目前的統計,解出25~49題的人數將近四萬(38392),而解出超過450題的人數僅有45。
各種賞心悅目的成就徽章。
可先上去看看有哪些題目,若想開始解題,先註冊個帳號,然後便可輸入答案。每道問題經過設計,都以某數字作為明確的答案,在題目頁面的欄位「Answer:」裡提交,便可檢查是否正確解題。
若答案正確,可進入該問題的論壇,看看別的語言的解法,觀看其他人的討論文章。當然啦,網路上也可找到他人公布分享的程式原始碼,但如果直接照抄就失去意義了。
接下來試試第1~10題吧,我將以Python語言來解題,原始碼放在GitHub(yehnan/project_euler_python),基本上我都是先採用暴力法,若運算量太大,無法在1分鐘內算出答案,才著手研究想出更快的演算法。當題目的數字很小時,隨便寫就能在1分鐘內算出答案,但是,也就只有幾題而已,大部分的題目都需要大量計算,可沒那麼好混。

注意:並提供詳細的說明與分析。

Problem 1: Multiples of 3 and 5,小於1000的正整數中,找出3或5的倍數,相加求總和。

底下函式直接使用暴力法,傳入參數1000可算出此題的答案。

def m35(n):
    result = 0
    for i in range(1, n):
        if i % 3 == 0 or i % 5 == 0:
            result += i
    return result


也可縮成一行,其中運用了產生器運算式(generator expression)。

def m35(n):
    return sum(i for i in range(1, n) if i%3==0 or i%5==0)


Problem 2: Even Fibonacci numbers,找出小於四百萬而且是偶數的費波那契數,相加求總和。

首先是計算費波那契數(fibonacci number)的函式,採取遞迴形式,但會記憶算過的值,免得重複計算。

memo = {0: 0, 1: 1}
def fib(n):
    if n not in memo:
        memo[n] = fib(n-1) + fib(n-2)
    return memo[n]

然後是題目本身,使用暴力法,從小到大算出費波那契數,小於四百萬而且是偶數的話就加起來。底下函式傳入參數4000000,即可得到答案。

def efn(max):
    result = 0
    n = 1
    fn = fib(n)
    while fn <= max:
        if fn % 2 == 0:
            result += fn
        n += 1
        fn = fib(n)
    return result

Problem 3: Largest prime factor,找出600851475143最大的質因數。

首先是能從小到大、不斷給出質數的函式,底下prime_gf是產生器函式,prime_f則是產生器。

def prime_gf():
    primes = [2]
    def is_prime(n):
        for i in primes:
            if n % i == 0:
                return False
        return True
    while True:
        p = primes[-1]
        yield p
        while True:
            p += 1
            if is_prime(p):
                primes.append(p)
                break
prime_f = prime_gf()


然後是題目本身,底下函式傳入參數600851475143,即可得到答案。

def lpf(n):
    max = 1
    for i in prime_f:
        if i <= n:
            while n % i == 0:
                n /= i
                max = i
        else:
            break
    return max


Problem 4: Largest palindrome product,迴文數,譬如9009,從頭端看過來、與從尾端看過去都一樣,而且9009是兩個2位數相乘後的迴文數中最大的(91 × 99),請問兩個3位數相乘後最大的迴文數是多少。

首先,底下函式可判斷是否為迴文數。

def is_palindrome(n):
    ns = str(n)
    a, b = 0, len(ns)-1
    while a < b:
        if ns[a] != ns[b]:
            return False
        else:
            a += 1
            b -= 1
    return True


然後,底下函式的參數d代表幾位數,以暴力法列出所有的兩個d位數組合,相乘後判斷是不是迴文數,通通記錄起來。

def palindromes(d):
    result = []
    for x in range(10**d-1, 10**(d-1)-1, -1):
        for y in range(x, 10**(d-1)-1, -1):
            if is_palindrome(x * y):
                result.append((x*y, x, y))
    return result


最後進行排序,便可得到最大的迴文數。底下函式傳入參數3可得到答案。

def palindromes_max(d):
    return sorted(palindromes(d), key=lambda x: x[0])[-1][0]


Problem 5: Smallest multiple,能夠被1到20所有數字整除的正整數裡,最小是多少。

首先是計算最大公因數與最小公倍數的函式。

def gcd(a, b):
    while b:
        a, b = b, a%b
    return a

def lcm(a, b):
    return a * b // gcd(a, b)


然後是題目本身,一次處理兩個數,不斷算出兩個數的最小公倍數,1跟2的最小公倍數是2,2跟3的最小公倍數是6,6跟4的最小公倍數是12,依此類推。底下函式傳入參數1與20,即可算出答案。

def sm(x, y):
    m = 1
    for n in range(x, y+1):
        m = lcm(m, n)
    return m


Problem 6: Sum square difference,前100個自然數(正整數)的平方和與和平方的差值。

呃,沒什麼好說的,數字小、計算量少,直接暴力法即可。底下函式傳入參數100可得到答案。

def ssd(n):
    sum_sq = sum(i**2 for i in range(1, n+1))
    sq_sum = sum(range(1, n+1)) ** 2
    return sq_sum - sum_sq


Problem 7: 10001st prime,求第10001個質數。

底下函式可判斷參數n是否為質數,檢查小於n的數之中,是否有可以整除n的數存在,但不需要檢查「小於n」的數,只需檢查「小於n開根號+1」的數即可。

def is_prime(n):
    if n % 2 == 0:
        return False
       
    p = 3
    n_sqrt = n**0.5 + 1
    while p < n_sqrt:
        if n % p == 0:
            return False
        p += 2
    return True

然後是問題本身,從3開始,逐一判斷下一個(+2)是否為質數,直到目標。底下函式傳入參數10001便可算出答案。

def nth_prime(n):
    prime = 2
    count = 1
    p_next = 3
    while count < n:
        if is_prime(p_next):
            prime = p_next
            count += 1
        p_next += 2
    return prime

Problem 8: Largest product in a series,給定一個有100位數的數字,找出其中相鄰的13個位數、位數的乘積最大。

暴力法,逐一算出乘積,比較哪個大。

def lpias(num, n):
    num_str = str(num)
    num_len = len(num_str)
   
    def prod(ns):
        p = 1
        for i in ns:
            p *= int(i)
        return p
       
    max = 1
    for i in range(0, num_len - n + 1):
        m = prod(num_str[i:i+n])
        if m > max:
            max = m
    return max


Problem 9: Special Pythagorean triplet,只有一組畢氏三元數組a、b、c能使得a+b+c等於1000,請算出a*b*c。

首先,判斷a、b、c是否為畢氏三元數組。

def is_pythagorean_triplet(a, b, c):
    return a**2 + b**2 == c**2


然後,底下產生器函式可產生出a + b + c = n而且a < b < c的三數組。注意,其中range的參數並非最佳化。

def abc_gf(n):
    for c in range(n//3, n):
        for b in range((n-c)//2, n-c):
            a = n - c - b
            if a < b < c:
                yield a, b, c


最後是問題本身,底下函式傳入參數1000,即可找出答案。

def spr(n):
    for a, b, c in abc_gf(n):
        if is_pythagorean_triplet(a, b, c):
            return a*b*c


Problem 10: Summation of primes,找出小於兩百萬的所有質數,加總求和。

使用埃拉托斯特尼篩法找出小於n的質數。

def prime_sieve(n):
    primes = set(range(2, n+1))
    for i in range(2, (n+1+1) // 2):
        if i in primes:
            m = 2
            while i*m <= n:
                primes.discard(i*m)
                m += 1
    return primes


然後加總即可。底下函式傳入參數2000000即可算出答案。

def sop(n):
    return sum(prime_sieve(n))


參考資料: