整数を対称行列のindexと対応付ける方法

Question

自分で実装してみたアルゴリズムがあるのですが，もっと良い方法があれば教えてほしいです．以下はC++で実装してますが，PythonやFortranでの回答でも構いません．

N×N の対称行列のインデックスを(i,j)と書きます．i<N, j<N です．
ここで，適当な整数の入力ele ( <= (N-1)(N+2)/2 ) に対して，対称行列のi <= j を満たすインデックス(i,j)を対応させます．

たとえば N = 4 で入力が ele=3 の場合は，

[0,,,*]

[1,2,,]
[3,4,5,*]
[6,7,8,9]

なので，出力は(2,0)，同じく，N = 4 で ele = 8 ならば，出力は(3,2) といった具体です．

自分が実装したコードは以下です．高校数学の群数列のアイディアを思い出して，n行目の最小がn*(n+3)/2 - n で最大が n*(n+3)/2 であることを使ってみました．皆さんだったらどう実装しますか？

tuple<int, int> get_mtx_idx(int ele, int N){

    /*
     * ele map to index (i, j) of symmetric matrix
     * [[ 0, *, *, *, *, *, *]
     *  [ 1, 2, *, *, *, *, *]
     *  [ 3, 4, 5, *, *, *, *]
     *  [ 6, 7, 8, 9, *, *, *]
     *  [10,11,12,13,14, *, *]
     *  [15,16,17,18,19,20, *]
     *  [21,22,23,24,25,26,27]
     *
     *  min in n-th row is n*(n+3)/2 - n
     *  max in n-th row is n*(n+3)/2
     *
     * */

    for (int n = 0 ; n <= N; n++){
        if ( n * ( n+3 ) /2 - n <= ele and ele <= n * ( n+3 ) /2 ){
            int i = n;
            int j = ele - n * ( n+3 )/2 + n;

            return{i, j};

        }
    }

}

Answer 1

こんにちは
Pythonを使いました。

import math

def get_mtx_idx(ele):
  i = int(math.sqrt(2*ele+0.25) - 0.5)
  j = ele - int(i*(i+1)/2)
  return i, j

以下は上記を使った動作確認です。

for ele in range(21):
    print(f'{ele}: {get_mtx_idx(ele)}')

出力結果：

0: (0, 0)

1: (1, 0)
2: (1, 1)
3: (2, 0)
4: (2, 1)
5: (2, 2)
6: (3, 0)
7: (3, 1)
8: (3, 2)
9: (3, 3)
10: (4, 0)
11: (4, 1)
12: (4, 2)
13: (4, 3)
14: (4, 4)
15: (5, 0)
16: (5, 1)
17: (5, 2)
18: (5, 3)
19: (5, 4)
20: (5, 5)

以下で確認できます。

• 動作確認用サンプル: https://repl.it/@jun68ykt/Q293964

参考になれば幸いです。

追記

コメントからいただきました

ただ，もしよろしければ，int( math.sqrt( 2ele+0.25 ) - 0.5 ) がどこからやってきたか教えてもらうことは可能でしょうか？

に回答します。

なお、以下の説明で、プログラムのコードではなく数式の中では、乗算記号を省略します。また、i^2 は i の二乗を表すものとします。

まず、i （≧0）行目の先頭（0列目)の要素は、

i (i+1) / 2

です。たとえば、 i = 3 の行の先頭要素は、上記の式のiに3を代入した 3 (3 + 1) / 2 から得られる 6です。上記の先頭要素の式を使うと、与えられた ele が配置される行インデクスi は、以下の不等式

i (i+1) / 2 ≦ ele

を満たす最大の非負整数です。この不等式を変形して、左辺が i になるようにしていきます。まず、両辺を２倍し、左辺のカッコを展開します。

i^2 + i ≦ 2ele

両辺に 0.25 を加えて、

i^2 + i + 0.25 ≦ 2ele + 0.25

としておいて、左辺を因数分解すると

(i + 0.5)^2 ≦ 2ele + 0.25

となります。両辺の値は0より大きく、i + 0.5も0より大きいので、両辺の平方根をとると以下が得られます。

i + 0.5  ≦ √(2ele + 0.25)  

左辺の 0.5 を右辺に移項すると、左辺がi の不等式

i ≦ √(2ele + 0.25) - 0.5

が得られました。これを満たす最大の非負整数が求めたいiになります。それは、右辺の値の小数部分を切り捨てれば得られますので、i を求める式をPythonで書くと

i = int(math.sqrt(2*ele+0.25) - 0.5)

となります。i が得られると、列インデクスj は、ele と 行の先頭要素i (i+1) / 2 との差なので

j = ele - int(i*(i+1)/2)

で得ることができます。

Comments

[physics303]

ありがとうございます！！
とてもシンプルで分かりやすいです！

ただ，もしよろしければ，int( math.sqrt( 2ele+0.25 ) - 0.5 ) がどこからやってきたか教えてもらうことは可能でしょうか？

[jun68ykt]

どういたしまして。

> int( math.sqrt( 2ele+0.25 ) - 0.5 )

の導き方を回答に追記しました。

[physics303]

おおおお，大変よく分かりました！ありがとうございます！

[jun68ykt]

????

Answer 2

やってみたけど、アルゴリズムのキモは同じやね。
※ 時間計算量はO(√ele)

#include <iostream>
#include <tuple>

std::tuple<int, int> get_mtx_idx(int ele, int N){
  int delta = 0;
  int sum = 0;
  while ( ele >= sum ) {
    sum += ++delta;
  }
  return { delta - 1, ele - sum + delta };
}

// おためし
int main() {
  const int N = 4;
  for ( int i = 0; i < 10; ++i ) {
    auto idx = get_mtx_idx(i, 4);
    std::cout << i << " (" << std::get<0>(idx) << ',' << std::get<1>(idx) << ")\n";
  }
}

...あ、引数Nは要らんわ。