/* sqmatrix.h
K.Shiota
2007.7 version
*/

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

#define RANDOMIZE()	srand(time(NULL))
#define RANDOM(x)	(rand()%(x))		// 0..x-1 間の整数乱数
#define REALRANDOM(x)	(((x)*rand())/RAND_MAX)	// 0..x 間の実数乱数

/* 行列の最大サイズを設定 */
#define MAXSIZE 128

/* 実行対象の行列サイズを設定
   （デフォルト値は 3, 主プログラムで指定せよ ） */
int SIZE=3; 

/* MATRIX, VECTOR, ROWVECTOR という型を定義 */
typedef struct { double ent[MAXSIZE]; } VECTOR;
typedef struct { double ent[MAXSIZE]; } ROWVECTOR;
typedef struct { double ent[MAXSIZE][MAXSIZE]; } MATRIX;

/********************************************************

VECTOR v;
と宣言すると v.ent[i] が列ベクトル v の 第i成分を表す。 

ROWVECTOR v;
と宣言すると v.ent[i] が行ベクトル v の 第i成分を表す。 

MATRIX a;
と宣言すると a.ent[i][j] が行列 a の (i,j)-成分を表す。 

*********************************************************/

/* 列ベクトルを表示する関数 */
void vecprint(VECTOR v);

/* 行ベクトルを表示する関数 */
void rowvecprint(ROWVECTOR v);

/* 行列を表示する関数 */
void matprint(MATRIX a);

/* 列ベクトル v のスカラー倍 w = c v を返す関数 */
VECTOR vecscmul(double c, VECTOR v);

/* 行ベクトル v のスカラー倍 w = c v を返す関数 */
ROWVECTOR rowvecscmul(double c, ROWVECTOR v);

/* 行列 A のスカラー倍 B = c A を返す関数 */
MATRIX matscmul(double c, MATRIX a);

/* ふたつの列ベクトル u, v の和 w = u + v を返す関数 */
VECTOR vecadd(VECTOR u, VECTOR v);

/* ふたつの列ベクトル u, v の差 w = u - v を返す関数 */
VECTOR vecsub(VECTOR u, VECTOR v);

/* ふたつの行ベクトル u, v の和 w = u + v を返す関数 */
ROWVECTOR rowvecadd(ROWVECTOR u, ROWVECTOR v);

/* ふたつの行ベクトル u, v の差 w = u - v を返す関数 */
ROWVECTOR rowvecsub(ROWVECTOR u, ROWVECTOR v);

/* ふたつの行列 A, B の和 C = A + B を返す関数 */
MATRIX matadd(MATRIX a, MATRIX b);

/* ふたつの行列 A, B の差 C = A - B を返す関数 */
MATRIX matsub(MATRIX a, MATRIX b);

/* 行ベクトル u と列ベクトルv の積 w = u v を返す関数 */
double rowvecvecmul(ROWVECTOR u, VECTOR v);

/* 列ベクトル u と行ベクトル v の積 A = u v を返す関数 */
MATRIX vecrowvecmul(VECTOR u, ROWVECTOR v);

/* 行列 A と列ベクトル v の積 w = A v を返す関数 */
VECTOR matvecmul(MATRIX a, VECTOR v);

/* 行ベクトル v と行列 A の積 w = v A を返す関数 */
ROWVECTOR rowvecmatmul(ROWVECTOR b, MATRIX a);

/* 行列 A, B の積 C = A B を返す関数 */
MATRIX matmul(MATRIX a, MATRIX b);

/* ゼロベクトル（列ベクトル）を返す関数 */
VECTOR zerovec();

/* ゼロベクトル（行ベクトル）を返す関数 */
ROWVECTOR zerorowvec();

/* ゼロ行列を返す関数 */
MATRIX zeromat();

/* 単位行列を返す関数 */
MATRIX unitmat();

/* 第 i,j 座標の角度 x の回転を表す行列を返す関数 */
MATRIX rotmat(int i, int j, double x);

/* 列ベクトル v の転置ベクトルを返す関数 */
ROWVECTOR transvec(VECTOR v);

/* 行ベクトル v の転置ベクトルを返す関数 */
VECTOR transrowvec(ROWVECTOR v);

