2009年3月28日土曜日

ErlangでFizzBuzz問題、それをC++のテンプレートで表現する

このブログ記事をはてなブックマークに追加

今更ながらFizzBuzz問題を解いてみる。今回はErlangで解いてみた。それだけだとつまらないので、それをC++のテンプレートで表現してみた。因みにErlangを使っているけど並列処理はしていない。

まずは普通に何の工夫もなく書いてみた。

-module(fizzbuzz). -export([fb/1]). fb(0) -> ok; fb(N) -> fb(N - 1), if N rem 15 =:= 0 -> io:put_chars("FizzBuzz"); N rem 3 =:= 0 -> io:put_chars("Fizz"); N rem 5 =:= 0 -> io:put_chars("Buzz"); true -> io:write(N) end, io:nl().

うーん、面白くない。どうせなら if をなくしてしまおう。というわけで、今度は以下のように書いた。

-module(fizzbuzz). -export([fb/3, fizzbuzz/1]). fb(_, 0, 0) -> io:put_chars("FizzBuzz\n"); fb(_, 0, _) -> io:put_chars("Fizz\n"); fb(_, _, 0) -> io:put_chars("Buzz\n"); fb(N, _, _) -> io:write(N), io:nl(). fizzbuzz(0) -> ok; fizzbuzz(N) -> fizzbuzz(N - 1), fb(N, N rem 3, N rem 5).

少しは関数型プログラミングっぽいかなぁ? これをC++のテンプレートで書いてみると以下のようになる。

#include <iostream> template<int N, int M3, int M5> struct FB { FB() { std::cout << N << std::endl; } }; template<int N, int M5> struct FB<N, 0, M5> { FB() { std::cout << "Fizz" << std::endl; } }; template<int N, int M3> struct FB<N, M3, 0> { FB() { std::cout << "Buzz" << std::endl; } }; template<int N> struct FB<N, 0, 0> { FB() { std::cout << "FizzBuzz" << std::endl; } }; template<int N> struct FizzBuzz { FizzBuzz() { FizzBuzz<N-1>(); FB<N, N%3, N%5>(); } }; template<> struct FizzBuzz<0> {}; int main() { FizzBuzz<100>(); return 0; }

どうせなら2つある関数(C++ならクラス)を1つにしてみよう。

-module(fizzbuzz). -export([fb/1]). fb(0, 0, 0) -> ok; fb(N, 0, 0) -> fb(N-1, (N-1) rem 3, (N-1) rem 5), io:put_chars("FizzBuzz\n"); fb(N, 0, _) -> fb(N-1, (N-1) rem 3, (N-1) rem 5), io:put_chars("Fizz\n"); fb(N, _, 0) -> fb(N-1, (N-1) rem 3, (N-1) rem 5), io:put_chars("Buzz\n"); fb(N, _, _) -> fb(N-1, (N-1) rem 3, (N-1) rem 5), io:write(N), io:nl(). fb(N) -> fb(N, N rem 3, N rem 5).

C++のテンプレートでは以下の通り。

#include <iostream> template<int N, int M3=N%3, int M5=N%5> struct FB { FB() { FB<N-1, (N-1)%3, (N-1)%5>(); std::cout << N << std::endl; } }; template<int N, int M5> struct FB<N, 0, M5> { FB() { FB<N-1, (N-1)%3, (N-1)%5>(); std::cout << "Fizz" << std::endl; } }; template<int N, int M3> struct FB<N, M3, 0> { FB() { FB<N-1, (N-1)%3, (N-1)%5>(); std::cout << "Buzz" << std::endl; } }; template<int N> struct FB<N, 0, 0> { FB() { FB<N-1, (N-1)%3, (N-1)%5>(); std::cout << "FizzBuzz" << std::endl; } }; template<> struct FB<0, 0, 0> {}; int main() { FB<100>(); return 0; }

何だかちょっとやり過ぎな気がするなぁ。

続きを読む...

2009年3月25日水曜日

Google App Engine: Googleはてブを作ってみた

このブログ記事をはてなブックマークに追加

Googleの検索結果にリンク付きのはてなブックマーク件数を表示できれば便利だと思い、Google App Engineを使ってウェブアプリケーションとして実装してみた。

最初、はてなブックマーク件数取得APIを使ってみたが、非常に遅い。これはXML-RPC APIによる実装だが、件数を取得するまでブラウザに表示できないのでかなり待たされる。酷いときにはGoogle App Engineの制限によりタイムアウトしてしまうこともあった。仕方がないので素直にhttp://b.hatena.ne.jp/entry/を使って画像表示を行うことにした。


