题目127：考察在给定大小以下的abc组的数量
abc-hitsThe radical of n, rad(n), is the product of distinct prime factors of n. For example, 504 = 23 × 32 × 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 于 2017-7-17 17:04 编辑
跑了20多分钟。。。
from numpy import ones
from math import gcd
from numba import jit
primes = *120000
primes[:2] = *2
for i, prime in enumerate(primes):
if prime:
for j in range(i*i, 120000, i):
primes = False
plist =
matrix = ones((120000, 1), dtype='int64')
for p in plist:
matrix *= p
@jit
def solve(limit):
s = 0
for c in range(3, limit):
if matrix < c:
for a in range(1, c//2+1):
if matrix*matrix*matrix < c:
if gcd(a, c) == 1 and gcd(a, c-a) == 1:
s += c
return s
print(solve(120000))
18407904
新人回陈贴，希望没违规。
跑了5s多，但比最初的780多秒还是要好多了。{:10_250:}
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 = *MAXC
rad = 1
for i in range(2, MAXC): # 先计算rad值
if rad > 0:
continue
for j in range(i, MAXC, i):
if rad == 0:
rad = i
else:
rad *= i
DicRADN = {}
countc, sumc = 0, 0
for i in range(1, MAXC): # 以rad值为key，原数字为value，生成字典
try:
DicRADN].append(i)
except KeyError:
DicRADN] =
for c in range(3, MAXC):
maxradab = c//rad # 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:
if a > (c-1)//2: # 因为a<b,所以a<c/2
break
b = c-a
if rad == b: # b为质数的一次方的乘积。易证不合题意。
continue
if rad*rad >= 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
页:
[1]