Алгоритм для индексных чисел треугольных матричных коэффициентов
Я думаю, что это должно быть просто, но я не могу понять это правильно.
У меня есть треугольная матрица MxM, коэффициенты которой хранятся в векторе, строка за строкой.
Например:
M = [ m00 m01 m02 m03 ]
[ m11 m12 m12 ]
[ 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
и т.д., если подсчет индекса основан на нуле.
EDIT: несколько ответов с использованием итерации. Кто-нибудь знает об алгебраических выражениях?
Ответы
Ответ 1
Здесь алгебраическое (в основном) решение:
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;
}
EDIT: fixed row_index()
Ответ 2
Одиночные вкладыши в конце этого ответа, объяснение следует: -)
В массиве коэффициентов есть: M элементов для первой строки (строка 0, в вашей индексации), (M-1) для второго (строка 1) и т.д. для всего M + (M-1) + & hellip; + 1 = M (M + 1)/2 элементов.
Немного легче работать с конца, потому что массив коэффициентов всегда имеет 1 элемент для последней строки (строка M-1), 2 элемента для второй строки (строка M-2), 3 элемента для строки М-3 и т.д. Последние K-строки занимают последние K (K + 1)/2 позиции массива коэффициентов.
Теперь предположим, что вам предоставляется индекс я в массиве коэффициентов. Есть элементы M (M + 1)/2-1-i в положениях, больших i. Назовите это число я '; вы хотите найти наибольшее k такое, что k (k + 1)/2 & le; я '. Форма этой проблемы является источником ваших страданий - насколько я вижу, вы не можете избежать квадратных корней: -)
Хорошо, пусть это все равно: k (k + 1) & le; 2i 'означает (k + 1/2) 2 - 1/4 & le; 2i 'или эквивалентно k & le; (SQRT (8i '+ 1) -1)/2. Позвольте мне назвать наибольший такой k, как K, то есть K строк, которые позже (из общего числа M строк), поэтому row_index (i, M) является M-1-K. Что касается индекса столбца, то 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 - (К + 1) (К + 2)/2]. Отрадно, что оба выражения дают один и тот же ответ!]
(Pseudo-) Code [используя ii вместо я ', поскольку некоторые языки не допускают переменные с именами, такими как 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. Возможно, вам придется быть осторожными с ошибками точности, но есть способы найти целочисленный квадратный корень, который не требует вычисления с плавающей запятой, Кроме того, процитировать Кнута: "Остерегайтесь ошибок в приведенном выше коде, я только доказал это правильно, не пробовал".
Если я могу рискнуть сделать еще одно замечание: просто сохраняя все значения в M & times; массив M будет занимать не более двух раз больше памяти, и в зависимости от вашей проблемы коэффициент 2 может быть незначителен по сравнению с алгоритмическими улучшениями, и, возможно, стоит торговать, возможно, дорогостоящим вычислением квадратного корня для более простых выражений, которые у вас будут.
[Edit: BTW, вы можете доказать, что floor ((sqrt (8ii + 1) -1)/2) = (isqrt (8ii + 1) -1)/2 где isqrt (x) = floor (sqrt ( x)) представляет собой целочисленный квадратный корень, а деление - целочисленное деление (усечение, значение по умолчанию в C/С++/Java и т.д.) - поэтому, если вы беспокоитесь о проблемах с точностью, вам нужно только беспокоиться о реализации правильного целочисленного квадрата корень.]
Ответ 3
Должно быть, что
i == col + row*(M-1)-row*(row-1)/2
Итак, один способ найти col и row - перебирать возможные значения строки:
for(row = 0; row < M; row++){
col = i - row*(M-1)-row*(row-1)/2
if (row <= col < M) return (row,column);
}
Это, по крайней мере, нерекурсивный, я не знаю, можете ли вы сделать это без итерации.
Как видно из этого и других ответов, вы в значительной степени должны вычислять строку, чтобы вычислить столбец, поэтому было бы разумно выполнять обе функции в одной функции.
Ответ 4
Объяснение ShreevatsaR превосходно и помогло мне решить мою проблему. Однако объяснение и код, предоставленные для индекса столбца, дают немного другой ответ, чем тот, который задает вопрос.
Чтобы повторить, в строке после я есть 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
Ответ ShreevatsaR, K - j ', является индексом столбца при начале отсчета (с нулем) по диагонали матрицы. Следовательно, его вычисление дает column_index (7,4) = 0 и не, как указано в вопросе, column_index (7,4) = 2.
Ответ 5
Для них может быть умный один вкладыш, но (минус любая проверка ошибок):
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;
}
Ответ 6
Я подумал немного, и я получил следующий результат. Обратите внимание, что вы получаете как строки, так и столбцы на один снимок.
Предположения: строки начинаются с 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 года действительно хороший код для этого был добавлен в GNU/Octave. Потенциально они будут расширять sub2ind и ind2sub. Код может быть найден на данный момент как частные функции ind2sub_tril и sub2ind_tril
Ответ 7
Пришло время понять, что вам нужно!:)
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.
}