今時のブラウザは遅延読み込みによる画像表示のため待たされずにサクサク表示できる。


一方、API版は画面に表示されるまでかなり待たされる。

コードについては、あまり時間をかけずにザザッと書いたので結構いい加減かも。日本語処理で少しはまって適当に処置した。最後にソースコードを示しておく。

画像表示版:

class GoogleHatenaImage(webapp.RequestHandler): def __init__(self): self.google_url = "http://www.google.co.jp/" self.google_search_url = "http://www.google.co.jp/search" def get(self): uri_data = urlparse.urlsplit(self.request.uri) baseuri = uri_data[0] + "://" + "".join(uri_data[1:3]) url = self.request.get("content").encode("utf-8") query = self.request.query_string.replace("Shift_JIS", "UTF-8") if not url: result = urlfetch.fetch("%s?%s" % (self.google_search_url, query)) else: result = urlfetch.fetch(query.replace("content=", "")) if result.status_code == 200: content = result.content.decode("cp932") content = re.sub("charset=Shift_JIS", "charset=UTF-8", content) p = re.compile("a href=[\'\"]\/(.+?)[\'\"]") p2 = re.compile("<a href=[\'\"]http(.+?)[\'\"](.+?)\/a>") p3 = re.compile("action=\"\/search") content = p.sub("a href=\"%s?content=%s%s\"" % (baseuri, self.google_url, "\g<1>"), content) content = p2.sub("<a href=\"http\g<1>\"\g<2>/a><a href=\'http://b.hatena.ne.jp/entry/http\g<1>\'><img src=\'http://b.hatena.ne.jp/entry/image/http\g<1>\' style=\'border: none;\'/></a>", content) content = p3.sub("action=\"%s" % baseuri, content) content = content.replace("</head>", "</head><base href=%s />" % self.google_url) self.response.out.write(content)

XML-RPC API版:

class GoogleHatenaXmlrpc(webapp.RequestHandler): def __init__(self): self.s = xmlrpclib.ServerProxy("http://b.hatena.ne.jp/xmlrpc") self.google_url = "http://www.google.co.jp/" self.google_search_url = "http://www.google.co.jp/search" def insert_hatena(self, matchobj): cnt = self.s.bookmark.getCount("http" + matchobj.group(1)).values()[0] if cnt > 0: us = "users" if cnt > 1 else "user" style_format = "style='font-size: 0.9em; text-decoration: underline; white-space: pre; background-color: #fff0f0; font-weight: bold; color: #f06060;'" if cnt >= 10: style_format = "style='font-size: 0.9em; text-decolation: underline; white-space: pre; background-color: #ffcccc; font-weight: bold; color: red;'" return "<a href=\"http%s\"%s/a><a %s href=\'http://b.hatena.ne.jp/entry/http%s\'>%d %s</a>" % (matchobj.group(1), matchobj.group(2), style_format, matchobj.group(1), cnt, us) else: return "<a href=\"http%s\"%s/a>" % (matchobj.group(1), matchobj.group(2)) def get(self): uri_data = urlparse.urlsplit(self.request.uri) baseuri = uri_data[0] + "://" + "".join(uri_data[1:3]) url = self.request.get("content").encode("utf-8") query = self.request.query_string.replace("Shift_JIS", "UTF-8") if not url: result = urlfetch.fetch("%s?%s" % (self.google_search_url, query)) else: result = urlfetch.fetch(query.replace("content=", "")) if result.status_code == 200: content = result.content.decode("cp932") content = re.sub("charset=Shift_JIS", "charset=UTF-8", content) p = re.compile("a href=[\'\"]\/(.+?)[\'\"]") p2 = re.compile("<a href=[\'\"]http(.+?)[\'\"](.+?)\/a>") p3 = re.compile("action=\"\/search") content = p.sub("a href=\"%s?content=%s%s\"" % (baseuri, self.google_url, "\g<1>"), content) content = p2.sub(self.insert_hatena, content) content = p3.sub("action=\"%s" % baseuri, content) content = content.replace("</head>", "</head><base href=%s />" % self.google_url) self.response.out.write(content)

追記(2009/3/26):

日本語処理を若干修正。

続きを読む...

2009年3月23日月曜日

はじめてのCUDAプログラミングで分子動力学計算

このブログ記事をはてなブックマークに追加

追記(2009/12/12): ここで示してるコードはかなりひどいので少し手直しして、CUDAで作成した分子動力学計算プログラムを書き直してみたという記事にまとめた。こちらも参照して欲しい。

最近、自分の周りではCUDAを使ったGPGPUが流行っているので、どんなものかこの週末にいじってみた。

