题目127：考察在给定大小以下的abc组的数量

欧拉计划

abc-hits

The radical of n, rad(n), is the product of distinct prime factors of n. For example, 504 = 2³ × 3² × 7, so rad(504) = 2 × 3 × 7 = 42.

We shall define the triplet of positive integers (a, b, c) to be an abc-hit if:

• GCD(a, b) = GCD(a, c) = GCD(b, c) = 1
• a < b
• a + b = c
• rad(abc) < c
  

For example, (5, 27, 32) is an abc-hit, because:

• GCD(5, 27) = GCD(5, 32) = GCD(27, 32) = 1
• 5 < 27
• 5 + 27 = 32
• rad(4320) = 30 < 32
  

It turns out that abc-hits are quite rare and there are only thirty-one abc-hits for c < 1000, with ∑c = 12523.

Find ∑c for c < 120000.

题目：

n 的根函数，rad(n)，定义为其所有不同质因子之积。例如，504 = 23 × 32 × 7， 所以 rad(504) = 2 × 3 × 7 = 42。

如果正整数三元组 (a, b, c) 满足以下条件，我们将其称为一个 abc 组：

• GCD(a, b) = GCD(a, c) = GCD(b, c) = 1
• a < b
• a + b = c
• rad(abc) < c
  

例如，(5, 27, 32) 是一个 abc 组，因为：

• GCD(5, 27) = GCD(5, 32) = GCD(27, 32) = 1
• 5 < 27
• 5 + 27 = 32
• rad(4320) = 30 < 32
  

事实上 abc 组非常的稀少，对于 c < 1000 只有 31 个 abc 组，∑c = 12523。

对于 c < 120000，求 ∑c 。

jerryxjr1220

跑了20多分钟。。。

from numpy import ones
from math import gcd
from numba import jit
primes = [True]*120000
primes[:2] = [False]*2
for i, prime in enumerate(primes):
        if prime:
                for j in range(i*i, 120000, i):
                        primes[j] = False
plist = [k for k,trueprime in enumerate(primes) if trueprime]
matrix = ones((120000, 1), dtype='int64')
for p in plist:
        matrix[p::p] *= p
@jit
def solve(limit):
        s = 0
        for c in range(3, limit):
                if matrix[c,0] < c:
                        for a in range(1, c//2+1):
                                if matrix[a,0]*matrix[c-a,0]*matrix[c,0] < c:
                                        if gcd(a, c) == 1 and gcd(a, c-a) == 1:
                                                s += c
        return s
print(solve(120000))

18407904
[Finished in 1433.2s]

l910607

新人回陈贴，希望没违规。
跑了5s多，但比最初的780多秒还是要好多了。

import time
t0 = time.time()

def iscoprime(b, a):    # 检查互质
    yushu = b % a
    while yushu != 0:
        b = a
        a = yushu
        yushu = b % a
    if a == 1:
        return True
    else:
        return False


MAXC = 120000
rad = [0]*MAXC
rad[1] = 1
for i in range(2, MAXC):    # 先计算rad值
    if rad[i] > 0:
        continue
    for j in range(i, MAXC, i):
        if rad[j] == 0:
            rad[j] = i
        else:
            rad[j] *= i
DicRADN = {}
countc, sumc = 0, 0
for i in range(1, MAXC):    # 以rad值为key，原数字为value，生成字典
    try:
        DicRADN[rad[i]].append(i)
    except KeyError:
        DicRADN[rad[i]] = [i]

for c in range(3, MAXC):
    maxradab = c//rad[c]    # rad(a)*rad(b)的最大值
    if maxradab == 1:       # 这种情况只能是rada=radb=1，即a=b=1，不合题意
        continue
    for rada in DicRADN.keys():
        if rada > maxradab//2:  # 因为b一定至少有一个质因子是2次方以上，所以radb≥2
            break
        for a in DicRADN[rada]:
            if a > (c-1)//2:    # 因为a<b,所以a<c/2
                break
            b = c-a
            if rad[b] == b:     # b为质数的一次方的乘积。易证不合题意。
                continue
            if rad[a]*rad[b] >= maxradab:
                continue
            if iscoprime(b, a):
                countc += 1
                sumc += c

print("Answer: [数量, 求和]:", countc, sumc)
print("Start time: ", time.asctime(time.localtime(t0)))
print("Running time: ", time.time()-t0, "s")

Answer: [数量, 求和]: 456 18407904
Running time:  5.613126516342163 s