#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <ctime>
#include <clocale>
#include <list>

void PrintTable(const wchar_t *message, double **x, double *a, wchar_t aWchar, int aSize, double *b, wchar_t bWchar, int bSize){
    wprintf(message);
    for(int i = 0; i < bSize; i++){
        wprintf(L"\t%c(%3.f)", bWchar, b[i]);
    }
    wprintf(L"\n");
    for(int i = 0; i < aSize; i++){
        wprintf(L"%c(%3.f)", aWchar, a[i]);
        for(int j = 0; j < bSize; j++){
            wprintf(L"\t%3.f", x[i][j]);
        }
        wprintf(L"\n");
    }
    wprintf(L"\n");
}

bool CheckTask(double *a, int &aSize, double *b, int &bSize, double **c, double **x){
    double sumA = 0., sumB = 0.;
    bool closed = true;
    for(int i = 0; i < aSize; i++){
        sumA += a[i];
    }
    for(int i = 0; i < bSize; i++){
        sumB += b[i];
    }
    if(fabs(sumA - sumB) < 0.0001){
        wprintf(L"Задача закрытая\n\n");
    }
    else{
        closed = false;
        wprintf(L"Задача открытая\n");
        if(sumA > sumB){
            wprintf(L"\nВведен фиктивный потребитель\n");
            b[bSize++] = sumA - sumB;
        }
        else{
            wprintf(L"\nВведен фиктивный поставщик\n");
            a[aSize++] = sumB - sumA;
        }
    }
    return closed;
}

double GetFullPrice(int aSize, int bSize, double** x, double** c){
    double sum = 0.;
    for(int i = 0; i < aSize; i++){
        for(int j = 0; j < bSize; j++){
            if(x[i][j] > 0.){
                sum += x[i][j]*c[i][j];
            }
        }
    }
    return sum;
}

void NorthwestCornerApproximation(double* a, int aSize, double* b, int bSize, double** x){
    double *newA = (double*)malloc(aSize*sizeof(double));
    for(int i = 0; i < aSize; i++){
        newA[i] = a[i];
    }
    double *newB = (double*)malloc(bSize*sizeof(double));
    for(int i = 0; i < bSize; i++){
        newB[i] = b[i];
    }
    int aIndex = 0, bIndex = 0;
    while(aIndex < aSize && bIndex < bSize){
        if(newA[aIndex] > newB[bIndex]){
            x[aIndex][bIndex] = newB[bIndex];
            newA[aIndex] -= newB[bIndex];
            bIndex++;
        }
        else{
            x[aIndex][bIndex] = newA[aIndex];
            newB[bIndex] -= newA[aIndex];
            aIndex++;
        }
    }
    free(newA);
    free(newB);
}

void MinimalPricesApproximation(double* a, int aSize, double* b, int bSize, double** x, double **c){
    double *newA = (double*)malloc(aSize*sizeof(double));
    for(int i = 0; i < aSize; i++){
        newA[i] = a[i];
    }
    double *newB = (double*)malloc(bSize*sizeof(double));
    for(int i = 0; i < bSize; i++){
        newB[i] = b[i];
    }
    int minI, minJ;
    bool result = true;
    bool *usedI = (bool*)calloc(aSize, sizeof(bool));
    bool *usedJ = (bool*)calloc(bSize, sizeof(bool));
    while(result){
        result = false;
        for(int i = 0; i < aSize; i++){
            if(!usedI[i]){
                minI = i;
                break;
            }
        }
        for(int i = 0; i < bSize; i++){
            if(!usedJ[i]){
                minJ = i;
                break;
            }
        }
        for(int i = 0; i < aSize; i++){
            for(int j = 0; j < bSize; j++){
                if(c[minI][minJ] >= c[i][j] && !usedI[i] && !usedJ[j] && c[i][j] != 0){
                    minI = i;
                    minJ = j;
                    result = true;
                }
            }
        }
        if(!result && !c[minI][minJ] && !usedI[minI] && !usedJ[minJ]){
            result = true;
        }
        if(result){
            if(newA[minI] > newB[minJ]){
                x[minI][minJ] = newB[minJ];
                newA[minI] -= newB[minJ];
                usedJ[minJ] = true;
            }
            else{
                x[minI][minJ] = newA[minI];
                newB[minJ] -= newA[minI];
                usedI[minI] = true;
            }
        }
    }
    free(newA);
    free(newB);
    free(usedI);
    free(usedJ);
}

void GetPotentialsInLine(double* u, double *v, int line, int aSize, int bSize, double **x, double **c, bool *found);

void GetPotentialsInColumn(double* u, double *v, int column, int aSize, int bSize, double **x, double **c, bool *found);

bool SearchInLine(double **x, std::list<std::pair<int, int>> &path, int aSize, int bSize, int line, int column, int startI, int startJ);

bool SearchInColumn(double **x, std::list<std::pair<int, int>> &path, int aSize, int bSize, int line, int column, int startI, int startJ);

