алгоритм для числовых индексов треугольных матричных коэффициентов - PullRequest
22 голосов
/ 28 октября 2008

Я думаю, это должно быть просто, но я не могу понять это правильно ...

У меня есть 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 и т. Д., Если подсчет индекса начинается с нуля.

РЕДАКТИРОВАТЬ: было дано несколько ответов с использованием итерации. Кто-нибудь знает алгебраические выражения?

Ответы [ 7 ]

20 голосов
/ 28 октября 2008

Однострочники в конце этого ответа, пояснение следующее: -)

Массив коэффициентов содержит: 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 позиции массива коэффициентов.

Теперь предположим, что вам дан индекс i в массиве коэффициентов. Есть элементы 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 элемента из i 'находятся в более поздних строках, поэтому в этой строке есть j' = i'-K (K + 1) / 2 более поздних элемента (который имеет K + 1 элементов), поэтому индекс столбца K-j '. [Или, что эквивалентно, эта строка начинается с (K + 1) (K + 2) / 2 с конца, поэтому нам просто нужно вычесть начальную позицию этой строки из i: i- [M (M + 1) / 2 - (К + 1) (К + 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. Возможно, вам следует быть осторожным с ошибками точности, но есть способы найти целочисленный квадратный корень, который не требует вычисления с плавающей запятой , Кроме того, процитирую Кнута: «Остерегайтесь ошибок в приведенном выше коде; я только доказал, что это правильно, а не пробовал».

Если я позволю себе сделать еще одно замечание: простое хранение всех значений в массиве M & times; M займет максимум в два раза больше памяти, и в зависимости от вашей проблемы коэффициент 2 может быть незначительным по сравнению с алгоритмическими улучшениями, и, возможно, стоит обменять, возможно, дорогостоящее вычисление квадратного корня на более простые выражения, которые у вас будут.

[Редактировать: Кстати, вы можете доказать, что floor ((sqrt (8ii + 1) -1) / 2) = (isqrt (8ii + 1) -1) / 2, где isqrt (x) = floor (sqrt ( x)) является целочисленным квадратным корнем, а деление является целочисленным делением (усечение; значение по умолчанию в C / C ++ / Java и т. д.), поэтому, если вас беспокоит проблема точности, вам нужно только беспокоиться о реализации правильного целочисленного квадрата корень.]

7 голосов
/ 28 октября 2008

Вот алгебраическое (в основном) решение:

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 голосов
/ 23 января 2015

Объяснение ShreevatsaR превосходно и помогло мне решить мою проблему. Тем не менее, объяснение и код, предоставленные для индекса столбца, дают немного другой ответ, чем вопрос.

Повторюсь, в строке после i находится j '= i' - K (K + 1) / 2 элемента. Но ряд, как и любой другой ряд, имеет М элементов. Следовательно, (основанный на нуле) индекс столбца равен 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.

2 голосов
/ 28 октября 2008

Может быть один умный лайнер для них, но (минус любая проверка ошибок):

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;
}
2 голосов
/ 28 октября 2008

Должно быть, что

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

Это как минимум нерекурсивно, я не знаю, можно ли сделать это без итерации.

Как видно из этого и других ответов, вам в значительной степени приходится вычислять строку, чтобы вычислить столбец, поэтому было бы разумно выполнить обе функции в одной функции.

1 голос
/ 30 июня 2010

Я немного подумал и получил следующий результат. Обратите внимание, что вы получаете как строки, так и столбцы за один выстрел.

Допущения: строки начинаются с 0. Столбцы начинаются с 0. Индекс начинается с 0

Обозначение

N = размер матрицы (был М в исходной задаче)

m = Индекс элемента

Код псевдо

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
}

И некоторый октавный / матричный код

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

0 голосов
/ 28 октября 2008

Мне потребовалось время, чтобы понять, что вам нужно! :)

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.
}
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...