iframe-proxy | Sunbelt Computer Software

History

219 lines (177 loc) · 6.28 KB

Raw

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

"""

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(matrix: Matrix, vector: Matrix) -> Matrix:

"""

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(matrix)

augmented: Matrix = [[0 for _ in range(size + 1)] for _ in range(size)]

row: int

row2: int

col: int

col2: int

pivot_row: int

ratio: float

for row in range(size):

for col in range(size):

augmented[row][col] = matrix[row][col]

augmented[row][size] = vector[row][0]

row = 0

col = 0

while row < size and col < size:

# pivoting

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[pivot_row] = augmented[pivot_row], augmented[row]

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

# back substitution

for col in range(1, size):

for row in range(col):

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 [

[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)

>>> interpolate([1, 8])(3)

>>> interpolate([1, 8, 27])(4)

>>> interpolate([1, 8, 27, 64])(6)

216

"""

size: int = len(y_list)

matrix: Matrix = [[0 for _ in range(size)] for _ in range(size)]

vector: Matrix = [[0] for _ in range(size)]

coeffs: Matrix

x_val: int

y_val: int

col: int

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

coeffs = solve(matrix, vector)

def interpolated_func(var: int) -> int:

"""

>>> interpolate([1])(3)

>>> interpolate([1, 8])(3)

>>> interpolate([1, 8, 27])(4)

>>> 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:

"""

The generating function u as specified in the question.

>>> question_function(0)

>>> question_function(1)

>>> question_function(5)

8138021

>>> question_function(10)

9090909091

"""

return (

- 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] = 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

not equal u(x).

>>> solution(lambda n: n ** 3, 3)

"""

data_points: List[int] = [func(x_val) for x_val in range(1, order + 1)]

polynomials: List[Callable[[int], int]] = [

interpolate(data_points[:max_coeff]) for max_coeff in range(1, order + 1)

]

ret: int = 0

poly: int

x_val: int

for poly in polynomials:

x_val = 1

while func(x_val) == poly(x_val):

x_val += 1

ret += poly(x_val)

return ret

if __name__ == "__main__":

print(f"{solution() = }")

Sunbelt Computer Software

PL/B Language Development and Support

Navigation Menu

Sunbelt Computer Software

PL/B Language Development and Support

FilesExpand file tree

sol1.py

Latest commit

History

sol1.py

File metadata and controls