半三角矩阵的压缩存取及矩阵乘法的实现

在数组的压缩存储中我们学过了对称矩阵的存储方式，半三角矩阵的存储方式是类似的。对于n阶方阵，总共需要存储1+2+3…+n = n(n+1)/2 个元素。若以行序为主序，
1.已知i,j，则对应的offest(i,j) = (i+1)*i/2 + j。
2.已知k=offest(i,j)，则i = (-1+sqrt(1 + 8k)) / 2 向下取整，j = k - (i+1)*i/2。

下面来推导下半角矩阵的乘法公式，以3x3矩阵为例，C=AxB则
c20 = a20 * b00 + a21 * b10 + a22 * b20
c21 = a21 * b11 + a22 * b21
c22 = a22 * b22

观察上面几个式子我们很容易得到cij的数学表达式：
若n阶下半三角矩阵采用一维数组压缩存储，矩阵C中第k个元素的值的求法如下：

//矩阵乘法
i = (sqrt(1 + 8 * k) - 1) / 2;  //获取位置i,j
j = k - i * (i + 1) / 2;
for (int l = j; l <= i; l++){
	c += a[i* (i + 1) / 2 + l] * b[l * (l + 1) / 2 + j];
}       

实现代码：

#include<iostream>
#include<math.h>
using namespace std;

void printMartix(int* martix, int n) {
	cout << "计算结果：" << endl;
	for (int i = 0; i < n; i++) {
		for (int j = 0; j <= i; j++)
		{
			cout << martix[i * (i + 1) / 2 + j] << " ";
		}
		cout << endl;
	}
}
void init(int *martix1, int *martix2, int *martix3, int n)
{
	int elem;

	cout << "矩阵1" << endl;
	for (int i = 0; i < (n + 1) * n / 2; i++) {
		cin >> elem;
		martix1[i] = elem;
	}

	cout << "矩阵2" << endl;
	for (int i = 0; i < (n + 1) * n / 2; i++) {
		cin >> elem;
		martix2[i] = elem;
	}
}


void multi(int* martix1, int* martix2, int* martix3, int n) {
	for (int k = 0; k < n * (n + 1) / 2; k++) {
		int sum = 0;
		int i, j;
		i = (sqrt(1 + 8 * k) - 1) / 2;
		j = k - i * (i + 1) / 2;

		for (int l = j; l <= i; l++) {
			sum += martix1[i * (i + 1) / 2 + l] * martix2[l * (l + 1) / 2 + j];
		}
		martix3[k] = sum;
	}
}



int main() {
	int n;
	cout << "请输入矩阵阶数：";
	cin >> n;
	int max_size = n * (n + 1) / 2;
	int *martix1 = new int[max_size];
	int *martix2 = new int[max_size];
	int *martix3 = new int[max_size];
	init(martix1, martix2, martix3, n);
	multi(martix1, martix2, martix3, n);
	printMartix(martix3, n);
}