#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "f.h"
#define LEN 5000
#define eps 1e-40
 int LU(double *a,int n)
{
	int i,j,k,t;
	double s;
	double *p,*q;
	for(i=0;i<n;i++)
	{
		p=a+i*n;
		for(j=i;j<n;j++)
		{
				s=0;
				q=a+i;
				t=0;
				for(k=0;k<i;k++)
				{
					s+=p[k]*q[t];
					t+=n;
				}
				p[i]-=s;
				p+=n;


		}
		p=a+i*n;
		q=a+i+1;
		for(j=i+1;j<n;j++)
		{
				t=0;
				s=0;
				for(k=0;k<i;k++)
				{
					s+=p[k]*q[t];
					t+=n;
				}
				p[j]-=s;
				if(fabs(p[i])< eps)	return 1;
				p[j]/=p[i];
				q++;
		}
	}
	return 0;
}

void obr_LU(double*a,int n)
{
	double e[LEN];
	double k;
	double *p,*q;
	int i,j,r;
	p=a+(n-1)*n;
	for(i=n-1;i>=0;i--) //U^(-1)
	{
		e[i]=1;
		q=a+n*(i+1);
		for(j=i+1;j<n;j++)
		{
			k=p[j];
			e[j]-=k;
			for(r=j+1;r<n;r++)
			{
				e[r]-=k*q[r];
			}
			q+=n;
		}
		for(j=i+1;j<n;j++)
		{
			p[j]=e[j];
			e[j]=0;
		}
		e[i]=0;
		p-=n;
	}
	p=a;
	for(i=0;i<n;i++) //L^(-1)
	{
		e[i]=1;
		q=a;
		for(j=0;j<i;j++)
		{
			k=p[j];
			for(r=0;r<=j;r++)
			{
				e[r]-=k*q[r];
			}
			q+=n;
		}
		for(j=0;j<=i;j++)
		{
			p[j]=e[j]/p[i];
			e[j]=0;
		}
		e[i]=0;
		p+=n;
	}
}

void L_r(double *a,int n)
{
	int i,j;
	double d;
	double *p=a,*b=a+n*n-1;
	for(i=0;i<n/2;i++)
	{
		for(j=0;j<=i;j++)
		{
			d=p[j];
			printf("p[j]=%lf\n",p[j]);
			p[j]=b[-j];
			b[-j]=d;
		}
		p+=n;
		b-=n;
	}
}

void mult_matr(double *a,double *d,int n)
{
	int i,j,k,t;
	double s;
	double *p=a,*q=a+n-1,*l=d;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{

			if(i>=j)
			{
				s=0;
				s=p[j];
				t=-j*n+i;
				for(k=i+1;k<n;k++)
				{
					s+=p[k]*q[t];
					t--;
				}
				l[j]=s;
			}
			else
			{
				s=0;
				t=0;
				for(k=j;k<n;k++)
				{
					s+=p[k]*q[t];
					t++;
				}
				l[j]=s;
			}
		q--;
		}
		p+=n;
		l+=n;
	}
}

double func(double *a,double *d,int n)
{
	int i,j,k;
	double sum,max=-1,s;
	for(i=0;i<n;i++)
	{
		sum=0;
		for(j=0;j<n;j++)
		{
			s=0;
			for(k=0;k<n;k++)
			{
				s+=a[i*n+k]*d[k*n+j];
			}
			if(i==j) s--;
			sum+=fabs(s);
		}
		if(sum>max) max=sum;
	}
	return max;
}