一日目はNVIDIA GeForce 8800 GTSを搭載したWindows XPマシンにCUDAをインストールをしてリファレンスやサンプルプログラムをパラパラと読んでみた。二日目になんちゃって分子動力学(MD)計算プログラムを書いた。擬似原子を扱っており、全原子のvdW力とクーロン力のみを考慮して分子内結合については考慮しない、液滴中のMD計算だ。取り敢えず、半径20Åの液滴中に1229原子を入れて1000ステップ計算させてみた。C++で書いたプログラムでは4分30秒かかった計算がCUDAでは32秒で計算できた。約9倍ほどの速度が出ているようだ。

まあ、C++のコードも併せて一日で書いたので最適化や高速化については考慮していない。後でスレッドでの分割がやりやすいかもと思ってiとjをそのまま全部計算しているし。気が向いたらちゃんとしたプログラムを書いてみるかも。

CUDAアーキテクチャについてはまだまだ理解が足りない。効率のよいメモリの扱い方を覚えないとなぁ。あと、sharedメモリをうまく扱えてないっぽい。さすがに一日二日じゃキツイか。

初心者が一日で書いたコードなので参考になるかわからないが、以下にソースコードとサンプルの入力データを示しておく。入力データについては、上記で用いた1229原子ではブログ記事するには大きすぎるので154原子の小さなデータを用意した。また、記事の表示量が多くなりすぎるのでC++のソースは割愛させてもらった。Windows上でのコンパイルは下記の通り。

nvcc -lcutil32 -lkernel32 --host-compilation c++ -Xcompiler /EHsc,/W3,/nologo,/O2,/Zi,/MT simple_md_gpu.cu -o simple_md_gpu.exe

入力ファイルをmd.datとして、100ステップ計算させる場合は、以下のように実行する。

simple_md_gpu md.dat 100

最初に初期構造が表示され、次いでステップごとのエネルギー、最後に最終構造が出力される。

simple_md_gpu.cu

