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

这是一个代数(大部分)解:

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;
}

编辑:固定行索引()

本回复末尾有一行,解释如下:—)

系数数组的第一行(索引中的第0行)有m个元素,第二行(第1行)有m-1个元素,依此类推,共有m+(m-1)+…+1=m(m+1)/2个元素。

从末端开始工作稍微容易一些,因为系数数组 总是 最后一行有1个元素(行M-1),第二行有2个元素(行M-2),行M-3有3个元素,依此类推。最后k行占据系数数组的最后k(k+1)/2个位置。

现在假设在系数数组中有一个索引i。在大于i的位置有m(m+1)/2-1-i元素。称这个数字为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行),所以行索引(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]中减去此行的起始位置即可。令人欣慰的是,这两个表达给出了相同的答案!]

(伪)代码[使用i i而不是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,您可以证明Lead((SqRT(8II+ 1)-1)/2)=(ISQRT(8II+1)-1)/2,其中ISQRT(x)=Lead(Sqt(x))是整数平方根,并且除法是整数除法(截断;C/C/+/Java等中的缺省值)-因此,如果您担心精度问题,则只需要担心实现正确的整数平方根。

可能有一个聪明的一行程序,但是(减去任何错误检查):

unsigned int row_index( unsigned int i, unsigned int M ){
    unsigned int row = 0;
    unsigned int delta = M - 1;
    for( unsigned int x = delta; x < i; x += delta-- ){
        row++;
    }
    return row;
}

unsigned int column_index( unsigned int i, unsigned int M ){
    unsigned int row = 0;
    unsigned int delta = M - 1;
    unsigned int x;
    for( x = delta; x < i; x += delta-- ){
        row++;
    }
    return M + i - x - 1;
}

一定是那样的

i == col + row*(M-1)-row*(row-1)/2

因此,查找列和行的一种方法是迭代行的可能值:

for(row = 0; row < M; row++){
  col = i - row*(M-1)-row*(row-1)/2
  if (row <= col < M) return (row,column);
}

这至少是非递归的,我不知道您是否可以在不进行迭代的情况下完成这项工作。

从这个和其他答案中可以看到,为了计算列,您几乎必须计算行,所以在一个函数中同时执行这两个操作是明智的。

史里瓦察的解释很好,帮助我解决了我的问题。但是,为列索引提供的解释和代码给出的答案与问题要求的答案稍有不同。

重申一下,在i之后的行中有j'=i'-k(k+1)/2个元素,但与其他行一样,行中也有m个元素。因此,(从零开始)列索引是y=m-1-j'。

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

Shrevatsar,k-j'给出的答案是当从矩阵的对角线开始计数(零)时的列索引。因此,他的计算给出了列_index(7,4)=0,而不是问题中规定的列_index(7,4)=2。

我想了一下,得到了下面的结果。请注意,在一个快照中同时获得行和列。

假设:行从0开始。列从0开始。索引从0开始

记数法

n=矩阵大小(原始问题中是m)

m=元素索引

psuedo代码是

function ind2subTriu(m,N)
{
  d = 0;
  i = -1;
  while d < m
  {
    i = i + 1
    d = i*(N-1) - i*(i-1)/2
  }
  i0 = i-1;
  j0 = m - i0*(N-1) + i0*(i0-1)/2 + i0 + 1;
  return i0,j0
}

和一些八度/matlab代码

function [i0 j0]= ind2subTriu(m,N)
 I = 0:N-2;
 d = I*(N-1)-I.*(I-1)/2;
 i0 = I(find (d < m,1,'last'));
 j0 = m - d(i0+1) + i0 + 1;

你怎么认为?

截至2011年12月,GNU/Octave增加了很好的代码。它们可能会扩展sub2ind和ind2sub。代码现在可以作为私有函数找到。 ind2sub_tril 和 sub2ind_tril

我花了些时间来理解你需要什么!:)

unsigned int row_index(int i, int m)
{
    int iCurrentRow = 0;
    int iTotalItems = 0;
    for(int j = m; j > 0; j--)
    {
        iTotalItems += j;

        if( (i+1) <= iTotalItems)
            return iCurrentRow;

        iCurrentRow ++;
    }

    return -1; // Not checking if "i" can be in a MxM matrix.
}

抱歉,忘记了其他功能…..

unsigned int column_index(int i, int m)
{
    int iTotalItems = 0;
    for(int j = m; j > 0; j--)
    {
        iTotalItems += j;

        if( (i+1) <= iTotalItems)
            return m - (iTotalItems - i);
    }

    return -1; // Not checking if "i" can be in a MxM matrix.
}