bool Supply(double **x, int aSize, int bSize, int startI, int startJ){
    std::list<std::pair<int, int>> path;
    path.push_back(std::make_pair(startI, startJ));
    if(!SearchInLine(x, path, aSize, bSize, startI, startJ, startI, startJ)){
        return false;
    }
    std::list<std::pair<int, int>>::iterator it;
    int sign = 1;
    double minimal = 0.;
    for(it = path.begin(); it != path.end(); ++it){
        if(++sign%2){
            if(!minimal){
                minimal = x[it->first][it->second];
            }
            else{
                minimal = std::min(minimal, x[it->first][it->second]);
            }
        }
    }
    for(it = path.begin(); it != path.end(); ++it){
        if(++sign%2){
            x[it->first][it->second] -= minimal;
        }
        else{
            x[it->first][it->second] += minimal;
        }
    }
    path.clear();
    return true;
}

void MakeSolution(double* a, int aSize, double* b, int bSize, double** x, double** c){
    double *u = (double*)malloc(aSize*sizeof(double));
    double *v = (double*)malloc(bSize*sizeof(double));
    bool *found = (bool*)malloc((aSize + bSize)*sizeof(bool));
    int minTI, minTJ;
    bool optimal, catchMin;
    double tMin;
    while(true){
        optimal = true;
        catchMin = false;
        for(int i = 0; i < aSize; i++){
            u[i] = 0.;
            found[i] = false;
        }
        for(int j = 0; j < bSize; j++){
            v[j] = 0.;
            found[aSize + j] = false;
        }
        GetPotentialsInLine(u, v, 0, aSize, bSize, x, c, found);
        wprintf(L"\n\tТаблица оценок опорного плана\n");
        for(int i = 0; i < bSize; i++){
            wprintf(L"\tv(%3.f)", v[i]);
        }
        wprintf(L"\n");
        for(int i = 0; i < aSize; i++){
            wprintf(L"u(%3.f)", u[i]);
            for(int j = 0; j < bSize; j++){
                wprintf(L"\t%3.f", c[i][j] - u[i] - v[j]);
                if(c[i][j] - u[i] - v[j] < 0){
                    optimal = false;
                    if(!catchMin){
                        tMin = c[i][j] - u[i] - v[j];
                        minTI = i;
                        minTJ = j;
                        catchMin = true;
                    }
                    else{
                        if(c[i][j] - u[i] - v[j] < tMin){
                            tMin = c[i][j] - u[i] - v[j];
                            minTI = i;
                            minTJ = j;
                        }
                    }
                }
            }
            wprintf(L"\n");
        }
        wprintf(L"\n");
        if(!optimal){
            wprintf(L"Решение не оптимально\n");
            if(Supply(x, aSize, bSize, minTI, minTJ)){
                PrintTable(L"\n\tНовый опорный план\n", x, a, 'a', aSize, b, 'b', bSize);
                wprintf(L"Общая стоимость: %.f\n", GetFullPrice(aSize, bSize, x, c));
            }
            else{
                wprintf(L"Ошибка построения пути\n");
                optimal = true;
            }
        }
        else{
            wprintf(L"Решение оптимально\n");
        }
        if(optimal){
            break;
        }
    }
    free(u);
    free(v);
    free(found);
}

void MakeRandomTask(int aSize, int bSize, int maxNum){
    double *a, *b, **x, **c;
    int lastASize = aSize;
    a = (double*)calloc(aSize + 1, sizeof(double));
    b = (double*)calloc(bSize + 1, sizeof(double));
    x = (double**)calloc(aSize + 1, sizeof(double*));
    c = (double**)calloc(aSize + 1, sizeof(double*));
    //srand(time(NULL));
    for(int i = 0; i < aSize; i++){
        if(maxNum > 20){
            a[i] = (rand()%(maxNum/10))*10 + 10;
        }
        else{
            a[i] = rand()%maxNum + 1;
        }
        x[i] = (double*)calloc(bSize + 1, sizeof(double));
        c[i] = (double*)calloc(bSize + 1, sizeof(double));
        for(int j = 0; j < bSize; j++){
            if(!i){
                if(maxNum > 20){
                    b[j] = (rand()%(maxNum/10))*10 + 10;
                }
                else{
                    b[j] = rand()%maxNum + 1;
                }
            }
            if(maxNum > 20){
                c[i][j] = (rand()%(maxNum/10))*10 + 10;
            }
            else{
                c[i][j] = rand()%maxNum + 1;
            }
        }
    }
    x[aSize] = (double*)calloc(bSize + 1, sizeof(double));
    c[aSize] = (double*)calloc(bSize + 1, sizeof(double));
    PrintTable(L"\tНачальное условие. Таблица стоимостей\n", c, a, 'a', aSize, b, 'b', bSize);
    if(!CheckTask(a, aSize, b, bSize, c, x)){
        PrintTable(L"\n\tНовое условие. Таблица стоимостей\n", c, a, 'a', aSize, b, 'b', bSize);
    }
    for(int i = 0; i < bSize; i++){
        a[0] += 0.00001;
        b[i] += 0.00001;
    }
    NorthwestCornerApproximation(a, aSize, b, bSize, x);
    PrintTable(L"\tНачальное приближение методом северо-западного угла\n", x, a, 'a', aSize, b, 'b', bSize);
    //MinimalPricesApproximation(a, aSize, b, bSize, x, c);
    //PrintTable(L"\tНачальное приближение методом минимальных стоимостей\n", x, a, 'a', aSize, b, 'b', bSize);
    wprintf(L"Общая стоимость: %.f\n", GetFullPrice(aSize, bSize, x, c));
    MakeSolution(a, aSize, b, bSize, x, c);
    for(int i = 0; i < lastASize + 1; i++){
        free(x[i]);
        free(c[i]);
    }
    free(a);
    free(b);
    free(x);
    free(c);
}