//////////////////////////////////////////////////////////////// // Simple Molecular Dynamics using GPU by nox, 2009.03.22. // // Perform molecular dynamics for pseudo-atoms without internal // interactions, such as bonds, angles, dihedrals. Only vdW and // Coulomb interactions. #include <iostream> #include <iomanip> #include <fstream> #include <sstream> #include <cmath> #include <cutil_inline.h> using namespace std; const int BLOCK = 256; const int WIDTH = 2048; __device__ float Distance(float* crd, int i, int j) { float dx = crd[i * 3] - crd[j * 3]; float dy = crd[i * 3 + 1] - crd[j * 3 + 1]; float dz = crd[i * 3 + 2] - crd[j * 3 + 2]; return sqrtf(dx * dx + dy * dy + dz * dz); } __device__ float Energy(float* ene, int num_atoms) { float sum = 0.0f; for (int i = 0; i < num_atoms; i++) sum += ene[i]; return sum / 2.0f; } __global__ void Center(float* crd,int num_atoms) { __shared__ float d[3]; for (int i = 0; i < 3; i++) d[i] = 0.0f; for (int i = 0; i < num_atoms; i++) for (int j = 0; j < 3; j++) d[j] += crd[i * 3 + j]; for (int i = 0; i < 3; i++) d[i] /= num_atoms; for (int i = 0; i < num_atoms; i++) for (int j = 0; j < 3; j++) crd[i * 3 + j] -= d[j]; } // Calculate energies and forces // Energy = vdW + Coulomb // vdW // U = eps * [(Rij / r[i])^12 - 2 * (Rij / r[i])^6] // F = -12 * eps / Rij * [(Rij / r[i])^13 - (Rij / r[i])^7] * r_xyz / r[i] // Coulomb // U = SUM_i>j qiqj / r[i] // F = SUM_j qiqj / r[i]^3 * r_xyz __global__ void Calc(float* crd, float* f, float* ene, float* chg, int num_atoms, float* tene) { int dtx = blockIdx.x * blockDim.x + threadIdx.x; const float eps = 0.2f; const float Rij = 2.5f; float r, r_xyz, df, r12, r6, q; for (int i = dtx; i < num_atoms; i += WIDTH) ene[i] = 0.0f; __syncthreads(); for (int i = dtx; i < num_atoms; i += WIDTH) { for (int j = 0; j < num_atoms; j++) { if (i == j) continue; r = Distance(crd, i, j); q = chg[i] * chg[j]; r12 = powf(Rij / r, 12); r6 = powf(Rij / r, 6); ene[i] += eps * (r12 - 2.0f * r6) + q / r; for (int k = 0; k < 3; k++) { r_xyz = crd[i * 3 + k] - crd[j * 3 + k]; df = (12 * eps / Rij) * (r12 * Rij / r - r6 * Rij / r) * r_xyz / r; df += q / powf(r, 3) * r_xyz; f[i * 3 + k] += df; } } } __syncthreads(); *tene = Energy(ene, num_atoms); } __global__ void Move(float* crd, float* f, int num_atoms, float cap_range) { int dtx = blockIdx.x * blockDim.x + threadIdx.x; float r, dr; for (int i = dtx; i < num_atoms; i += WIDTH) { for (int j = 0; j < 3; j++) { crd[i * 3 + j] += f[i * 3 + j] * 0.01f; f[i * 3 + j] = 0.0f; } r = 0.0f; for (int j = 0; j < 3; j++) r += crd[i * 3 + j] * crd[i * 3 + j]; dr = r - cap_range * cap_range; if (dr > 0.0f) for (int j = 0; j < 3; j++) f[i * 3 + j] -= crd[i * 3 + j] / cap_range * dr * 0.01f; } } class SimpleMD { private: int num_steps; int num_atoms; float cap_range; float *h_crd, *h_f, *h_ene, *h_chg; float *d_crd, *d_f, *d_ene, *d_chg; float tene; public: SimpleMD(int n, char*); ~SimpleMD(); void ReadCrd(char*); void PrintCrd(); void Run(); }; SimpleMD::SimpleMD(int n, char* fname) : num_steps(n) { ReadCrd(fname); h_f = new float[num_atoms * 3]; h_ene = new float[num_atoms]; for (int i = 0; i < num_atoms; i++) { for (int j = 0; j < 3; j++) h_f[i * 3 + j] = 0.0f; h_ene[i] = 0.0f; } cudaMalloc((void**)&d_crd, sizeof(float) * num_atoms * 3); cudaMalloc((void**)&d_f, sizeof(float) * num_atoms * 3); cudaMalloc((void**)&d_ene, sizeof(float) * num_atoms); cudaMalloc((void**)&d_chg, sizeof(float) * num_atoms); } SimpleMD::~SimpleMD() { cudaFree(d_chg); cudaFree(d_ene); cudaFree(d_f); cudaFree(d_crd); delete[] h_ene; delete[] h_f; delete[] h_chg; delete[] h_crd; } void SimpleMD::ReadCrd(char* fname) { fstream fs(fname, ios_base::in); string line; stringstream ss; if (!fs.is_open()) { cerr << "File open error: " << fname << endl; exit(1); } getline(fs, line); cout << line << endl; getline(fs, line); ss.str(line); ss >> num_atoms; ss.clear(); getline(fs, line); ss.str(line); ss >> cap_range; ss.clear(); h_crd = new float[num_atoms * 3]; h_chg = new float[num_atoms]; for (int i = 0; i < num_atoms; i++) { getline(fs, line); ss.str(line); ss >> h_crd[i * 3] >> h_crd[i * 3 + 1] >> h_crd[i * 3 + 2] >> h_chg[i]; ss.clear(); ss.str(""); } fs.close(); } void SimpleMD::PrintCrd() { cout << endl; cout << num_atoms << endl; cout << cap_range << endl; for (int i = 0; i < num_atoms; i++) { for (int j = 0; j < 3; j++) cout << " " << fixed << setw(10) << setprecision(6) << h_crd[i * 3 + j]; cout << " " << fixed << setw(8) << setprecision(4) << h_chg[i]; cout << endl; } } void SimpleMD::Run() { float *d_tene; cudaMalloc((void**)&d_tene, sizeof(float)); cudaMemcpy(d_crd, h_crd, sizeof(float) * num_atoms * 3, cudaMemcpyHostToDevice); cudaMemcpy(d_f, h_f, sizeof(float) * num_atoms * 3, cudaMemcpyHostToDevice); cudaMemcpy(d_ene, h_ene, sizeof(float) * num_atoms, cudaMemcpyHostToDevice); cudaMemcpy(d_chg, h_chg, sizeof(float) * num_atoms, cudaMemcpyHostToDevice); dim3 grid(WIDTH / BLOCK, 1, 1); dim3 threads(BLOCK, 1, 1); Center<<<1, 1>>>(d_crd, num_atoms); for (int i = 0; i < num_steps; i++) { Calc<<<grid, threads>>>(d_crd, d_f, d_ene, d_chg, num_atoms, d_tene); Move<<<grid, threads>>>(d_crd, d_f, num_atoms, cap_range); Center<<<1, 1>>>(d_crd, num_atoms); cudaMemcpy(&tene, d_tene, sizeof(float), cudaMemcpyDeviceToHost); cout << "Energy (" << setw(7) << i + 1 << "): "; cout << fixed << setw(15) << setprecision(5) << tene << endl; } cudaMemcpy(h_crd, d_crd, sizeof(float) * num_atoms * 3, cudaMemcpyDeviceToHost); cudaMemcpy(h_f, d_f, sizeof(float) * num_atoms * 3, cudaMemcpyDeviceToHost); cudaMemcpy(h_ene, d_ene, sizeof(float) * num_atoms, cudaMemcpyDeviceToHost); cudaFree(d_tene); } int main(int argc, char** argv) { if (argc != 3) { cerr << "Usage: " << argv[0] << " input_file number_of_steps" << endl; cerr << "Input: line 1 : title" << endl; cerr << " line 2 : number of atoms" << endl; cerr << " line 3 : radius of droplet" << endl; cerr << " line 4-: x-crd y-crd z-crd charge" << endl; exit(1); } stringstream ss; int nstep; cutilDeviceInit(1, argv); ss.str(argv[2]); ss >> nstep; SimpleMD md(nstep, argv[1]); md.PrintCrd(); md.Run(); md.PrintCrd(); return 0; }

