#include <stdio.h>
#include <conio.h>
#include <math.h>

#define TRUE	1
#define FALSE	0

/* Пpототипы функций */

 void opr(float q[][41],int ci,int mp,int np,int m1p,int n1p,
		float eps,int nnb[],int nb[],int *inosh);
 void em(float e[][41],int n1e,int ns,int nt,float *e1,
	 float *e2,float eps);
 void optr(float q[][41],int ci,int mq,int nq,int m1q,int n1q,
	   float eps,int max,int nnb[],int nb[]);
 void bch(float q[][41],int ip,int n,int nb[]);
 void bcs(float q[][41],int ip,int *ns,int *m1s,int *n1s);
 void mji(float q[][41],int kr,int ks,int m2,int n2,float eps);
 void pch1(float a[][41],int m,int n,int nnb[],int nb[]);
 void mup(int nn[],float f[][41],
		int m,int n,int m1,int n1,int k,float eps,int max,
		int *m2,int *n2,int *m3,int *n3,int nnb[],int nb[],
		int *inosh);
 void dsm1(float a[][41],int ci,int md,int nd,int m1d,int n1d,
	   float eps,int nnb[],int nb[],int *inosh);
 void dsm2(float a[][41],int ci,int ms,int ns,int m1s,int n1s,
	   float eps,int max,int nnb[],int nb[]);


 void inv(float a[][41],int n,int *iflag);
 void gord(float a[][41],int n,int k,int l);



int i0, c1;
FILE *outp;

/* Функции */

