スキップしてメイン コンテンツに移動

テトレーション・フラクタル

5年ほど前、このブログで虚数のテトレーションという記事を書いたことがある。テトレーション(tetration)とは、自らのべき乗を指定された回数反復する演算のことで、na と表現する。35 の場合、555 = 1.911×102185 となる。Pythonの関数で表現すれば以下のようになる。

def tetration(a, n):
b = 1.0
for i in range(n): b = a**b
return b
view raw gistfile1.py hosted with ❤ by GitHub

以前の記事では、∞i が 0.43828+0.36059i に収束することを見つけたのだが、今回はそれに関連したフラクタルについて紹介したい。

テトレーション na の a には複素数を指定することができる。このとき、a = x + yi として、それを複素平面に置く。ここで正の整数nを大きくしていき、発散と判定されたnに対応した色で平面を色分けする。発散しない場合、予め決めておいたnの上限値を使う。これで得られる図をテトレーション・フラクタル(tetration fractal)と呼ぶ。

n(x + yi)=a + biとしたとき、n+1(x + yi)=a' + b'i は以下のように計算できる。




Pythonでの実装はActiveSate CodeにTetration Fractalとして載っている。今回はこれを少し修正したコード、C++で書いたコード、更にTBBで並列化したコードを用意して、それらの実行速度のベンチマークを取ってみた。

Python C++ C++ with TBB
258.732秒 16.999秒 0.976秒

ここで、複素平面の領域は(-1.5, 0)-(-0.75, 0.75)であり、最大繰り返し数は256としている。画像の大きさは1,024×1,024とした。実行環境はIntel Xeon X5650の2CPUで、12コア・24スレッドとなっている。

以下に、ソースコードを示す。詳細についてはコードを読んでいただきたい。

Pythonによる実装:

% ./tetration.py tetration.png -1.5 0.0 0.75 0.75 1024 1024

#!/usr/bin/env python
import sys, os, math, time
from PIL import Image
def tetration(l, t, w, h, sw, sh, mit):
data = []
for ky in range(sh):
for kx in range(sw):
x = l + w * kx / (sw - 1)
y = t + h * ky / (sh - 1)
for i in range(mit):
try:
e = math.exp(-0.5 * math.pi * y)
p = math.pi * x / 2
x = e * math.cos(p)
y = e * math.sin(p)
except:
break
if math.hypot(x, y) > 1000000:
break
data.append((i % 4 * 64, i % 8 * 32, i % 16 * 16))
return data
def main(args):
if len(args) < 6:
print >>sys.stderr, "Usage: %s imagefile left top width height [screen_width=512 screen_height=512 max_iters=256]" % os.path.basename(args[0])
print >>sys.stderr, " Ex. %s tetration.png -1.5 0.0 0.75 0.75" % os.path.basename(args[0])
sys.exit(1)
imagefile = args[1]
l, t, w, h = map(float, args[2:6])
sw, sh, mit = 512, 512, 256
if len(args) >= 7: sw = int(args[6])
if len(args) >= 8: sh = int(args[7])
if len(args) >= 9: mit = int(args[8])
print "Range: (%f, %f) - (%f, %f)" % (l, t, l + w, t + h)
print "Image: %d x %d" % (sw, sh)
print "Max number of iterations: %d" % mit
start_time = time.time()
data = tetration(l, t, w, h, sw, sh, mit)
print "Time: %f" % (time.time() - start_time)
im = Image.new("RGB", (sw, sh))
for i, rgb in enumerate(data):
im.putpixel((i % sw, i / sh), rgb)
im.show()
im.save(imagefile)
if __name__ == "__main__": main(sys.argv)
view raw tetration.py hosted with ❤ by GitHub


C++による実装:

% g++ tetration.cpp -o tetration -O3
% ./tetration tetration.out -1.5 0.0 0.75 0.75 1024 1024
% ./image_tetration.py tetration.out tetration.png

#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
#include <cstdlib>
#include <cmath>
#ifdef _MSC_VER
inline double get_time()
{
return static_cast<double>(std::clock()) / CLOCKS_PER_SEC;
}
#else
#include <sys/time.h>
inline double get_time()
{
timeval tv;
gettimeofday(&tv, 0);
return tv.tv_sec + 1e-6 * tv.tv_usec;
}
#endif
using namespace std;
const double PI = 3.1415926535897931;
// left, top, width , height, screen width, screen height, max iterations, data[sw*sh]
void tetration(float l, float t, float w, float h, int sw, int sh, int mit, vector<int>& data)
{
for (int ky = 0; ky < sh; ky++) {
for (int kx = 0; kx < sw; kx++) {
float x = l + w * kx / (sw - 1);
float y = t + h * ky / (sh - 1);
int i = 0;
for (; i < mit - 1; i++) {
float e = exp(-0.5 * PI * y);
if (e != e) break;
float p = PI * x * 0.5;
x = e * cos(p);
y = e * sin(p);
if (x * x + y * y > 1e12) break;
}
data.push_back(i);
}
}
}
int main(int argc, char* argv[])
{
if (argc < 6) {
cerr << "Usage: " << argv[0] << " outfile left top width height [screen_width=512 screen_height=512 max_iters=256]" << endl;
cerr << " Ex. " << argv[0] << " tetration.out -1.5 0.0 0.75 0.75" << endl;
exit(1);
}
fstream fs(argv[1], ios_base::out);
float l = atof(argv[2]);
float t = atof(argv[3]);
float w = atof(argv[4]);
float h = atof(argv[5]);
int sw = 512;
int sh = 512;
int mit = 256;
if (argc >= 7) sw = atoi(argv[6]);
if (argc >= 8) sh = atoi(argv[7]);
if (argc >= 9) mit = atoi(argv[8]);
vector<int> data;
cout << "Range: (" << l << ", " << t << ") - (" << l + w << ", " << t + h << ")" << endl;
cout << "Image: " << sw << " x " << sh << endl;
cout << "Max number of iterations: " << mit << endl;
double start_time = get_time();
tetration(l, t, w, h, sw, sh, mit, data);
cout << "Time: " << get_time() - start_time << endl;
int cnt = 0;
fs << sw << " " << sh << endl;
for (vector<int>::iterator p = data.begin(); p != data.end(); ++p, cnt++) {
fs << *p << " ";
if ((cnt + 1) % sw == 0) fs << endl;
}
return 0;
}
view raw tetration.cpp hosted with ❤ by GitHub


