题目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多分钟。。。

1. from numpy import ones
   
2. from math import gcd
   
3. from numba import jit
   
4. primes = [True]*120000
   
5. primes[:2] = [False]*2
   
6. for i, prime in enumerate(primes):
   
7.         if prime:
   
8.                 for j in range(i*i, 120000, i):
   
9.                         primes[j] = False
   
10. plist = [k for k,trueprime in enumerate(primes) if trueprime]
    
11. matrix = ones((120000, 1), dtype='int64')
    
12. for p in plist:
    
13.         matrix[p::p] *= p
    
14. @jit
    
15. def solve(limit):
    
16.         s = 0
    
17.         for c in range(3, limit):
    
18.                 if matrix[c,0] < c:
    
19.                         for a in range(1, c//2+1):
    
20.                                 if matrix[a,0]*matrix[c-a,0]*matrix[c,0] < c:
    
21.                                         if gcd(a, c) == 1 and gcd(a, c-a) == 1:
    
22.                                                 s += c
    
23.         return s
    
24. print(solve(120000))

复制代码

18407904
[Finished in 1433.2s]

l910607

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

1. import time
   
2. t0 = time.time()
   
3. 
   
4. def iscoprime(b, a):    # 检查互质
   
5.     yushu = b % a
   
6.     while yushu != 0:
   
7.         b = a
   
8.         a = yushu
   
9.         yushu = b % a
   
10.     if a == 1:
    
11.         return True
    
12.     else:
    
13.         return False
    
14. 
    
15. 
    
16. MAXC = 120000
    
17. rad = [0]*MAXC
    
18. rad[1] = 1
    
19. for i in range(2, MAXC):    # 先计算rad值
    
20.     if rad[i] > 0:
    
21.         continue
    
22.     for j in range(i, MAXC, i):
    
23.         if rad[j] == 0:
    
24.             rad[j] = i
    
25.         else:
    
26.             rad[j] *= i
    
27. DicRADN = {}
    
28. countc, sumc = 0, 0
    
29. for i in range(1, MAXC):    # 以rad值为key，原数字为value，生成字典
    
30.     try:
    
31.         DicRADN[rad[i]].append(i)
    
32.     except KeyError:
    
33.         DicRADN[rad[i]] = [i]
    
34. 
    
35. for c in range(3, MAXC):
    
36.     maxradab = c//rad[c]    # rad(a)*rad(b)的最大值
    
37.     if maxradab == 1:       # 这种情况只能是rada=radb=1，即a=b=1，不合题意
    
38.         continue
    
39.     for rada in DicRADN.keys():
    
40.         if rada > maxradab//2:  # 因为b一定至少有一个质因子是2次方以上，所以radb≥2
    
41.             break
    
42.         for a in DicRADN[rada]:
    
43.             if a > (c-1)//2:    # 因为a<b,所以a<c/2
    
44.                 break
    
45.             b = c-a
    
46.             if rad[b] == b:     # b为质数的一次方的乘积。易证不合题意。
    
47.                 continue
    
48.             if rad[a]*rad[b] >= maxradab:
    
49.                 continue
    
50.             if iscoprime(b, a):
    
51.                 countc += 1
    
52.                 sumc += c
    
53. 
    
54. print("Answer: [数量, 求和]:", countc, sumc)
    
55. print("Start time: ", time.asctime(time.localtime(t0)))
    
56. print("Running time: ", time.time()-t0, "s")
    

复制代码

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