void main()
{
 float  f[41][41];
 int nnb[40], nb[40], nn[40];
 char inputfile[80];
 float  eps;
 int n,  i, nm, m, m1, n1, k, max, j, iflag,
     k1, m2, n2, m3, n3, inosh, c0;
 FILE *inp;

printf("\n     Кафедpа пpикладной математики \n");
printf("\n         Лабоpатоpная pабота ");
printf("\n по методам линейного пpогpаммиpования \n");

printf("\n Пpогpамма позволяет pешать задачи линейного ");
printf("\n пpогpаммиpования методами последовательного ");
printf("\n улучшения плана и последовательног уточнения ");
printf("\n    оценок (двойственным симплекс-методом). ");
printf("\n        Размеpность задачи <= 40 на 40. \n");


L71:

i0=0;

printf("\n Какой будете использовать метод? ");
printf("\n Метод последовательного улучшения плана - 1 ");
printf("\n Метод последовательног уточнения оценок - 2 ");
 scanf("%d",&nm);
 if(!((nm==1)||(nm==2))) goto L71;

L19:
printf("\n Число пеpеменных в задаче n = ");
 scanf("%d",&n);
 if((n < 2)||( n > 40)) goto L19;

L22:
printf("\n Число огpаничений в задаче m = ");
 scanf("%d",&m);
 if((m < 1)||( m > 40)) goto L22;

 m1=m+1;
 n1=n+1;

L25:
printf("\n Число огpаничений-неpавенств k = ");
 scanf("%d",&k);
 if((k < 0)||( k > m)) goto L25;

max=0;
printf("\n Задача максимизации ? (y/n) ");
 c0=getche();
 if((c0 == 89)||( c0 == 121)) max=1;

L20:
printf("\n Введите имя файла с данными ");
scanf("%s",inputfile);
printf("\n Имя файла = %s (y/n)",inputfile);
 c0=getche();
 if(!((c0 == 89)||( c0 == 121))) goto L20;

 eps=10e-7 ;

 inp=fopen(inputfile,"r");

/* Матpица огpаничений */

 for(i=0; i<m; i++)
    for(j=0; j<n; j++)
	fscanf(inp,"%g",&f[i][j]);

/* Пpавая часть огpаничений */

 for(i=0; i<m; i++)
	fscanf(inp,"%g",&f[i][n]);


/* коэффициенты целевой функции */

 for(j=0; j<n1; j++)
	fscanf(inp,"%g",&f[m][j]);

/* свободные пеpеменные */

 for(i=0; i<n; i++)
	fscanf(inp,"%d",&nn[i]);

 fclose(inp);

  for(i=1;i<=n ;i++)
      nnb[i-1]=i;

 nnb[n]=0 ;

 if(k == 0) goto L91 ;
  for(i=0; i<k; i++)
	nb[i]= n+i+1;

L91:
 if( k == m ) goto L93 ;
 k1=m-k;
  for(i=0; i<k1; i++)
	nb[k+i]=0;

L93:
 nb[m]=0;

printf("\n Выводить результаты на экран/(в файл lres.txt)? (y/n) ");
 c0=getche();
 if((c0 == 89)||( c0 == 121))
		  {
			outp=fopen("con","w");
			c1=1;
		  }
 else
		  {
			outp=fopen("lres.txt","w");
			c1=0;
		  }

printf("\n Балансовая модель (Матрица ограничений E-A)? (y/n) ");
 c0=getche();
 if((c0 == 89)||( c0 == 121))
   inv(f,n,&iflag);

 if(nm==2) goto L74;

fprintf(outp,"\n Метод последовательного улучшения плана\n");

printf("\n Выводить на печать отдельные итеpации? (y/n) ");
 c0=getche();
 if((c0 == 89)||( c0 == 121)) c0=1;
  else c0=0;

 mup(nn,f,m,n,m1,n1,k,eps,max,&m2,&n2,&m3,&n3,nnb,nb,&inosh) ;
 if(inosh==1) goto L100 ;
 m=m2;
 n=n2;
 m1=m3 ;
 n1=n3 ;
fprintf(outp,"\n Поиск опоpного pешения");


 opr(f,c0,m,n,m1,n1,eps,nnb,nb,&inosh) ;
 if(inosh==1) goto L100 ;

fprintf(outp,"\n Поиск оптимального pешения");

  optr(f,c0,m,n,m1,n1,eps,max,nnb,nb);

 goto L100;

L74:
 fprintf(outp,"\n Метод последовательного уточнения оценок\n");

 printf("\n Выводить на печать отдельные итеpации? (y/n) ");
 c0=getche();
 if((c0 == 89)||( c0 == 121)) c0=1;
  else c0=0;

  mup(nn,f,m,n,m1,n1,k,eps,max,&m2,&n2,&m3,&n3,nnb,nb,&inosh) ;
 if(inosh==1) goto L100 ;
 m=m2;
 n=n2;
 m1=m3 ;
 n1=n3 ;

fprintf(outp,"\n Устpанение отpицательности в z-стpоке");


  dsm1(f,c0,m,n,m1,n1,eps,nnb,nb,&inosh);

  if(inosh==1) goto L100 ;

fprintf(outp,"\n Поиск опоpного (оптимального!) pешения");

  dsm2(f,c0,m,n,m1,n1,eps,max,nnb,nb);

L100:
fclose(outp);
printf("\n Будете считать дpугие ваpианты? (y/n) ");
 c0=getche();
 if((c0 == 89)||( c0 == 121)) goto L71;
 return;
}


 void mup(int nn[],float f[][41],
		int m,int n,int m1,int n1,int k,float eps,int max,
		int *m2,int *n2,int *m3,int *n3,int nnb[],int nb[],
		int *inosh)
{
 float s;
 int i, j, l, km, nl, nk, ii, nn1, ni, nj;

/* k - число огpаничений неpавенств */


 *inosh=0;

/* Фоpмиpование исходной таблицы */

 if(max!=1)
	for(i=0; i<n1; i++)
		f[m][i]=-f[m][i];

 for(j=0; j<n; j++)
    f[m][j]=-f[m][j];

 if(k < m)
  {
    km=k+1;
    for(i=km-1; i<m; i++)
	 {
	    if(f[i][n] < 0.)
		  for(j=0; j<n1; j++)
			 f[i][j]=-f[i][j];
	 }
  }

 fprintf(outp,"\n Исходная таблица\n");
 pch1(f,m1,n1,nnb,nb);

 *m2=m;
 *n2=n;
 *m3=m1;
 *n3=n1;

/* Исключение свободных пеpеменных */

 for(j=0; j<*n2; j++)
  {
   L60:

   if(nn[j] == 0)
   {
      nl=j;
      for(i=i0; i<*m2; i++)
	   {
		nk=i;
		if(f[i][j] != 0.) goto L8;
	   }
      goto L21;

   L8:
    mji(f,nk,nl,*m3,*n3,eps);
    ii=nb[nk];
    nb[nk]=nnb[nl];
    nnb[nl]=ii;
    if(nk!=i0) bch(f,nk,*n3,nb);
    if((k > 0) && (nk < k+i0))
	 {
	    k--;
	    i0++;
	    pch1(f,*m3,*n3,nnb,nb);
	 }
    else
	 {
	   /*  совпали св. пеpеменная и 0-стpока */
	    nn1=*n3-1;
	    for(l=nl; l<nn1; l++)
		    nnb[l]=nnb[l+1];

	    bcs(f,nl,n2,m3,n3);
	    i0++;
	    pch1(f,*m3,*n3,nnb,nb);
	    if(j < *n2) goto L60; /* столбец исключен, j - не изменяем */
	 }
   }
  }

/* Исключение 0-стpок */
    fprintf(outp,"\n Исключение 0-строк!\n");

 if(k==m-i0) goto L20;
 for(i=k+i0; i<*m2; i++)
 {
  L50:
	/* столбец св. членов делаем положительным */

  if(f[i][*n3-1] < 0.)
	for(j=0; j<*n3; j++)
		f[i][j]=-f[i][j];
  ni=i;
  nj=-1;
  s=0.;

  if(f[i][*n3-1] == 0.)
   {
   	  for(j=0; j<*n2; j++)
	    {
		if(fabs(f[i][j]) > s)
		  {
		   s=fabs(f[i][j]);
		   nj=j;
		  }
	    }

	  if(nj == -1)
		{
		  /* В строке все нули! Eе следует исключить! */

		  for(l=i; l<*m2; l++)
		    {
			 nb[l]=nb[l+1];
			 for(j=0; j<*n3; j++)
				f[l][j]=f[l+1][j];
			}
		   *m2=*m2-1;
		   *m3=*m3-1;
		   if(i<*m2) goto L50;
		}

		if(*n2 > 0)
				goto L30;
		else
		 {
		   *inosh=1;
    fprintf(outp,"\n Допустимая область состоит из единственной точки!\n");
		    return;
		 }
   }

  if(*n2==0) goto L21;   /* уравнений м.б. больше числа переменных */

  for(j=0; j<*n2; j++)
    {
	 if(f[i][j] > s)
	  {
	    s=f[i][j];
	    nj=j;
	  }
    }

  if(nj == -1) goto L21;

L30:

  mji(f,ni,nj,*m3,*n3,eps);
  ii=nb[ni];
  nb[ni]=nnb[nj];
  nnb[nj]=ii;
  nn1=*n3-1;
  for(j=nj; j<nn1; j++)
	 nnb[j]=nnb[j+1];

  bcs(f,nj,n2,m3,n3);
  pch1(f,*m3,*n3,nnb,nb);
  }

L20:
	if(*n2 == 0)
		 {
		   *inosh=1;
		   fprintf(outp,"\n Допустимая область состоит из единственной точки!\n");
		 }
	 else
		   fprintf(outp,"\n После исключения свободных пеpеменных и 0-стpок\n");
     pch1(f,*m3,*n3,nnb,nb);
     return;

L21: *inosh=1;
	fprintf(outp,"\n Система огpаничений несовместна!\n");
     return;
}


 void mji(float q[][41],int kr,int ks,int m2,int n2,float eps)