サンプル入力ファイル: md.dat

simple_md_gpu sample data 154 10.0 -7.00000 -7.00000 -1.00000 -1.000 -7.00000 -4.00000 -4.00000 1.000 -7.00000 -4.00000 -1.00000 -1.000 -7.00000 -4.00000 2.00000 1.000 -7.00000 -4.00000 5.00000 -1.000 -7.00000 -1.00000 -7.00000 1.000 -7.00000 -1.00000 -4.00000 -1.000 -7.00000 -1.00000 -1.00000 1.000 -7.00000 -1.00000 2.00000 -1.000 -7.00000 -1.00000 5.00000 1.000 -7.00000 2.00000 -4.00000 -1.000 -7.00000 2.00000 -1.00000 1.000 -7.00000 2.00000 2.00000 -1.000 -7.00000 2.00000 5.00000 1.000 -7.00000 5.00000 -4.00000 -1.000 -7.00000 5.00000 -1.00000 1.000 -7.00000 5.00000 2.00000 -1.000 -7.00000 5.00000 5.00000 1.000 -4.00000 -7.00000 -4.00000 -1.000 -4.00000 -7.00000 -1.00000 1.000 -4.00000 -7.00000 2.00000 -1.000 -4.00000 -7.00000 5.00000 1.000 -4.00000 -4.00000 -7.00000 -1.000 -4.00000 -4.00000 -4.00000 1.000 -4.00000 -4.00000 -1.00000 -1.000 -4.00000 -4.00000 2.00000 1.000 -4.00000 -4.00000 5.00000 -1.000 -4.00000 -4.00000 8.00000 1.000 -4.00000 -1.00000 -7.00000 -1.000 -4.00000 -1.00000 -4.00000 1.000 -4.00000 -1.00000 -1.00000 -1.000 -4.00000 -1.00000 2.00000 1.000 -4.00000 -1.00000 5.00000 -1.000 -4.00000 -1.00000 8.00000 1.000 -4.00000 2.00000 -7.00000 -1.000 -4.00000 2.00000 -4.00000 1.000 -4.00000 2.00000 -1.00000 -1.000 -4.00000 2.00000 2.00000 1.000 -4.00000 2.00000 5.00000 -1.000 -4.00000 2.00000 8.00000 1.000 -4.00000 5.00000 -7.00000 -1.000 -4.00000 5.00000 -4.00000 1.000 -4.00000 5.00000 -1.00000 -1.000 -4.00000 5.00000 2.00000 1.000 -4.00000 5.00000 5.00000 -1.000 -4.00000 8.00000 -4.00000 1.000 -4.00000 8.00000 -1.00000 -1.000 -4.00000 8.00000 2.00000 1.000 -1.00000 -7.00000 -7.00000 -1.000 -1.00000 -7.00000 -4.00000 1.000 -1.00000 -7.00000 -1.00000 -1.000 -1.00000 -7.00000 2.00000 1.000 -1.00000 -7.00000 5.00000 -1.000 -1.00000 -4.00000 -7.00000 1.000 -1.00000 -4.00000 -4.00000 -1.000 -1.00000 -4.00000 -1.00000 1.000 -1.00000 -4.00000 2.00000 -1.000 -1.00000 -4.00000 5.00000 1.000 -1.00000 -4.00000 8.00000 -1.000 -1.00000 -1.00000 -7.00000 1.000 -1.00000 -1.00000 -4.00000 -1.000 -1.00000 -1.00000 -1.00000 1.000 -1.00000 -1.00000 2.00000 -1.000 -1.00000 -1.00000 5.00000 1.000 -1.00000 -1.00000 8.00000 -1.000 -1.00000 2.00000 -7.00000 1.000 -1.00000 2.00000 -4.00000 -1.000 -1.00000 2.00000 -1.00000 1.000 -1.00000 2.00000 2.00000 -1.000 -1.00000 2.00000 5.00000 1.000 -1.00000 2.00000 8.00000 -1.000 -1.00000 5.00000 -7.00000 1.000 -1.00000 5.00000 -4.00000 -1.000 -1.00000 5.00000 -1.00000 1.000 -1.00000 5.00000 2.00000 -1.000 -1.00000 5.00000 5.00000 1.000 -1.00000 5.00000 8.00000 -1.000 -1.00000 8.00000 -4.00000 1.000 -1.00000 8.00000 -1.00000 -1.000 -1.00000 8.00000 2.00000 1.000 -1.00000 8.00000 5.00000 -1.000 2.00000 -7.00000 -4.00000 1.000 2.00000 -7.00000 -1.00000 -1.000 2.00000 -7.00000 2.00000 1.000 2.00000 -7.00000 5.00000 -1.000 2.00000 -4.00000 -7.00000 1.000 2.00000 -4.00000 -4.00000 -1.000 2.00000 -4.00000 -1.00000 1.000 2.00000 -4.00000 2.00000 -1.000 2.00000 -4.00000 5.00000 1.000 2.00000 -4.00000 8.00000 -1.000 2.00000 -1.00000 -7.00000 1.000 2.00000 -1.00000 -4.00000 -1.000 2.00000 -1.00000 -1.00000 1.000 2.00000 -1.00000 2.00000 -1.000 2.00000 -1.00000 5.00000 1.000 2.00000 -1.00000 8.00000 -1.000 2.00000 2.00000 -7.00000 1.000 2.00000 2.00000 -4.00000 -1.000 2.00000 2.00000 -1.00000 1.000 2.00000 2.00000 2.00000 -1.000 2.00000 2.00000 5.00000 1.000 2.00000 2.00000 8.00000 -1.000 2.00000 5.00000 -7.00000 1.000 2.00000 5.00000 -4.00000 -1.000 2.00000 5.00000 -1.00000 1.000 2.00000 5.00000 2.00000 -1.000 2.00000 5.00000 5.00000 1.000 2.00000 5.00000 8.00000 -1.000 2.00000 8.00000 -4.00000 1.000 2.00000 8.00000 -1.00000 -1.000 2.00000 8.00000 2.00000 1.000 2.00000 8.00000 5.00000 -1.000 5.00000 -7.00000 -4.00000 1.000 5.00000 -7.00000 -1.00000 -1.000 5.00000 -7.00000 2.00000 1.000 5.00000 -7.00000 5.00000 -1.000 5.00000 -4.00000 -7.00000 1.000 5.00000 -4.00000 -4.00000 -1.000 5.00000 -4.00000 -1.00000 1.000 5.00000 -4.00000 2.00000 -1.000 5.00000 -4.00000 5.00000 1.000 5.00000 -1.00000 -7.00000 -1.000 5.00000 -1.00000 -4.00000 1.000 5.00000 -1.00000 -1.00000 -1.000 5.00000 -1.00000 2.00000 1.000 5.00000 -1.00000 5.00000 -1.000 5.00000 -1.00000 8.00000 1.000 5.00000 2.00000 -7.00000 -1.000 5.00000 2.00000 -4.00000 1.000 5.00000 2.00000 -1.00000 -1.000 5.00000 2.00000 2.00000 1.000 5.00000 2.00000 5.00000 -1.000 5.00000 2.00000 8.00000 1.000 5.00000 5.00000 -7.00000 -1.000 5.00000 5.00000 -4.00000 1.000 5.00000 5.00000 -1.00000 -1.000 5.00000 5.00000 2.00000 1.000 5.00000 5.00000 5.00000 -1.000 5.00000 8.00000 -1.00000 1.000 5.00000 8.00000 2.00000 -1.000 8.00000 -4.00000 -4.00000 1.000 8.00000 -4.00000 -1.00000 -1.000 8.00000 -4.00000 2.00000 1.000 8.00000 -1.00000 -4.00000 -1.000 8.00000 -1.00000 -1.00000 1.000 8.00000 -1.00000 2.00000 -1.000 8.00000 -1.00000 5.00000 1.000 8.00000 2.00000 -4.00000 -1.000 8.00000 2.00000 -1.00000 1.000 8.00000 2.00000 2.00000 -1.000 8.00000 2.00000 5.00000 1.000 8.00000 5.00000 -1.00000 -1.000 8.00000 5.00000 2.00000 1.000

