다른 분들이 올려주신 문제만 풀다가 저도 한 번 올려봅니다.
길가의 풀님이 소개해주신 프로젝트 오일러의 101번 문제입니다. 문제 원문
영어인데 100번까지 한국어로 번역된 사이트가 있습니다. 다시 보니 113번까지 갱신된 듯 하네요. 번역된 문제
사이냅소프트에서 제공하고 있습니다.
어떤 수열의 첫 k개의 항이 주어졌다고 해도, 그 다음 항을 확신할 수는 없습니다. k개의 항을 만족하는 다항식은 그야말로 무한정 만들어낼 수 있기 때문입니다.
예를 들어, 세제곱 수의 수열을 생각해 봅시다.
un = n^3 : 1, 8, 27, 64, 125, 216...
여기서 앞의 숫자 두 항만을 알고 있다고 가정합시다. "단순한 게 장땡"이라는 원칙에 따라, 우리는 이 수열이 공차가 7인 등차수열이라고 가정하고 다음 항이 15일 것이라고 추측할 수 있습니다. 첫 3개항이 주어지더라도, 2차 다항식으로 오해될 여지가 있습니다.
처음 k개의 항이 주어지고, 이 수열을 바탕으로 만들 수 있는 최소차 다항식으로 구할 수 있는 n번째 항을 OP(k, n)이라고 합시다. n≤k일 때는 무조건 정답이 나올 수밖에 없습니다. 하지만 OP(k, k+1)은 첫 오답(first incorrect term : 이하 FIT이라고 통칭합니다)이 될 가능성이 있습니다. 이런 경우 OP(k, k+1)은 bad OP(BOP)라고 부르기로 합니다.
위 정의대로 한다면, 수열의 첫 항만이 주어진 경우에는 OP(1, n) =u1이라고 가정하게 될 것입니다.
그러므로 우리는 다음과 같은 순서로 OP를 구할 수 있습니다.
OP(1, n) = 1 1, 1, 1, 1, ...
OP(2, n) = 7n-6 1, 8, 15, ...
OP(3, n) = 6n^2-11n+6 1, 8, 27, 58, ...
OP(4, n) = n^3 1, 8, 27, 64, 125, ...
k≥4일 때 BOP가 존재하지 않으리라는 사실은 명백합니다.
BOP에서 나온 FIT을 모두 더하면(예제에서 빨간색으로 표시해 두었습니다), 그 값은 1 + 15 + 58 = 74가 됩니다.
자 그럼 다음 10차 다항함수를 잘 보세요.
un = 1 - n + n^2 - n^3 + n^4 - n^5 + n^6 - n^7 + n^8 - n^9 + n^10
이 함수의 BOP에서 나온 FIT의 총합은 얼마입니까?
11개의 풀이가 있습니다.
def get_diff(arr)
# [1, 2, 5] -> [1, 3]
arr2 = arr.each_cons(2).map{|a,b| b - a}
end
def get_fit(terms)
if terms.length == 1
return terms[0]
end
# [1, 8, 27]
# [7, 19]
# [12]
diff = []
diff[1] = get_diff(terms)
(2...terms.length).each{|i| diff[i] = get_diff(diff[i-1]) }
# [1, 8, 27, 58] 58 = 27 + 31
# [7, 19, 31] 31 = 19 + 12
# [12, 12] <- same difference at lowest degree
(terms.length-2).downto(1){|i|
t = diff[i+1].last
diff[i] << diff[i].last + t
}
terms.last + diff[1].last
end
def u(n)
[1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1].each_with_index.map{ |e, i| e * (n ** i) }.inject(:+)
end
result = 0
arr = []
(1..10).each do |i|
arr << u(i)
result += get_fit(arr)
end
puts result
f=function(x) 1-x+x^2-x^3+x^4-x^5+x^6-x^7+x^8-x^9+x^10
u.n=f(1:10)
op=function(k){
if(k==1) return(rep(u.n[1], 2))
a=solve(outer(1:k, (k-1):0, "^"), u.n[1:k])
c(u.n[1:k], (k+1)^((k-1):0)%*%a)
}
list.answers=sapply(1:10, op)
sum(sapply(list.answers, function(x) tail(x, n=1)))
R입니다
op(k, n) = un은 미지수가 k개고 식도 k개이며 n의 최대 차수가 k-1인 연립방정식입니다. 역행렬을 통해 연립방정식의 해를 구하고, fit을 구하기 우해 n=k+1일 때의 값을 구합니다. k=11일 때 op(k, n) = un이 되므로 k가 1부터 10일 때 op(k, n)은 bop입니다. 따라서 모든 bop의 fit 합을 구하면
> sum(sapply(list.answers, function(x) tail(x, n=1)))
[1] 37076114526
파이썬입니다. 풀이방식은 위의 분과 거의 같아요
import numpy as np
def u(n):
return 1 - n + n**2 - n**3 + n**4 - n**5 + n**6 - n**7 + n**8 - n**9 + n**10
l = [u(x) for x in range(11)[1:]]
l1 = []
for i in range(10):
a = np.array([[x**y for y in range(i+1)[::-1]] for x in range(i+2)[:0:-1]][::-1])
b = np.array([[x] for x in l][:i+1])
x = np.linalg.solve(a,b)
c = np.array([(i+2)**y for y in range(i+1)[::-1]])
l1.append(np.dot(c,x))
print "%10.f" % sum(l1)
import numpy as np
def f(n):
return sum((-n)**x for x in range(11))
def arr(i):
li = []
for n in range(i):
li.append([])
for p in range(i):
li[n].append((-(n+1))**p)
return np.array(li)
def gen(func):
x = 1
while 1:
yield func(x)
x += 1
if __name__ == '__main__':
g = gen(f); result = []
full = np.array([list(next(g) for x in range(11))])
for i in range(1, 11):
X = arr(i).T
B = full[0][:i]
A = [email protected](X)
X = arr(i+1).T[:i][:]
k = (A@X)[i]
result.append(k)
print(sum(result))
@ 연산자는 파이썬 3.5에서 추가되었습니다. 근데 왜 파이썬은 역행렬이 기본 라이브러리에 없을까 ㅠ
파이썬 3.5.1
Ruby
Kim Jaeju님 답을 보고 문제가 이해되었네요. 알고리즘 그대로 차용했습니다. 함수 fn_u가 주어졌을 때 1~10을 넣서 만들어진 수열 op[n] 각각에 대해서 fit를 구합니다. op[n]의 fit는 차이수열을 만들어낸 뒤 이를 이용해 op[n]의 증분을 알아냅니다. 주어진 op[n]의 마지막 값과 증분의 마지막 값을 더해서 fit를 구합니다.
df = ->arr { arr.each_cons(2).map {|a,b| b-a } }
get_fit = ->op do
return op[0] if op.size == 1
dfs = (2...op.size).reduce([df[op]]) {|a,_| a.unshift df[a[0]] }
inc = (1..dfs.size).reduce([0,*dfs]) {|a,k| a[k] = a[k][-1] + a[k-1]; a }
op.last + inc.last
end
fit_sum = ->fn { op=[]; (1..10).reduce(0) {|sum,n| sum + get_fit[op << fn[n]] } }
fn_u = ->n { (0..10).reduce(0) {|value,i| value += (-1)**i * (n**i) } }
Test
expect(fit_sum[fn_u]).to eq 37076114526
한번도 답을 보고 푼적이 없는데... 도저히 모르겟드라구여 ㅠㅠ 이번풀이 Kim Jaeju 님 것은 배낀것과 같은데 정말 좋은 알고리즘이고 많이 배워서 답을 보고 푼것을 잘했다고 생각합니다
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double* getDifArr(double* arr, int size);
void main() {
double un[12]; // 5
int size = 11; // 4
for(int i=1;i<=size;i++) {
//un[i] = pow((double)i, 3);
un[i] = 1 -i + pow((double)i, 2) - pow((double)i, 3) + pow((double)i, 4) - pow((double)i, 5)
+ pow((double)i, 6) - pow((double)i, 7) + pow((double)i, 8) - pow((double)i, 9) + pow((double)i, 10);
}
for(int i=1;i<=size;i++) {
printf("%f\n", un[i]);
}
double result = 0;
for(int i=1;i<size;i++) {
if(i==1) result = result + un[1];
else {
int size = i;
double* arr = un;
double temp = 0;
while(size > 1) {
temp = temp + arr[size];
size--;
arr = getDifArr(arr, size);
}
result = result + temp + arr[1];
}
}
printf("%f", result);
}
double* getDifArr(double arr[], int size) {
double* rarr = (double*) malloc (sizeof(double) * size+1);
for(int i=1;i<=size;i++) {
rarr[i] = arr[i+1] - arr[i];
}
return rarr;
}
Python 수열의 차이를 반복해서 구해낸뒤 마지막항을 다 더하면 FIT이 됩니다.
ans = []
for i in range(1, 12):
tmp_sum = 0
for j in range(11):
k = 1 if j%2 == 0 else -1
tmp_sum += k * i**j
ans.append(tmp_sum)
print(ans)
fit = []
#ans = [1, 8, 27, 64]
for i in range(len(ans)-1):
if i == 0:
fit.append(ans[i])
else:
cha = ans[:i+1]
fit_var = ans[i]
while len(cha) != 1:
tmp = []
for j in range(len(cha)-1):
tmp.append(cha[j+1]-cha[j])
cha = tmp
fit_var += cha[-1]
fit.append(fit_var)
print(sum(fit))
import numpy as np
un = lambda n: 1 - n + n**2 - n**3 + n**4 - n**5 + n**6 - n**7 + n**8 - n**9 + n**10
li = [un(i) for i in range(1,12)]
r=0
for m in range(1,11):
mat = [[(i+1)**(j) for j in range(len(li[:m]))] for i in range(len(li[:m]))]
a = np.array(mat)
y = np.array(li[:m])
r += np.linalg.inv(a)@y.T @ np.array([(m+1)**(i) for i in range(m)])
mat.clear()
print(r)
37076114526.04962
import sympy
import re
def lagrange(inp_LIST) : # 라그랑주 보간
if len(inp_LIST) == 1 :
return inp_LIST[0][1]
else :
x, inter, kou = sympy.Symbol('x'), [], []
for i in range(0, len(inp_LIST)) :
for m in range(1, len(inp_LIST)) :
inter.append('((x-%s)/(%s-%s))'% (inp_LIST[(i+m)%len(inp_LIST)][0], inp_LIST[i][0], inp_LIST[(i+m)%len(inp_LIST)][0] ))
kou.append("*".join(inter)+'*'+str(inp_LIST[i][1]))
inter = []
return sympy.expand(('+'.join(kou)))
def opt_poly(count, un) :
y_list, result, x = [], [], sympy.Symbol('x')
while True :
y_list.append([str(count), str(input("%s번째 항 : " %count))])
count += 1
lag_func = sympy.simplify(lagrange(y_list))
print(lag_func)
result.append(lag_func.subs(x, count)) # FIT 을 result에 추가
if sympy.simplify(lag_func) == sympy.simplify(un) :
print(sum(result[:-1]))
#최종 합을 출력
break
opt_poly(1, '1-x+x**2-x**3+x**4-x**5+x**6-x**7+x**8-x**9+x**10')
# opt_poly(1, 'x**3')
158번 문제에서 힌트를 얻어, 라그랑주 보간법을 활용해봤습니다. 사용자 입력으로 각 수열의 값을 입력받아, 누적된 수열에 기초해 일반항을 계산해 내고, FIT을 계산해서 result에 추가합니다. 애초에 의도되었던 일반항과 라그랑주 보간법으로 계산된 일반항이 일치할 경우, result의 요소들을 합해서 출력하고 루프에서 빠져나옵니다.
1번째 항 : 1
1
2번째 항 : 683
682*x - 681
3번째 항 : 44287
21461*x**2 - 63701*x + 42241
4번째 항 : 838861
118008*x**3 - 686587*x**2 + 1234387*x - 665807
5번째 항 : 8138021
210232*x**4 - 1984312*x**3 + 6671533*x**2 - 9277213*x + 4379761
6번째 항 : 51828151
159060*x**5 - 2175668*x**4 + 11535788*x**3 - 29116967*x**2 + 34305227*x - 14707439
7번째 항 : 247165843
58542*x**6 - 1070322*x**5 + 8069182*x**4 - 31492582*x**3 + 65955241*x**2 - 68962861*x + 27442801
8번째 항 : 954437177
11165*x**7 - 254078*x**6 + 2524808*x**5 - 13814218*x**4 + 44083303*x**3 - 80663539*x**2 + 76941359*x - 28828799
9번째 항 : 3138105961
1111*x**8 - 28831*x**7 + 352528*x**6 - 2514688*x**5 + 11126621*x**4 - 30669221*x**3 + 50572225*x**2 - 44806465*x + 15966721
10번째 항 : 9090909091
54*x**9 - 1319*x**8 + 18149*x**7 - 157772*x**6 + 902054*x**5 - 3416929*x**4 + 8409499*x**3 - 12753575*x**2 + 10628639*x - 3628799
11번째 항 : 23775972551
x**10 - x**9 + x**8 - x**7 + x**6 - x**5 + x**4 - x**3 + x**2 - x + 1
37076114526
"""
op(4, n) = n^3 :
1. 8, 27, 64, 125
7 19 37 61
12 18 24
15 6
BOP = 1 + 15 + 58
= 1 + (7+8) + (12+19+27) = (1+8+27) + (7+19) + 12
""
def un():
lst = []
for n in range(1, 11):
v = 1
tmp=1
for _ in range(1, 11):
tmp *= (-n)
v += tmp
lst.append(v)
return lst
steplst = []
steplst.append(un())
for n in range(9):
templst = []
for i in range(1, len(steplst[n])):
templst.append(steplst[n][i]- steplst[n][i-1])
steplst.append(templst)
sumlst = 0
for n in range(10):
sumlst += sum(steplst[n][:10-n])
print(sumlst)
JAVA에서는 다루기 쉬운 행렬 라이브러리가 없어 Python으로 작성했습니다.
def getPolynomialValue(n, Un):
sum = 0
for i in range(0, len(Un)):
sum += (Un[i])*(n**i)
return sum
import numpy as np
Un = np.array([1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1])
bopSum = 0
for i in range(1, len(Un)):
OP = []
polyValue = []
for j in range(0, i):
OP.append([])
polyValue.append(getPolynomialValue(j+1, Un))
for k in range(0, i):
OP[j].append((j+1)**k)
OP = np.array(OP)
polyValue = np.array(polyValue)
ans = np.linalg.solve(OP, polyValue)
if(getPolynomialValue(i+1, ans) != getPolynomialValue(i+1, Un)):
bopSum += getPolynomialValue(i+1, ans)
print(bopSum)