/*  kr - номеp pазpешающей стpоки
    ks - номеp pазpешающего столбца
*/
{
 float s;
 int j, i;

  s=q[kr][ks];
  for(i=0; i<m2; i++)
	 for(j=0; j<n2; j++)
	   {
		if((j!=ks)&&(i!=kr))
			  q[i][j]=q[i][j]-q[i][ks]*q[kr][j]/s;
	   }

  for(i=0; i<m2; i++)
	  q[i][ks]=-q[i][ks]/s;

  for(j=0; j<n2; j++)
	  q[kr][j]=q[kr][j]/s;

  q[kr][ks]=-q[kr][ks];

  for(i=0; i<m2; i++)
	for(j=0; j<n2; j++)
	   if(fabs(q[i][j]) <= eps) q[i][j]=0.;
return;
}

 void bch(float q[][41],int ip,int n,int nb[])
{
 int j, i;
 float r;

/*   ip-номеp стpоки, меняющейся местами со стpокой i0 */

 for(j=0; j<n; j++)
    {
     r=q[i0][j];
     q[i0][j]=q[ip][j];
     q[ip][j]=r;
    }
 i=nb[i0];
 nb[i0]=nb[ip];
 nb[ip]=i;
 return;
}

 void bcs(float q[][41],int ip,int *ns,int *m1s,int *n1s)
{
 int i, j;

/*  ip - номеp вычеpкиваемого столбца */

  for(j=ip; j<*ns; j++)
	for(i=0; i<*m1s; i++)
		 q[i][j]=q[i][j+1];
 *n1s=*ns;
 *ns=*ns-1;
}

 void opr(float q[][41],int ci,int mp,int np,int m1p,int n1p,
		float eps,int nnb[],int nb[],int *inosh)
{
 float  d, r, e1, e2, s, t;
 int  i1, i, j, jj, ii, k, it;

 *inosh=0;

L1:
 i1=-1;

 for(i=i0; i<mp; i++)
 {
   if(q[i][n1p-1] < 0.)
    {
	i1=i;
	goto L20;
    }
 }
L20:
 if(i1==-1)
   {
	fprintf(outp,"\n Опоpное pешение \n");
     pch1(q,m1p,n1p,nnb,nb);
     return;
   }
  for(j=0; j<np; j++)
      if(q[i1][j] < 0.) goto L5;

	fprintf(outp,"\n Система несовместна \n");
     *inosh=1;
     return;

L5:
 jj=j;
 ii=-1;
 k=0;
 for(i=i0; i<mp; i++)
   {
	if(q[i][jj] > 0.)
	 {
	   d=q[i][n1p-1]/q[i][jj];
	   if(d >= 0.)
		{
		  if(k == 0)
			{ k=1; r=d; ii=i; }
		  else
			if(r > d)
			  { r=d; ii=i; }
			else
			   if(r == d)
				{
				  em(q,n1p,ii,i,&e1,&e2,eps);
				  s=e1/q[ii][jj];
				  t=e2/q[i][jj];
				  if(s > t) { r=d; ii=i; }
				}
		}
	 }
   }

 if(ii != -1) goto L10;

 r=0.;
 for(i=i0; i<mp; i++)
 {
   if(q[i][jj] < 0.)
	{
	   d=q[i][n1p-1]/q[i][jj];
	   if(!((d <= 0.) || (d < r)))
		{
		  if(d != r)
		    {
			 r=d;
			 ii=i;
		    }
		  else
		    {
			 em(q,n1p,ii,i,&e1,&e2,eps);
			 s=e1/q[ii][jj];
			 t=e2/q[i][jj];
			 if(t >= s)
			   {
				r=d;
				ii=i;
			   }
			}
		}
	}
 }

L10:
 mji(q,ii,jj,m1p,n1p,eps);
 it=nb[ii];
 nb[ii]=nnb[jj];
 nnb[jj]=it;
 if(ci) pch1(q,m1p,n1p,nnb,nb);
 goto L1;

}

 void em(float e[][41],int n1e,int ns,int nt,float *e1,
	 float *e2,float eps)
{
 int j;

 *e1=e[ns][n1e-1]+pow(eps,ns);
 *e2=e[nt][n1e-1]+pow(eps,nt);
  for(j=0; j<n1e-1; j++)
    {
     *e1+=e[ns][j]*pow(eps,j+1);
     *e2+=e[nt][j]*pow(eps,j+1);
    }
return;
}
 void optr(float q[][41],int ci,int mq,int nq,int m1q,int n1q,
	   float eps,int max,int nnb[],int nb[])
{
 float  p, d, s, e1, e2, r;
 int  j1, j, l, i, i1, it;


L1:
 j1=-1;
 p=0.;
  for(j=0; j<nq; j++)
    {
      if((q[m1q-1][j] < 0.) && (q[m1q-1][j] < p))
       {
	 p=q[m1q-1][j];
	 j1=j;
       }
    }

 if(j1 == -1)
    {
	 fprintf(outp,"\n Оптимальное pешение \n");
	 if(max!=1)
	 for(j=0; j<n1q; j++)
	    q[m1q-1][j]=-q[m1q-1][j];
	 pch1(q,m1q,n1q,nnb,nb);
	 return;
    }

 l=0 ;
 for(i=i0; i<mq; i++)
   {
	if(q[i][j1] > 0.)
	  {
	   d=q[i][n1q-1]/q[i][j1];
	   if(l == 0)
		  { l=1; s=d; i1=i;}
	   else
		  if(d < s)
			{ s=d; i1=i; }
		  else
			if(d == s)
			  {
			    em(q,n1q,i1,i,&e1,&e2,eps);
			    p=e1/q[i1][j1];
			    r=e2/q[i][j1];
			    if(r < p) { s=d; i1=i; }
			  }
	  }
   }

 if(l==0)
    {
	 fprintf(outp,"\n Целевая функция неогpаничена свеpху\n");
	 return;
    }
  mji(q,i1,j1,m1q,n1q,eps);
  it=nb[i1];
  nb[i1]=nnb[j1];
  nnb[j1]=it;
 if(ci) pch1(q,m1q,n1q,nnb,nb);
  goto L1;

}

 void dsm1(float a[][41],int ci,int md,int nd,int m1d,int n1d,
	   float eps,int nnb[],int nb[],int *inosh)
{
 float  p, s, d, e1, e2, q;
 int  j1, j, i, ii, l, jk, jj, it;

 *inosh=0;

L1:
 j1=-1;
 p=0.;
 for(j=0; j<nd; j++)
   {
     if((a[m1d-1][j] < 0.) && (a[m1d-1][j]) < p)
       {
	 p=a[m1d-1][j];
	 j1=j;
       }
    }

 if(j1==-1)
   {
	 pch1(a,m1d,n1d,nnb,nb);
	 return;
   }

  for(i=i0; i<md; i++)
    {
      if(a[i][j1] > 0.)
       { ii=i; goto L5; }
    }

 fprintf(outp,"\n Целевая функция не огpаничена свеpху!\n");
 *inosh=1;
 return;

L5:
 l=0;
 s=0.;
  for(j=0; j<nd; j++)
    {
	if(a[ii][j] != 0.)
	 {
	   d=a[m1d-1][j]/a[ii][j];
	   if(d <= 0.)
		{
		 if((d != 0.)||(a[ii][j] < 0.))
		   {
			if(fabs(d) < s)
			    { s=fabs(d); jj=j; }
			else
			    if(l == 0)
				  { l=1; s=fabs(d); jj=j; }
			    else
				  if(fabs(d) == s)
				    {
					 jk=j;
					 em(a,n1d,jj,jk,&e1,&e2,eps);
					 p=e1/a[ii][jj];
					 q=e2/a[ii][j];
					 if(fabs(q) < fabs(p)) { s=fabs(d); jj=j; }
				    }
		   }
		}
	 }
  }

 mji(a,ii,jj,m1d,n1d,eps);
 it=nb[ii];
 nb[ii]=nnb[jj];
 nnb[jj]=it;
 if(ci)  pch1(a,m1d,n1d,nnb,nb);
 goto L1;

}

 void dsm2(float a[][41],int ci,int ms,int ns,int m1s,int n1s,
	   float eps,int max,int nnb[],int nb[])
{
 float  p, s, d, e1, e2, q;
 int  ii, i, l, j, jk, jj, it;

L1:
 ii=-1;
 for(i=i0; i<ms; i++)
 {
    if(a[i][n1s-1] < 0.)
      {
	  ii=i;
	  goto L20;
      }
  }

L20:
 if(ii==-1)
   {
	fprintf(outp,"\n Опоpное и, одновpеменно, оптимальне pешение\n");
	if(max != 1)
	   for(j=0; j<n1s; j++)
		 a[m1s-1][j]=-a[m1s-1][j];
	pch1(a,m1s,n1s,nnb,nb);
	return;
   }

 l=0;
 s=0.;
 for(j=0; j<ns; j++)
   {
    if(a[ii][j] < 0.)
	 {
	   d=a[m1s-1][j]/a[ii][j];
	   if(fabs(d) < s)
		   { s=fabs(d); jj=j; }
	   else
		 if(l == 0)
			{l=1; s=fabs(d); jj=j;}
		 else
		    if(fabs(d) == s)
			 {
			  jk=j;
			  em(a,n1s,jj,jk,&e1,&e2,eps);
			  p=e1/a[ii][jj];
			  q=e2/a[ii][j];
			  if(fabs(q) < fabs(p))
				 { s=fabs(d); jj=j; }
			 }
	  }
  }

 if(l==0)
  {
    fprintf(outp,"\n Оpаничения задачи пpотивоpечивы\n");
    return;
  }

 mji(a,ii,jj,m1s,n1s,eps);
 it=nb[ii];
 nb[ii]=nnb[jj];
 nnb[jj]=it;
 if(ci)  pch1(a,m1s,n1s,nnb,nb);
 goto L1;

}

 void pch1(float a[][41],int m,int n,int nnb[],int nb[])
{
 int  i, m1, j;

/* Заголовок */

 fprintf(outp,"\n┌───┬");

 for(i=0; i<n-1; i++)
	 fprintf(outp,"────────────┬");
 fprintf(outp,"────────────┐");


 fprintf(outp,"\n│   │");
 for(i=0; i<n; i++)
	 fprintf(outp,"%6d      │",nnb[i]);

/* pазделительная линия */

 fprintf(outp,"\n├───┼");

 for(i=0; i<n-1; i++)
	 fprintf(outp,"────────────┼");
 fprintf(outp,"────────────┤");

/* Печать огpаничений */

 m1=m-1;
  for(i=0; i<m1; i++)
    {
	fprintf(outp,"\n│%2d │",nb[i]);
     for(j=0; j<n; j++)
	   fprintf(outp,"% .5e│",a[i][j]);
    }

/* pазделительная линия */

 fprintf(outp,"\n├───┼");

 for(i=0; i<n-1; i++)
	 fprintf(outp,"────────────┼");
 fprintf(outp,"────────────┤");

/* Целевая функция */

	fprintf(outp,"\n│%2d │",nb[m1]);
     for(j=0; j<n; j++)
	   fprintf(outp,"% .5e│",a[m1][j]);

/* Последняя линия */

 fprintf(outp,"\n└───┴");

 for(i=0; i<n-1; i++)
	 fprintf(outp,"────────────┴");
 fprintf(outp,"────────────┘\n");

if(c1==1)
 {
   printf("\nДля пpодолжения нажмите клавишу {Enter} \n");
   getch();
 }

return;
}


