[Euler Project 203] 소수제곱없는 이항계수의 합 구하기

0C0 부터 50C50 까지의 이항계수들 가운데에 소인수분해하여, 제곱이상의 인자가 없는 것들의 합을 구하는 것. 트릭은 nCr = n x ((n-1)C(r-1) / r 이라는 관계식을 이용하는 것, 수를 50이하의 소수 (15개)의 거듭제곱수의 리스트로 표현하여 곱과 나누기의 오버헤드를 줄이고, 제곱 이상의 인자가 없는지 확인하기 편하게 하는 것.

아래 코드는 좀 지져분하다.

---

#!/usr/bin/env python

def MUL(l1, l2, l3):
       n = len(primes)
       #l3 = [ 0 ] * n
       for i in range(n):
               l3[i] = l1[i]+l2[i]
       return

def DIV(l1, l2, l3):
       n = len(primes)
       #l3 = [ 0 ] * n
       for i in range(n):
               l3[i] = l1[i]-l2[i]
       return

def is_sqare_free(l):
       n = len(primes)
       for i in range(n):
               if l[i] >= 2:
                       return False
       return True

#------------ get primes under 50
primes = [ 2 ]

n = 3
while n <= 50:
       for p in primes:
               if p*p > n:
                       primes.append(n)
                       break
               if n%p == 0:
                       break
       n+=1

print primes

#------------ get primes compositions of n under 50

N = [ [ 0 for i in range(len(primes)) ] for j in range(52) ]

for n in range(1, 51):
       nn = n
       print n
       for i in range(len(primes)):
               k = 0
               while nn%primes[i] == 0:
                       nn /= primes[i]
                       k += 1
               N[n][i] = k
       print n, N[n]

for n in range(1, 51):
       print n, N[n]
#------------- get primes composition of nCr from 0C0 ... 50C50
NCR = [ [ [ 0 for i in range(len(primes)) ] for j in range(52) ] for k in range(52) ]
ncr = [ [ 1 for i in range(51) ] for j in range(51) ]
#T = [ 0 ] *len(primes)
s = [ ]
print N
for n in range(51):
       print "-------", n
       print N
       for r in range(n+1):
               if r == 0:
                       #NCR[n][r] = [0]*len(primes)
                       ncr[n][r] = 1
               else:
                       #MUL(N[n], NCR[n-1][r-1], T)
                       for i in range(len(primes)):
                               NCR[n][r][i] = N[n][i] + NCR[n-1][r-1][i] - N[r][i]
                       ncr[n][r] = n*ncr[n-1][r-1]/r
               print n, r, ncr[n][r], NCR[n][r]
               if is_sqare_free(NCR[n][r]) and not ncr[n][r] in s:
                       s.append(ncr[n][r])
print sum(s)

---