サポートしているSIMD拡張を調べるプログラム
Intel AVX について勉強しているが、使っているCPUがどのSIMD拡張をサポートしているのかを知りたくなった。なので色々調べてインラインアセンブラを使ったコードを書いてみた。Sandy Bridge (Core i7 2600K)でAVXが使えることは確認できた。
// sup_ext.c #include <stdio.h> int main() { int REG_edx, REG_ecx, REG_eax=-1; __asm__("mov eax, 1\n\t" "cpuid\n\t" "mov %0, edx\n\t" "mov %1, ecx\n\t" : "=r" (REG_edx), "=r" (REG_ecx) ); if(REG_edx && 0x800000) printf("MMX "); // EDX bit 23 if(REG_edx && 0x2000000) printf("SSE "); // EDX bit 25 if(REG_edx && 0x4000000) printf("SSE2 "); // EDX bit 26 if(REG_ecx && 0x00000001) printf("SSE3 "); // ECX bit 0 if(REG_ecx && 0x00000200) printf("SSSE3 "); // ECX bit 9 if(REG_ecx && 0x00080000) printf("SSE4.1 "); // ECX bit 19 if(REG_ecx && 0x00100000) printf("SSE4.2 "); // ECX bit 20 if(REG_ecx && 0x00200000) printf("AESNI "); // ECX bit 25 if(REG_ecx && 0x00000002) printf("PCLMULQDQ "); // ECX bit 1 __asm__("mov eax, 1\n\t" "cpuid\n\t" "and ecx, 402653184\n\t" // 0x01800000 = 402653184 "cmp ecx, 402653184\n\t" "jne not_supported\n\t" "mov ecx, 0\n\t" "xgetbv\n\t" "and eax, 6\n\t" "cmp eax, 6\n\t" "mov %0, eax\n" "not_supported:\n\t" : "=r" (REG_eax) ); if (REG_eax == 0x6) printf("AVX "); // EAX bit 1&2 printf("\n"); return 0; }
最短経路の距離行列から経路行列を構築するアルゴリズム
全点対最短経路問題(All-Pairs Shortest Paths Problem)において、Floyd-Warshall法などで最短経路の距離行列が求められている場合に、その最短距離行列から対応する経路行列を構築するアルゴリズムについてのメモ。アルゴリズムはO(n^3)のアルゴリズムで単純な実装ではFloyd-Warshall法と同程度の計算処理時間がかかる。以下はCによる実装例。
void construct_aps_path(const int n, const int **A, const int **D, int **P) { for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { for (int k = 0; k < n; k++) { if (D[i][j] == D[i][k]+A[k][j]) { P[i*n+j] = k; break; } if (k != n-1) { P[i*n+j] = NIL; } } } } }
nが問題サイズ、Aが入力距離行列、Dが最短経路の距離行列、Pは求める最短経路のpredecessor行列。初期条件等はCLRSの"Introduction to Algorithms"の"Constructing a shortest path"の項を参照(この問題はCLRSの演習問題にもなっている)。入力グラフによっては、Floyd-Warshall法の計算途上で経路行列も記録させた場合と結果が違うものになることもある。
Sun Ultra24 Workstation での GSL の CBLAS の SGEMM の性能
Sun Ultra 24 Workstation (プロセッサは Intel Core2Duo E8600 (3.3GHz)) で GSL (GNU Scientific Library) の CBLAS の SGEMM の性能を測った。
結果は以下図(行列サイズは16の倍数)。僕が実装したスカラーとブロックの行列乗算アルゴリズムと比較した(ブロックのは、8x8でブロック化)。時間があったら他のBLASも試す。
中置記法を後置記法に変換するプログラム
https://www.spoj.pl/problems/ONP/を解くために書いたRubyのコードを、後々使えそうだからポストする。
$priority = [[['+',0],['-',1],['*',2],['/',3],['^',4]]] def recur(s, idx) op_stack = Array.new p = idx while p < s.length c = s[p].chr case c when 'a'..'z' print c when '(' p = recur(s, p+1) when ')' until op_stack.empty? print op_stack.pop end return p else while !op_stack.empty? && $priority[op_stack.last] > $priority[c] print op_stack.pop end op_stack.push c end p += 1 end until op_stack.empty? print op_stack.pop end end T = gets.to_i T.times do s = gets.chomp recur(s, 0) puts end
Google Code Jam 2009 - Round 1
Round 1 落ちです。悔しいな。去年よりは実力がついたとは思うが、それでもまだまだ対応できる問題の種類が少なく、レベルが低い。
Google Code Jam 2009 - Qualification Round
Code Jam - Google’s Coding Competitions
参戦記録として書いておく。
結果は約1時間半で全問解けて、65位。それぞれ30分ぐらいずつかかった。出だし好調と言える結果だが、全問解いた方は約2400人もいる。Round 1の合計通過人数が3000人なので、Round 1を通過するだけでも大変そうだ。
以下、適当な解説と書いたプログラムのソースコード。
Problem A. Alien Language
宇宙人のメッセージを解読。簡単な正規表現。実際、正規表現をサポートしている言語を使って問題を解いた参加者もいたみたいだ。僕は、メッセージの一文字ごとの候補を作り、それをwordの文字と比較して、表現できるかどうかを調べた。
#include <iostream> #include <set> #include <vector> #include <string> #include <sstream> #include <cstdio> using namespace std; int main(void) { int L, D, N; string line; getline(cin, line); istringstream iss(line); iss >> L >> D >> N; set<string> words; while (D-- > 0) { getline(cin, line); words.insert(line); } for (int caseNo = 1; caseNo <= N; caseNo++) { getline(cin, line); istringstream iss2(line); vector<string> token(L, ""); bool isInPara = false; char c; int i = 0; while (iss2 >> c) { if (c == '(') { isInPara = true; } else if (c == ')') { isInPara = false; i++; } else if (isInPara) { token[i] += string(1,c); } else { token[i] = string(1,c); i++; } } int K = 0; for (set<string>::const_iterator wd = words.begin(); wd != words.end(); wd++) { for (int i = 0; i < L; i++) { if (token[i].find((*wd)[i]) == string::npos) break; if (i == L-1) { K++; break; } } } printf("Case #%d: %d\n", caseNo, K); } return 0; }
Problem B. Watersheds
どこに水がたまるか。シミュレーション。再帰で解いた。1回調べたマスはもう調べないようにして計算量を削減した。
#include <iostream> #include <vector> #include <cstdio> using namespace std; int H, W; vector<vector<int> > basin; vector<vector<char> > sink; char go(const int y, const int x, const char label) { static int d[4][2] = { {-1, 0}, { 0,-1}, { 0, 1}, { 1, 0} }; if (sink[y][x] != '?') return sink[y][x]; int lowest = -1; int lowestApptitude = INT_MAX; for (int i = 0; i < 4; i++) { const int dy = y + d[i][0]; const int dx = x + d[i][1]; if (0<=dy&&dy<H && 0<=dx&&dx<W && basin[dy][dx] < lowestApptitude && basin[dy][dx] < basin[y][x]) { lowest = i; lowestApptitude = basin[dy][dx]; } } if (lowest == -1) { sink[y][x] = label; } else { sink[y][x] = go(y+d[lowest][0],x+d[lowest][1],label); } return sink[y][x]; } int main(void) { int T; cin >> T; for (int caseNo = 1; caseNo <= T; caseNo++) { cin >> H >> W; basin = vector<vector<int> >(H, vector<int>(W)); sink = vector<vector<char> >(H, vector<char>(W,'?')); for (int y = 0; y < H; y++) for (int x = 0; x < W; x++) cin >> basin[y][x]; char label = 'a'; for (int y = 0; y < H; y++) for (int x = 0; x < W; x++) if (sink[y][x] == '?') { if (go(y, x, label) == label) label++; } printf("Case #%d:\n", caseNo); for (int y = 0; y < H; y++) for (int x = 0; x < W; x++) { cout << sink[y][x]; if (x == W-1) cout << endl; else cout << " "; } } return 0; }
Problem C. Welcome to Code Jam
文字列中に"welcome to code jam"が何個隠れているかを数える。入力例があまり親切とは言えなく、どういう問題だと解釈すればいいのか迷った。通ったので解釈はあっていたのだと思う。だけど、このプログラム、余計な処理を行っている気がする。気が向いたら削れるかどうかを考える。
#include <iostream> #include <string> #include <sstream> #include <vector> using namespace std; int main(void) { string line; getline(cin, line); int N; istringstream(line) >> N; const string message = "welcome to code jam"; const int mesLen = message.length(); for (int caseNo = 1; caseNo <= N; caseNo++) { int counter = 0; getline(cin, line); const int lineLen = line.length(); vector<vector<int> > memo(lineLen,vector<int>(mesLen,0)); for (int i = 0; i < lineLen; i++) if (line[i] == message[0]) memo[i][0] = 1; for (int i = 0; i < lineLen; i++) { for (int j = mesLen-1; j > 0; j--) { if (line[i] != message[j]) continue; for (int k = 0; k < i; k++) memo[i][j] += memo[k][j-1]; memo[i][j] %= 10000; } } for (int i = 0; i < lineLen; i++) { counter += memo[i][mesLen-1]; counter %= 10000; } printf("Case #%d: %04d\n", caseNo, counter); } return 0; }