TBBを利用したC++による並列化実装:

% g++ tetration_tbb.cpp -o tetration_tbb -ltbb -O3
% ./tetration_tbb tetration.out -1.5 0.0 0.75 0.75 1024 1024
% ./image_tetration.py tetration.out tetration.png

#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
#include <cstdlib>
#include <cmath>
#include <sys/time.h>
#include "tbb/task_scheduler_init.h"
#include "tbb/parallel_for.h"
#include "tbb/blocked_range.h"
#include "tbb/tick_count.h"
using namespace std;
using namespace tbb;
const double PI = 3.1415926535897931;
// Tetration fractal
class Tetration {
private:
float l, t, w, h;
int sw, sh, mit;
vector<int>& data;
public:
// left, top, width , height, screen width, screen height, max iterations, data[sw*sh]
Tetration(float l, float t, float w, float h, int sw, int sh, int mit, vector<int>& data)
: l(l), t(t), w(w), h(h), sw(sw), sh(sh), mit(mit), data(data) { }
void operator()(const blocked_range<int>& r) const {
for (int n = r.begin(); n != r.end(); n++) {
float x = l + w * (n % sw) / (sw - 1);
float y = t + h * (n / sw) / (sh - 1);
int i = 0;
for (; i < mit - 1; i++) {
float e = exp(-0.5 * PI * y);
if (e != e) break;
float p = PI * x * 0.5;
x = e * cos(p);
y = e * sin(p);
if (x * x + y * y > 1e12) break;
}
data[n] = i;
}
}
};
int main(int argc, char* argv[])
{
if (argc < 6) {
cerr << "Usage: " << argv[0] << " outfile left top width height [screen_width=512 screen_height=512 max_iters=256]" << endl;
cerr << " Ex. " << argv[0] << " tetration.out -1.5 0.0 0.75 0.75" << endl;
exit(1);
}
task_scheduler_init init;
fstream fs(argv[1], ios_base::out);
float l = atof(argv[2]);
float t = atof(argv[3]);
float w = atof(argv[4]);
float h = atof(argv[5]);
int sw = 512;
int sh = 512;
int mit = 256;
if (argc >= 7) sw = atoi(argv[6]);
if (argc >= 8) sh = atoi(argv[7]);
if (argc >= 9) mit = atoi(argv[8]);
vector<int> data(sw * sh);
cout << "Range: (" << l << ", " << t << ") - (" << l + w << ", " << t + h << ")" << endl;
cout << "Image: " << sw << " x " << sh << endl;
cout << "Max number of iterations: " << mit << endl;
tick_count start_time = tick_count::now();
parallel_for(blocked_range<int>(0, sw * sh), Tetration(l, t, w, h, sw, sh, mit, data), auto_partitioner());
cout << "Time: " << (tick_count::now() - start_time).seconds() << endl;
// output data
int cnt = 0;
fs << sw << " " << sh << endl;
for (vector<int>::iterator p = data.begin(); p != data.end(); ++p, cnt++) {
fs << *p << " ";
if ((cnt + 1) % sw == 0) fs << endl;
}
return 0;
}
view raw tetration_tbb.cpp hosted with ❤ by GitHub


C++の実行ファイルで生成された出力ファイルを画像ファイルに変換するスクリプトを以下に示す。

#!/usr/bin/env python
import sys, os
from PIL import Image
def main(args):
if len(args) < 3:
print >>sys.stderr, "Usage: %s datafile imagefile" % os.path.basename(args[0])
sys.exit(1)
datafile, imagefile = args[1:3]
f = file(datafile)
sw, sh = map(int, f.readline().strip().split())
im = Image.new("RGB", (sw, sh))
for y, l in enumerate(f):
for x, data in enumerate(map(int, l.strip().split())):
rgb = (data % 4 * 64, data % 8 * 32, data % 16 * 16)
im.putpixel((x, y), rgb)
#im.show()
im.save(imagefile)
if __name__ == "__main__": main(sys.argv)
view raw image_tetration.py hosted with ❤ by GitHub

ラベル

ラベル:

コメント

hideta さんのコメント…
初めまして。
Cython を使うと、 C++ のコードと同程度の速度で計算が行なえました。手軽なので試す価値はあるように思います。
14/5/12 10:56
nox さんの投稿…
hidetaさん、こんにちは。

CythonはPythonの高速化の他にも、Cへのラッパーとしても使えて便利ですね。
NumPyモジュールとの組み合わせなども有用だと思います。
14/5/12 18:40
コメントを投稿