#!/usr/local/bin/tcc -run
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define MAX_ENT 100 /* Max entry size */


int solve_leq(double **a,double *b,int n);

double *zero_vect(int n)
{
  return (double *)calloc(n,sizeof(double));
}

double **zero_mat(int n)
{
  int i,j;
  double **a;

  a = (double **)malloc(n*sizeof(double *));
  for ( i = 0; i < n; i++ )
    a[i] = zero_vect(n);
  return a;
}


double **random_mat(int n)
{
  int i,j;
  double **a;

  a = zero_mat(n);
  for ( i = 0; i < n; i++ )
    for ( j = 0; j < n; j++ )
      a[i][j] = (double)(random()%MAX_ENT);
  return a;
}

double *random_vect(int n)
{
  int i;
  double *v;

  v = zero_vect(n);
  for ( i = 0; i < n; i++ )
    v[i] = (double)(random()%MAX_ENT);
  return v;
}

double **copy_mat(double **a,int n)
{
  int i,j;
  double **c;

  c = zero_mat(n);
  for ( i = 0; i < n; i++ )
    for ( j = 0; j < n; j++ )
      c[i][j] = a[i][j];
  return c;
}

double *copy_vect(double *a,int n)
{
  int i;
  double *c;

  c = zero_vect(n);
  for ( i = 0; i < n; i++ )
    c[i] = a[i];
  return c;
}

void print_mat(double **a,int n)
{
  int i,j; 
  for ( i = 0; i < n; i++ ) {
    for ( j = 0; j < n; j++ ) printf("%lf ",a[i][j]);
    printf("\n");
  }
}

void print_vect(double *a,int n)
{
  int i; 
  for ( i = 0; i < n; i++ ) printf("%lf ",a[i]);
  printf("\n");
}


double check_sol(double **a,double *x,double *b,int n)
{
  double t,s;
  int i,j;

  s = 0;
  for ( i = 0; i < n; i++ ) {
    t = -b[i];
    for ( j = 0; j < n; j++ )
     t += a[i][j]*x[j];
    s += t*t;
  }
  return s;
}

int check_isnan(double *a,int n)
{
  int i;

  for ( i = 0; i < n; i++ )
    if ( isnan(a[i]) ) {
      fprintf(stderr,"x[%d] is NaN\n",i);
        return 1;
    }
  return 0;
}


/* Main Function */
int main(int argc,char **argv)
{
	int n,n1,ret;
	double s;
	double **a,**a1;
	double *b,*x;
	clock_t t0;

	srandom((int)time(0));
	//if ( argc == 2 ) {
		// n = atoi(argv[1]);  /* Input size */
		printf("行�Eサイズを�E力してください\n");
		scanf("%d", &n); 
		printf("行�Eサイズ n = %d\n", n);
		
		a = random_mat(n);  /* Random matrix */
		b = random_vect(n); /* Random vector */
	//} else {
	//	printf("Please input the size\n");
	//	exit(1); 
	//}
	
	
	// printf("A=\n"); print_mat(a,n);
	// printf("b=\n"); print_vect(b,n);
	a1 = copy_mat(a,n);
	x = copy_vect(b,n);
	t0 = clock();
	ret = solve_leq(a1,x,n);
	t0 = clock()-t0;
	printf("%lfsec\n",t0/(double)CLOCKS_PER_SEC);
	
	if ( ret == 0 ) {
		printf("A is singular\n");
	} else {
		if ( check_isnan(x,n) ) {
			fprintf(stderr,"solution contains NaN\n");
			exit(1);
		}
		printf("x=\n"); print_vect(x,n);
		printf("||Ax-b||^2=%lf\n",s = check_sol(a,x,b,n));
	}
}

int solve_leq(double **a,double *b,int n)
{
	/* Write your code here  */
	int i,j,k;
	double m,temp;
	for (j = 0; j < n; j++) {
		i = -1;
		for (k = j; k < n; k++) {
			if ((i == -1) && (a[k][j] != 0)) {
				i = k;
			}
		}
		if (i == -1) {
			return 0;
		}
		if (i != j) {
			for (k = 1; k < n; k++) {
				temp = a[i][k];
				a[i][k] = a[j][k];
				a[j][k] = temp;
			}
			temp = b[i];
			b[i] = b[j];
			b[j] = temp;
		}
		for (i = j + 1; i < n; i++) {
			m = a[i][j] / a[j][j];
			for (k = j; k < n; k++) {
				a[i][k] = a[i][k] - m * a[j][k];
			}
			b[i] = b[i] - m * b[j];
		}
	}
	for (i = n - 1; i >= 0; i--) {
		for (j = i + 1; n - 1; j++) {
			b[i] = b[i] - a[i][j] * b[j];
		}
		b[i] = b[i] / a[i][j];
	}
	return 1; /* success */
}