三角矩阵系数的索引数算法

我认为这一定很简单,但我说得不对劲......

我有一个MxM三角矩阵,其系数逐行存储在矢量中.例如:

M =   [ m00 m01 m02 m03 ] 
      [     m11 m12 m13 ]
      [         m22 m23 ]
      [             m33 ]

存储为

coef[ m00 m01 m02 m03 m11 m12 m13 m22 m23 m33 ]

现在我正在寻找一种非递归算法,它给出了矩阵大小M和系数数组索引i

unsigned int row_index(i,M)

和

unsigned int column_index(i,M)

它所指的矩阵元素.所以 row_index(9,4) == 3,column_index(7,4) == 2等等.如果指数计数是从零开始.

编辑:已经给出了几个使用迭代的回复.有谁知道代数表达式？

1> ShreevatsaR..：

---

本回复结尾的单行,解释如下:-)

系数数组具有:第一行的M个元素(索引中的第0行),第二行(第1行)的(M-1),依此类推,总共M +(M-1)+ ... + 1 = M(M + 1)/ 2个元素.

从最后开始工作会稍微容易一些,因为系数数组总是有1个元素用于最后一行(行M-1),2个元素用于倒数第二行(行M-2),3个元素用于行M-3 , 等等.最后的K行占据系数阵列的最后K(K + 1)/ 2个位置.

现在假设在系数数组中给出了索引i.在大于i的位置处存在M(M + 1)/ 2-1-i个元素.拨打这个号码我'; 你想找到最大的k,使得k(k + 1)/2≤i'.这个问题的形式是你苦难的根源 - 据我所知,你不能避免采取平根:-)

好吧,我们无论如何都要这样做:k(k + 1)≤2i'表示(k + 1/2)² - 1/4≤2i',或等效地k≤(sqrt(8i'+ 1)-1)/ 2.让我将最大的k称为K,然后有K行(稍后是总共M行),因此row_index(i,M)是M-1-K.对于列索引,i'中的K(K + 1)/ 2个元素在后面的行中,因此该行中有j'= i'-K(K + 1)/ 2个后面的元素(具有总共K + 1个元素),因此列索引为K-j'.[或等效地,此行从结尾开始于(K + 1)(K + 2)/ 2,因此我们只需要从i中减去该行的起始位置:i- [M(M + 1)/ 2 - (K + 1)(K + 2)/ 2].令人鼓舞的是,两个表达都给出了相同的答案!]

(伪)代码[使用ii代替i',因为某些语言不允许使用名称为i';-)的变量:

row_index(i, M):
    ii = M(M+1)/2-1-i
    K = floor((sqrt(8ii+1)-1)/2)
    return M-1-K

column_index(i, M):
    ii = M(M+1)/2-1-i
    K = floor((sqrt(8ii+1)-1)/2)
    return i - M(M+1)/2 + (K+1)(K+2)/2

当然,您可以通过将表达式替换为ii和K来将它们转换为单行.您可能必须小心精度错误,但是有一些方法可以找到不需要浮点计算的整数平方根.另外,引用Knuth:"注意上面代码中的错误;我只是证明它是正确的,没有尝试过."

如果我可以冒昧地进一步评论:简单地将所有值保存在M×M阵列中将最多占用内存的两倍,并且根据您的问题,与算法改进相比,因子2可能是微不足道的,并且可能是值得交易平方根的可能昂贵的计算,以获得更简单的表达式.

[编辑:BTW,你可以证明楼层((sqrt(8ii + 1)-1)/ 2)=(isqrt(8ii + 1)-1)/ 2其中isqrt(x)= floor(sqrt(x))是整数平方根,除法是整数除法(截断; C/C++/Java等中的默认值) - 所以如果你担心精度问题,你只需要担心实现一个正确的整数平方根.

---

@frank当然,它很相似.最后的K行一起具有K(K-1)/ 2个元素.因此,给定i(基于0的索引),我们想要最大K,使得K(K-1)/2≤M(M-1)/ 2-1-i.这可以达到K = floor(1/2*(1 + sqrt(4M*(M-1)-8*i + 7))).利用该值K,行索引(从0开始)是M-1-K.

2> Matt Lewis..：

---

这是一个代数(大多数)解决方案:

unsigned int row_index( unsigned int i, unsigned int M ){
    double m = M;
    double row = (-2*m - 1 + sqrt( (4*m*(m+1) - 8*(double)i - 7) )) / -2;
    if( row == (double)(int) row ) row -= 1;
    return (unsigned int) row;
}


unsigned int column_index( unsigned int i, unsigned int M ){
    unsigned int row = row_index( i, M);
    return  i - M * row + row*(row+1) / 2;
}

编辑:修复row_index()

---

什么是神奇的数字2,1,4,8和7？