結果
| 問題 |
No.453 製薬会社
|
| コンテスト | |
| ユーザー |
 toyuzuko
|
| 提出日時 | 2022-07-25 00:37:53 |
| 言語 | PyPy3 (7.3.15) |
| 結果 |
AC
|
| 実行時間 | 40 ms / 2,000 ms |
| コード長 | 13,971 bytes |
| コンパイル時間 | 438 ms |
| コンパイル使用メモリ | 82,416 KB |
| 実行使用メモリ | 57,280 KB |
| 最終ジャッジ日時 | 2024-07-06 20:25:20 |
| 合計ジャッジ時間 | 1,291 ms |
|
ジャッジサーバーID （参考情報） |
judge2 / judge3 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 4 |
| other | AC * 9 |
ソースコード
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)