最適化（線形計画法）

<?php
/****************************/
/* 線形計画法               */
/*      coded by Y.Suganuma */
/****************************/
#include 

/*************/
/* クラス LP */
/*************/
class LP
{
		private $n;   // 制約条件の数
		private $m;   // 変数の数
		private $m_s;   // スラック変数の数
		private $m_a;   // 人為変数の数
		private $mm;   // m + m_s + m_a
		private $a = array();   // 人為変数があるか否か
		private $cp = array();   // 比較演算子（-1: 左辺 < 右辺, 0: 左辺 = 右辺, 1: 左辺 > 右辺）
		private $row = array();   // 各行の基底変数の番号
		private $z = array();   // 目的関数の係数
		private $s = array();   // 単体表
		private $eps;   // 許容誤差
		private $err;   // エラーコード （0:正常終了, 1:解無し）

		/******************************/
		/* コンストラクタ             */
		/*      n1 : 制約条件の数     */
		/*      m1 : 変数の数         */
		/*      z1 : 目的関数の係数   */
		/*      eq_l : 制約条件の左辺 */
		/*      eq_r : 制約条件の右辺 */
		/*      cp1 : 比較演算子      */
		/******************************/
		function LP($n1, $m1, $z1, $eq_l, $eq_r, $cp1)
		{
					// 初期設定
			$this->eps = 1.0e-10;
			$this->err = 0;
			$this->n   = $n1;
			$this->m   = $m1;
			for ($i1 = 0; $i1 < $this->n; $i1++) {
				$this->a[$i1]  = 0;
				$this->cp[$i1] = $cp1[$i1];
			}
			for ($i1 = 0; $i1 < $this->m; $i1++)
				$this->z[$i1] = $z1[$i1];
					// スラック変数と人為変数の数を数える
			$this->m_s = 0;
			$this->m_a = 0;
			for ($i1 = 0; $i1 < $this->n; $i1++) {
				if ($this->cp[$i1] == 0) {
					$this->m_a++;
					if ($eq_r[$i1] < 0.0) {
						$eq_r[$i1] = -$eq_r[$i1];
						for ($i2 = 0; $i2 < $this->m; $i2++)
							$eq_l[$i1][$i2] = -$eq_l[$i1][$i2];
					}
				}
				else {
					$this->m_s++;
					if ($eq_r[$i1] < 0.0) {
						$this->cp[$i1]   = -$this->cp[$i1];
						$eq_r[$i1] = -$eq_r[$i1];
						for ($i2 = 0; $i2 < $this->m; $i2++)
							$eq_l[$i1][$i2] = -$eq_l[$i1][$i2];
					}
					if ($this->cp[$i1] > 0)
						$this->m_a++;
				}
			}
					// 単体表の作成
							// 初期設定
			$this->mm  = $this->m + $this->m_s + $this->m_a;
			for ($i1 = 0; $i1 <= $this->n; $i1++) {
				$this->s[$i1] = array($this->mm+1);
				if ($i1 < $this->n) {
					$this->s[$i1][0] = $eq_r[$i1];
					for ($i2 = 0; $i2 < $this->m; $i2++)
						$this->s[$i1][$i2+1] = $eq_l[$i1][$i2];
					for ($i2 = $this->m+1; $i2 <= $this->mm; $i2++)
						$this->s[$i1][$i2] = 0.0;
				}
				else {
					for ($i2 = 0; $i2 <= $this->mm; $i2++)
						$this->s[$i1][$i2] = 0.0;
				}
			}
							// スラック変数
			$k = $this->m + 1;
			for ($i1 = 0; $i1 < $this->n; $i1++) {
				if ($this->cp[$i1] != 0) {
					if ($this->cp[$i1] < 0) {
						$this->s[$i1][$k] = 1.0;
						$this->row[$i1]  = $k - 1;
					}
					else
						$this->s[$i1][$k] = -1.0;
					$k++;
				}
			}
							// 人為変数
			for ($i1 = 0; $i1 < $this->n; $i1++) {
				if ($this->cp[$i1] >= 0) {
					$this->s[$i1][$k] = 1.0;
					$this->row[$i1]  = $k - 1;
					$this->a[$i1]    = 1;
					$k++;
				}
			}
							// 目的関数
			if ($this->m_a == 0) {
				for ($i1 = 0; $i1 < $this->m; $i1++)
					$this->s[$this->n][$i1+1] = -$this->z[$i1];
			}
			else {
				for ($i1 = 0; $i1 <= $this->m+$this->m_s; $i1++) {
					for ($i2 = 0; $i2 < $this->n; $i2++) {
						if ($this->a[$i2] > 0)
							$this->s[$this->n][$i1] -= $this->s[$i2][$i1];
					}
				}
			}
		}

