線形計画法

#***************************/
# 線形計画法               */
#      coded by Y.Suganuma */
#***************************/

#************/
# クラス LP */
#************/
class LP
	#*****************************/
	# コンストラクタ             */
	#      n1 : 制約条件の数     */
	#      m1 : 変数の数         */
	#      z1 : 目的関数の係数   */
	#      eq_l : 制約条件の左辺 */
	#      eq_r : 制約条件の右辺 */
	#      cp1 : 比較演算子      */
	#*****************************/
	def initialize(n1, m1, z1, eq_l, eq_r, cp1)
				# 初期設定
		@_eps = 1.0e-10   # 許容誤差
		@_err = 0   # エラーコード （0:正常終了, 1:解無し）
		@_n   = n1   # 制約条件の数
		@_m   = m1   # 変数の数
		@_a   = Array.new(@_n)   # 人為変数があるか否か
		@_cp  = Array.new(@_n)   # 比較演算子（-1: 左辺 < 右辺, 0: 左辺 = 右辺, 1: 左辺 > 右辺）
		for i1 in 0 ... @_n
			@_a[i1]  = 0
			@_cp[i1] = cp1[i1]
		end
		@_z = Array.new(@_m)   # 目的関数の係数
		for i1 in 0 ... @_m
			@_z[i1] = z1[i1]
		end
				# スラック変数と人為変数の数を数える
		@_m_s = 0   # スラック変数の数
		@_m_a = 0   # 人為変数の数
		for i1 in 0 ... @_n
			if @_cp[i1] == 0
				@_m_a += 1
				if eq_r[i1] < 0.0
					eq_r[i1] = -eq_r[i1]
					for i2 in 0 ... @_m
						eq_l[i1][i2] = -eq_l[i1][i2]
					end
				end
			else
				@_m_s += 1
				if eq_r[i1] < 0.0
					@_cp[i1] = -@_cp[i1]
					eq_r[i1] = -eq_r[i1]
					for i2 in 0 ... @_m
						eq_l[i1][i2] = -eq_l[i1][i2]
					end
				end
				if @_cp[i1] > 0
					@_m_a += 1
				end
			end
		end
				# 単体表の作成
						# 初期設定
		@_mm  = @_m + @_m_s + @_m_a   # m + m_s + m_a
		@_row = Array.new(@_n)   # 各行の基底変数の番号
		@_s   = Array.new(@_n+1)   # 単体表
		for i1 in 0 ... @_n+1
			@_s[i1] = Array.new(@_mm+1)
			if i1 < @_n
				@_s[i1][0] = eq_r[i1]
				for i2 in 0 ... @_m
					@_s[i1][i2+1] = eq_l[i1][i2]
				end
				for i2 in @_m+1 ... @_mm+1
					@_s[i1][i2] = 0.0
				end
			else
				for i2 in 0 ... @_mm+1
					@_s[i1][i2] = 0.0
				end
			end
		end
						# スラック変数
		k = @_m + 1
		for i1 in 0 ... @_n
			if @_cp[i1] != 0
				if @_cp[i1] < 0
					@_s[i1][k] = 1.0
					@_row[i1]  = k - 1
				else
					@_s[i1][k] = -1.0
				end
				k += 1
			end
		end
						# 人為変数
		for i1 in 0 ... @_n
			if @_cp[i1] >= 0
				@_s[i1][k] = 1.0
				@_row[i1]  = k - 1
				@_a[i1]    = 1
				k += 1
			end
		end
						# 目的関数
		if @_m_a == 0
			for i1 in 0 ... @_m
				@_s[@_n][i1+1] = -@_z[i1]
			end
		else
			for i1 in 0 ... @_m+@_m_s+1
				for i2 in 0 ... @_n
					if @_a[i2] > 0
						@_s[@_n][i1] -= @_s[i2][i1]
					end
				end
			end
		end
	end

	#******************************/
	# 最適化                      */
	#      sw : =0 : 中間結果無し */
	#           =1 : 中間結果あり */
	#      return : =0 : 正常終了 */
	#             : =1 : 解無し   */
	#******************************/
	def optimize(sw)
				# フェーズ１
		if sw > 0
			if @_m_a > 0
				printf("\nphase 1\n")
			else
				printf("\nphase 2\n")
			end
		end
		opt_run(sw)
				# フェーズ２
		if @_err == 0 && @_m_a > 0
						# 目的関数の変更
			@_mm -= @_m_a
			for i1 in 0 ... @_mm+1
				@_s[@_n][i1] = 0.0
			end
			for i1 in 0 ... @_n
				k = @_row[i1]
				if k < @_m
					@_s[@_n][0] += @_z[k] * @_s[i1][0]
				end
			end
			for i1 in 0 ... @_mm
				for i2 in 0 ... @_n
					k = @_row[i2]
					if k < @_m
						@_s[@_n][i1+1] += @_z[k] * @_s[i2][i1+1]
					end
				end
				if i1 < @_m
					@_s[@_n][i1+1] -= @_z[i1]
				end
			end
						# 最適化
			if sw > 0
				printf("\nphase 2\n")
			end
			opt_run(sw)
		end

		return @_err
	end

	#******************************/
	# 最適化（単体表の変形）      */
	#      sw : =0 : 中間結果無し */
	#           =1 : 中間結果あり */
	#******************************/
	def opt_run(sw)
		@_err = -1
		while @_err < 0
				# 中間結果
			if sw > 0
				printf("\n")
				output()
			end
				# 列の選択（巡回を防ぐため必ずしも最小値を選択しない，Bland の規則）
			q = -1
			for i1 in 1 ... @_mm+1
				if @_s[@_n][i1] < -@_eps
					q = i1 - 1
				end
				if q >= 0
					break
				end
			end
				# 終了（最適解）
			if q < 0
				@_err = 0
				# 行の選択（ Bland の規則を採用）
			else
				p   = -1
				k   = -1
				min = 0.0
				for i1 in 0 ... @_n
					if @_s[i1][q+1] > @_eps
						x = @_s[i1][0] / @_s[i1][q+1]
						if p < 0 || x < min || x == min && @_row[i1] < k
							min = x
							p   = i1
							k   = @_row[i1]
						end
					end
				end
						# 解無し
				if p < 0
					@_err = 1
						# 変形
				else
					x        = @_s[p][q+1]
					@_row[p] = q
					for i1 in 0 ... @_mm+1
						@_s[p][i1] /= x
					end
					for i1 in 0 ... @_n+1
						if i1 != p
							x = @_s[i1][q+1]
							for i2 in 0 ... @_mm+1
								@_s[i1][i2] -= x * @_s[p][i2]
							end
						end
					end
				end
			end
		end
	end

	#***************/
	# 単体表の出力 */
	#***************/
	def output()
		for i1 in 0 ... @_n+1
			if i1 < @_n
				printf("x%d", @_row[i1]+1)
			else
				printf(" z")
			end
			for i2 in 0 ... @_mm+1
				printf(" %f", @_s[i1][i2])
			end
			printf("\n")
		end
	end

	#*************/
	# 結果の出力 */
	#*************/
	def result()
		if @_err > 0
			printf("\n解が存在しません\n")
		else
			printf("\n(")
			for i1 in 0 ... @_m
				x = 0.0
				for i2 in 0 ... @_n
					if @_row[i2] == i1
						x = @_s[i2][0]
						break
					end
				end
				if i1 == 0
					printf("%f", x)
				else
					printf(", %f", x)
				end
			end
			printf(") のとき，最大値 %f\n", @_s[@_n][0])
		end
	end
end

				# 入力
						# 変数の数と式の数
str  = gets();
a    = str.split(" ");
m    = Integer(a[0]);
n    = Integer(a[1]);
cp   = Array.new(n)
z    = Array.new(m)
eq_r = Array.new(n)
eq_l = Array.new(n)
for i1 in 0 ... n
	eq_l[i1] = Array.new(m)
end
						# 目的関数の係数
str = gets();
a   = str.split(" ");
for i1 in 0 ... m
	z[i1] = Float(a[i1])
end
						# 制約条件
for i1 in 0 ... n
	str = gets();
	a   = str.split(" ");
	for i2 in 0 ... m
		eq_l[i1][i2] = Float(a[i2])
	end
	if a[m] == '<'
		cp[i1] = -1
	elsif a[m] == '>'
		cp[i1] = 1
	else
		cp[i1] = 0
	end
	eq_r[i1] = Float(a[m+1])
end
				# 実行
lp = LP.new(n, m, z, eq_l, eq_r, cp)
sw = lp.optimize(1)
				# 結果の出力
lp.result()