追記:

特に記述しなかったが、出力エネルギーはポテンシャルエネルギーなので、念のため。また、ここでは時間積分について考慮していない。NVEアンサンブル(定エネルギー)であれば、時間積分は速度ベルレ法がもっとも実装しやすいように思う。dt/2で速度を求め、そこからdtの座標を計算し、その座標からさらにdt/2で速度を求める。

続きを読む...

2009年3月13日金曜日

C++のテンプレートで素数を求める

このブログ記事をはてなブックマークに追加

ふと、C++のテンプレートで素数の計算をさせてみようと思い立ち、コードを書いてみた。書いた後でググってみたら同じことを考えた人がやっぱりいたけど気にしない。これは車輪の再発明ではなくてパズル遊びだから。コードの書き方も違うしね。因みに、このテンプレートによる素数計算ではコンパイル完了時に既に素数が求まっており、実行ではそれを表示するだけなので、ある意味、最速の素数計算プログラムと言えるかも。

短いコードだから読めばわかると思うけど、一応説明しておく。Primesに渡したテンプレートパラメータnを1ずつ減らしながら順番にPrimeに渡して、それが素数なら出力する。素数の判定では、Primeのテンプレートパラメータdの初期値を2とし、dを増やしながらnとの剰余を求め、途中で剰余0となれば特殊化されたPrime<n, d, 0, false>が呼ばれて何もしない。dについてはテンプレートの入れ子を減らしてコンパイル時間を短くするために、2および3で割り切れない値を次のdとしている。dの二乗がnより大きくなった時点(つまりnがnの平方根以下の整数で割りきれなかった場合)、特殊化されたPrime<n, d, r, true>が呼ばれて値を出力する。

