組み合せで表現可能な数字の和

30分プログラム、その293。昨日の続きで、『組み合わせで表現可能な数の和』を計算するコード。ハッシュを使ったりして高速化をはかるも挫折。
明日は方針を変えて、Xが過剰数のときN-Xが過剰数であるかをチェックする方法でやってみよう。

使い方

>>> combination(xrange(1,100))
19700

ソースコード

#! /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/28 22:10:13
#
# 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

## 素因数分解を利用した約数の計算
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,primes):
    xs = factor(n,primes)
    return product(map(lambda (p,a):(p**(a+1)-1)/(p-1),xs))

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

def combination(xs):
    dic = {}
    sum = 0
    for x in xs:
        for y in xs:
            if x+y > 28123:
                break
            elif dic.get(x+y) == None:
                sum += x+y
                dic[x+y] = True
    return sum


def main():
    primes = sieve(28123)
    xs = xrange(1,28123)
    y = sum(filter(lambda x:is_abudant(x,primes),xs))
    return sum(xs) - y

参考

• 過去の30分プログラム