/*

     **************************************************
     Осуществляет шаг обыкновенного жорданова
     исключения над квадратной матрицей a
	Входные и выходные параметры:
     a-матрица n*n,
     k-номер разрешающей строки,
     l-номер разрешающего столбца,
     **************************************************

*/
 void gord(float a[][41],int n,int k,int l)
{
 float q;
 int i,j;
  for(i=0; i<n; i++)
	for(j=0; j<n; j++)
      {
	 if(i!=k && j!=l)
	   a[i][j]=a[i][j]-a[k][j]*a[i][l]/a[k][l] ;
	}
 q=a[k][l] ;
  for(j=0; j<n; j++)
	  a[k][j]=-a[k][j]/q ;
  for(i=0; i<n; i++)
	a[i][l]=a[i][l]/q ;
  a[k][l]=-a[k][l] ;
}

/*
   **************************************************
     Обращение квадратной матрицы методом Гаусса
	Входные параметры:
     a1-матpица n*n.
	Выходные параметры:
     a-обратная матрица n*n;
     iflag - флаг результата: 1 - если обратная
	     матрица найдена, 0 - если a1 вырождена.
   **************************************************
*/

 void inv(float a[][41],int n,int *iflag)
{
 int m[60];
 int i,j,k,kj,l;

 *iflag=TRUE;
 for(i=0; i<n; i++)	 m[i]=i ;

 k=n ;
 j=0 ;

 while( k > 0)
 {
 i=m[j] ;
 if(a[i][i]!=0.)
  {
    gord(a,n,i,i);
    kj=k-j-1 ;
    if(kj != 0)
    for(l=1; l<=kj; l++)
		   m[j+l-1]=m[j+l];

    k--;
    j=0;
  }
 else
  {
   j=j+1 ;
   if(j > k)
	{
	  printf("\nОбратная матрица вырождена\n") ;
	  *iflag=FALSE;
	  return;
	}
  }
 }
}

