// MTask.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include <iostream>
#include <iomanip>
#include <vector>
#include <sstream>
#include <string>
#include <math.h>
// максимальное количество итераций метода
#define DEBUG_MAX_ITER_NUM 10
/****************************************************************************
* Обертка для вещественных чисел. Нужна для M-задачи.
* Такое число представляет две компоненты - вещественную часть
* и коэффициент при некотором числе M.
*/
class number {
private:
// вещественная часть
double real;
// коэффициет при M
double coef;
// флаг "не число"
bool is_nan;
public:
/**
* Конструктор.
* r - вещественная компонента
* c - коэффициент при M
* num - если этот аргумент будет false, то число считается как NaN
*/
number(double r = 0, double c = 0, bool num = true) {
real = r;
coef = c;
is_nan = ! num;
}
double get_real() {
return real;
}
double get_coef() {
return coef;
}
/**
* Далее идут методы перегрузки арифметических операторов.
*/
number operator+(number num) {
number res;
res.real = real + num.real;
res.coef = coef + num.coef;
return res;
}
void operator+=(number num) {
real += num.real;
coef += num.coef;
}
number operator-(number num) {
number res;
res.real = real - num.real;
res.coef = coef - num.coef;
return res;
}
void operator-=(number num) {
real -= num.real;
coef -= num.coef;
}
number operator*(number num) {
number res;
res.real = real * num.real;
if (real != 0 && coef == 0 && num.coef != 0) {
res.coef = real * num.coef;
} else if (num.real != 0 && num.coef == 0 && coef != 0) {
res.coef = coef * num.real;
}
return res;
}
void operator*=(number num) {
real *= num.real;
if (real != 0 && coef == 0 && num.coef != 0) {
coef = real * num.coef;
} else if (num.real != 0 && num.coef == 0 && coef != 0) {
coef = coef * num.real;
}
}
number operator/(number num) {
number res;
res.real = real / num.real;
return res;
}
void operator/=(number num) {
real /= num.real;
}
bool operator!=(number num) {
return (real != num.real && coef != num.coef);
}
bool operator!=(double num) {
return (real != num);
}
bool operator<(number num) {
if (coef != num.coef) {
return (coef < num.coef);
}
return (real < num.real);
}
bool operator>(number num) {
if (coef != num.coef) {
return (coef > num.coef);
}
return (real > num.real);
}
bool operator==(number num) {
return (real == num.real && coef == num.coef);
}
bool operator>=(number num) {
return this->operator>(num) || this->operator==(num);
}
bool operator<=(number num) {
return this->operator<(num) || this->operator==(num);
}
/**
* Вывод в виде строки в формате
* [R][+,-][C]M
* где R - вещественная часть,
* C - коэффициент при M.
*/
std::string str() {
if (is_nan) {
return std::string("NaN");
}
if (real == 0 && coef == 0) {
return std::string("0");
}
std::stringstream ss;
ss << std::setprecision(2);
if (real != 0) {
ss << real;
}
if (coef != 0) {
if (real != 0 && coef > 0) {
ss << '+';
} else if (coef < 0) {
ss << '-';
}
if (coef != 1 && coef != -1) {
ss << abs(coef);
}
ss << 'M';
}
return ss.str();
}
/**
* Возврашает дробную часть вещественной компоненты
*/
number fract() {
return (real > 0) ? number(real - floor(real)) : number(abs(floor(real)) + (real));
}
/**
* Проверяет является ли число (почти) целым
*/
bool is_integer() {
return this->fract() < 0.00001;
}
/**
* Модуль числа
*/
number absolute() {
return number(fabs(real));
}
};
/****************************************************************************
* Несколько тайпдефов
*/
typedef std::vector<number> vector_nums;
typedef std::vector<int> vector_ints;
typedef std::vector<std::vector<number> > matrix_nums;
/****************************************************************************
* Несколько функций для вывода в консоль
*/
void print_line() {
std::cout << std::endl;
for (int i = 0; i < 80; ++i) {
std::cout << (char) 0xCD;
}
std::cout << std::endl;
}
/**
* Печать строки
*/
void print_str(char* msg) {
std::cout << std::endl << msg << std::endl;
}
/**
* Печатает число
*/
void print_num(char* msg, number n) {
std::cout << std::endl << msg << n.str() << std::endl;
}
/**
* Печатает вещественное число
*/
void print_num(char* msg, double n) {
std::cout << std::endl << msg << n << std::endl;
}
/**
* Печатает целое число
*/
void print_num(char* msg, int n) {
std::cout << std::endl << msg << n << std::endl;
}
/**
* Печает вектор
*/
void print_vector(char* msg, vector_nums v) {
int i, j;
int width = 8;
std::cout << std::endl << msg << std::endl << (char) 0xDA;
for (i = 0; i < v.size(); ++i) {
for (j = 0; j < width - 1; ++j) {
std::cout << (char) 0xC4;
}
if (i < v.size() - 1) {
std::cout << (char) 0xC2;
}
}
std::cout << (char) 0xBF << std::endl << (char) 0xB3;
for (i = 0; i < v.size(); ++i) {
std::cout << std::setw(width - 1) << v[i].str() << (char) 0xB3;
}
std::cout << std::endl << (char) 0xC0;
for (i = 0; i < v.size(); ++i) {
for (j = 0; j < width - 1; ++j) {
std::cout << (char) 0xC4;
}
if (i < v.size() - 1) {
std::cout << (char) 0xC1;
}
}
std::cout << (char) 0xD9 << std::endl;
}
/**
* Печатает вектор целых чисел
*/
void print_vector(char* msg, std::vector<int> v) {
std::cout << std::endl << msg << std::endl << (char) 0xDA;
int i, j;
int width = 8;
for (i = 0; i < v.size(); ++i) {
for (j = 0; j < width - 1; ++j) {
std::cout << (char) 0xC4;
}
if (i < v.size() - 1) {
std::cout << (char) 0xC2;
}
}
std::cout << (char) 0xBF << std::endl << (char) 0xB3;
for (i = 0; i < v.size(); ++i) {
std::cout << std::setw(width - 1) << v[i] << (char) 0xB3;
}
std::cout << std::endl << (char) 0xC0;
for (i = 0; i < v.size(); ++i) {
for (j = 0; j < width - 1; ++j) {
std::cout << (char) 0xC4;
}
if (i < v.size() - 1) {
std::cout << (char) 0xC1;
}
}
std::cout << (char) 0xD9 << std::endl;
}
/**
* Печатает матрицу
*/
void print_matrix(char* msg, matrix_nums m) {
int i, j, k;
int width = 8;
std::cout << std::endl << msg << std::endl << (char) 0xDA;
for (i = 0; i < m[0].size(); ++i) {
for (j = 0; j < width - 1; ++j) {
std::cout << (char) 0xC4;
}
if (i < m[0].size() - 1) {
std::cout << (char) 0xC2;
}
}
std::cout << (char) 0xBF << std::endl;
for (i = 0; i < m.size(); ++i) {
std::cout << (char) 0xB3;
for (j = 0; j < m[i].size(); ++j) {
std::cout << std::setw(width - 1) << std::right << m[i][j].str() << (char) 0xB3;
}
std::cout << std::endl;
if (i < m.size() - 1) {
std::cout << (char) 0xC3;
for (j = 0; j < m[i].size(); ++j) {
for (k = 0; k < width - 1; ++k) {
std::cout << (char) 0xC4;
}
if (j < m[i].size() - 1) {
std::cout << (char) 0xC5;
} else {
std::cout << (char) 0xB4;
}
}
std::cout << std::endl;
}
}
std::cout << (char) 0xC0;
for (i = 0; i < m[0].size(); ++i) {
for (j = 0; j < width - 1; ++j) {
std::cout << (char) 0xC4;
}
if (i < m[0].size() - 1) {
std::cout << (char) 0xC1;
}
}
std::cout << (char) 0xD9 << std::endl;
std::cout << std::endl;
}
/****************************************************************************
* Исключение
*/
class simplex_error {
public:
std::string msg;
simplex_error(char* m) {
msg = std::string(m);
}
};
/****************************************************************************
* структура таблицы симплекс-метода
*/
struct opt_table {
vector_nums cib;
vector_ints bp;
vector_nums br;
vector_nums c;
matrix_nums a;
vector_nums z;
vector_nums delta;
};
/****************************************************************************
* Проверка целочисленности компонент вектора
*/
bool vector_is_integer(vector_nums v) {
for (int i = 0; i < v.size(); ++i) {
if (! v[i].is_integer()) {
return false;
}
}
return true;
}
/**
* Тип экстремума
*/
enum extremum { EXT_MIN, EXT_MAX };
/****************************************************************************
* Симплекс метод поиска максимума
*/
opt_table simplex(extremum type, opt_table task, bool negative = false, bool int_break = false) {
int i, j;
int n = task.c.size();
int m = task.bp.size();
vector_nums cib(m);
vector_nums z(n);
vector_nums delta(n);
vector_nums min(m);
matrix_nums tmp_a(m, vector_nums(n));
vector_nums tmp_br(m);
for (i = 0; i < m; ++i) {
if (task.bp[i] == -1) {
continue;
}
cib[i] = task.c[ task.bp[i] ];
}
for (int k = 0; k < DEBUG_MAX_ITER_NUM; ++k) {
print_line();
// определяем разрешающий столбец
number deltaExtr;
if (negative) {
if (type == EXT_MAX) {
deltaExtr = number(1000000);
} else {
deltaExtr = number(-1000000);
}
} else {
deltaExtr = delta[0];
}
int r = -1;
for (j = 0; j < n; ++j) {
z[j] = number();
for (i = 0; i < m; ++i) {
if (task.bp[i] == -1) {
continue;
}
z[j] += cib[i] * task.a[i][j];
}
delta[j] = task.c[j] - z[j];
if (type == EXT_MAX) {
// наименьшая по модулю отрицательная оценка
if (negative) {
if (delta[j] < number(0) && delta[j].absolute() < deltaExtr) {
deltaExtr = delta[j].absolute();
r = j;
}
// наибольшая оценка
} else {
if (deltaExtr < delta[j]) {
deltaExtr = delta[j];
r = j;
}
}
} else {
// наибольшая по модулю положительная оценка
if (negative) {
if (delta[j] >= number(0) && delta[j].absolute() > deltaExtr) {
deltaExtr = delta[j].absolute();
r = j;
}
// наименьшая оценка
} else {
if (deltaExtr > delta[j]) {
deltaExtr = delta[j];
r = j;
}
}
}
}
print_vector("Coefs of variables (c)", task.c);
print_vector("Coefs of basic variables (cib)", cib);
print_vector("Basic variables numbers (bp)", task.bp);
print_vector("Basic variables values (br)", task.br);
print_matrix("Coefs of system (a)", task.a);
print_vector("Values (z)", z);
print_vector("Relative valuations (delta)", delta);
print_num("Permitted column number (r): ", r);
if (r == -1) {
break;
}
// определяем разрешающую строку
number min_row = number(100000);
int s = -1;
for (i = 0; i < m; ++i) {
if (task.bp[i] == -1) {
s = i;
break;
}
if (task.a[i][r] != 0) {
min[i] = task.br[i] / task.a[i][r];
} else {
min[i] = number(0, 0, false);
continue;
}
if (min[i] < 0) {
min[i] = number(0, 0, false);
continue;
}
if (min_row > min[i]) {
min_row = min[i];
s = i;
}
}
//print_vector("Relations x/a (br[i]/a[i][r], min): ", min);
if (s == -1) {
throw simplex_error("Permitted row not found");
}
print_num("Permitted row number (s): ", s);
number element = task.a[s][r];
print_num("Permitted value (element): ", element);
for (i = 0; i < m; ++i) {
for (j = 0; j < n; ++j) {
tmp_a[i][j] = task.a[i][j];
}
tmp_br[i] = task.br[i];
}
// вносим переменную в базис
task.bp[s] = r;
cib[s] = task.c[r];
for (j = 0; j < n; ++j) {
task.a[s][j] /= element;
}
task.br[s] /= element;
for (i = 0; i < m; ++i) {
if (i == s) {
continue;
}
number air = tmp_a[i][r];
for (j = 0; j < n; ++j) {
task.a[i][j] -= (air * tmp_a[s][j]) / element;
}
task.br[i] -= (air * tmp_br[s]) / element;
}
if (int_break && vector_is_integer(task.br)) {
break;
}
}
print_line();
print_vector("Result: ", task.br);
opt_table tab;
tab.cib = cib;
tab.bp = task.bp;
tab.br = task.br;
tab.c = task.c;
tab.a = task.a;
tab.z = z;
tab.delta = delta;
return tab;
}
/****************************************************************************
* Сравнение по модулю 1
*/
bool cmp_mod_one(number a, number b) {
return (a - b).is_integer();
}
/****************************************************************************
* Проверка значения элемента в массиве
*/
bool in_vector(int a, vector_ints v) {
for (int i = 0; i < v.size(); ++i) {
if (a == v[i]) {
return true;
}
}
return false;
}
/****************************************************************************
* Метод Гомори
*/
opt_table gomory(extremum type, opt_table tab) {
if (type == EXT_MIN) {
print_str("Minimum");
} else {
print_str("Maximum");
}
tab = simplex(type, tab, false, false);
int i, j;
// общее количество переменных
int n = tab.c.size();
// количество базисных переменных
int m = tab.bp.size();
for (int k = 0; k < DEBUG_MAX_ITER_NUM; ++k) {
print_line();
std::cout << "Gomory iteration #" << k << std::endl;
print_line();
// условие выхода
if (vector_is_integer(tab.br)) {
return tab;
}
// максимальная дробная часть
number max_fract = 0;
// номер строки с максимальной дробной частью
int r = -1;
// поиск переменной с максимальной дробной частью
for (i = 0; i < m; ++i) {
number fract = tab.br[i].fract();
if (fract >= max_fract) {
max_fract = fract;
r = i;
}
}
// составляем новую таблицу, добавляя один столбец и одну строку
++n;
++m;
opt_table ext_tab;
// вектор коэффициентов всех переменных
ext_tab.c = vector_nums(n);
for (j = 0; j < n - 1; ++j) {
ext_tab.c[j] = tab.c[j];
}
ext_tab.c[n - 1] = 0;
// матрица a
ext_tab.a = matrix_nums(m, vector_nums(n));
for (i = 0; i < m - 1; ++i) {
for (j = 0; j < n - 1; ++j) {
ext_tab.a[i][j] = tab.a[i][j];
}
ext_tab.a[i][n - 1] = 0;
}
for (j = 0; j < n - 1; ++j) {
if (in_vector(j, tab.bp)) {
ext_tab.a[m - 1][j] = 0;
} else {
ext_tab.a[m - 1][j] = tab.a[r][j].fract();
}
}
ext_tab.a[m - 1][n - 1] = -1;
// вектор-решение
ext_tab.bp = vector_ints(m);
ext_tab.br = vector_nums(m);
for (i = 0; i < m - 1; ++i) {
//ext_tab.cib[i] = tab.cib[i];
ext_tab.bp[i] = tab.bp[i];
ext_tab.br[i] = tab.br[i];
}
// ext_tab.cib[m - 1] = ?
ext_tab.bp[m - 1] = -1; // ?
/*
if (type == EXT_MAX) {
ext_tab.bp[m - 1] = -1;
} else {
ext_tab.bp[m - 1] = 1;
}
*/
ext_tab.br[m - 1] = tab.br[r].fract();
tab = simplex(type, ext_tab, true, true);
}
return tab;
}
/****************************************************************************
* Мейн, ага
*/
int main(int argc, char* argv[]) {
int n = 5;
int m = 2;
opt_table task;
task.c = vector_nums(n);
task.a = matrix_nums(m, vector_nums(n));
task.bp = vector_ints(m);
task.br = vector_nums(m);
task.c[0] = number(2);
task.c[1] = number(1);
task.c[2] = number(0);
task.a[0][0] = number(15);
task.a[0][1] = number(30);
task.a[0][2] = number(1);
task.bp[0] = 2;
task.br[0] = number(96);
/*task.c[0] = number(1);
task.c[1] = number(-1);
task.c[2] = number(0);
task.c[3] = number(0);
task.c[4] = number(0, -1);
task.a[0][0] = number(-1);
task.a[0][1] = number(2);
task.a[0][2] = number(-1);
task.a[0][3] = number(0);
task.a[0][4] = number(1);
task.a[1][0] = number(3);
task.a[1][1] = number(2);
task.a[1][2] = number(0);
task.a[1][3] = number(1);
task.a[1][4] = number(0);
task.bp[0] = 4;
task.bp[1] = 3;
task.br[0] = number(4);
task.br[1] = number(14);*/
/**/
/*task.c[0] = number(-1);
task.c[1] = number(2);
task.c[2] = number(-1);
task.c[3] = number(-1);
task.a[0][0] = number(-1);
task.a[0][1] = number(1);
task.a[0][2] = number(1);
task.a[0][3] = number(0);
task.a[1][0] = number(1);
task.a[1][1] = number(1);
task.a[1][2] = number(0);
task.a[1][3] = number(1);
task.bp[0] = 2;
task.bp[1] = 3;
task.br[0] = number(2);
task.br[1] = number(4);*/
gomory(EXT_MIN, task);
return 0;
}