/* 行列 A の転置行列を返す関数 */
MATRIX transmat(MATRIX a);

/* ふたつのベクトル u, v の標準内積を返す関数 */
double innerprod(VECTOR u, VECTOR v);

/* 行ベクトル v の長さ（ユークリッドノルム）を返す関数 */
double vecnorm(VECTOR v);

/* 列ベクトル v の長さ（ユークリッドノルム）を返す関数 */
double rowvecnorm(ROWVECTOR v);

/* ランダムな列ベクトルを返す関数（成分の最大絶対値 < x） */
VECTOR randomvec(double x);

/* ランダムな整数成分列ベクトルを返す関数（成分の最大絶対値 < x） */
VECTOR randomintvec(int x);

/* ランダムな行ベクトルを返す関数（成分の最大絶対値 < x） */
ROWVECTOR randomrowvec(double x);

/* ランダムな整数成分行ベクトルを返す関数（成分の最大絶対値 < x） */
ROWVECTOR randomintrowvec(int x);

/* ランダムな行列を返す関数（成分の最大絶対値 < x）*/
MATRIX randommat(double x);

/* ランダムな整数成分行列を返す関数（成分の最大絶対値 < x）*/
MATRIX randomintmat(int x);

/********** 以下関数定義 **********/

void vecprint(VECTOR v)
{
	int i;
	for(i=0;i<SIZE;i++) printf("%12.7lf\n",v.ent[i]);
}

void rowvecprint(ROWVECTOR v)
{
	int i;
	for(i=0;i<SIZE;i++) printf("%12.7lf  ",v.ent[i]);
	printf("\n");
}

void matprint(MATRIX a)
{
	int i,j;
	for(i=0;i<SIZE;i++){
		for(j=0;j<SIZE;j++) printf("%12.7lf",a.ent[i][j]);
		printf("\n");
	}
}

VECTOR vecscmul(double c, VECTOR v)
{
	VECTOR w;
	int i;
	for(i=0;i<SIZE;i++) w.ent[i]=c*v.ent[i];
	return w;
}

ROWVECTOR rowvecscmul(double c, ROWVECTOR v)
{
	ROWVECTOR w;
	int i;
	for(i=0;i<SIZE;i++) w.ent[i]=c*v.ent[i];
	return w;
}

MATRIX matscmul(double c, MATRIX a)
{
	MATRIX b;
	int i,j;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			b.ent[i][j]=c*a.ent[i][j];
	return b;
}

VECTOR vecadd(VECTOR u, VECTOR v)
{
	VECTOR w;
	int i;
	for(i=0;i<SIZE;i++) w.ent[i]=u.ent[i]+v.ent[i];
	return w;
}

VECTOR vecsub(VECTOR u, VECTOR v)
{
	VECTOR w;
	int i;
	for(i=0;i<SIZE;i++) w.ent[i]=u.ent[i]-v.ent[i];
	return w;
}

ROWVECTOR rowvecadd(ROWVECTOR u, ROWVECTOR v)
{
	ROWVECTOR w;
	int i;
	for(i=0;i<SIZE;i++) w.ent[i]=u.ent[i]+v.ent[i];
	return w;
}

ROWVECTOR rowvecsub(ROWVECTOR u, ROWVECTOR v)
{
	ROWVECTOR w;
	int i;
	for(i=0;i<SIZE;i++) w.ent[i]=u.ent[i]-v.ent[i];
	return w;
}

MATRIX matadd(MATRIX a, MATRIX b)
{
	MATRIX c;
	int i,j;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			c.ent[i][j]=a.ent[i][j]+b.ent[i][j];
	return c;
}

MATRIX matsub(MATRIX a, MATRIX b)
{
	MATRIX c;
	int i,j;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			c.ent[i][j]=a.ent[i][j]-b.ent[i][j];
	return c;
}

double rowvecvecmul(ROWVECTOR u, VECTOR v)
{
	double x=0.0;
	int k;
	for(k=0;k<SIZE;k++) x+=u.ent[k]*v.ent[k];
	return x;
}

MATRIX vecrowvecmul(VECTOR u, ROWVECTOR v)
{
	MATRIX a;
	int i,j;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			a.ent[i][j]=u.ent[i]*v.ent[j];
	return a;
}

