summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--lab04/Makefile4
-rw-r--r--lab04/README.md7
-rw-r--r--lab04/matrix.h381
-rw-r--r--lab04/tests/Makefile2
-rw-r--r--lab04/tests/test.c271
-rw-r--r--lab04/ui.c248
6 files changed, 907 insertions, 6 deletions
diff --git a/lab04/Makefile b/lab04/Makefile
index dc9f8ab..15c60c0 100644
--- a/lab04/Makefile
+++ b/lab04/Makefile
@@ -16,8 +16,10 @@ test:
$(TARGET): $(OBJECTS)
$(CC) $(CFLAGS) $(OBJECTS) $(LDLIBS) -o $@
+ui.o: matrix.h
+
%.o: %.c
$(CC) $(CFLAGS) -c $< -o $@
clean:
- rm -f $(TARGET) $(OBJECTS)
+ rm -f $(TARGET) $(OBJECTS) tests/test_runner
diff --git a/lab04/README.md b/lab04/README.md
index c2b753b..b5dd0a5 100644
--- a/lab04/README.md
+++ b/lab04/README.md
@@ -7,8 +7,9 @@
- 线性表的添加、插入、删除、修改、查找和遍历；
- 向量加法、减法和夹角余弦计算；
- 顺序表多项式的加减、乘法和求导；
-- 链表多项式的加减、乘法和求导。
-- 四则运算表达式求值。
+- 链表多项式的加减、乘法和求导；
+- 四则运算表达式求值；
+- 实矩阵加减乘、转置、行列式和特征值计算（特征值可为复数）。
底层实现使用自定义的 `vec` 动态数组和循环双向链表。
@@ -38,7 +39,7 @@ make clean
./calculator
```
-在主菜单中选择线性表、向量、顺序表多项式，链表多项式或者四则运算表达式功能。菜单支持方向键或 `j`/`k` 移动，按 Enter 确认，按 `Esc` 或 `q` 返回或退出。
+在主菜单中选择线性表、向量、顺序表多项式、链表多项式、表达式计算或矩阵运算。矩阵运算入口为菜单项 6，先输入矩阵 A、B 的行列数及各行元素，再选择运算。菜单支持方向键或 `j`/`k` 移动，按 Enter 确认，按 `Esc` 或 `q` 返回或退出。
运行测试：
diff --git a/lab04/matrix.h b/lab04/matrix.h
new file mode 100644
index 0000000..5d3abf8
--- /dev/null
+++ b/lab04/matrix.h
@@ -0,0 +1,381 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+/* Copyright (C) 2026 Richard Qin */
+
+#ifndef MATRIX_H
+#define MATRIX_H
+
+#include <float.h>
+#include <math.h>
+
+#include "vec.h"
+
+/* Real elements in row-major order: elements[row * cols + col]. */
+struct matrix {
+ size_t rows;
+ size_t cols;
+ struct vec elements;
+};
+
+struct matrix_eigenvalue {
+ double real;
+ double imag;
+};
+
+static inline void matrix_init(struct matrix *matrix)
+{
+ matrix->rows = 0;
+ matrix->cols = 0;
+ vec_init(&matrix->elements, double);
+}
+
+static inline void matrix_destroy(struct matrix *matrix)
+{
+ vec_destroy(&matrix->elements);
+ matrix->rows = 0;
+ matrix->cols = 0;
+}
+
+/* Retains the leading elements of the flat storage; new elements are zeroed. */
+static inline bool matrix_resize(struct matrix *matrix, size_t rows, size_t cols)
+{
+ if (rows == 0 || cols == 0 || rows > SIZE_MAX / cols) return false;
+ if (!vec_resize_zero(&matrix->elements, rows * cols)) return false;
+ matrix->rows = rows;
+ matrix->cols = cols;
+ return true;
+}
+
+/* Indices must be within the matrix dimensions. */
+static inline double *matrix_at(struct matrix *matrix, size_t row, size_t col)
+{
+ return (double *)matrix->elements.data + row * matrix->cols + col;
+}
+
+static inline const double *matrix_at_const(const struct matrix *matrix,
+ size_t row, size_t col)
+{
+ return (const double *)matrix->elements.data + row * matrix->cols + col;
+}
+
+static inline void __matrix_replace(struct matrix *dst, struct matrix *src)
+{
+ struct matrix old = *dst;
+ *dst = *src;
+ *src = old;
+ matrix_destroy(src);
+}
+
+/* All matrix arguments must be initialized. Destinations may alias operands.
+ * On failure, destinations and output values are unchanged.
+ */
+static inline bool matrix_assign(struct matrix *dst, size_t rows, size_t cols,
+ const double *elements)
+{
+ struct matrix temporary;
+ matrix_init(&temporary);
+ if (!matrix_resize(&temporary, rows, cols)) return false;
+ memcpy(temporary.elements.data, elements, temporary.elements.size * sizeof(double));
+ __matrix_replace(dst, &temporary);
+ return true;
+}
+
+static inline bool matrix_copy(struct matrix *dst, const struct matrix *src)
+{
+ struct matrix temporary;
+ matrix_init(&temporary);
+ if (!vec_append(&temporary.elements, src->elements.data, src->elements.size)) {
+ matrix_destroy(&temporary);
+ return false;
+ }
+ temporary.rows = src->rows;
+ temporary.cols = src->cols;
+ __matrix_replace(dst, &temporary);
+ return true;
+}
+
+static inline bool matrix_add(struct matrix *dst, const struct matrix *lhs,
+ const struct matrix *rhs)
+{
+ if (lhs->rows != rhs->rows || lhs->cols != rhs->cols) return false;
+ struct matrix temporary;
+ matrix_init(&temporary);
+ if (!matrix_resize(&temporary, lhs->rows, lhs->cols)) return false;
+ const double *a = lhs->elements.data;
+ const double *b = rhs->elements.data;
+ double *result = temporary.elements.data;
+ for (size_t i = 0; i < lhs->elements.size; ++i) result[i] = a[i] + b[i];
+ __matrix_replace(dst, &temporary);
+ return true;
+}
+
+static inline bool matrix_sub(struct matrix *dst, const struct matrix *lhs,
+ const struct matrix *rhs)
+{
+ if (lhs->rows != rhs->rows || lhs->cols != rhs->cols) return false;
+ struct matrix temporary;
+ matrix_init(&temporary);
+ if (!matrix_resize(&temporary, lhs->rows, lhs->cols)) return false;
+ const double *a = lhs->elements.data;
+ const double *b = rhs->elements.data;
+ double *result = temporary.elements.data;
+ for (size_t i = 0; i < lhs->elements.size; ++i) result[i] = a[i] - b[i];
+ __matrix_replace(dst, &temporary);
+ return true;
+}
+
+static inline bool matrix_mul(struct matrix *dst, const struct matrix *lhs,
+ const struct matrix *rhs)
+{
+ if (lhs->cols != rhs->rows) return false;
+ struct matrix temporary;
+ matrix_init(&temporary);
+ if (!matrix_resize(&temporary, lhs->rows, rhs->cols)) return false;
+ const double *a = lhs->elements.data;
+ const double *b = rhs->elements.data;
+ double *result = temporary.elements.data;
+ for (size_t i = 0; i < lhs->rows; ++i) {
+ for (size_t k = 0; k < lhs->cols; ++k) {
+ double value = a[i * lhs->cols + k];
+ for (size_t j = 0; j < rhs->cols; ++j) {
+ result[i * rhs->cols + j] += value * b[k * rhs->cols + j];
+ }
+ }
+ }
+ __matrix_replace(dst, &temporary);
+ return true;
+}
+
+static inline bool matrix_transpose(struct matrix *dst, const struct matrix *src)
+{
+ struct matrix temporary;
+ matrix_init(&temporary);
+ if (!matrix_resize(&temporary, src->cols, src->rows)) return false;
+ const double *a = src->elements.data;
+ double *result = temporary.elements.data;
+ for (size_t i = 0; i < src->rows; ++i) {
+ for (size_t j = 0; j < src->cols; ++j) {
+ result[j * src->rows + i] = a[i * src->cols + j];
+ }
+ }
+ __matrix_replace(dst, &temporary);
+ return true;
+}
+
+/* Gaussian elimination with partial pivoting. A singular matrix has det = 0. */
+static inline bool matrix_det(const struct matrix *matrix, double *result)
+{
+ if (matrix->rows == 0 || matrix->rows != matrix->cols) return false;
+ size_t n = matrix->rows;
+ VEC(work, long double);
+ if (!vec_resize(&work, matrix->elements.size)) return false;
+ long double *a = work.data;
+ const double *source = matrix->elements.data;
+ for (size_t i = 0; i < work.size; ++i) a[i] = source[i];
+ long double determinant = 1.0L;
+ for (size_t k = 0; k < n; ++k) {
+ size_t pivot = k;
+ for (size_t i = k + 1; i < n; ++i) {
+ if (fabsl(a[i * n + k]) > fabsl(a[pivot * n + k])) pivot = i;
+ }
+ if (a[pivot * n + k] == 0.0L) {
+ determinant = 0.0L;
+ break;
+ }
+ if (pivot != k) {
+ for (size_t j = k; j < n; ++j) {
+ long double value = a[k * n + j];
+ a[k * n + j] = a[pivot * n + j];
+ a[pivot * n + j] = value;
+ }
+ determinant = -determinant;
+ }
+ determinant *= a[k * n + k];
+ for (size_t i = k + 1; i < n; ++i) {
+ long double factor = a[i * n + k] / a[k * n + k];
+ a[i * n + k] = 0.0L;
+ for (size_t j = k + 1; j < n; ++j) {
+ a[i * n + j] -= factor * a[k * n + j];
+ }
+ }
+ }
+ *result = (double)determinant;
+ vec_destroy(&work);
+ return true;
+}
+
+/* Build a unit vector v for the reflector I - 2*v*v^T. Return the new first
+ * component of the reflected column, or zero for an already zero column.
+ */
+static inline long double __matrix_householder(long double *v, size_t count)
+{
+ long double norm = 0.0L;
+ for (size_t i = 0; i < count; ++i) norm = hypotl(norm, v[i]);
+ if (norm == 0.0L) return 0.0L;
+ long double first = -copysignl(norm, v[0]);
+ v[0] -= first;
+ norm = 0.0L;
+ for (size_t i = 0; i < count; ++i) norm = hypotl(norm, v[i]);
+ for (size_t i = 0; i < count; ++i) v[i] /= norm;
+ return first;
+}
+
+/* Apply the same reflector on the left and right, preserving eigenvalues. */
+static inline void __matrix_reflect(long double *a, size_t n,
+ size_t start, size_t count,
+ const long double *v,
+ size_t first_col, size_t rows)
+{
+ for (size_t j = first_col; j < n; ++j) {
+ long double dot = 0.0L;
+ for (size_t k = 0; k < count; ++k) dot += v[k] * a[(start + k) * n + j];
+ for (size_t k = 0; k < count; ++k) a[(start + k) * n + j] -= 2.0L * v[k] * dot;
+ }
+ for (size_t i = 0; i < rows; ++i) {
+ long double dot = 0.0L;
+ for (size_t k = 0; k < count; ++k) dot += a[i * n + start + k] * v[k];
+ for (size_t k = 0; k < count; ++k) a[i * n + start + k] -= 2.0L * dot * v[k];
+ }
+}
+
+static inline void __matrix_hessenberg(long double *a, size_t n, long double *v)
+{
+ for (size_t k = 0; k + 2 < n; ++k) {
+ size_t count = n - k - 1;
+ for (size_t i = 0; i < count; ++i) v[i] = a[(k + 1 + i) * n + k];
+ long double first = __matrix_householder(v, count);
+ if (first == 0.0L) continue;
+ __matrix_reflect(a, n, k + 1, count, v, k, n);
+ a[(k + 1) * n + k] = first;
+ for (size_t i = k + 2; i < n; ++i) a[i * n + k] = 0.0L;
+ }
+}
+
+/* Solve the trailing 2-by-2 block, including complex conjugate roots. */
+static inline void __matrix_eigen_pair(long double a, long double b,
+ long double c, long double d,
+ struct matrix_eigenvalue *first,
+ struct matrix_eigenvalue *second)
+{
+ long double half_difference = (a - d) / 2.0L;
+ long double product = b * c;
+ long double discriminant = half_difference * half_difference + product;
+ if (discriminant < 0.0L) {
+ first->real = second->real = (double)(d + half_difference);
+ first->imag = (double)sqrtl(-discriminant);
+ second->imag = -first->imag;
+ } else {
+ long double root = half_difference + copysignl(sqrtl(discriminant), half_difference);
+ first->real = (double)(d + root);
+ second->real = (double)(root == 0.0L ? d : d - product / root);
+ first->imag = second->imag = 0.0;
+ }
+}
+
+/* Compute all eigenvalues of a finite real square matrix using Hessenberg
+ * reduction and implicit double-shift QR. Output order is unspecified.
+ * values must be initialized for struct matrix_eigenvalue elements.
+ * Returns false for a non-square/empty matrix, allocation failure or failure
+ * to converge. The source and previous output remain unchanged on failure.
+ */
+static inline bool matrix_eigenvalues(const struct matrix *matrix, struct vec *values)
+{
+ if (matrix->rows == 0 || matrix->rows != matrix->cols) return false;
+ size_t n = matrix->rows;
+ VEC(work, long double);
+ VEC(vector, long double);
+ VEC(roots, struct matrix_eigenvalue);
+ bool success = false;
+ if (!vec_resize(&work, matrix->elements.size) ||
+ !vec_resize(&vector, n) || !vec_resize(&roots, n)) goto done;
+ long double *a = work.data;
+ long double *v = vector.data;
+ struct matrix_eigenvalue *eigenvalues = roots.data;
+ const double *source = matrix->elements.data;
+ for (size_t i = 0; i < work.size; ++i) a[i] = source[i];
+ __matrix_hessenberg(a, n, v);
+
+ size_t high = n;
+ size_t iterations = 0;
+ while (high > 0) {
+ size_t last = high - 1;
+ size_t low = last;
+ while (low > 0) {
+ long double lower = fabsl(a[low * n + low - 1]);
+ if (lower == 0.0L) break;
+ long double diagonal = fabsl(a[(low - 1) * n + low - 1]) + fabsl(a[low * n + low]);
+ if (lower <= LDBL_EPSILON * diagonal) {
+ /* Keep a tiny subdiagonal when its coupling to the upper
+ * entry still changes the eigenvalues appreciably.
+ */
+ long double upper = fabsl(a[(low - 1) * n + low]);
+ long double ab = fmaxl(lower, upper);
+ long double ba = fminl(lower, upper);
+ long double difference = fabsl(a[(low - 1) * n + low - 1] - a[low * n + low]);
+ long double aa = fmaxl(fabsl(a[low * n + low]), difference);
+ long double bb = fminl(fabsl(a[low * n + low]), difference);
+ long double scale = aa + ab;
+ if (ba * (ab / scale) <= LDBL_EPSILON * bb * (aa / scale)) {
+ a[low * n + low - 1] = 0.0L;
+ break;
+ }
+ }
+ --low;
+ }
+ if (low == last) {
+ eigenvalues[last] = (struct matrix_eigenvalue){ .real = (double)a[last * n + last] };
+ --high;
+ iterations = 0;
+ continue;
+ }
+ if (low + 1 == last) {
+ __matrix_eigen_pair(a[low * n + low], a[low * n + last],
+ a[last * n + low], a[last * n + last],
+ &eigenvalues[low], &eigenvalues[last]);
+ high -= 2;
+ iterations = 0;
+ continue;
+ }
+ if (++iterations > 1000 * (high - low)) goto done;
+
+ long double x = a[last * n + last];
+ long double y = a[(last - 1) * n + last - 1];
+ long double product = a[last * n + last - 1] * a[(last - 1) * n + last];
+ if (iterations % 20 == 0) {
+ /* Exceptional shifts break stalled QR cycles. */
+ long double scale = fabsl(a[last * n + last - 1]) +
+ fabsl(a[(last - 1) * n + last - 2]);
+ x += 0.75L * scale;
+ y = x;
+ product = -0.4375L * scale * scale;
+ }
+ for (size_t k = low; k < last; ++k) {
+ size_t count = k + 2 <= last ? 3 : 2;
+ if (k == low) {
+ long double diagonal = a[k * n + k];
+ long double lower = a[(k + 1) * n + k];
+ v[0] = (diagonal - x) * (diagonal - y) - product + a[k * n + k + 1] * lower;
+ v[1] = lower * (diagonal + a[(k + 1) * n + k + 1] - x - y);
+ v[2] = lower * a[(k + 2) * n + k + 1];
+ } else {
+ for (size_t i = 0; i < count; ++i) v[i] = a[(k + i) * n + k - 1];
+ }
+ long double first = __matrix_householder(v, count);
+ if (first == 0.0L) continue;
+ __matrix_reflect(a, n, k, count, v, k == low ? k : k - 1, high);
+ if (k != low) {
+ a[k * n + k - 1] = first;
+ for (size_t i = 1; i < count; ++i) a[(k + i) * n + k - 1] = 0.0L;
+ }
+ }
+ }
+ vec_swap(values, &roots);
+ success = true;
+
+done:
+ vec_destroy(&work);
+ vec_destroy(&vector);
+ vec_destroy(&roots);
+ return success;
+}
+
+#endif /* MATRIX_H */
diff --git a/lab04/tests/Makefile b/lab04/tests/Makefile
index 5724f6f..93e1d65 100644
--- a/lab04/tests/Makefile
+++ b/lab04/tests/Makefile
@@ -8,7 +8,7 @@ TARGET := test_runner
all: $(TARGET)
-$(TARGET): test.c ../container_of.h ../list.h ../vec.h ../polynomial.h ../expression.h ../stack.h
+$(TARGET): test.c ../container_of.h ../list.h ../vec.h ../polynomial.h ../expression.h ../stack.h ../matrix.h
$(CC) $(CFLAGS) -I.. test.c $(LDLIBS) -o $@
run: $(TARGET)
diff --git a/lab04/tests/test.c b/lab04/tests/test.c
index 6ce4c7e..870d717 100644
--- a/lab04/tests/test.c
+++ b/lab04/tests/test.c
@@ -4,6 +4,7 @@
#include "container_of.h"
#include "expression.h"
#include "list.h"
+#include "matrix.h"
#include "polynomial.h"
#include "vec.h"
@@ -750,6 +751,269 @@ static bool test_expression_conditional_and_call_errors(void)
return true;
}
+static bool test_matrix_storage(void)
+{
+ struct matrix a, b;
+ matrix_init(&a);
+ matrix_init(&b);
+ CHECK(a.rows == 0 && a.cols == 0 && a.elements.size == 0);
+ double data[] = { 1, 2, 3, 4, 5, 6 };
+ CHECK(matrix_assign(&a, 2, 3, data));
+ CHECK(a.rows == 2 && a.cols == 3 && a.elements.size == 6);
+ CHECK(*matrix_at_const(&a, 1, 2) == 6);
+ CHECK(matrix_assign(&a, 2, 3, a.elements.data));
+ CHECK(matrix_copy(&b, &a));
+ *matrix_at(&b, 1, 2) = 9;
+ CHECK(*matrix_at_const(&a, 1, 2) == 6);
+ CHECK(matrix_copy(&b, &b));
+ CHECK(*matrix_at_const(&b, 1, 2) == 9);
+ CHECK(matrix_resize(&a, 3, 3));
+ CHECK(memcmp(a.elements.data, data, sizeof(data)) == 0);
+ for (size_t i = 6; i < 9; ++i) CHECK(((double *)a.elements.data)[i] == 0);
+ CHECK(matrix_resize(&a, 1, 2));
+ CHECK(a.rows == 1 && a.cols == 2 && a.elements.size == 2);
+ CHECK(*matrix_at_const(&a, 0, 1) == 2);
+ matrix_destroy(&a);
+ CHECK(a.rows == 0 && a.cols == 0 && a.elements.data == NULL);
+ matrix_init(&a);
+ CHECK(matrix_copy(&b, &a));
+ CHECK(b.rows == 0 && b.cols == 0 && b.elements.size == 0);
+ matrix_destroy(&a);
+ matrix_destroy(&b);
+ return true;
+}
+
+static bool test_matrix_add_sub_and_alias(void)
+{
+ struct matrix a, b, result;
+ matrix_init(&a);
+ matrix_init(&b);
+ matrix_init(&result);
+ double left[] = { 1, 2, 3, 4, 5, 6 };
+ double right[] = { 6, 5, 4, 3, 2, 1 };
+ double sum[] = { 7, 7, 7, 7, 7, 7 };
+ double difference[] = { -5, -3, -1, 1, 3, 5 };
+ CHECK(matrix_assign(&a, 2, 3, left));
+ CHECK(matrix_assign(&b, 2, 3, right));
+ CHECK(matrix_add(&result, &a, &b));
+ CHECK(result.rows == 2 && result.cols == 3);
+ CHECK(memcmp(result.elements.data, sum, sizeof(sum)) == 0);
+ CHECK(matrix_sub(&result, &a, &b));
+ CHECK(memcmp(result.elements.data, difference, sizeof(difference)) == 0);
+ CHECK(memcmp(a.elements.data, left, sizeof(left)) == 0);
+ CHECK(memcmp(b.elements.data, right, sizeof(right)) == 0);
+ CHECK(matrix_add(&a, &a, &b));
+ CHECK(memcmp(a.elements.data, sum, sizeof(sum)) == 0);
+ CHECK(matrix_sub(&b, &a, &b));
+ CHECK(memcmp(b.elements.data, left, sizeof(left)) == 0);
+ CHECK(matrix_sub(&a, &a, &a));
+ for (size_t i = 0; i < a.elements.size; ++i) {
+ CHECK(((double *)a.elements.data)[i] == 0);
+ }
+ matrix_destroy(&a);
+ matrix_destroy(&b);
+ matrix_destroy(&result);
+ return true;
+}
+
+static bool test_matrix_multiply_and_alias(void)
+{
+ struct matrix a, b, result;
+ matrix_init(&a);
+ matrix_init(&b);
+ matrix_init(&result);
+ double left[] = { 1, 2, 3, 4, 5, 6 };
+ double right[] = { 7, 8, 9, 10, 11, 12 };
+ double product[] = { 58, 64, 139, 154 };
+ CHECK(matrix_assign(&a, 2, 3, left));
+ CHECK(matrix_assign(&b, 3, 2, right));
+ CHECK(matrix_mul(&result, &a, &b));
+ CHECK(result.rows == 2 && result.cols == 2);
+ CHECK(memcmp(result.elements.data, product, sizeof(product)) == 0);
+ CHECK(memcmp(a.elements.data, left, sizeof(left)) == 0);
+ CHECK(memcmp(b.elements.data, right, sizeof(right)) == 0);
+ CHECK(matrix_mul(&a, &a, &b));
+ CHECK(a.rows == 2 && a.cols == 2);
+ CHECK(memcmp(a.elements.data, product, sizeof(product)) == 0);
+ CHECK(matrix_assign(&a, 2, 3, left));
+ CHECK(matrix_mul(&b, &a, &b));
+ CHECK(b.rows == 2 && b.cols == 2);
+ CHECK(memcmp(b.elements.data, product, sizeof(product)) == 0);
+ double square[] = { 1, 2, 3, 4 };
+ double squared[] = { 7, 10, 15, 22 };
+ CHECK(matrix_assign(&a, 2, 2, square));
+ CHECK(matrix_mul(&a, &a, &a));
+ CHECK(memcmp(a.elements.data, squared, sizeof(squared)) == 0);
+ matrix_destroy(&a);
+ matrix_destroy(&b);
+ matrix_destroy(&result);
+ return true;
+}
+
+static bool test_matrix_transpose(void)
+{
+ struct matrix a, result;
+ matrix_init(&a);
+ matrix_init(&result);
+ double data[] = { 1, 2, 3, 4, 5, 6 };
+ double transposed[] = { 1, 4, 2, 5, 3, 6 };
+ CHECK(matrix_assign(&a, 2, 3, data));
+ CHECK(matrix_transpose(&result, &a));
+ CHECK(result.rows == 3 && result.cols == 2);
+ CHECK(memcmp(result.elements.data, transposed, sizeof(transposed)) == 0);
+ CHECK(memcmp(a.elements.data, data, sizeof(data)) == 0);
+ CHECK(matrix_transpose(&a, &a));
+ CHECK(a.rows == 3 && a.cols == 2);
+ CHECK(memcmp(a.elements.data, transposed, sizeof(transposed)) == 0);
+ CHECK(matrix_transpose(&a, &a));
+ CHECK(a.rows == 2 && a.cols == 3);
+ CHECK(memcmp(a.elements.data, data, sizeof(data)) == 0);
+ CHECK(matrix_assign(&a, 1, 6, data));
+ CHECK(matrix_transpose(&a, &a));
+ CHECK(a.rows == 6 && a.cols == 1);
+ CHECK(memcmp(a.elements.data, data, sizeof(data)) == 0);
+ matrix_destroy(&a);
+ matrix_destroy(&result);
+ return true;
+}
+
+static bool test_matrix_determinant(void)
+{
+ struct matrix matrix;
+ matrix_init(&matrix);
+ const struct {
+ size_t size;
+ double data[9];
+ double expected;
+ } cases[] = {
+ { 1, { -7 }, -7 },
+ { 2, { 0, 2, 3, 4 }, -6 },
+ { 2, { 1, 2, 2, 4 }, 0 },
+ { 2, { 0, 0, 0, 0 }, 0 },
+ { 3, { 6, 1, 1, 4, -2, 5, 2, 8, 7 }, -306 },
+ { 2, { 0, 1e200, 1e-200, 0 }, -1 },
+ };
+ for (size_t i = 0; i < sizeof(cases) / sizeof(cases[0]); ++i) {
+ size_t n = cases[i].size;
+ double result = 123;
+ CHECK(matrix_assign(&matrix, n, n, cases[i].data));
+ CHECK(matrix_det(&matrix, &result));
+ CHECK(fabs(result - cases[i].expected) < 1e-12);
+ CHECK(memcmp(matrix.elements.data, cases[i].data, n * n * sizeof(double)) == 0);
+ }
+ matrix_destroy(&matrix);
+ return true;
+}
+
+static bool test_matrix_eigenvalues(void)
+{
+ struct matrix matrix;
+ matrix_init(&matrix);
+ VEC(values, struct matrix_eigenvalue);
+ const struct {
+ size_t size;
+ double data[16];
+ struct matrix_eigenvalue expected[4];
+ double tolerance;
+ } cases[] = {
+ { 1, { -3 }, { { -3, 0 } }, 1e-12 },
+ { 2, { 2, 1, 1, 2 }, { { 1, 0 }, { 3, 0 } }, 1e-12 },
+ { 2, { 0, -1, 1, 0 }, { { 0, 1 }, { 0, -1 } }, 1e-12 },
+ { 2, { 0, 0, 0, 0 }, { { 0, 0 }, { 0, 0 } }, 1e-12 },
+ { 3, { 0, 0, 6, 1, 0, -11, 0, 1, 6 },
+ { { 1, 0 }, { 2, 0 }, { 3, 0 } }, 1e-10 },
+ { 3, { 0, 0, 2, 1, 0, -1, 0, 1, 2 },
+ { { 2, 0 }, { 0, 1 }, { 0, -1 } }, 1e-10 },
+ { 3, { 2, 4, -1, 0, -3, 2, 0, 0, 5 },
+ { { 2, 0 }, { -3, 0 }, { 5, 0 } }, 1e-12 },
+ { 3, { 2, 0, 0, 0, 2, 0, 0, 0, -1 },
+ { { 2, 0 }, { 2, 0 }, { -1, 0 } }, 1e-12 },
+ /* Similar to diag(6, 18) combined with [0, -12; 12, 0]. */
+ { 4, { 9, 17, -5, 1, 3, 9, -9, 3, -21, 27, -3, 15, -3, 15, -3, 9 },
+ { { 6, 0 }, { 18, 0 }, { 0, 12 }, { 0, -12 } }, 1e-10 },
+ /* Similar to diag(6, 12, 18, 24). */
+ { 4, { 21, -7, 7, -11, -3, 17, 1, 1, 12, -8, 26, -16, 12, -8, 8, -4 },
+ { { 6, 0 }, { 12, 0 }, { 18, 0 }, { 24, 0 } }, 1e-10 },
+ { 2, { 0, -1e-200, 1e-200, 0 },
+ { { 0, 1e-200 }, { 0, -1e-200 } }, 1e-212 },
+ { 2, { 0, -1e200, 1e200, 0 },
+ { { 0, 1e200 }, { 0, -1e200 } }, 1e188 },
+ { 2, { 1, 1e300, 1e-300, 1 }, { { 0, 0 }, { 2, 0 } }, 1e-12 },
+ };
+ for (size_t i = 0; i < sizeof(cases) / sizeof(cases[0]); ++i) {
+ size_t n = cases[i].size;
+ CHECK(matrix_assign(&matrix, n, n, cases[i].data));
+ CHECK(matrix_eigenvalues(&matrix, &values));
+ CHECK(values.size == n);
+ bool matched[4] = { false };
+ for (size_t j = 0; j < n; ++j) {
+ bool found = false;
+ const struct matrix_eigenvalue *value = vec_at_const(&values, j);
+ for (size_t k = 0; k < n; ++k) {
+ if (!matched[k] &&
+ fabs(value->real - cases[i].expected[k].real) < cases[i].tolerance &&
+ fabs(value->imag - cases[i].expected[k].imag) < cases[i].tolerance) {
+ matched[k] = true;
+ found = true;
+ break;
+ }
+ }
+ CHECK(found);
+ }
+ CHECK(memcmp(matrix.elements.data, cases[i].data, n * n * sizeof(double)) == 0);
+ }
+ matrix_destroy(&matrix);
+ vec_destroy(&values);
+ return true;
+}
+
+static bool test_matrix_error_outputs(void)
+{
+ struct matrix a, b, output, empty;
+ matrix_init(&a);
+ matrix_init(&b);
+ matrix_init(&output);
+ matrix_init(&empty);
+ double data[] = { 1, 2, 3, 4, 5, 6 };
+ double sentinel[] = { 91, 92 };
+ CHECK(matrix_assign(&a, 2, 3, data));
+ CHECK(matrix_assign(&b, 2, 2, data));
+ CHECK(matrix_assign(&output, 1, 2, sentinel));
+ CHECK(!matrix_resize(&a, 0, 2));
+ CHECK(!matrix_resize(&a, 2, 0));
+ CHECK(!matrix_resize(&a, SIZE_MAX, 2));
+ CHECK(!matrix_resize(&a, SIZE_MAX / sizeof(double) + 1, 1));
+ CHECK(!matrix_assign(&a, 0, 2, data));
+ CHECK(a.rows == 2 && a.cols == 3);
+ CHECK(memcmp(a.elements.data, data, sizeof(data)) == 0);
+ CHECK(!matrix_add(&output, &a, &b));
+ CHECK(!matrix_sub(&output, &a, &b));
+ CHECK(!matrix_mul(&output, &a, &b));
+ CHECK(!matrix_transpose(&output, &empty));
+ CHECK(!matrix_add(&output, &empty, &empty));
+ CHECK(!matrix_mul(&output, &empty, &empty));
+ CHECK(output.rows == 1 && output.cols == 2);
+ CHECK(memcmp(output.elements.data, sentinel, sizeof(sentinel)) == 0);
+ double determinant = 123;
+ CHECK(!matrix_det(&a, &determinant));
+ CHECK(!matrix_det(&empty, &determinant));
+ CHECK(determinant == 123);
+ VEC(values, struct matrix_eigenvalue);
+ struct matrix_eigenvalue root = { 4, 5 };
+ CHECK(vec_push_back(&values, &root));
+ CHECK(!matrix_eigenvalues(&a, &values));
+ CHECK(!matrix_eigenvalues(&empty, &values));
+ CHECK(values.size == 1);
+ const struct matrix_eigenvalue *retained = vec_at_const(&values, 0);
+ CHECK(retained->real == 4 && retained->imag == 5);
+ matrix_destroy(&a);
+ matrix_destroy(&b);
+ matrix_destroy(&output);
+ matrix_destroy(&empty);
+ vec_destroy(&values);
+ return true;
+}
+
struct test_case {
const char *name;
bool (*function)(void);
@@ -792,6 +1056,13 @@ int main(void)
{ "expression definition errors", test_expression_definition_errors },
{ "expression conditional and call errors",
test_expression_conditional_and_call_errors },
+ { "matrix storage", test_matrix_storage },
+ { "matrix add, sub and alias", test_matrix_add_sub_and_alias },
+ { "matrix multiply and alias", test_matrix_multiply_and_alias },
+ { "matrix transpose", test_matrix_transpose },
+ { "matrix determinant", test_matrix_determinant },
+ { "matrix eigenvalues", test_matrix_eigenvalues },
+ { "matrix error outputs", test_matrix_error_outputs },
};
size_t failures = 0;
diff --git a/lab04/ui.c b/lab04/ui.c
index e6c96ba..7c36bd7 100644
--- a/lab04/ui.c
+++ b/lab04/ui.c
@@ -5,6 +5,7 @@
#include "polynomial.h"
#include "expression.h"
+#include "matrix.h"
#include "tui.h"
#include <ctype.h>
@@ -35,6 +36,13 @@ enum ui_vector_result_kind {
UI_VECTOR_RESULT_COSINE,
};
+enum ui_matrix_result_kind {
+ UI_MATRIX_RESULT_NONE,
+ UI_MATRIX_RESULT_VALUE,
+ UI_MATRIX_RESULT_DETERMINANT,
+ UI_MATRIX_RESULT_EIGENVALUES,
+};
+
struct ui_state {
struct vec linear;
@@ -44,6 +52,13 @@ struct ui_state {
enum ui_vector_result_kind vector_result_kind;
double vector_cosine;
+ struct matrix matrix_a;
+ struct matrix matrix_b;
+ struct matrix matrix_result;
+ struct vec matrix_eigenvalues;
+ enum ui_matrix_result_kind matrix_result_kind;
+ double matrix_determinant;
+
struct poly_vec poly_vec_a;
struct poly_vec poly_vec_b;
struct poly_vec poly_vec_result;
@@ -531,6 +546,10 @@ static bool ui_state_init(struct ui_state *state)
vec_init(&state->vector_a, double);
vec_init(&state->vector_b, double);
vec_init(&state->vector_result, double);
+ matrix_init(&state->matrix_a);
+ matrix_init(&state->matrix_b);
+ matrix_init(&state->matrix_result);
+ vec_init(&state->matrix_eigenvalues, struct matrix_eigenvalue);
poly_vec_init(&state->poly_vec_a);
poly_vec_init(&state->poly_vec_b);
poly_vec_init(&state->poly_vec_result);
@@ -579,6 +598,10 @@ static void ui_state_destroy(struct ui_state *state)
vec_destroy(&state->vector_a);
vec_destroy(&state->vector_b);
vec_destroy(&state->vector_result);
+ matrix_destroy(&state->matrix_a);
+ matrix_destroy(&state->matrix_b);
+ matrix_destroy(&state->matrix_result);
+ vec_destroy(&state->matrix_eigenvalues);
poly_vec_destroy(&state->poly_vec_a);
poly_vec_destroy(&state->poly_vec_b);
poly_vec_destroy(&state->poly_vec_result);
@@ -963,6 +986,225 @@ static void ui_vector_page(struct ui_state *state)
}
}
+static void ui_print_matrix(const struct matrix *matrix)
+{
+ if (matrix->rows == 0) {
+ printf("（未输入）\n");
+ return;
+ }
+ printf("%zu × %zu\n", matrix->rows, matrix->cols);
+ for (size_t row = 0; row < matrix->rows; ++row) {
+ printf(" [");
+ for (size_t col = 0; col < matrix->cols; ++col) {
+ printf(" %12.6g", *matrix_at_const(matrix, row, col));
+ }
+ printf(" ]\n");
+ }
+}
+
+static void ui_draw_matrix(const struct ui_state *state, size_t selected)
+{
+ static const char *const items[] = {
+ "输入矩阵 A",
+ "输入矩阵 B",
+ "A + B",
+ "A - B",
+ "A × B",
+ "转置矩阵",
+ "计算行列式",
+ "计算特征值",
+ };
+
+ tui_clear();
+ ui_print_header("矩阵运算");
+ printf(COLOR_CYAN "A = " STYLE_RESET);
+ ui_print_matrix(&state->matrix_a);
+ printf(COLOR_CYAN "B = " STYLE_RESET);
+ ui_print_matrix(&state->matrix_b);
+ putchar('\n');
+ ui_print_menu(items, sizeof(items) / sizeof(items[0]), selected, "返回");
+
+ switch (state->matrix_result_kind) {
+ case UI_MATRIX_RESULT_VALUE:
+ printf("\n" COLOR_GREEN "结果：" STYLE_RESET);
+ ui_print_matrix(&state->matrix_result);
+ break;
+ case UI_MATRIX_RESULT_DETERMINANT:
+ printf("\n" COLOR_GREEN "行列式 = %.12g" STYLE_RESET "\n",
+ state->matrix_determinant);
+ break;
+ case UI_MATRIX_RESULT_EIGENVALUES:
+ printf("\n" COLOR_GREEN "特征值：" STYLE_RESET "\n");
+ for (size_t i = 0; i < state->matrix_eigenvalues.size; ++i) {
+ const struct matrix_eigenvalue *value =
+ vec_at_const(&state->matrix_eigenvalues, i);
+ printf(" λ%zu = %.12g", i + 1, value->real);
+ if (value->imag != 0.0) printf("%+.12gi", value->imag);
+ putchar('\n');
+ }
+ break;
+ case UI_MATRIX_RESULT_NONE:
+ break;
+ }
+ ui_print_status(state);
+ ui_print_footer("↑/↓ 或 j/k 选择 Enter 确认 q/Esc 返回");
+ fflush(stdout);
+}
+
+static void ui_matrix_input(struct matrix *destination, const char *name,
+ struct ui_state *state)
+{
+ size_t rows, cols;
+ if (!ui_read_size_value(name, "行数（正整数）：", &rows, state) ||
+ !ui_read_size_value(name, "列数（正整数）：", &cols, state)) {
+ return;
+ }
+ if (rows == 0 || cols == 0) {
+ ui_set_status(state, "错误：矩阵行数和列数必须为正整数");
+ return;
+ }
+
+ struct matrix temporary;
+ matrix_init(&temporary);
+ if (!matrix_resize(&temporary, rows, cols)) {
+ ui_set_status(state, "错误：矩阵尺寸过大或内存分配失败");
+ return;
+ }
+ VEC(row, double);
+ for (size_t i = 0; i < rows; ++i) {
+ char input[512];
+ char description[128];
+ snprintf(description, sizeof(description),
+ "第 %zu/%zu 行：输入 %zu 个实数，以空格分隔", i + 1, rows, cols);
+ if (!ui_prompt_inline(name, description, input, sizeof(input))) {
+ ui_set_status(state, "已取消");
+ goto done;
+ }
+ if (!ui_parse_vector(input, &row) || row.size != cols) {
+ ui_set_status(state, "错误：每行必须输入 %zu 个实数", cols);
+ goto done;
+ }
+ for (size_t j = 0; j < cols; ++j) {
+ double value = *(const double *)vec_at_const(&row, j);
+ if (!isfinite(value)) {
+ ui_set_status(state, "错误：矩阵元素必须是有限实数");
+ goto done;
+ }
+ *matrix_at(&temporary, i, j) = value;
+ }
+ }
+ struct matrix old = *destination;
+ *destination = temporary;
+ temporary = old;
+ ui_set_status(state, "已更新%s", name);
+done:
+ vec_destroy(&row);
+ matrix_destroy(&temporary);
+}
+
+static void ui_matrix_binary(struct ui_state *state, char operation)
+{
+ const struct matrix *a = &state->matrix_a;
+ const struct matrix *b = &state->matrix_b;
+ if (a->rows == 0 || b->rows == 0) {
+ ui_set_status(state, "错误：请先输入矩阵 A 和 B");
+ return;
+ }
+ if (operation == '*' ? a->cols != b->rows :
+ (a->rows != b->rows || a->cols != b->cols)) {
+ ui_set_status(state, operation == '*' ?
+ "错误：A 的列数必须等于 B 的行数" :
+ "错误：两个矩阵的行数和列数必须相同");
+ return;
+ }
+ bool success;
+ switch (operation) {
+ case '+': success = matrix_add(&state->matrix_result, a, b); break;
+ case '-': success = matrix_sub(&state->matrix_result, a, b); break;
+ default: success = matrix_mul(&state->matrix_result, a, b); break;
+ }
+ if (!success) {
+ ui_set_status(state, "错误：结果尺寸过大或内存分配失败");
+ return;
+ }
+ state->matrix_result_kind = UI_MATRIX_RESULT_VALUE;
+ ui_set_status(state, "已计算 A %c B", operation);
+}
+
+static void ui_matrix_unary(struct ui_state *state, int choice)
+{
+ const char *title = choice == 5 ? "转置矩阵" :
+ choice == 6 ? "计算行列式" : "计算特征值";
+ size_t selected;
+ if (!ui_read_size_value(title, "选择矩阵：1 = A，2 = B", &selected, state)) {
+ return;
+ }
+ if (selected != 1 && selected != 2) {
+ ui_set_status(state, "错误：请输入 1 或 2");
+ return;
+ }
+ const struct matrix *matrix = selected == 1 ? &state->matrix_a : &state->matrix_b;
+ if (matrix->rows == 0) {
+ ui_set_status(state, "错误：请先输入所选矩阵");
+ return;
+ }
+ if (choice != 5 && matrix->rows != matrix->cols) {
+ ui_set_status(state, "错误：行列式和特征值计算要求方阵");
+ return;
+ }
+ bool success;
+ enum ui_matrix_result_kind kind;
+ switch (choice) {
+ case 5:
+ success = matrix_transpose(&state->matrix_result, matrix);
+ kind = UI_MATRIX_RESULT_VALUE;
+ break;
+ case 6:
+ success = matrix_det(matrix, &state->matrix_determinant);
+ kind = UI_MATRIX_RESULT_DETERMINANT;
+ break;
+ default:
+ success = matrix_eigenvalues(matrix, &state->matrix_eigenvalues);
+ kind = UI_MATRIX_RESULT_EIGENVALUES;
+ break;
+ }
+ if (!success) {
+ ui_set_status(state, choice == 7 ?
+ "错误：特征值计算失败（内存不足或迭代未收敛）" :
+ "错误：结果尺寸过大或内存分配失败");
+ return;
+ }
+ state->matrix_result_kind = kind;
+ ui_set_status(state, "%s完成（矩阵 %c）", title, selected == 1 ? 'A' : 'B');
+}
+
+static void ui_matrix_action(struct ui_state *state, int choice)
+{
+ ui_clear_status(state);
+ state->matrix_result_kind = UI_MATRIX_RESULT_NONE;
+ switch (choice) {
+ case 0: ui_matrix_input(&state->matrix_a, "矩阵 A", state); break;
+ case 1: ui_matrix_input(&state->matrix_b, "矩阵 B", state); break;
+ case 2: ui_matrix_binary(state, '+'); break;
+ case 3: ui_matrix_binary(state, '-'); break;
+ case 4: ui_matrix_binary(state, '*'); break;
+ case 5:
+ case 6:
+ case 7: ui_matrix_unary(state, choice); break;
+ }
+}
+
+static void ui_matrix_page(struct ui_state *state)
+{
+ size_t selected = 0;
+ for (;;) {
+ ui_draw_matrix(state, selected);
+ int choice = ui_menu_input(&selected, 8);
+ if (choice == -1) return;
+ if (choice >= 0) ui_matrix_action(state, choice);
+ }
+}
+
static void ui_draw_poly_vec(const struct ui_state *state, size_t selected)
{
static const char *const items[] = {
@@ -1563,6 +1805,7 @@ static void ui_draw_main(const struct ui_state *state, size_t selected)
"顺序表多项式",
"链表多项式",
"表达式计算",
+ "矩阵运算",
};
(void)state;
@@ -1598,7 +1841,7 @@ int ui_run(void)
size_t selected = 0;
for (;;) {
ui_draw_main(&state, selected);
- int choice = ui_menu_input(&selected, 5);
+ int choice = ui_menu_input(&selected, 6);
if (choice == -1) {
break;
}
@@ -1618,6 +1861,9 @@ int ui_run(void)
case 4:
ui_expression_page(&state);
break;
+ case 5:
+ ui_matrix_page(&state);
+ break;
default:
break;
}
generated by cgit v1.2.3 (git 2.46.0) at 2026-10-28 22:30:55 +0000