void MakeKnownTask(){
    double *a, *b, **x, **c;
    int aSize = 3, bSize = 5;
    int lastASize = aSize;
    a = (double*)calloc(aSize + 1, sizeof(double));
    b = (double*)calloc(bSize + 1, sizeof(double));
    x = (double**)calloc(aSize + 1, sizeof(double*));
    c = (double**)calloc(aSize + 1, sizeof(double*));
    for(int i = 0; i < aSize + 1; i++){
        x[i] = (double*)calloc(bSize + 1, sizeof(double));
        c[i] = (double*)calloc(bSize + 1, sizeof(double));
    }
    a[0] = 300;
    a[1] = 150;
    a[2] = 250;
    b[0] = 170;
    b[1] = 110;
    b[2] = 100;
    b[3] = 120;
    b[4] = 200;
    c[0][0] = 70;
    c[0][1] = 50;
    c[0][2] = 15;
    c[0][3] = 80;
    c[0][4] = 70;
    c[1][0] = 80;
    c[1][1] = 90;
    c[1][2] = 40;
    c[1][3] = 60;
    c[1][4] = 85;
    c[2][0] = 50;
    c[2][1] = 10;
    c[2][2] = 90;
    c[2][3] = 11;
    c[2][4] = 25;
    PrintTable(L"\tНачальное условие. Таблица стоимостей\n", c, a, 'a', aSize, b, 'b', bSize);
    if(!CheckTask(a, aSize, b, bSize, c, x)){
        PrintTable(L"\n\tНовое условие. Таблица стоимостей\n", c, a, 'a', aSize, b, 'b', bSize);
    }
    //NorthwestCornerApproximation(a, aSize, b, bSize, x);
    //PrintTable(L"\tНачальное приближение методом северо-западного угла\n", x, a, 'a', aSize, b, 'b', bSize);
    MinimalPricesApproximation(a, aSize, b, bSize, x, c);
    PrintTable(L"\tНачальное приближение методом минимальных стоимостей\n", x, a, 'a', aSize, b, 'b', bSize);
    wprintf(L"Общая стоимость: %.f\n", GetFullPrice(aSize, bSize, x, c));
    MakeSolution(a, aSize, b, bSize, x, c);
    for(int i = 0; i < lastASize + 1; i++){
        free(x[i]);
        free(c[i]);
    }
    free(a);
    free(b);
    free(x);
    free(c);
}

void GetPotentialsInLine(double* u, double *v, int line, int aSize, int bSize, double **x, double **c, bool *found){
    for(int i = 0; i < bSize; i++){
        if(x[line][i] > 0. && !v[i] && !found[aSize + i]){
            v[i] = c[line][i] - u[line];
            found[aSize + i] = true;
            GetPotentialsInColumn(u, v, i, aSize, bSize, x, c, found);
        }
    }
}

void GetPotentialsInColumn(double* u, double *v, int column, int aSize, int bSize, double **x, double **c, bool *found){
    for(int i = 1; i < aSize; i++){
        if(x[i][column] > 0. && !u[i] && !found[i]){
            u[i] = c[i][column] - v[column];
            found[i] = true;
            GetPotentialsInLine(u, v, i, aSize, bSize, x, c, found);
        }
    }
}

bool SearchInLine(double **x, std::list<std::pair<int, int>> &path, int aSize, int bSize, int line, int column, int startI, int startJ){
    for(int j = 0; j < bSize; j++){
        if(j != column && x[line][j] > 0.){
            if(j == startJ){
                path.push_back(std::make_pair(line, j));
                return true;
            }
            if(SearchInColumn(x, path, aSize, bSize, line, j, startI, startJ)){
                path.push_back(std::make_pair(line, j));
                return true;
            }
        }
    }
    return false;
}

bool SearchInColumn(double **x, std::list<std::pair<int, int>> &path, int aSize, int bSize, int line, int column, int startI, int startJ){
    for(int i = 0; i < aSize; i++){
        if(i != line && x[i][column] > 0.){
            if(SearchInLine(x, path, aSize, bSize, i, column, startI, startJ)){
                path.push_back(std::make_pair(i, column));
                return true;
            }
        }
    }
    return false;
}

int main()
{
    setlocale(LC_ALL, ".866");
    //MakeKnownTask();
    MakeRandomTask(3, 5, 250);
    return 0;
}