From 0df84c808dfdfccb0629ede6a46d22ee4a104ab4 Mon Sep 17 00:00:00 2001 From: Freddy Pringle Date: Mon, 7 Dec 2020 17:31:18 +0100 Subject: [PATCH 1/7] Added solution for Project Euler problem 101 --- project_euler/problem_101/__init__.py | 0 project_euler/problem_101/sol1.py | 198 ++++++++++++++++++++++++++ 2 files changed, 198 insertions(+) create mode 100644 project_euler/problem_101/__init__.py create mode 100644 project_euler/problem_101/sol1.py diff --git a/project_euler/problem_101/__init__.py b/project_euler/problem_101/__init__.py new file mode 100644 index 000000000000..e69de29bb2d1 diff --git a/project_euler/problem_101/sol1.py b/project_euler/problem_101/sol1.py new file mode 100644 index 000000000000..2c959fdeeb38 --- /dev/null +++ b/project_euler/problem_101/sol1.py @@ -0,0 +1,198 @@ +""" +If we are presented with the first k terms of a sequence it is impossible to say with +certainty the value of the next term, as there are infinitely many polynomial functions +that can model the sequence. + +As an example, let us consider the sequence of cube +numbers. This is defined by the generating function, +u(n) = n3: 1, 8, 27, 64, 125, 216, ... + +Suppose we were only given the first two terms of this sequence. Working on the +principle that "simple is best" we should assume a linear relationship and predict the +next term to be 15 (common difference 7). Even if we were presented with the first three +terms, by the same principle of simplicity, a quadratic relationship should be +assumed. + +We shall define OP(k, n) to be the nth term of the optimum polynomial +generating function for the first k terms of a sequence. It should be clear that +OP(k, n) will accurately generate the terms of the sequence for n ≤ k, and potentially +the first incorrect term (FIT) will be OP(k, k+1); in which case we shall call it a +bad OP (BOP). + +As a basis, if we were only given the first term of sequence, it would be most +sensible to assume constancy; that is, for n ≥ 2, OP(1, n) = u(1). + +Hence we obtain the +following OPs for the cubic sequence: + +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, ... + +Clearly no BOPs exist for k ≥ 4. + +By considering the sum of FITs generated by the BOPs (indicated in red above), we +obtain 1 + 15 + 58 = 74. + +Consider the following tenth degree polynomial generating function: + +1 - n + n^2 - n^3 + n^4 - n^5 + n^6 - n^7 + n^8 - n^9 + n^10 + +Find the sum of FITs for the BOPs. +""" + + +from typing import Callable, List, Union + +Matrix = List[List[Union[float, int]]] + + +def solve(A: Matrix, b: Matrix) -> Matrix: + """ + Solve the linear system of equations Ax = b for x using Gaussian elimination + and back substitution. We assume that A is an invertible square matrix and + that b is a column vector of the same height. + >>> solve([[1, 0], [0, 1]], [[1],[2]]) + [[1.0], [2.0]] + >>> solve([[2, 1, -1],[-3, -1, 2],[-2, 1, 2]],[[8], [-11],[-3]]) + [[2.0], [3.0], [-1.0]] + """ + size: int = len(A) + augmented: Matrix = [[0 for _ in range(size + 1)] for _ in range(size)] + row: int + col: int + i_max: int + f: float + j: int + + for row in range(size): + for col in range(size): + augmented[row][col] = A[row][col] + + augmented[row][size] = b[row][0] + + row = 0 + col = 0 + while row < size and col < size: + # pivoting + i_max = max([(abs(augmented[i][col]), i) for i in range(col, size)])[1] + if augmented[i_max][col] == 0: + col += 1 + continue + else: + augmented[row], augmented[i_max] = augmented[i_max], augmented[row] + + for i in range(row + 1, size): + f = augmented[i][col] / augmented[row][col] + augmented[i][col] = 0 + for j in range(col + 1, size + 1): + augmented[i][j] -= augmented[row][j] * f + + row += 1 + col += 1 + + # back substitution + for col in range(1, size): + for row in range(col): + f = augmented[row][col] / augmented[col][col] + for j in range(col, size + 1): + augmented[row][j] -= augmented[col][j] * f + + # round to get rid of numbers like 2.000000000000004 + return [ + [round(augmented[row][size] / augmented[row][row], 10)] for row in range(size) + ] + + +def interpolate(y_list: List[int]) -> Callable[[int], int]: + """ + Given a list of data points (1,y0),(2,y1), ..., return a function that + interpolates the data points. We find the coefficients of the interpolating + polynomial by solving a system of linear equations corresponding to + x = 1, 2, 3... + + >>> interpolate([1])(3) + 1 + >>> interpolate([1, 8])(3) + 15 + >>> interpolate([1, 8, 27])(4) + 58 + >>> interpolate([1, 8, 27, 64])(6) + 216 + """ + + size: int = len(y_list) + A: Matrix = [[0 for _ in range(size)] for _ in range(size)] + b: Matrix = [[0] for _ in range(size)] + a: Matrix + i: int + y: int + + for i, y in enumerate(y_list): + for j in range(size): + A[i][j] = (i + 1) ** (size - j - 1) + b[i][0] = y + + a = solve(A, b) + + return lambda x: sum(round(a[i][0]) * (x ** (size - i - 1)) for i in range(size)) + + +def u(n: int) -> int: + """ + The generating function u as specified in the question. + >>> u(0) + 1 + >>> u(1) + 1 + >>> u(5) + 8138021 + >>> u(10) + 9090909091 + """ + return ( + 1 + - n + + n ** 2 + - n ** 3 + + n ** 4 + - n ** 5 + + n ** 6 + - n ** 7 + + n ** 8 + - n ** 9 + + n ** 10 + ) + + +def solution(func: Callable[[int], int] = u, order: int = 10) -> int: + """ + Find the sum of the FITs of the BOPS. For each interpolating polynomial of order + 1, 2, ... , 10, find the first x such that the value of the polynomial at x does + not equal u(x). + >>> solution(lambda n: n ** 3, 3) + 74 + """ + data_points: List[int] = list(map(func, range(1, order + 1))) + + polynomials: List[Callable[[int], int]] = [ + interpolate(data_points[:i]) for i in range(1, order + 1) + ] + + ret: int = 0 + i: int + x: int + + for i in range(order): + x = 1 + while func(x) == polynomials[i](x): + x += 1 + + ret += polynomials[i](x) + + return ret + + +if __name__ == "__main__": + print(f"{solution() = }") From 89782777f8924a2fc6ef8b6df64db9925ba3f89e Mon Sep 17 00:00:00 2001 From: Freddy Pringle Date: Sat, 12 Dec 2020 13:17:49 +0100 Subject: [PATCH 2/7] Got rid of map functions --- project_euler/problem_101/sol1.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/project_euler/problem_101/sol1.py b/project_euler/problem_101/sol1.py index 2c959fdeeb38..4aabd0096cae 100644 --- a/project_euler/problem_101/sol1.py +++ b/project_euler/problem_101/sol1.py @@ -174,7 +174,7 @@ def solution(func: Callable[[int], int] = u, order: int = 10) -> int: >>> solution(lambda n: n ** 3, 3) 74 """ - data_points: List[int] = list(map(func, range(1, order + 1))) + data_points: List[int] = [func(x) for x in range(1, order + 1)] polynomials: List[Callable[[int], int]] = [ interpolate(data_points[:i]) for i in range(1, order + 1) From 14629250e5e14b98ff7d6e01e7407cda04c89cba Mon Sep 17 00:00:00 2001 From: github-actions <${GITHUB_ACTOR}@users.noreply.github.com> Date: Sun, 13 Dec 2020 11:18:04 +0000 Subject: [PATCH 3/7] updating DIRECTORY.md --- DIRECTORY.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/DIRECTORY.md b/DIRECTORY.md index 929a986b0f3b..1f1bb9907e52 100644 --- a/DIRECTORY.md +++ b/DIRECTORY.md @@ -744,6 +744,8 @@ * [Sol1](https://github.com/TheAlgorithms/Python/blob/master/project_euler/problem_097/sol1.py) * Problem 099 * [Sol1](https://github.com/TheAlgorithms/Python/blob/master/project_euler/problem_099/sol1.py) + * Problem 101 + * [Sol1](https://github.com/TheAlgorithms/Python/blob/master/project_euler/problem_101/sol1.py) * Problem 112 * [Sol1](https://github.com/TheAlgorithms/Python/blob/master/project_euler/problem_112/sol1.py) * Problem 113 From dd70ff408d141ea6b4ad147da037217667f79f3f Mon Sep 17 00:00:00 2001 From: Freddy Pringle Date: Sun, 13 Dec 2020 12:46:11 +0100 Subject: [PATCH 4/7] Better function/variable names --- project_euler/problem_101/sol1.py | 57 ++++++++++++++++--------------- 1 file changed, 29 insertions(+), 28 deletions(-) diff --git a/project_euler/problem_101/sol1.py b/project_euler/problem_101/sol1.py index 4aabd0096cae..d0c73d217ea2 100644 --- a/project_euler/problem_101/sol1.py +++ b/project_euler/problem_101/sol1.py @@ -48,17 +48,18 @@ Matrix = List[List[Union[float, int]]] -def solve(A: Matrix, b: Matrix) -> Matrix: +def solve(matrix: Matrix, vector: Matrix) -> Matrix: """ - Solve the linear system of equations Ax = b for x using Gaussian elimination - and back substitution. We assume that A is an invertible square matrix and - that b is a column vector of the same height. + Solve the linear system of equations Ax = b (A = "matrix", b = "vector") + for x using Gaussian elimination and back substitution. We assume that A + is an invertible square matrix and that b is a column vector of the + same height. >>> solve([[1, 0], [0, 1]], [[1],[2]]) [[1.0], [2.0]] >>> solve([[2, 1, -1],[-3, -1, 2],[-2, 1, 2]],[[8], [-11],[-3]]) [[2.0], [3.0], [-1.0]] """ - size: int = len(A) + size: int = len(matrix) augmented: Matrix = [[0 for _ in range(size + 1)] for _ in range(size)] row: int col: int @@ -68,9 +69,9 @@ def solve(A: Matrix, b: Matrix) -> Matrix: for row in range(size): for col in range(size): - augmented[row][col] = A[row][col] + augmented[row][col] = matrix[row][col] - augmented[row][size] = b[row][0] + augmented[row][size] = vector[row][0] row = 0 col = 0 @@ -123,50 +124,50 @@ def interpolate(y_list: List[int]) -> Callable[[int], int]: """ size: int = len(y_list) - A: Matrix = [[0 for _ in range(size)] for _ in range(size)] - b: Matrix = [[0] for _ in range(size)] + matrix: Matrix = [[0 for _ in range(size)] for _ in range(size)] + vector: Matrix = [[0] for _ in range(size)] a: Matrix i: int y: int for i, y in enumerate(y_list): for j in range(size): - A[i][j] = (i + 1) ** (size - j - 1) - b[i][0] = y + matrix[i][j] = (i + 1) ** (size - j - 1) + vector[i][0] = y - a = solve(A, b) + a = solve(matrix, vector) return lambda x: sum(round(a[i][0]) * (x ** (size - i - 1)) for i in range(size)) -def u(n: int) -> int: +def question_function(variable: int) -> int: """ The generating function u as specified in the question. - >>> u(0) + >>> question_function(0) 1 - >>> u(1) + >>> question_function(1) 1 - >>> u(5) + >>> question_function(5) 8138021 - >>> u(10) + >>> question_function(10) 9090909091 """ return ( 1 - - n - + n ** 2 - - n ** 3 - + n ** 4 - - n ** 5 - + n ** 6 - - n ** 7 - + n ** 8 - - n ** 9 - + n ** 10 + - variable + + variable ** 2 + - variable ** 3 + + variable ** 4 + - variable ** 5 + + variable ** 6 + - variable ** 7 + + variable ** 8 + - variable ** 9 + + variable ** 10 ) -def solution(func: Callable[[int], int] = u, order: int = 10) -> int: +def solution(func: Callable[[int], int] = question_function, order: int = 10) -> int: """ Find the sum of the FITs of the BOPS. For each interpolating polynomial of order 1, 2, ... , 10, find the first x such that the value of the polynomial at x does From 73c2aa128d034f3cc70b23c8a72c985761bb7ddb Mon Sep 17 00:00:00 2001 From: Freddy Pringle Date: Fri, 18 Dec 2020 11:32:28 +0100 Subject: [PATCH 5/7] Better variable names --- project_euler/problem_101/sol1.py | 70 +++++++++++++++++-------------- 1 file changed, 38 insertions(+), 32 deletions(-) diff --git a/project_euler/problem_101/sol1.py b/project_euler/problem_101/sol1.py index d0c73d217ea2..9f90ac29522a 100644 --- a/project_euler/problem_101/sol1.py +++ b/project_euler/problem_101/sol1.py @@ -62,10 +62,11 @@ def solve(matrix: Matrix, vector: Matrix) -> Matrix: size: int = len(matrix) augmented: Matrix = [[0 for _ in range(size + 1)] for _ in range(size)] row: int + row2: int col: int - i_max: int - f: float - j: int + col2: int + pivot_row: int + ratio: float for row in range(size): for col in range(size): @@ -77,18 +78,20 @@ def solve(matrix: Matrix, vector: Matrix) -> Matrix: col = 0 while row < size and col < size: # pivoting - i_max = max([(abs(augmented[i][col]), i) for i in range(col, size)])[1] - if augmented[i_max][col] == 0: + pivot_row = max( + [(abs(augmented[row2][col]), row2) for row2 in range(col, size)] + )[1] + if augmented[pivot_row][col] == 0: col += 1 continue else: - augmented[row], augmented[i_max] = augmented[i_max], augmented[row] + augmented[row], augmented[pivot_row] = augmented[pivot_row], augmented[row] - for i in range(row + 1, size): - f = augmented[i][col] / augmented[row][col] - augmented[i][col] = 0 - for j in range(col + 1, size + 1): - augmented[i][j] -= augmented[row][j] * f + for row2 in range(row + 1, size): + ratio = augmented[row2][col] / augmented[row][col] + augmented[row2][col] = 0 + for col2 in range(col + 1, size + 1): + augmented[row2][col2] -= augmented[row][col2] * ratio row += 1 col += 1 @@ -96,9 +99,9 @@ def solve(matrix: Matrix, vector: Matrix) -> Matrix: # back substitution for col in range(1, size): for row in range(col): - f = augmented[row][col] / augmented[col][col] - for j in range(col, size + 1): - augmented[row][j] -= augmented[col][j] * f + ratio = augmented[row][col] / augmented[col][col] + for col2 in range(col, size + 1): + augmented[row][col2] -= augmented[col][col2] * ratio # round to get rid of numbers like 2.000000000000004 return [ @@ -126,18 +129,21 @@ def interpolate(y_list: List[int]) -> Callable[[int], int]: size: int = len(y_list) matrix: Matrix = [[0 for _ in range(size)] for _ in range(size)] vector: Matrix = [[0] for _ in range(size)] - a: Matrix - i: int - y: int + coeffs: Matrix + x_val: int + y_val: int + col: int - for i, y in enumerate(y_list): - for j in range(size): - matrix[i][j] = (i + 1) ** (size - j - 1) - vector[i][0] = y + for x_val, y_val in enumerate(y_list): + for col in range(size): + matrix[x_val][col] = (x_val + 1) ** (size - col - 1) + vector[x_val][0] = y_val - a = solve(matrix, vector) + coeffs = solve(matrix, vector) - return lambda x: sum(round(a[i][0]) * (x ** (size - i - 1)) for i in range(size)) + return lambda var: sum( + round(coeffs[x_val][0]) * (var ** (size - x_val - 1)) for x_val in range(size) + ) def question_function(variable: int) -> int: @@ -175,22 +181,22 @@ def solution(func: Callable[[int], int] = question_function, order: int = 10) -> >>> solution(lambda n: n ** 3, 3) 74 """ - data_points: List[int] = [func(x) for x in range(1, order + 1)] + data_points: List[int] = [func(x_val) for x_val in range(1, order + 1)] polynomials: List[Callable[[int], int]] = [ - interpolate(data_points[:i]) for i in range(1, order + 1) + interpolate(data_points[:max_coeff]) for max_coeff in range(1, order + 1) ] ret: int = 0 - i: int - x: int + poly: int + x_val: int - for i in range(order): - x = 1 - while func(x) == polynomials[i](x): - x += 1 + for poly in polynomials: + x_val = 1 + while func(x_val) == poly(x_val): + x_val += 1 - ret += polynomials[i](x) + ret += poly(x_val) return ret From 5a16c20378f46414c4089b4888819f5218cff7e5 Mon Sep 17 00:00:00 2001 From: Freddy Pringle Date: Fri, 18 Dec 2020 11:34:53 +0100 Subject: [PATCH 6/7] Type hints --- project_euler/problem_101/sol1.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/project_euler/problem_101/sol1.py b/project_euler/problem_101/sol1.py index 9f90ac29522a..786c38f098e3 100644 --- a/project_euler/problem_101/sol1.py +++ b/project_euler/problem_101/sol1.py @@ -141,9 +141,11 @@ def interpolate(y_list: List[int]) -> Callable[[int], int]: coeffs = solve(matrix, vector) - return lambda var: sum( - round(coeffs[x_val][0]) * (var ** (size - x_val - 1)) for x_val in range(size) - ) + def interpolated_func(var: int) -> int: + return sum( + round(coeffs[x_val][0]) * (var ** (size - x_val - 1)) + for x_val in range(size) + ) def question_function(variable: int) -> int: From 7f5c25a386c9186ab8c98059b76cf2ed40ac5273 Mon Sep 17 00:00:00 2001 From: Freddy Pringle Date: Fri, 18 Dec 2020 11:39:02 +0100 Subject: [PATCH 7/7] Doctest for nested function --- project_euler/problem_101/sol1.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/project_euler/problem_101/sol1.py b/project_euler/problem_101/sol1.py index 786c38f098e3..e66316090fb2 100644 --- a/project_euler/problem_101/sol1.py +++ b/project_euler/problem_101/sol1.py @@ -142,11 +142,23 @@ def interpolate(y_list: List[int]) -> Callable[[int], int]: coeffs = solve(matrix, vector) def interpolated_func(var: int) -> int: + """ + >>> interpolate([1])(3) + 1 + >>> interpolate([1, 8])(3) + 15 + >>> interpolate([1, 8, 27])(4) + 58 + >>> interpolate([1, 8, 27, 64])(6) + 216 + """ return sum( round(coeffs[x_val][0]) * (var ** (size - x_val - 1)) for x_val in range(size) ) + return interpolated_func + def question_function(variable: int) -> int: """