以下に100までの素数を出力するコードを示す。main関数のPrimes<100>();の数値を変更すれば任意の数までの素数を出力できる。

#include <iostream> template<int n, int d = 2, int r = 1, bool is_prime = false> struct Prime { Prime() { Prime<n, (d > 2) ? ((d + 2) % 3 ? d + 2 : d + 4) : d + 1, n % d, (n < d * d)>(); } }; template<int n, int d> struct Prime<n, d, 0, false> {}; template<int n, int d, int r> struct Prime<n, d, r, true> { Prime() { std::cout << n << std::endl; } }; template<int n> struct Primes { Primes() { Primes<n - 1>(); Prime<n>(); } }; template<> struct Primes<1> {}; int main() { Primes<100>(); return 0; }

続きを読む...

2009年3月6日金曜日

C++: コンパイラは何を知りたがっている?

このブログ記事をはてなブックマークに追加

template<class T> void print_vector(const std::vector<T>& v) { for (std::vector<T>::const_iterator p = v.begin(); p != v.end(); ++p) std::cout << *p << std::endl; }

問題: 上に書いたコードを含むソースをコンパイルしたらエラーが出た。コンパイラにとって何かの情報が足りないようだ。どこが間違っているのだろうか?

答えは以下の通りで、「typenameが抜けている」。

template<class T> void print_vector(const std::vector<T>& v) { for (typename std::vector<T>::const_iterator p = v.begin(); p != v.end(); ++p) std::cout << *p << std::endl; }

コンパイラからはテンプレートパラメータTを含むstd::vector<T>::const_iteratorが型かどうか判断できないので、型であることを明示的に教える必要がある。

まあ、コンパイラによってはtypenameがなくても通っちゃったりするんだけどね。

続きを読む...

2009年3月3日火曜日

Bloggerの記事が表示されない不具合

このブログ記事をはてなブックマークに追加

Blogger(blogspot)の記事(エントリ)が正常に表示されない不具合が確認されている。こちらでもその不具合を確認した。

ブログのタイトルやサイドバーは表示されるのだが、記事そのものがブランクとなり表示されない。表示されない記事はランダムで決まるようだ。Googleのディスカッションボードでもかなりの被害が報告されている

取り敢えずの対策として、新しい記事をドラフトで保存しキャッシュをクリアする方法がある。ただし、これは一時的な対策であり、しばらくするとまた表示されなくなったりする。

<script type='text/javascript'> google.friendconnect.container.renderUrlCanvasGadget({'id': 'gadget-canvas'}); </script>

表示されないページのソースを読むと上記のようなコードが吐き出されていることから、Google Friend Connect関連のエラーかと思い、ダッシュボードのプロフィールの設定で「読者になっているサイトを表示」のチェックをはずしたところ、今のところ不具合は出ていない。もっとも今後再発するかもしれないがそしたらそのとき考えよう。

