計算某個範圍內的質數和的辦法

lt發表於2017-06-08

在解尤拉計劃606題時，遇到一個求質數立方和的問題，pe論壇中有人給出了一個不求質數而求立方和的演算法，感覺很神奇，查閱他提到的維基百科，有點明白了。

為簡單起見，先看第10題:質數求和，在此基礎上，立方和就好理解了。 enter image description here
我修改了原始碼，增加了中間結果輸出，基本思路是先求範圍內所有數的和，然後減去合數的和，用容斥定理去重複。

def P10(n):
    r = int(n**0.5)
    assert r*r <= n and (r+1)**2 > n
    V = [n//i for i in range(1,r+1)]
    print(V)
    V += list(range(V[-1]-1,0,-1))
    print(V)
    S = {i:i*(i+1)//2-1 for i in V} #等差數列求和
    print(S)
    for p in range(2,r+1):
        if S[p] > S[p-1]:  # p is prime
            sp = S[p-1]  # sum of primes smaller than p
            p2 = p*p
            for v in V:
                if v < p2: break
                S[v] -= p*(S[v//p] - sp)
            print(p,S)
    return S[n]

#輸出
>>>> P10(90)
[90, 45, 30, 22, 18, 15, 12, 11, 10]
[90, 45, 30, 22, 18, 15, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1]
{90: 4094, 45: 1034, 30: 464, 22: 252, 18: 170, 15: 119, 12: 77, 11: 65, 10: 54, 9: 44, 8: 35, 7: 27, 6: 20, 5: 14, 4: 9, 3: 5, 2: 2, 1: 0}
(2, {90: 2026, 45: 530, 30: 226, 22: 122, 18: 82, 15: 65, 12: 37, 11: 37, 10: 26, 9: 26, 8: 17, 7: 17, 6: 10, 5: 10, 4: 5, 3: 5, 2: 2, 1: 0})
(3, {90: 1354, 45: 341, 30: 154, 22: 77, 18: 58, 15: 41, 12: 28, 11: 28, 10: 17, 9: 17, 8: 17, 7: 17, 6: 10, 5: 10, 4: 5, 3: 5, 2: 2, 1: 0})
(5, {90: 1089, 45: 281, 30: 129, 22: 77, 18: 58, 15: 41, 12: 28, 11: 28, 10: 17, 9: 17, 8: 17, 7: 17, 6: 10, 5: 10, 4: 5, 3: 5, 2: 2, 1: 0})
(7, {90: 963, 45: 281, 30: 129, 22: 77, 18: 58, 15: 41, 12: 28, 11: 28, 10: 17, 9: 17, 8: 17, 7: 17, 6: 10, 5: 10, 4: 5, 3: 5, 2: 2, 1: 0})
963

參考資料 https://en.wikipedia.org/wiki/Prime-counting_function

相關文章