再度プログラマを目指す弁理士

Project Euler 9

ピタゴラス数(ピタゴラスの定理を満たす自然数)とは a < b < c で以下の式を満たす数の組である.

a2 + b2 = c2
例えば, 32 + 42 = 9 + 16 = 25 = 52 である.

a + b + c = 1000 となるピタゴラスの三つ組が一つだけ存在する.
これらの積 abc を計算しなさい.
project euler 9

MyCode

プロジェクトオイラー1問目

プロジェクトオイラーの1問目を解いてみた。

■ワンライナー+できるだけ文字数少なく
print sum([x for x in range(1000)if x%3*x%5==0])

説明
if の部分は、((( x % 3 ) * x ) % 5 )==0の順で評価される。
一つ一つ分解して説明

(x%3)

 3の倍数なら0、倍数でないなら0以外(1か2)

((x%3)*x)

 3の倍数なら0、倍数でないならxまたは、2x

　※xが5の倍数なら2xも5の倍数、逆も成立

(((x%3) * x ) % 5 

　3の倍数もしくは5の倍数なら0

(((x%3) * x ) % 5)==0 

　左辺の結果が0(=xが3の倍数もしくは5の倍数)なら真

 →リストに追加



我ながらもはや意味わからず。

■if文なし繰り返し
p=[0]*1000
p[3:1000:3] = range(3,1000,3)
p[5:1000:5] = range(5,1000,5)
print sum(p)

プロジェクトオイラー 35問目

35問目
前回の別解をつくってみた。

get_listは、mapに不特定の数のリストを指定する方法がわからなかったので、作ってみた。

import math

class Prime:

  def __init__(self):

    pass

  def is_prime(self,number):

    if number % 2 :

      if number > 2:

        return self.prime[(number-3)//2]

      else:

        return False

    else:

      if number == 2:

        return True

      else:

        return False

  def get_prime(self,max):

    return p.prime

  def ert(self,max):

    self.prime = range(3,max+1,2)

    sqrt_max = int(math.sqrt(max))+1

    l = (max+1)/2 - 1

    for  i in xrange(sqrt_max/2):

      m = 2*(i)+3

      if self.prime[i]:

        r = (m*m-3)/2

        self.prime[r:l:m] = [0] * ((l-r-1)//m+1)



def get_str(x,y):

  return  "%s%s" % (x,y)

def mapping(f,l1,l2):

  return [ f(x,y) for x in l1 for y in l2 ]

def get_list(f,list,number=2):

  if number == 2:

    return mapping(f,list,list)

  if number > 2:

    return mapping(f,list,get_list(f,list,number - 1))

  else :

    return list



def is_junkai_prime(p,number1,number2=None):

  if number1 == number2:

    return True

  elif not number2:

    number2 = number1

  if p.is_prime(int(number2)):

    if int(number2) > 9:

      return is_junkai_prime(p,number1,get_str(number2[-1],number2[:-1]))

    else:

      return True

  else:

    return False



if __name__ == '__main__':

  number_list=[1,3,7,9]

  max_len = 6

  max = 10**max_len

  p = Prime()

  p.ert(max)

  #2,5 is not count prime

  count = 2

  for i in range(1,max_len+1):

    for number in get_list(get_str,number_list,i):

      if is_junkai_prime(p,number):

        print number

        count += 1

  print count

実行結果：

real    0m0.161s

user    0m0.128s

sys     0m0.019s



意外に速い

エラトステネスの篩はなぜ速いのか

素数判定の有力なアルゴリズムといわれるエラトステネスの篩であるが、なぜ速いのか。

Web上には、単純な誤解からエラトステネスは遅いとの文献も散見される。
そういった誤解を解くためにもわかりやすい図を書いて説明する。

なお、エラトステネスの篩が遅いといわれる理由は、下記2点による誤解のようである。（後述の剰余計算のアルゴリズムの改悪版をエラトステネスの篩と誤解しているようである）
1.素数のリストを作るアルゴリズムだと思われていること
2.素数を判定するための計算では、剰余計算が一般的であること

例えば30が素数かどうかを判定するとする。
エラトステネスの篩では、平方根以下の2,3,5の倍数を計算するため、30に関する計算だけでも、3回の計算が必要とされる。
一方、剰余で計算する場合には、2で割り算する（＝割り切れる）だけなので、1回の計算で済む。
一つの数字に絞った場合は、剰余で計算したほうが効率よさそうに見える。

しかし、集合を相手にした場合は、エラトステネスの篩のほうがものすごく効率よくなる。
エラトステネスの篩のアルゴリズムは、ある集合から倍数を除いていき、最後に残った数が素数であると判定するものである。
200万以下でいうと、2の倍数の計算は100万回である。1000付近の素数の計算でも、2000回ちょっとの倍数計算で済む。（しかも足し算で倍数を計算していけるので、計算コストが低くなる）

剰余で200万以下の素数を洗い出すためには、2で割り算する回数は200万回である。1000付近ですら、素数の数(14万数千個）の回数以上の剰余計算をする必要がある。

どう見てもダンチでエラトステネスの篩のほうが速い。

答え：
計算回数が少ないから

プロジェクトオイラー 32問目

32問目

a * b = c
としたとき、cを素因数分解し、その後素因数たちを掛け算してaを求める。cをaで割ってbを得る。
そうして得られたa,b,cの組が、条件[1から9が一回ずつ出現する]を満たすかどうかを検討すればよい。

a,b < c からcは4桁未満にならないことが容易に分かる。
また、99*99<10000(5桁)より、cは5桁以上にもならない。
よって、探索する範囲は1000-10000でよい。

なお、本問を解くに当たって、素因数分解するメソッド、分解された素因数から因数を求めるメソッド、数値を一文字ごとに分離するメソッドを作成した。これがかなりめんどくさかった。

処理が遅いのは、このメソッドたちがあまり効率的でないからだと思う。後ほど最適化を行いたい。とくにgetFactorsとかソースがカオス過ぎて自分でも読めない。

def splitInteger(number):

  ret=[]

  strings = str(number)

  for i in range(len(strings)):

    ret.append(int(strings[i]))

  return ret



def getFactors(factors,min=1,p=None):

  i=0

  while i < len(factors):

    fact = copy.deepcopy(factors)

    f = fact.pop(i)

    f_tmp = f

    if min < f:

      yield f

      min = f

      for pd in getFactors(fact,min):

        yield pd*f_tmp

      while len(fact) > i and f == fact[i]:

        fact.pop(i)

        f_tmp *= f

        yield f_tmp

        for pd in getFactors(fact,min):

          yield pd*f_tmp

    i+=1



def primeFactorization(number,p=None):

  ret = []

  if number == 1:

    return ret

  if not p:

    p = Primes(number)

  p.eratosthen(number)

  for i in p.get_primes():

    if number < i ** 2:

      if number != 1:

        ret.append(number)

      break

    while not number % i:

      ret.append(i)

      number //= i

  return ret

MIN=1000

MAX=10000

TARGET=[1,2,3,4,5,6,7,8,9]



def condition(n1,n2,n3):

  if len(str(n1)) + len(str(n2)) + len(str(n3)) == len(TARGET):

    arr = splitInteger(n1)

    arr += splitInteger(n2)

    arr += splitInteger(n3)

    sorted_arr = sorted(set(arr),key=arr.index)

    sorted_arr.sort()

    if TARGET == sorted_arr:

      return True

  return False



def get_tri(n,prime):

  for i in getFactors(primeFactorization(n,prime)):

    d = n // i

    if d < i:

      continue

    if condition(i,d,n):

      print i,d,n

      return True

  return False



sum=0

prime=Primes(MAX)

for i in range(MIN,MAX):

  if get_tri(i,prime): sum += i

  #print i,sum

#print get_tri(7254,prime)

print sum

プロジェクトオイラー31問目

31問目
一度カウントに入れた通貨以下もののしか使用しないように注意して再帰呼び出しで終了。
一度1ペンスが使われたら1通りしかないため、1を返すようにした。

結果
73682

real 0m0.232s
user 0m0.221s
sys 0m0.010s

若干かかり過ぎなきはす。

pen = [200,100,50,20,10,5,2,1]



def get_pens(max, use=200):

  if use == 1:

    return 1

  count = 0

  for i in pen:

    if i > max or i > use:

      continue

    elif i == max:

      count +=1

    else:

      if i == use:

        count += get_pens(max-i,use)

      else:

        count += get_pens(max-i,i)



  #print "get_pens max:%d use:%d ,count:%d" % (max,use,count)

  return count



if __name__ == '__main__':

  print get_pens(200)

プロジェクトオイラー 45問目

45問目

一番大きそうな6角数を元に探索。そのほかも条件を満たせば終了。

class TriangularNumber:

  def get(self,n):

    return n*(n+1)//2

  def is_number(self,n):

    return not (-1+math.sqrt(1+8*n))%2



class PentagonalNumber:

  def get(self,n):

    return n*(3*n-1)//2

  def is_number(self,n):

    return not (1+math.sqrt(1+24*n)) % 6



class HexagonalNumber:

  def get(self,n):

    return n*(2*n-1)

  def is_number(self,n):

    return not (1+math.sqrt(1+8*n))%4

tn = TriangularNumber()

pn = PentagonalNumber()

hn = HexagonalNumber()



def is_tri_pen_hex(n):

  if tn.is_number(n) and pn.is_number(n) and hn.is_number(n):

    return True

  else:

    return False



i=0

while 1:

  n = hn.get(i)

  if is_tri_pen_hex(n):

    print n

    if n > 40755:

      break

  i+=1

プロジェクトオイラー 10問目 素数判定を速くする その2

その1

以前のコードの致命的結果に気がついた。
エラトステネスの篩とか言いながら、そもそも剰余で素数判定をしているところがおかしい。

1.剰余で素数判定する->素因数分解の困難性により、計算量が膨大になる
2.掛け算で合成数を除く->重複が（剰余で素数判定に比べ）少なくなり、計算量が削減できる

1のほうで素数の集合を作り、その中に素数が含まれているかを計算していたのだから計算量が膨大になって、通常の剰余判定と変わらないか、余計な偶数まで判定するため遅くなるわけである。

原因が分かったので、修正した。
また、すべての数を持つならハッシュテーブル（辞書）より配列のほうがいいということに気がつき、配列にしてみた。

結果：2秒

今まで何やっていたのか。

class Primes():

  def __init__(self,interval=10**4):

    self.prime = [False,True]

    self.max = 2

    self.interval = interval



  def eratosthen(self,max):

    if max < self.max:

      return

    self.prime += (max - self.max) * [True]

    self.max = max

    for p in range(2,int(math.sqrt(max))+1):

      if self.prime[p]:

        for q in range(p*2,max,p):

          self.prime[q] = False



  def is_prime(self,q):

    return self.prime[q]



  def is_prime_with_eratosthen(self,q):

    while q > self.max:

      self.eratosthen(self.max+self.interval)

    return self.is_prime(q)

プロジェクトオイラー 48問目

48問目

単純に足して終了

def sum(fanction,start,end):

  sum = 0

  for i in range(start,end+1):

    sum += fanction(i)

  return sum



def f(n,digit=10):

  return n ** n



print str(sum(f,1,1000))[-10:]

プロジェクトオイラー 44問目

プロジェクトオイラー44問目
---
五角数は Pn = n(3n-1)/2で生成される. 最初の10項は

1, 5, 12, 22, 35, 51, 70, 92, 117, 145, ...
である.

P4 + P7 = 22 + 70 = 92 = P8である. しかし差 70 - 22 = 48は五角数ではない.

五角数のペア PjとPkについて, 差と和が五角数になるものを考える. 差を D = |Pk - Pj| と書く. 差 D の最小値を求めよ.
---

解答方針
5角数を配列に記憶し、その配列に存在するかどうかをチェックすることで、容易に五角数かどうかチェックできるようにする。
また、Dの最小値が最小値と確定するのは、あるnでPn-Pn-1>Dとなったときである。ちなみに、Pn-Pn-1=3n-2なので、その数と得られた差分を比較し、差分を超えたときに終了とする。

・・・つもりだったが、そのままだといつまでたっても終わらんかった。
とりあえず最初に見つかった値を答えと入力してみてみたところ正解だった。

class PentagonalNumber:

  def __init__(self):

    self.pn = [1]



  def get(self,n):

    if len(self.pn) < n:

      self.set(self._get_pentagonal(n))

    return self.pn[n-1]



  def set(self,n):

    while self.pn[-1] < n:

      self.pn.append(self._get_pentagonal(len(self.pn)))



  def _get_pentagonal(self,n):

    return n*(3*n-1)//2



  def is_pentagonal(self,n):

    self.set(n)

    return n in self.pn







if __name__ == '__main__' :

  d=0

  i=2

  pn = PentagonalNumber()



  def get_threthold(n):

    return 3*n-2



  while 1:

  #print pn

    if d and get_threthold(i) >= d:

      break

    pn1 = pn.get(i)

    for j in range(1,i):

      pn2 = pn.get(i-j)

      if d and pn1 - pn2 >= d:

        break

      if pn.is_pentagonal(pn1 - pn2) and pn.is_pentagonal(pn1 + pn2):

        d = pn1 - pn2

        print d

    #print pn1, pn2 ,d

    i+=1



  print d