UOJ的这道题

Test#32错误地输出了Infeasible

Test#47与答案相差0.000015

救救蒟蒻吧

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
namespace L
{
	void swap(int &a,int &b) {a^=b,b^=a,a^=b;}
	const int M_Size=1e4+10,N_Size=1e3+10;
	const double eps=1e-14,INF=1e14;
	bool deq(double a,double b){return fabs(a-b)<=eps;}
	int n,m,B[M_Size],N[N_Size];
	double a[M_Size][N_Size],b[M_Size],c[N_Size],v,x[N_Size];
//	void print()
//	{
//		printf("z=%g",v);
//		for(int i=1;i<=n;i++)
//		if(c[i])
//		{
//			if(deq(c[i],-1)) printf("-x_%d",N[i]);
//			else if(deq(c[i],+1)) printf("+x_%d",N[i]);
//			else if(c[i]>0) printf("+%gx_%d",c[i],N[i]);
//			else printf("-%gx_%d",-c[i],N[i]);
//		}
//		printf("\n");
//		for(int i=1;i<=m;i++)
//		{
//			printf("x_%d=%g",B[i],b[i]);
//			for(int j=1;j<=n;j++)
//			if(!deq(a[i][j],0))
//			{
//				if(deq(a[i][j],1)) printf("-x_%d",N[j]);
//				else if(deq(a[i][j],-1)) printf("+x_%d",N[j]);
//				else if(a[i][j]>0) printf("-%gx_%d",a[i][j],N[j]);
//				else printf("+%gx_%d",-a[i][j],N[j]);
//			}
//			printf("\n");
////			printf("%g\n",a[1][1],-a[1][1]);
//		}
//	}
	void pivot(int l,int e)
	{
		swap(B[l],N[e]);
		b[l]/=a[l][e];
		for(int j=1;j<=n;j++)
		if(j!=e) a[l][j]/=a[l][e];
		a[l][e]=1/a[l][e];
//		printf("PIVOTING:\n");
//		print();
		for(int i=1;i<=m;i++)
		if(i!=l)
		{
			b[i]-=a[i][e]*b[l];
			for(int j=1;j<=n;j++)
			if(j!=e) a[i][j]-=a[i][e]*a[l][j];
			a[i][e]=-a[i][e]*a[l][e];
//			printf("when %d end: %g\n",i,a[1][1]);
		}
//		printf("ALSO PIVOTING(%d,%g):\n",e,a[1][1]);print();
		v+=c[e]*b[l];
		for(int j=1;j<=n;j++)
		if(j!=e) c[j]-=c[e]*a[l][j];
		c[e]=-c[e]*a[l][e];
	}
	double simplex()
	{
//		printf("SIMPLEX:\n");
//		print();
		while(1)
		{
//			printf("!");
			int e=1,l=0;
			for(;e<=n;e++)
			if(c[e]>eps) break;
			if(e==n+1) break;
			double delta=INF;
			for(int i=1;i<=m;i++)
			if(a[i][e]>eps)
			{
				if(b[i]/a[i][e]<delta)
				{
					delta=b[i]/a[i][e];
					l=i;
				}
			}
			if(l==0) return INF;
//			print();
//			printf("pivot %d(%d) %d(%d)\n",l,B[l],e,N[e]);
			pivot(l,e);
//			print();
		}
		memset(x,0,sizeof x);
		for(int i=1;i<=m;i++)
		if(B[i]<=n) x[B[i]]=b[i];
		return v;
	}
	double old_c[N_Size];
//	bool bz[N];
	bool initialize()
	{
//		printf("INITIALIZE:\n");print();
		int k=1;
		for(int i=2;i<=m;i++)
		if(b[i]<b[k]) k=i;
		if(b[k]>=-eps) return true;
		memcpy(old_c,c,sizeof c);
		memset(c,0,sizeof c);
		c[++n]=-1;N[n]=n+m+1;v=0;
		for(int i=1;i<=m;i++) a[i][n]=-1;
//		printf("BEFORE PIVOT:\n");print();
		pivot(k,n);
//		printf("AFTER PIVOT:\n");print();
		if(!deq(simplex(),0)) return false;
		for(int l=1,e;l<=m;l++)
		if(B[l]==n+m+1)
		{
			for(e=1;e<=n;e++)
			if(!deq(a[l][e],0)) break;
			if(e==n+1) printf("Oop! Seems I forget this case...\n");
			pivot(l,e);
			break;
		}
//		print();
		for(int j=1;j<=n;j++)
		if(N[j]==n+m+1)
		{
//			printf("%d\n",j);
			swap(N[j],N[n]);
			for(int i=1;i<=m;i++)
			std::swap(a[i][j],a[i][n]);
			break;
		}
		n--;
		memset(c,0,sizeof c);
		for(int j=1;j<=n;j++)
		if(N[j]<=n) c[j]=old_c[N[j]];
		for(int i=1;i<=m;i++)
		if(B[i]<=n)
		{
			v+=old_c[B[i]]*b[i];
			for(int j=1;j<=n;j++)
			c[j]-=old_c[B[i]]*a[i][j];
		}
		return true;
	}
}
int main()
{
	int t;scanf("%d%d%d",&L::n,&L::m,&t);
	for(int i=1;i<=L::n;i++) scanf("%lf",&L::c[i]),L::N[i]=i;
	for(int i=1;i<=L::m;i++)
	{
		for(int j=1;j<=L::n;j++)
		scanf("%lf",&L::a[i][j]);
		scanf("%lf",&L::b[i]);
		L::B[i]=L::n+i;
	}
	if(!L::initialize()) printf("Infeasible");
	else if(L::simplex()==L::INF) printf("Unbounded");
	else
	{
		printf("%lf",L::v);
		if(t)
		{
			printf("\n");
			for(int i=1;i<=L::n;i++)
			printf("%lf ",L::x[i]);
		}
	}
}