プロジェクトオイラー
http://projecteuler.net/
Q233.
原点が中心で半径Nの円周上の格子点の数をf(N)とする。f(N) = 420となるNの総和を求めよ。
大変だった。
Nを素因数分解したとき、4で割って1余る素因数しか出てこないとき、この円周は第1象限に互いに素な格子点を持ち、素因数の重複を除いた個数をmとすると、対称の位置を重複しないで数えると、2m-1個格子点を持つ。互いに素でない場合も数えると、例えば、4で割って1余る素因数が3つ出てくるとき、
N = kp1e1p2e2p3e3
第1象限の右下半分に
4e1e2e3 + 2(e2e3 + e3e1 + e1e2) + e1 + e2 + e3
個ある。だから、これが52と等しくなるようなe1,e2,e3の組を求めればよい。これはしらみつぶしに求められる。ただ、しらみつぶしの方法も工夫が必要。
それぞれの指数の組について、Nを超えない4で割って1余る素因数の組を求める。そうして求まったベースの数でNを割って上のkの組を求める。ただし、kは4で割って1余る素因数を持っていてはいけない。こうして和を求める。実際には素因数の組は求めず、和だけ求めていく。
from itertools import count, product, ifilter, permutations
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 num_points(t):
n = len(t)
s = 0
for t2 in product( (0, 1), repeat = n):
m = 0
p = 1
for e, b in zip(t, t2):
if b:
m += 1
p *= e
if m:
s += (1 << (m - 1)) * p
return sdef gen_div(m, n, k = 1):
if m == 1:
if n >= k:
yield (n, )
else:
for e in range(k, n / m + 1):
for t in gen_div(m - 1, n - e, e):
yield (e, ) + tdef gen_solution(n):
for m in count(1):
if m == 1:
yield (n, )
else:
b = False
for k in count(m):
b1 = False
for t in gen_div(m, k):
s = num_points(t)
if s <= n:
b1 = True
if s == n:
yield t
if b1:
b = True
else:
break
if not b:
breakdef is_prime(n):
for p in primes:
if p * p > n:
return True
if n % p == 0:
return False
return Truedef gen_primes(m = 1):
for n in count(m + 1):
if is_prime(n):
yield ndef calc_sum(te, n, k = 0, l = 0):
e = te[k]
s = 0
if k == len(te) - 1:
for i in count(l):
p = prime_mod_1(i)
p_pow = p ** e
if n >= p_pow:
s += sum_not_have_factor_mod_1(n / p_pow) * p_pow
else:
return s
else:
for i in count(l):
p = prime_mod_1(i)
p_pow = p ** e
if n >= p_pow:
m = calc_sum(te, n / p_pow, k + 1, i + 1)
s += m * p_pow
if m == 0:
return s
else:
return sdef prime_mod_1(n):
if n >= len(primes_mod_1):
if len(primes_mod_1) == 0:
p0 = 2
else:
p0 = primes_mod_1[-1]
g = ifilter(lambda p: p % 4 == 1, gen_primes(p0))
while n >= len(primes_mod_1):
p = g.next()
primes_mod_1.append(p)
return primes_mod_1[n]def value(t2, t):
return reduce(lambda x, y: x * y,
map(lambda e: e[0] ** e[1], zip(t2, t)))def gen_order(t):
s = set()
for t1 in permutations(t):
s.add(t1)
for t1 in s:
yield t1def has_no_factor_mod_1(n):
for p in primes:
if p * p > n:
break
if n % p == 0:
if p % 4 == 1:
return False
n /= p
while n % p == 0:
n /= p
if n == 1:
return True
else:
return n % 4 != 1def sum_not_have_factor_mod_1(n):
if n < 5:
return n * (n + 1) / 2
else:
m = len(ary_sum_not_have_factor_mod_1)
if n >= m:
s = ary_sum_not_have_factor_mod_1[m-1]
for k in xrange(m, n + 1):
if has_no_factor_mod_1(k):
s += k
ary_sum_not_have_factor_mod_1.append(s)
return ary_sum_not_have_factor_mod_1[n]N = 420
M = (N - 4) / 8
L = 10 ** 11
s = 0
primes = [ 2 ]
make_primes(int(sqrt(L + 0.5)))
primes_mod_1 = [ ]
ary_sum_not_have_factor_mod_1 = [ k * (k + 1) / 2 for k in range(5) ]
for t in gen_solution(M):
for t3 in gen_order(t):
s += calc_sum(t3, L)
print s