		/*******************************/
		/* 最適化                      */
		/*      sw : =0 : 中間結果無し */
		/*           =1 : 中間結果あり */
		/*      return : =0 : 正常終了 */
		/*             : =1 : 解無し   */
		/*******************************/
		function optimize($sw)
		{
					// フェーズ１
			if ($sw > 0) {
				if ($this->m_a > 0)
					printf("\nphase 1\n");
				else
					printf("\nphase 2\n");
			}
			$this->opt_run($sw);
					// フェーズ２
			if ($this->err == 0 && $this->m_a > 0) {
							// 目的関数の変更
				$this->mm -= $this->m_a;
				for ($i1 = 0; $i1 <= $this->mm; $i1++)
					$this->s[$this->n][$i1] = 0.0;
				for ($i1 = 0; $i1 < $this->n; $i1++) {
					$k = $this->row[$i1];
					if ($k < $this->m)
						$this->s[$this->n][0] += $this->z[$k] * $this->s[$i1][0];
				}
				for ($i1 = 0; $i1 < $this->mm; $i1++) {
					for ($i2 = 0; $i2 < $this->n; $i2++) {
						$k = $this->row[$i2];
						if ($k < $this->m)
							$this->s[$this->n][$i1+1] += $this->z[$k] * $this->s[$i2][$i1+1];
					}
					if ($i1 < $this->m)
						$this->s[$this->n][$i1+1] -= $this->z[$i1];
				}
							// 最適化
				if ($sw > 0)
					printf("\nphase 2\n");
				$this->opt_run($sw);
			}

			return $this->err;
		}

		/*******************************/
		/* 最適化（単体表の変形）      */
		/*      sw : =0 : 中間結果無し */
		/*           =1 : 中間結果あり */
		/*******************************/
		function opt_run($sw)
		{
			$this->err = -1;
			while ($this->err < 0) {
					// 中間結果
				if ($sw > 0) {
					printf("\n");
					$this->output();
				}
					// 列の選択（巡回を防ぐため必ずしも最小値を選択しない，Bland の規則）
				$q = -1;
				for ($i1 = 1; $i1 <= $this->mm && $q < 0; $i1++) {
					if ($this->s[$this->n][$i1] < -$this->eps)
						$q = $i1 - 1;
				}
					// 終了（最適解）
				if ($q < 0)
					$this->err = 0;
					// 行の選択（ Bland の規則を採用）
				else {
					$p   = -1;
					$k   = -1;
					$min = 0.0;
					for ($i1 = 0; $i1 < $this->n; $i1++) {
						if ($this->s[$i1][$q+1] > $this->eps) {
							$x = $this->s[$i1][0] / $this->s[$i1][$q+1];
							if ($p < 0 || $x < $min || $x == $min && $this->row[$i1] < $k) {
								$min = $x;
								$p   = $i1;
								$k   = $this->row[$i1];
							}
						}
					}
							// 解無し
					if ($p < 0)
						$this->err = 1;
							// 変形
					else {
						$x             = $this->s[$p][$q+1];
						$this->row[$p] = $q;
						for ($i1 = 0; $i1 <= $this->mm; $i1++)
							$this->s[$p][$i1] /= $x;
						for ($i1 = 0; $i1 <= $this->n; $i1++) {
							if ($i1 != $p) {
								$x = $this->s[$i1][$q+1];
								for ($i2 = 0; $i2 <= $this->mm; $i2++)
									$this->s[$i1][$i2] -= $x * $this->s[$p][$i2];
							}
						}
					}
				}
			}
		}

		/****************/
		/* 単体表の出力 */
		/****************/
		function output()
		{
			for ($i1 = 0; $i1 <= $this->n; $i1++) {
				if ($i1 < $this->n)
					printf("x%d", $this->row[$i1]+1);
				else
					printf(" z");
				for ($i2 = 0; $i2 <= $this->mm; $i2++)
					printf(" %f", $this->s[$i1][$i2]);
				printf("\n");
			}
		}

		/**************/
		/* 結果の出力 */
		/**************/
		function result()
		{
			if ($this->err > 0)
				printf("\n解が存在しません\n");
			else {
				printf("\n(");
				for ($i1 = 0; $i1 < $this->m; $i1++) {
					$x = 0.0;
					for ($i2 = 0; $i2 < $this->n; $i2++) {
						if ($this->row[$i2] == $i1) {
							$x = $this->s[$i2][0];
							break;
						}
					}
					if ($i1 == 0)
						printf("%f", $x);
					else
						printf(", %f", $x);
				}
				printf(") のとき，最大値 %f\n", $this->s[$this->n][0]);
			}
		}
}

					// 入力
							// 変数の数と式の数
	fscanf(STDIN, "%d %d", $m, $n);
	$cp   = array($n);
	$z    = array($m);
	$eq_l = array($n);
	$eq_r = array($n);
	for ($i1 = 0; $i1 < $n; $i1++)
		$eq_l[$i1] = array($m);
							// 目的関数の係数
	$str  = trim(fgets(STDIN));
	$z[0] = floatval(strtok($str, " "));
	for ($i1 = 1; $i1 < $m; $i1++)
		$z[$i1] = floatval(strtok(" "));
							// 制約条件
	for ($i1 = 0; $i1 < $n; $i1++) {
		$str          = trim(fgets(STDIN));
		$eq_l[$i1][0] = floatval(strtok($str, " "));
		for ($i2 = 1; $i2 < $m; $i2++)
			$eq_l[$i1][$i2] = floatval(strtok(" "));
		$c = strtok(" ");
		if ($c == '<')
			$cp[$i1] = -1;
		else if ($c == '>')
			$cp[$i1] = 1;
		else
			$cp[$i1] = 0;
		$eq_r[$i1] = floatval(strtok(" "));
	}
					// 実行
	$lp = new LP($n, $m, $z, $eq_l, $eq_r, $cp);
	$sw = $lp->optimize(1);
					// 結果の出力
	$lp->result();

?>