続きを読む...

2009年3月1日日曜日

C言語: ポインタと配列について

このブログ記事をはてなブックマークに追加

C言語の入門書を読み始めたばかりの初心者である学生が、先生に質問をしているようです。どんな内容かちょっと聞いてみましょう。

ポインタの使い方

学生: 先生、ポインタの宣言ってこれで合っていますよね。

int *p;

でもp[0] = 5;のように使うとエラーが出るんです。

先生: それは当たり前だ。まだメモリを確保していないだろう? メモリを確保しないと使えないぞ。

学生: でも配列の場合、int a[10];のように宣言したら使えましたよ。ポインタって全然ダメですね。

先生: ポインタと配列はまったく違うものだから、どちらがダメでどちらが良いと云うことはないよ。配列はその場でメモリを確保して使えるが(上記の場合はintを10個分)、プログラム内で決め打ちした数値になってしまう。つまり後から大きさの変更ができないんだ。

一方、ポインタは宣言しただけでは使えないが、後から自由にメモリを割り当ててその大きさを変えられるから、少ないメモリで済む場合は小さく、たくさんのメモリが必要な場合は大きくメモリを自由に確保することができる。

学生: へぇ。ところで、どうやってメモリを確保すればいいんですか?

先生: それはこのようにすればいい。

p = (int*)malloc(10 * sizeof(int));

これは、intのポインタにintの大きさで10個分のメモリを割り当てるという意味だ。確保したメモリが不要になったなら

free(p);

で開放する。この辺についてはその入門書にも載っているだろう。

学生: ふ~ん。なんだか難しそうですね。

先生: 慣れれば難しいことなどないよ。C言語ならmallocでメモリを割り当て、freeで開放と覚えていればいい。

学生: 分かりました。早速使ってみます。

ポインタと配列の違い

学生: ポインタを使ってみましたが、結局、中の数値を見たり変更したりする場合は、p[2]のように配列の形にしないといけないんですよね。なんだか混乱するなぁ。

先生: 別にp[2]ではなくても*(p+2)としても良いぞ。どちらもまったく同じ意味だ。配列型と同じように見えるからといって配列ではないからその点は気をつけることだな。因みに、*(p+2)*(2+p)とも書けるので、2[p]とも書けることになり、実際にこれは文法上正しいのだ。

学生: えええ! そうだったんだ! 先生、これからは2[p]のように書くことにします!

先生: おいおい…、文法上正しくてもコーディングスタイルとしては正しくないぞ。文法も大事だが、それと共に分かりやすく書くことも大事なことなのだ。2[p]のような書き方はせいぜいコーディングの汚さを競うコンテスト、IOCCCぐらいでしか役に立たん。

学生: そうですかぁ…。じゃあ、使わないことにします。ところでさっき、先生はポインタと配列は違うようなことおっしゃってましたけど、配列からポインタに代入できますよね。それって結局同じことなんじゃないですか?

先生: それは違うぞ。例えばint a[10];という配列があるとして、int *p = a;とすることはできる。これはpに配列aの先頭の位置情報(アドレス)を入れるという意味だ。だからポインタpを使ってaの配列にアクセスできるようになる。

この場合、a[2]aの配列の先頭から2つ先の数値という意味だが、p[2]はポインタpに入っている位置情報(aの配列の先頭位置)に2を加えた場所に入っている数値という意味になる。結局は同じ数値にたどり着くんだが、そこにたどり着くまでの経路が違うんだ。

学生: …ややこしいですね。だったら、char *p = "string";char a[] = "string";も違うってことですか? 同じように見えますけど…。

先生: もちろん違う意味だ。aでは配列にそれぞれ's', 't', 'r', 'i', 'n', 'g', '\0'の文字が入るわけだが、pではどこか遠くのメモリ領域(場合によっては書き換え不可)に"string"という文字列があって、その先頭の位置情報(アドレス)が入るわけだ。だからa[0] = 'a'として配列の内容を変更することは問題ないが、p[0] = 'p'とするとメモリ違反のエラーになることもある。

学生: やっぱり配列とポインタは違うんですねぇ。あれ? でも関数の仮引数、例えば、

void func(int a[]) { ...; }

では配列で取ったりしますよね。あれはどういうからくりなんですか?

先生: これは便宜上このような記法で書けるだけで意味はただのポインタだ。だから、仮引数のint a[]int *aと読み替えれば間違いはない。

学生: なるほど、そうだったんですか~。ありがとうございました!

続きを読む...