Welcome to OStack Knowledge Sharing Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
419 views
in Technique[技术] by (71.8m points)

algorithm for index numbers of triangular matrix coefficients

I think this must be simple but I can't get it right...

I have an MxM triangular matrix, the coefficients of which are stored in a vector, row by row. For example:

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

is stored as

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

Now I'm looking for a non-recursive algorithm that gives me for matrix size M and coefficient array index i

unsigned int row_index(i,M)

and

unsigned int column_index(i,M)

of the matrix element that it refers to. So, row_index(9,4) == 3, column_index(7,4) == 2 etc. if the index counting is zero-based.

EDIT: Several replies using an iteration have been given. Does anyone know of algebraic expressions?

See Question&Answers more detail:os

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Answer

0 votes
by (71.8m points)

One-liners at the end of this reply, explanation follows :-)

The coefficient array has: M elements for the first row (row 0, in your indexing), (M-1) for the second (row 1), and so on, for a total of M+(M-1)+…+1=M(M+1)/2 elements.

It's slightly easier to work from the end, because the coefficient array always has 1 element for the last row (row M-1), 2 elements for the second-last row (row M-2), 3 elements for row M-3, and so on. The last K rows take up the last K(K+1)/2 positions of the coefficient array.

Now suppose you are given an index i in the coefficient array. There are M(M+1)/2-1-i elements at positions greater than i. Call this number i'; you want to find the largest k such that k(k+1)/2 ≤ i'. The form of this problem is the very source of your misery -- as far as I can see, you cannot avoid taking square roots :-)

Well let's do it anyway: k(k+1) ≤ 2i' means (k+1/2)2 - 1/4 ≤ 2i', or equivalently k ≤ (sqrt(8i'+1)-1)/2. Let me call the largest such k as K, then there are K rows that are later (out of a total of M rows), so the row_index(i,M) is M-1-K. As for the column index, K(K+1)/2 elements out of the i' are in later rows, so there are j'=i'-K(K+1)/2 later elements in this row (which has K+1 elements in total), so the column index is K-j'. [Or equivalently, this row starts at (K+1)(K+2)/2 from the end, so we just need to subtract the starting position of this row from i: i-[M(M+1)/2-(K+1)(K+2)/2]. It is heartening that both expressions give the same answer!]

(Pseudo-)Code [using ii instead of i' as some languages don't allow variables with names like 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

Of course you can turn them into one-liners by substituting the expressions for ii and K. You may have to be careful about precision errors, but there are ways of finding the integer square root which do not require floating-point computation. Also, to quote Knuth: "Beware of bugs in the above code; I have only proved it correct, not tried it."

If I may venture to make a further remark: simply keeping all the values in an M×M array will take at most twice the memory, and depending on your problem, the factor of 2 might be insignificant compared to algorithmic improvements, and might be worth trading the possibly expensive computation of the square root for the simpler expressions you'll have.

[Edit: BTW, you can prove that floor((sqrt(8ii+1)-1)/2) = (isqrt(8ii+1)-1)/2 where isqrt(x)=floor(sqrt(x)) is integer square root, and the division is integer division (truncation; the default in C/C++/Java etc.) -- so if you're worried about precision issues you only need to worry about implementing a correct integer square root.]


与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome to OStack Knowledge Sharing Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...