VECTOR matvecmul(MATRIX a, VECTOR v)
{
	VECTOR w;
	int i,k;
	double x;
	for(i=0;i<SIZE;i++){
		x=0.0;
		for(k=0;k<SIZE;k++) x+=a.ent[i][k]*v.ent[k];
		w.ent[i]=x;
	}
	return w;
}

ROWVECTOR rowmatvecmul(ROWVECTOR v, MATRIX a)
{
	ROWVECTOR w;
	int j,k;
	double x;
	for(j=0;j<SIZE;j++){
		x=0.0;
		for(k=0;k<SIZE;k++) x+=v.ent[k]*a.ent[k][j];
		w.ent[j]=x;
	}
	return w;
}

MATRIX matmul(MATRIX a, MATRIX b)
{
	MATRIX c;
	int i,j,k;
	double x;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++){
			x=0.0;
			for(k=0;k<SIZE;k++) x+=a.ent[i][k]*b.ent[k][j];
			c.ent[i][j]=x;
		}
	return c;
}

VECTOR zerovec()
{
	VECTOR v;
	int i;
	for(i=0;i<SIZE;i++) v.ent[i]=0;
	return v;
}

ROWVECTOR zerorowvec()
{
	ROWVECTOR v;
	int i;
	for(i=0;i<SIZE;i++) v.ent[i]=0;
	return v;
}

MATRIX zeromat()
{
	MATRIX a;
	int i,j;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			a.ent[i][j]=0.0;
	return a;
}

MATRIX unitmat()
{
	int i;
	MATRIX a=zeromat();
	for(i=0;i<SIZE;i++) a.ent[i][i]=1.0;
	return a;
}

MATRIX rotmat(int i, int j, double x)
{
	MATRIX a=unitmat();
	a.ent[i][i]=cos(x);
	a.ent[i][j]=-sin(x);
	a.ent[j][i]=sin(x);
	a.ent[j][j]=cos(x);
	return a;
}

ROWVECTOR transvec(VECTOR v)
{
	ROWVECTOR w;
	int k;
	for(k=0;k<SIZE;k++) w.ent[k]=v.ent[k];
	return w;
}

VECTOR transrowvec(ROWVECTOR v)
{
	VECTOR w;
	int k;
	for(k=0;k<SIZE;k++) w.ent[k]=v.ent[k];
	return w;
}

MATRIX transmat(MATRIX a)
{
	int i,j;
	MATRIX b;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			b.ent[j][i]=a.ent[i][j];
	return b;
}

double innerprod(VECTOR u, VECTOR v)
{
	double x=0.0;
	int k;
	for(k=0;k<SIZE;k++) x+=u.ent[k]*v.ent[k];
	return x;
}

double vecnorm(VECTOR v)
{
	return sqrt(innerprod(v,v));
}

double rowvecnorm(ROWVECTOR v)
{
	return vecnorm(transrowvec(v));
}

VECTOR randomvec(double x)
{
	VECTOR v;
	int i;
	for(i=0;i<SIZE;i++) v.ent[i]=2*REALRANDOM(x)-x;
	return v;
}

VECTOR randomintvec(int x)
{
	VECTOR v;
	int i;
	for(i=0;i<SIZE;i++) v.ent[i]=RANDOM(2*x-1)-x+1;
	return v;
}

ROWVECTOR randomrowvec(double x)
{
	ROWVECTOR v;
	int i;
	for(i=0;i<SIZE;i++) v.ent[i]=2*REALRANDOM(x)-x;
	return v;
}

ROWVECTOR randomintrowvec(int x)
{
	ROWVECTOR v;
	int i;
	for(i=0;i<SIZE;i++) v.ent[i]=RANDOM(2*x-1)-x+1;
	return v;
}

MATRIX randommat(double x)
{
	MATRIX a;
	int i,j;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			a.ent[i][j]=2*REALRANDOM(x)-x;
	return a;
}

MATRIX randomintmat(int x)
{
	MATRIX a;
	int i,j;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			a.ent[i][j]=RANDOM(2*x-1)-x+1;
	return a;
}

MATRIX array2mat(int a[SIZE][SIZE])
{
	int i,j;
	MATRIX b;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			b.ent[i][j]=a[i][j];
	return b;
}