EPS = 1e-10 def is_same(x, y): if abs(x - y) < EPS: return True return False def is_zero(x): if abs(x) < EPS: return True return False class Matrix(): def __init__(self, n, m, mat=None): self.n = n self.m = m self.mat = [[0.] * self.m for _ in range(self.n)] if mat: assert len(mat) == n and len(mat[0]) == m for i in range(self.n): self.mat[i] = mat[i].copy() def is_square(self): return self.n == self.m def __getitem__(self, key): if isinstance(key, slice): return self.mat[key] else: assert key >= 0 return self.mat[key] def id(n): res = Matrix(n, n) for i in range(n): res[i][i] = 1. return res def __len__(self): return len(self.mat) def __str__(self): return '\n'.join(' '.join(map(str, self[i])) for i in range(self.n)) def times(self, k): res = [[0.] * self.m for _ in range(self.n)] for i in range(self.n): res_i, self_i = res[i], self[i] for j in range(self.m): res_i[j] = k * self_i[j] return Matrix(self.n, self.m, res) def __pos__(self): return self def __neg__(self): return self.times(-1) def __add__(self, other): assert self.n == other.n and self.m == other.m res = [[0.] * self.m for _ in range(self.n)] for i in range(self.n): res_i, self_i, other_i = res[i], self[i], other[i] for j in range(self.m): res_i[j] = (self_i[j] + other_i[j]) return Matrix(self.n, self.m, res) def __sub__(self, other): assert self.n == other.n and self.m == other.m res = [[0.] * self.m for _ in range(self.n)] for i in range(self.n): res_i, self_i, other_i = res[i], self[i], other[i] for j in range(self.m): res_i[j] = (self_i[j] - other_i[j]) return Matrix(self.n, self.m, res) def __mul__(self, other): if other.__class__ == Matrix: assert self.m == other.n res = [[0.] * other.m for _ in range(self.n)] for i in range(self.n): res_i, self_i = res[i], self[i] for k in range(self.m): self_ik, other_k = self_i[k], other[k] for j in range(other.m): res_i[j] += self_ik * other_k[j] return Matrix(self.n, other.m, res) elif other.__class__ == Vector: assert self.m == other.n res = [0.] * self.n for i in range(self.n): for j in range(self.m): res[i] += self[i][j] * other[j] return Vector(self.n, res) else: return self.times(other) def __rmul__(self, other): if other.__class__ == Vector: #転置をとってかける assert other.n == self.n res = [0.] * self.m for i in range(self.m): for j in range(self.n): res[i] += self[j][i] * other[j] return Vector(self.m, res) return self.times(other) def __pow__(self, k): assert self.is_square() tmp = Matrix(self.n, self.n, self.mat) res = Matrix.id(self.n) while k: if k & 1: res *= tmp tmp *= tmp k >>= 1 return res def determinant(self): assert self.is_square() res = 1 tmp = Matrix(self.n, self.n, self.mat) for j in range(self.n): if is_zero(tmp[j][j]): for i in range(j + 1, self.n): if not is_zero(tmp[i][j]): break else: return 0 tmp.mat[j], tmp.mat[i] = tmp.mat[i], tmp.mat[j] res *= -1 tmp_j = tmp[j] inv = 1 / tmp_j[j] for i in range(j + 1, self.n): tmp_i = tmp[i] c = -inv * tmp_i[j] for k in range(self.n): tmp_i[k] += c * tmp_j[k] for i in range(self.n): res *= tmp[i][i] return res def inverse(self): assert self.is_square() res = Matrix.id(self.n) tmp = Matrix(self.n, self.n, self.mat) for j in range(self.n): if is_zero(tmp[j][j]): for i in range(j + 1, self.n): if not is_zero(tmp[i][j]): break else: return -1 tmp.mat[j], tmp.mat[i] = tmp.mat[i], tmp.mat[j] res.mat[j], res.mat[i] = res.mat[i], res.mat[j] tmp_j, res_j = tmp[j], res[j] inv = 1 / tmp_j[j] for k in range(self.n): tmp_j[k] *= inv res_j[k] *= inv for i in range(self.n): if i == j: continue c = tmp[i][j] tmp_i, res_i = tmp[i], res[i] for k in range(self.n): tmp_i[k] -= tmp_j[k] * c res_i[k] -= res_j[k] * c return res def linear_equations(self, vec): assert self.n == len(vec) aug = [self[i] + [vec[i]] for i in range(self.n)] rank = 0 p = [] q = [] for j in range(self.m + 1): for i in range(rank, self.n): if not is_zero(aug[i][j]): break else: q.append(j) continue if j == self.m: return -1, [], [] p.append(j) aug[rank], aug[i] = aug[i], aug[rank] inv = 1 / aug[rank][j] aug_rank = aug[rank] for k in range(self.m + 1): aug_rank[k] *= inv for i in range(self.n): if i == rank: continue aug_i = aug[i] c = -aug_i[j] for k in range(self.m + 1): aug_i[k] += c * aug_rank[k] rank += 1 dim = self.m - rank sol = [0.] * self.m for i in range(rank): sol[p[i]] = aug[i][-1] vecs = [[0.] * self.m for _ in range(dim)] for i in range(dim): vecs[i][q[i]] = 1. for i in range(dim): vecs_i = vecs[i] for j in range(rank): vecs_i[p[j]] = -aug[j][q[i]] return dim, sol, vecs def transpose(self): res = [[0.] * self.n for _ in range(self.m)] for i in range(self.n): for j in range(self.m): res[j][i] = self[i][j] return Matrix(self.m, self.n, res) def column(self, k): res = [0.] * self.n for i in range(self.n): res[i] = self[i][k] return Vector(self.n, res) def row(self, k): return Vector(self.m, self[k]) def is_regular(self): if is_zero(self.determinant()): return False return True class Vector(): def __init__(self, n, vec=None): self.n = n self.vec = [0.] * n if vec: assert len(vec) == n self.vec = vec.copy() def __getitem__(self, key): if isinstance(key, slice): return self.vec[key] else: assert key >= 0 return self.vec[key] def __setitem__(self, key, val): assert key >= 0 self.vec[key] = val def __len__(self): return len(self.vec) def __str__(self): return ' '.join(map(str, self.vec)) def __pos__(self): return self def __neg__(self): return self.times(-1) def __add__(self, other): assert self.n == other.n res = [0.] * self.n for i in range(self.n): res[i] = self[i] + other[i] return Vector(self.n, res) def __sub__(self, other): assert self.n == other.n res = [0.] * self.n for i in range(self.n): res[i] = self[i] - other[i] return Vector(self.n, res) def __mul__(self, other): if other.__class__ == Vector: #内積 assert self.n == other.n res = 0 for i in range(self.n): res += self[i] * other[i] return res else: return self.times(other) def __rmul__(self, other): return self.times(other) def times(self, k): res = [0.] * self.n for i in range(self.n): res[i] = self[i] * k return Vector(self.n, res) class SimplexMethod(): # Ax <= b, x >= 0 の条件下で c * x を最大化する def __init__(self, c, A, b): self.n = len(c) + len(b) #変数の個数(スラック変数を含む)もとの変数は n - m 個 self.m = len(b) #制約の個数 self.A = A self.b = b self.c_B = Vector(self.m) self.c_N = Vector(self.n - self.m, c.vec) self.B = Matrix.id(self.m) self.N = Matrix(self.m, self.n - self.m, A.mat) self.B_idx = list(range(self.n - self.m, self.n)) self.N_idx = list(range(self.n - self.m)) self.xB = None self.x_opt = None def simplex(self): while True: invB = self.B.inverse() b_ = invB * self.b c_N_ = self.c_N - self.N.transpose() * (invB.transpose() * self.c_B) for i in range(self.n - self.m): if c_N_[i] > 0: k = i # c_N_ が正になる一番最初の添字 break else: self.xB = b_ return 1 # c_N_ <= 0 のとき最適解となる a_k = invB * self.N.column(k) for i in range(self.m): if a_k[i] > 0: break else: return 0 # 非有界 theta_min = float('inf') for i in range(self.m): if a_k[i] > 0: theta = b_[i] / a_k[i] if theta_min > theta: theta_min = theta l = i # theta が最小になる一番最初の添字 self.swapNB(l, k) def swapNB(self, l, k): # 基底変数 l と非基底変数 k を入れ替える self.N_idx[k], self.B_idx[l] = self.B_idx[l], self.N_idx[k] self.c_N[k], self.c_B[l] = self.c_B[l], self.c_N[k] N_k = self.N.column(k) B_l = self.B.column(l) for i in range(self.m): self.N[i][k] = B_l[i] self.B[i][l] = N_k[i] def find_feasible_solution(self): b_min = min(self.b) if b_min >= 0: return 1 c_aux = Vector(self.n - self.m + 1, [-1] + [0] * (self.n - self.m)) A_aux_mat = [[0.] * (self.n - self.m + 1) for _ in range(self.m)] for i in range(self.m): A_aux_mat[i][0] = -1 for j in range(self.n - self.m): A_aux_mat[i][j + 1] = self.A[i][j] A_aux = Matrix(self.m, self.n - self.m + 1, A_aux_mat) S_aux = SimplexMethod(c_aux, A_aux, b) k = self.b.vec.index(b_min) S_aux.swapNB(k, 0) assert S_aux.simplex() # 非有界なことはあり得る? z_aux = S_aux.c_B * S_aux.xB if z_aux < 0: return 0 # 実行可能解がない if not 0 in S_aux.N_idx: # 退化して 0 が基底変数となっていた場合の処理は? r = S_aux.B_idx.index(0) for i in range(self.n - self.m + 1): S_aux.swapNB(r, i) if S_aux.B.is_regular(): break S_aux.swapNB(r, i) N_idx_aux = [S_aux.N_idx[i] - 1 for i in range(self.n - self.m + 1) if S_aux.N_idx[i] != 0] B_idx_aux = [S_aux.B_idx[i] - 1 for i in range(self.m)] idx_aux = N_idx_aux + B_idx_aux c_NB = self.c_N.vec + self.c_B.vec c_NB_sorted = [c_NB[idx_aux[i]] for i in range(self.n)] self.c_N = Vector(self.n - self.m, c_NB_sorted[:self.n - self.m]) self.c_B = Vector(self.m, c_NB_sorted[self.n - self.m:]) NB_mat = [[0.] * self.n for _ in range(self.m)] for i in range(self.m): for j in range(self.n - self.m): NB_mat[i][j] = self.N[i][j] for j in range(self.m): NB_mat[i][j + self.n - self.m] = self.B[i][j] NB_mat_sorted = [[0.] * self.n for _ in range(self.m)] for i in range(self.m): for j in range(self.n): NB_mat_sorted[i][j] = NB_mat[i][idx_aux[j]] N_mat = [[0.] * (self.n - self.m) for _ in range(self.m)] B_mat = [[0.] * self.m for _ in range(self.m)] for i in range(self.m): for j in range(self.n - self.m): N_mat[i][j] = NB_mat_sorted[i][j] for j in range(self.m): B_mat[i][j] = NB_mat_sorted[i][j + self.n - self.m] self.N = Matrix(self.m, self.n - self.m, N_mat) self.B = Matrix(self.m, self.m, B_mat) self.N_idx = N_idx_aux self.B_idx = B_idx_aux return 1 def solve(self): if not self.find_feasible_solution(): raise Exception('Infeasible') if not self.simplex(): raise Exception('Unbounded') z = self.c_B * self.xB x = Vector(self.n, self.xB.vec + [0] * (self.n - self.m)) x_sorted = [0.] * (self.n) idx = self.B_idx + self.N_idx for i in range(self.n): x_sorted[idx[i]] = x[i] self.x_opt = Vector(self.n - self.m, x_sorted[:self.n - self.m]) return z, self.x_opt C, D = map(int, input().split()) c = Vector(2, [1000, 2000]) A = Matrix(2, 2, [[3 / 4, 2 / 7], [1 / 4, 5 / 7]]) b = Vector(2, [C, D]) S = SimplexMethod(c, A, b) z, x = S.solve() print(z)