#include "config.h"
#include "system.h"
#include "coretypes.h"
#include "tm.h"
#include "ggc.h"
#include "varray.h"
#include "tree.h"
#include "lambda.h"
static void lambda_matrix_get_column (lambda_matrix, int, int,
lambda_vector);
lambda_matrix
lambda_matrix_new (int m, int n)
{
lambda_matrix mat;
int i;
mat = ggc_alloc (m * sizeof (lambda_vector));
for (i = 0; i < m; i++)
mat[i] = lambda_vector_new (n);
return mat;
}
void
lambda_matrix_copy (lambda_matrix mat1, lambda_matrix mat2,
int m, int n)
{
int i;
for (i = 0; i < m; i++)
lambda_vector_copy (mat1[i], mat2[i], n);
}
void
lambda_matrix_id (lambda_matrix mat, int size)
{
int i, j;
for (i = 0; i < size; i++)
for (j = 0; j < size; j++)
mat[i][j] = (i == j) ? 1 : 0;
}
bool
lambda_matrix_id_p (lambda_matrix mat, int size)
{
int i, j;
for (i = 0; i < size; i++)
for (j = 0; j < size; j++)
{
if (i == j)
{
if (mat[i][j] != 1)
return false;
}
else
{
if (mat[i][j] != 0)
return false;
}
}
return true;
}
void
lambda_matrix_negate (lambda_matrix mat1, lambda_matrix mat2, int m, int n)
{
int i;
for (i = 0; i < m; i++)
lambda_vector_negate (mat1[i], mat2[i], n);
}
void
lambda_matrix_transpose (lambda_matrix mat1, lambda_matrix mat2, int m, int n)
{
int i, j;
for (i = 0; i < n; i++)
for (j = 0; j < m; j++)
mat2[i][j] = mat1[j][i];
}
void
lambda_matrix_add (lambda_matrix mat1, lambda_matrix mat2,
lambda_matrix mat3, int m, int n)
{
int i;
for (i = 0; i < m; i++)
lambda_vector_add (mat1[i], mat2[i], mat3[i], n);
}
void
lambda_matrix_add_mc (lambda_matrix mat1, int const1,
lambda_matrix mat2, int const2,
lambda_matrix mat3, int m, int n)
{
int i;
for (i = 0; i < m; i++)
lambda_vector_add_mc (mat1[i], const1, mat2[i], const2, mat3[i], n);
}
void
lambda_matrix_mult (lambda_matrix mat1, lambda_matrix mat2,
lambda_matrix mat3, int m, int r, int n)
{
int i, j, k;
for (i = 0; i < m; i++)
{
for (j = 0; j < n; j++)
{
mat3[i][j] = 0;
for (k = 0; k < r; k++)
mat3[i][j] += mat1[i][k] * mat2[k][j];
}
}
}
static void
lambda_matrix_get_column (lambda_matrix mat, int n, int col,
lambda_vector vec)
{
int i;
for (i = 0; i < n; i++)
vec[i] = mat[i][col];
}
void
lambda_matrix_delete_rows (lambda_matrix mat, int rows, int from, int to)
{
int i;
int dist;
dist = to - from;
for (i = to; i < rows; i++)
mat[i - dist] = mat[i];
for (i = rows - dist; i < rows; i++)
mat[i] = NULL;
}
void
lambda_matrix_row_exchange (lambda_matrix mat, int r1, int r2)
{
lambda_vector row;
row = mat[r1];
mat[r1] = mat[r2];
mat[r2] = row;
}
void
lambda_matrix_row_add (lambda_matrix mat, int n, int r1, int r2, int const1)
{
int i;
if (const1 == 0)
return;
for (i = 0; i < n; i++)
mat[r2][i] += const1 * mat[r1][i];
}
void
lambda_matrix_row_negate (lambda_matrix mat, int n, int r1)
{
lambda_vector_negate (mat[r1], mat[r1], n);
}
void
lambda_matrix_row_mc (lambda_matrix mat, int n, int r1, int const1)
{
int i;
for (i = 0; i < n; i++)
mat[r1][i] *= const1;
}
void
lambda_matrix_col_exchange (lambda_matrix mat, int m, int col1, int col2)
{
int i;
int tmp;
for (i = 0; i < m; i++)
{
tmp = mat[i][col1];
mat[i][col1] = mat[i][col2];
mat[i][col2] = tmp;
}
}
void
lambda_matrix_col_add (lambda_matrix mat, int m, int c1, int c2, int const1)
{
int i;
if (const1 == 0)
return;
for (i = 0; i < m; i++)
mat[i][c2] += const1 * mat[i][c1];
}
void
lambda_matrix_col_negate (lambda_matrix mat, int m, int c1)
{
int i;
for (i = 0; i < m; i++)
mat[i][c1] *= -1;
}
void
lambda_matrix_col_mc (lambda_matrix mat, int m, int c1, int const1)
{
int i;
for (i = 0; i < m; i++)
mat[i][c1] *= const1;
}
static int lambda_matrix_inverse_hard (lambda_matrix, lambda_matrix, int);
int
lambda_matrix_inverse (lambda_matrix mat, lambda_matrix inv, int n)
{
if (n == 2)
{
int a, b, c, d, det;
a = mat[0][0];
b = mat[1][0];
c = mat[0][1];
d = mat[1][1];
inv[0][0] = d;
inv[0][1] = -c;
inv[1][0] = -b;
inv[1][1] = a;
det = (a * d - b * c);
if (det < 0)
{
det *= -1;
inv[0][0] *= -1;
inv[1][0] *= -1;
inv[0][1] *= -1;
inv[1][1] *= -1;
}
return det;
}
else
return lambda_matrix_inverse_hard (mat, inv, n);
}
static int
lambda_matrix_inverse_hard (lambda_matrix mat, lambda_matrix inv, int n)
{
lambda_vector row;
lambda_matrix temp;
int i, j;
int determinant;
temp = lambda_matrix_new (n, n);
lambda_matrix_copy (mat, temp, n, n);
lambda_matrix_id (inv, n);
for (j = 0; j < n; j++)
{
row = temp[j];
for (i = j; i < n; i++)
if (row[i] < 0)
{
lambda_matrix_col_negate (temp, n, i);
lambda_matrix_col_negate (inv, n, i);
}
while (lambda_vector_first_nz (row, n, j + 1) < n)
{
int min_col = lambda_vector_min_nz (row, n, j);
lambda_matrix_col_exchange (temp, n, j, min_col);
lambda_matrix_col_exchange (inv, n, j, min_col);
for (i = j + 1; i < n; i++)
{
int factor;
factor = -1 * row[i];
if (row[j] != 1)
factor /= row[j];
lambda_matrix_col_add (temp, n, j, i, factor);
lambda_matrix_col_add (inv, n, j, i, factor);
}
}
}
determinant = 1;
for (j = n - 1; j >= 0; j--)
{
int diagonal;
row = temp[j];
diagonal = row[j];
if (diagonal == 0)
abort ();
determinant = determinant * diagonal;
if (diagonal != 1)
{
for (i = 0; i < j; i++)
lambda_matrix_col_mc (inv, n, i, diagonal);
for (i = j + 1; i < n; i++)
lambda_matrix_col_mc (inv, n, i, diagonal);
row[j] = diagonal = 1;
}
for (i = j - 1; i >= 0; i--)
{
if (row[i])
{
int factor = -row[i];
lambda_matrix_col_add (temp, n, j, i, factor);
lambda_matrix_col_add (inv, n, j, i, factor);
}
}
}
return determinant;
}
void
lambda_matrix_hermite (lambda_matrix mat, int n,
lambda_matrix H, lambda_matrix U)
{
lambda_vector row;
int i, j, factor, minimum_col;
lambda_matrix_copy (mat, H, n, n);
lambda_matrix_id (U, n);
for (j = 0; j < n; j++)
{
row = H[j];
for (i = j; i < n; i++)
{
if (row[i] < 0)
{
lambda_matrix_col_negate (H, n, i);
lambda_vector_negate (U[i], U[i], n);
}
}
while (lambda_vector_first_nz (row, n, j + 1) < n)
{
minimum_col = lambda_vector_min_nz (row, n, j);
lambda_matrix_col_exchange (H, n, j, minimum_col);
lambda_matrix_row_exchange (U, j, minimum_col);
for (i = j + 1; i < n; i++)
{
factor = row[i] / row[j];
lambda_matrix_col_add (H, n, j, i, -1 * factor);
lambda_matrix_row_add (U, n, i, j, factor);
}
}
}
}
void
lambda_matrix_right_hermite (lambda_matrix A, int m, int n,
lambda_matrix S, lambda_matrix U)
{
int i, j, i0 = 0;
lambda_matrix_copy (A, S, m, n);
lambda_matrix_id (U, m);
for (j = 0; j < n; j++)
{
if (lambda_vector_first_nz (S[j], m, i0) < m)
{
++i0;
for (i = m - 1; i >= i0; i--)
{
while (S[i][j] != 0)
{
int sigma, factor, a, b;
a = S[i-1][j];
b = S[i][j];
sigma = (a * b < 0) ? -1: 1;
a = abs (a);
b = abs (b);
factor = sigma * (a / b);
lambda_matrix_row_add (S, n, i, i-1, -factor);
lambda_matrix_row_exchange (S, i, i-1);
lambda_matrix_row_add (U, m, i, i-1, -factor);
lambda_matrix_row_exchange (U, i, i-1);
}
}
}
}
}
void
lambda_matrix_left_hermite (lambda_matrix A, int m, int n,
lambda_matrix S, lambda_matrix V)
{
int i, j, i0 = 0;
lambda_matrix_copy (A, S, m, n);
lambda_matrix_id (V, m);
for (j = 0; j < n; j++)
{
if (lambda_vector_first_nz (S[j], m, i0) < m)
{
++i0;
for (i = m - 1; i >= i0; i--)
{
while (S[i][j] != 0)
{
int sigma, factor, a, b;
a = S[i-1][j];
b = S[i][j];
sigma = (a * b < 0) ? -1: 1;
a = abs (a);
b = abs (b);
factor = sigma * (a / b);
lambda_matrix_row_add (S, n, i, i-1, -factor);
lambda_matrix_row_exchange (S, i, i-1);
lambda_matrix_col_add (V, m, i-1, i, factor);
lambda_matrix_col_exchange (V, m, i, i-1);
}
}
}
}
}
int
lambda_matrix_first_nz_vec (lambda_matrix mat, int rowsize, int colsize,
int startrow)
{
int j;
bool found = false;
for (j = startrow; (j < rowsize) && !found; j++)
{
if ((mat[j] != NULL)
&& (lambda_vector_first_nz (mat[j], colsize, startrow) < colsize))
found = true;
}
if (found)
return j - 1;
return rowsize;
}
void
lambda_matrix_project_to_null (lambda_matrix B, int rowsize,
int colsize, int k, lambda_vector x)
{
lambda_matrix M1, M2, M3, I;
int determinant;
M1 = lambda_matrix_new (colsize, colsize);
lambda_matrix_transpose (B, M1, rowsize, colsize);
M2 = lambda_matrix_new (colsize, colsize);
lambda_matrix_mult (B, M1, M2, rowsize, colsize, rowsize);
M3 = lambda_matrix_new (colsize, colsize);
determinant = lambda_matrix_inverse (M2, M3, rowsize);
lambda_matrix_mult (M1, M3, M2, colsize, rowsize, rowsize);
lambda_matrix_mult (M2, B, M1, colsize, rowsize, colsize);
lambda_matrix_negate (M1, M1, colsize, colsize);
I = lambda_matrix_new (colsize, colsize);
lambda_matrix_id (I, colsize);
lambda_matrix_add_mc (I, determinant, M1, 1, M2, colsize, colsize);
lambda_matrix_get_column (M2, colsize, k - 1, x);
}
void
lambda_matrix_vector_mult (lambda_matrix matrix, int m, int n,
lambda_vector vec, lambda_vector dest)
{
int i, j;
lambda_vector_clear (dest, m);
for (i = 0; i < m; i++)
for (j = 0; j < n; j++)
dest[i] += matrix[i][j] * vec[j];
}
void
print_lambda_matrix (FILE * outfile, lambda_matrix matrix, int m, int n)
{
int i;
for (i = 0; i < m; i++)
print_lambda_vector (outfile, matrix[i], n);
fprintf (outfile, "\n");
}