http://projecteuler.net/index.php?section=problems&id=110
例えば解の個数が100より大きい最小の自然数を考えるとすると、重複を入れて201個以上が必要です。最も大きな素数を使う場合を考えると、各指数が1なので、35 = 243だから、
21・31・51・71・111
が考えられます。このとき、素数の数は5個です。素数の数を1〜5個として、指数を決めます。指数が実数なら、Lagrangeの未定定数法を使って、
g = Σk ek log pk - λ(Σk log(2ek+1) - log M)
M = 201
∂g/∂ek = log pk - 2λ / (2ek + 1) = 0
2ek + 1 = 2λ / log pk
これを、
Πk (2ek+1) = M
に代入して、
2mλm / Πk log pk = M
λ = (MΠk log pk)1/m / 2
ek = λ / log pk - 1 / 2
これでekを概算できます。そしてこの周りで計算を行います。指数が整数の範囲内ならこの概算より元の数が小さくなることはないので、それを手がかりにして指数の範囲を絞り込みます。
解の個数1億で計算すると、5秒ほどかかりました。きちんと書けばもっと速くなると思うのですが。
from itertools import * from math import log, exp N = 10 ** 8 M = N * 2 + 1 eps = 1e-10 def fst(x): return x[0] def snd(x): return x[1] def is_prime(n): return all(n % m != 0 for m in range(2, n)) def value_f(es): return reduce(lambda x, y: x * y[0] ** y[1], izip(primes, es), 1) def div_ceil(n, d): return (n + d - 1) / d def calc_in_float(m, n): ps = log_primes[:m] la = (n * reduce(lambda x, y: x * y, ps)) ** (1. / m) / 2 es = [ la / lp - 0.5 for lp in ps ] return sum(imap(lambda x: x[0] * x[1], izip(es, ps))) def calc_e(m, n): ps = log_primes[:m] la = (n * reduce(lambda x, y: x * y, ps)) ** (1. / m) / 2 return int(la / ps[-1] - 0.5) def calc_min(m, n = M, e_min = 1): def f(e): es, v = calc_min(m - 1, div_ceil(n, e * 2 + 1), e) return es + (e,), e * log_primes[m-1] + v def min_es(gen_e, min_x = 0): if min_x == 0: for e in gen_e: min_x = f(e) break min_es, min_v = min_x min_val = value_f(min_es) for e in gen_e: n1 = div_ceil(n, e * 2 + 1) v = e * log_primes[m-1] + calc_in_float(m - 1, n1) if v < min_v + eps: es, v = f(e) val = value_f(es) if val < min_val: min_es, min_v = es, v min_val = val else: break return min_es, min_v if m == 1: e = n / 2 return (e,), e * log_primes[0] p = primes[m-1] e_max = int((n ** (1. / m) - 1) / 2 + eps) if e_max == 1: es, v = f(1) else: e = calc_e(m, n) es, v = min_es(xrange(e, e_max + 1)) es, v = min_es(xrange(e - 1, e_min - 1, -1), (es, v)) return es, v nprimes = int(log(M) / log(3) + 1 - eps) primes = map(snd, takewhile(lambda x: fst(x) <= nprimes, izip(count(1), ifilter(is_prime, count(2))))) log_primes = map(log, primes) print primes a = [ calc_min(m, M) for m in range(1, nprimes + 1) ] a.sort(key = lambda x: x[1]) print min([ value_f(x[0]) for x in a if x[1] < a[0][1] + eps ])
