最初の方法のベクトル化

Cell

昨日の謎の表は56行と余り8ワードあるので、56行を4分割してベクトル化する

ベクトルの最初のワードに、最初の14行を割り当てる

  //    0:  509-282- 55-452-225-622-395-168-565-338-111
  //    1:  508-281- 54-451-224-621-394-167-564-337-110
  //    2:  507-280- 53-450-223-620-393-166-563-336-109
  //   ...
  //   11:  498-271- 44-441-214-611-384-157-554-327-100
  //   12:  497-270- 43-440-213-610-383-156-553-326- 99
  //   13:  496-269- 42-439-212-609-382-155-552-325- 98

次のワードには、次の14行

  //   14:  495-268- 41-438-211-608-381-154-551-324- 97
  //   15:  494-267- 40-437-210-607-380-153-550-323- 96
  //   16:  493-266- 39-436-209-606-379-152-549-322- 95
  //   ...
  //   25:  484-257- 30-427-200-597-370-143-540-313- 86
  //   26:  483-256- 29-426-199-596-369-142-539-312- 85
  //   27:  482-255- 28-425-198-595-368-141-538-311- 84

というように分けると、こんな感じにベクトルにおさまる

vector x_509 = {mt[509],mt[495],mt[481],mt[467]};
vector x_508 = {mt[508],mt[494],mt[480],mt[466]};
vector x_281 = {mt[281],mt[267],mt[253],mt[239]};

で、シャッフルとかしなくても正しい場所に計算相手がはいってるから

y = spu_selb(x_509, x_508, mask);

みたいにロードしたらすぐ取り掛かれるわけだ。

even*1+odd*3

Cell 
(y >> 1) ^ mag01[y & 1]

をeven*1+odd*3でやる話。XORはevenでやらないといけないから右シフトは当然qword右シフトになる。
mag01[y & 1] は gb してテーブルルックアップすることにすると

z = si_lqx(spu_slqw(spu_gather(y), 4), mag_lut);
r = spu_xor(spu_rlmaskqw(y,-1), z);

qword右シフトだと上からよけいなビットが降ってくるけど、
gatherした4ビットの中に同じビットがあるからmag_lutを工夫すれば打ち消せる。

これで、even*1+odd*4になったんだけど、これからoddを一命令減らすとか無理っぽいけど
続きは次回。

擬似コード

Cell

謎の表の順番で更新すれば1ワード1ロードにできるのを擬似コードにしてみた。
ベクトル化はなし。

#define N 624
#define N 397

word mt[N];
word mtx[N];

void transform(word dst[N], word src[N])
{
  for (int i = 0; i < N; i++) {
    // 謎の表の通りに並び替える
    // dst[0] = src[509]
    // dst[1] = src[282]
    // ...
    dst[i] = src[(509 + i * M) % N];
  }
}

word f(word x0, word x1, word xM)
{
  word y = (x0&UPPER_MASK)|(x1&LOWER_MASK);
  return xM ^ (y >> 1) ^ mag01[y & 0x1UL];  
}

unsigned int
genrand_mine(int num_rand)
{
  transform(mtx, mt);

  for (num_rand) {
    // 最初の行の一つ上の行をロード
    word mt510 = mtx[(607 + 510 * 613) % N]; // mt[510]
    word mt283 = mtx[(607 + 283 * 613) % N]; // mt[283]
    word mt56  = mtx[(607 +  56 * 613) % N]; // mt[ 56]
    word mt453 = mtx[(607 + 453 * 613) % N]; // mt[453]
    word mt226 = mtx[(607 + 226 * 613) % N]; // mt[226]
    word mt623 = mtx[(607 + 623 * 613) % N]; // mt[623]
    word mt396 = mtx[(607 + 396 * 613) % N]; // mt[396]
    word mt169 = mtx[(607 + 169 * 613) % N]; // mt[169]
    word mt566 = mtx[(607 + 566 * 613) % N]; // mt[566]
    word mt339 = mtx[(607 + 339 * 613) % N]; // mt[339]
    word mt112 = mtx[(607 + 112 * 613) % N]; // mt[112]
    word mt509 = mtx[(607 + 509 * 613) % N]; // mt[509]
    for (int i = 0; i < 56; i++) {
      // i行目をロード
      word mt282 = mtx[11 * i + 1];
      word mt55  = mtx[11 * i + 2];
      word mt452 = mtx[11 * i + 3];
      word mt225 = mtx[11 * i + 4];
      word mt622 = mtx[11 * i + 5];
      word mt395 = mtx[11 * i + 6];
      word mt168 = mtx[11 * i + 7];
      word mt565 = mtx[11 * i + 8];
      word mt338 = mtx[11 * i + 9];
      word mt111 = mtx[11 * i + 10];
      word mt508 = mtx[11 * i + 11];

      // i行目を更新。順番に注意
      word mt55n  = f(mt55,  mt56,  mt452);
      word mt225n = f(mt225, mt226, mt622);
      word mt168n = f(mt168, mt169, mt565);
      word mt111n = f(mt111, mt112, mt508);
      
      word mt282n = f(mt282, mt283, mt55n);
      word mt452n = f(mt452, mt453, mt225n);
      word mt395n = f(mt395, mt396, mt168n);
      word mt338n = f(mt338, mt339, mt111n);

      word mt509n = f(mt509, mt510, mt282n);
      word mt622n = f(mt622, mt623, mt395n);
      word mt565n = f(mt565, mt566, mt338n);

      // 更新した結果を保存
      mtx[11 * i + 0] = mt509n;
      mtx[11 * i + 1] = mt282n;
      mtx[11 * i + 2] = mt55n;
      // 略
      mtx[11 * i + 10] = mt111n;

      // 一行ずらして、次の行の更新へ
      mt510 = mt509;
      mt283 = mt282;
      mt56  = mt55;
      // 略
      mt112 = mt111;
      mt509 = mt508;
    }
    // 残り8ワードの処理
    // ...
    // temperingと和
    // ...
  }
}

で、i=0..13, 14..27, 28..41, 42..56 は並列に実行できるからベクトル化は簡単にできるーね。