プロジェクトオイラー
http://projecteuler.net/
Q245.
(n - φ(n)) / (n - 1)を約分して分子が1になる、合成数の1 < n ≤ 2×1011の総和を求めよ。
最難問題だそうだが、さほどてこずらなかった。
まず、nを素因数分解してべき乗の項があるとダメなのはすぐにわかる。また、偶数もダメ。nが2つの素数の積であるとき、
n = p1p2
(p1 + p2 - 1) / (p1p2 - 1) = 1 / m
(p1 - m)(p2 - m) = m2 - m + 1
m2 - m + 1の素因数分解は、n2 + 1のときと同じふるい的な手法がそのまま使えるので簡単。
しかし、3つ以上の素数から成る合成数ではこれが使えない。そこで、例えばp1 = 3、p2 = 5とすると、
(15p3 - 8(p3 - 1)) / (15p3 - 1) = (7p3 + 8) / (15p3 - 1) = 1 / m
となるから、mは2しか取りえず、そのときp3 = 17となる。
これを3個以上の素数から成る合成数で使った。
from itertools import count
from math import sqrtdef is_prime(n):
for p in primes:
if p * p > n:
return True
elif n % p == 0:
return False
return Truedef make_primes(n):
m = int(sqrt(n + 0.5))
for p in xrange(3, m + 1, 2):
if is_prime(p):
primes.append(p)
a = [ 0 ] * ((n + 15) / 16)
for p in primes:
if p * p > n:
break
if p == 2:
continue
for k in xrange(3 * p, n, 2 * p):
a[k>>4] |= 1 << (k & 15)
for k in xrange((m + 1) / 2 * 2 + 1, n + 1, 2):
if (a[k>>4] & (1 << (k & 15))) == 0:
primes.append(k)def is_prime2(n):
if n <= primes[-1]:
return n in set_primes
else:
return is_prime(n)def gen_divisor(f):
n = 1
for k in xrange(1, len(f), 2):
n *= f[k] + 1
for m in xrange( (n + 1) / 2):
d = 1
for k in xrange(0, len(f), 2):
e = m % (f[k+1] + 1)
m /= f[k+1] + 1
d *= f[k] ** e
yield ddef gen_phi(n0, l, k0 = 1):
for k in xrange(k0, len(primes)):
p = primes[k]
if l == 1:
if p > n0:
break
yield (p, (p - 1))
else:
if p ** l > n0:
break
n = n0 / p
if n > primes[k]:
for t in gen_phi(n, l - 1, k + 1):
yield (p * t[0], (p - 1) * t[1])def calc_num_primes_max(n):
for k in count():
n /= primes[k]
if n == 0:
return kdef sum_unit_fractions(n):
if n == 2:
return sum_unit_fractions_2()
else:
return sum_unit_fractions_other(n)def sum_unit_fractions_2():
M = int(sqrt(N)) / 2
a = [ m * (m - 1) + 1 for m in xrange(M + 1) ]
b = [ () for m in xrange(M + 1) ]
for m in xrange(1, M / 2 + 1):
p = a[m]
if p > 1:
if p > M:
g = xrange(m, m + 1)
else:
g = xrange(m, M + 1, p)
for k in g:
e = 1
while True:
a[k] /= p
if a[k] % p != 0:
break
e += 1
b[k] += (p, e)
for m in xrange(M / 2 + 1, M + 1):
p = a[m]
if p > 1:
b[m] += (p, 1)
s = 0
for m in xrange(2, M + 1, 2):
for d1 in gen_divisor(b[m]):
d2 = (m * (m - 1) + 1) / d1
p1 = d1 + m
p2 = d2 + m
n = p1 * p2
if n <= N and is_prime2(p1) and is_prime2(p2):
s += n
return sdef sum_unit_fractions_other(nprimes, n = 1, phi = 1, k0 = 1):
s = 0
if nprimes > 1:
for k in xrange(k0, len(primes)):
p = primes[k]
if n * p ** nprimes > N:
break
s += sum_unit_fractions_other(nprimes - 1,
n * p, phi * (p - 1), k + 1)
else:
p_min = primes[k0]
n_min = n * p_min
if n_min > N:
return 0
# (ap + b) / (cp - 1) = 1 / m
a = n - phi
b = phi
c = n
den = n_min - b
m_min = (n_min - 1 + den - 1) / den
m_max = (c - 1) / a
for m in xrange(m_min, m_max + 1):
if (b * m + 1) % (c - a * m) == 0:
p = (b * m + 1) / (c - a * m)
if p_min <= p and n * p <= N and is_prime2(p):
s += n * p
return sN = 10 ** 11 * 2
primes = [ 2 ]
make_primes(int(N ** (2./3)) / 2)
set_primes = set(primes)
print sum(map(sum_unit_fractions, range(2, calc_num_primes_max(N) + 1)))
