#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>

int GJ(int,double **);
double **TwoArrayAlloc(int,int);
void TwoArrayFree(double **);

void main()
{
	int i,j,n;
	double **a;
	n=4;
	a=TwoArrayAlloc(n,n);

	a[0][0]=5;	a[0][1]=7;	a[0][2]=6;	a[0][3]=5;
	a[1][0]=7;	a[1][1]=10;	a[1][2]=8;	a[1][3]=7;
	a[2][0]=6;	a[2][1]=8;	a[2][2]=10;	a[2][3]=9;
	a[3][0]=5;	a[3][1]=7;	a[3][2]=9;	a[3][3]=10;
	if(!GJ(n,a))
	{
		printf("矩陣求逆失敗\n");
		exit(1);
	}
	printf("該矩陣的逆為：\n");
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
			printf("%.2f\t",a[i][j]);
		printf("\n");
	}
}

int GJ(int n,double **a)
{
	int i,j,k;
	double p,q,*h;
	h=(double *)calloc(n,sizeof(double));
	if(h == NULL)
	{
		printf("內(nèi)存分配失敗\n");
		exit(1);
	}
	for(k=n;k>=1;k--)
	{
		p=a[0][0];
		if(p<=0)
		{
			free(h);
			return (0);
		}
		for(i=2;i<=n;i++)
		{
			q=a[i-1][0];
			if(i>k)
				h[i-1]=q/p;
			else
				h[i-1]=-q/p;
			for(j=2;j<=i;j++)
				a[i-2][j-2]=a[i-1][j-1]+q*h[j-1];
		}
		a[n-1][n-1]=1/p;
		for(i=2;i<=n;i++)
			a[n-1][i-2]=h[i-1];
	}
	free(h);
	return(1);
}

double **TwoArrayAlloc(int r,int c)
{
	double *x,**y;
	int n;
	x=(double *)calloc(r*c,sizeof(double));
	y=(double **)calloc(r,sizeof(double*));
	for(n=0;n<=r-1;++n)
		y[n]=&x[c*n];
	return (y);
}

void TwoArrayFree(double **x)
{
	free(x[0]);
	free(x);
}