Problem23が解けた

30分プログラム、その294。やっと、Problem23が解けた。

完全数とは, その数の真の約数の和がそれ自身と一致する数のことである. たとえば, 28の真の約数の和は, 1 + 2 + 4 + 7 + 14 = 28であるので, 28は完全数である.

真の約数の和がその数よりも少ないものを不足数といい, 真の約数の和がその数よりも大きいものを過剰数と呼ぶ.

12は, 1+2+3+4+6=16となるので, 最小の過剰数である. よって2つの過剰数の和で書ける最少の数は24である. 数学的な解析により, 28123より大きい任意の整数は2つの過剰数の和で書けることが知られている. 2つの過剰数の和で表せない最大の数がこの上限よりも小さいことは分かっているのだが, この上限を減らすことが出来ていない.

2つの過剰数の和で書き表せない正の整数の総和を求めよ.

ふう、やっと解けた。mapやfilterを使ったリスト操作だとうまくいかなかったので、ディクショナリやwhileを組み合せて解いてみた。

使い方

$ python problem23.py
...
28121 1976 26145
28122 12 28110
28123 88 28035
4179871
python problem23.py  74.50s user 5.47s system 69% cpu 1:54.54 total

ソースコード

#! /usr/bin/python
# -*- mode:python; coding:utf-8 -*-
#
# problem23.py -
#
# Copyright(C) 2008 by mzp
# Author: MIZUNO Hiroki / mzpppp at gmail dot com
# http://howdyworld.org
#
# Timestamp: 2008/04/29 21:48:24
#
# This program is free software; you can redistribute it and/or
# modify it under MIT Lincence.
#

def sieve(max):
    (n,xs) = (2,[])
    while n < max:
        if all(map(lambda x: n % x != 0,xs)):
            xs.append(n)
        n += 1
    return xs
primes = []

## 素因数分解を利用した約数の計算
def factor(n,primes):
    def f(x,y,n):
        if x % y == 0:
            return f(x/y,y,n+1)
        else:
            return n
    xs = []
    for prime in primes:
        if n < prime:
            return xs
        elif n % prime == 0:
            xs.append((prime,f(n,prime,0)))
    return xs

def product(xs):
    return reduce(lambda x,y: x*y,xs,1)

def divisor_sum(n):
    xs = factor(n,primes)
    return product(map(lambda (p,a):(p**(a+1)-1)/(p-1),xs))

## 過剰数の計算
def is_abundant(n):
    return divisor_sum(n) > n*2

def main(n):
    global primes
    primes = sieve(n)
    d = {}
    for x in xrange(12,n):
        if is_abundant(x):
            d[x] = True
    print d
    sum = 0
    for target in xrange(1,n):
        for x in d.keys():