2005年06月

覆面算の計算で各桁をテンプレートで生成してやることにした．これでソースは見やすくなった．しかも，全部を展開してあった前のコードより若干速いらしい．ちなみに，関数テンプレートの部分特化(Partial Specialization)は禁止されているようなので，部分的にインスタンス化するため仕方なくクラスを使う羽目になった．ま，次は答えをテンプレートで生成だな．

#include <iostream>
using namespace std;
template<const int m, const int l>
class Calc{
public:
    inline void calc(unsigned int &x, unsigned int &y, unsigned int (&digits)[m])
    {
        const int n = m/4*3 - l*3;
        unsigned int d0 = digits[n];
        for(int i0 = n; i0 < m; i0++){
            digits[n] = digits[i0];
            digits[i0] = d0;
            x = (digits[n]<<(4*(n/3))) | (x&((1<<(4*(n/3)))-1));
            unsigned int d1 = digits[n+1];
            for(int j0 = n+1; j0 < m; j0++){
                digits[n+1] = digits[j0];
                digits[j0] = d1;
                y = (digits[n+1]<<(4*(n/3))) | (y&((1<<(4*(n/3)))-1));
                if(n==0 && y >= x) {            // omit symmetric pattern
                    digits[j0] = digits[n+1];
                    continue;
                }
                unsigned int z1 = ((x * y) >> (4*(n/3)))& 0xf;
                int k1 = n+2;
                for(; k1 < m; k1++) if (digits[k1]==z1) break;
                if(k1 >= m) {
                    digits[j0] = digits[n+1];
                    continue;
                }
                digits[k1] = digits[n+2];
                digits[n+2] = z1;
                Calc<m,l-1>().calc(x, y, digits);
                digits[n+2] = digits[k1];
                digits[k1] = z1;
                digits[j0] = digits[n+1];
            }
            digits[n+1] = d1;
            digits[i0] = digits[n];
        }
        x &=((1<<(4*(n/3)))-1);
        digits[n] = d0;
    }
};
template<const int m>
class Calc<m,0>{
public:
    inline void calc(unsigned int &x, unsigned int &y, unsigned int (&digits)[m])
    {
        const int n = m/4*3;
        unsigned int z4 = ((x * y) >> (4*(n/3)));
        unsigned int ds[m-n];
        for(int i = 0; i < (m-n); i++){
            ds[i] = digits[i+n];
        }
        int i = 0;
        for(; i < m-n; i++){
            int k4 = i;
            unsigned int z5 = z4&0xf;
            z4>>=4;
            for(; k4 < m-n; k4++) if (ds[k4]==z5) break;
            if(k4 >= m-n) break;
            ds[k4] = ds[i];
        }
        if(i >=m-n && ds[m-n-1]!=0){
            cout << x << " * " << y << " = " << x*y << endl;
            cout << y << " * " << x << " = " << x*y << endl;
        }
    }
};
int main(int argc, char *argv[])
{
    const int m = 16;
    cout << hex;
    unsigned int digits[m];
    for(int i = 0; i < m; i++){
        digits[i] = i;
    }
    unsigned int x = 0;
    unsigned int y = 0;
    Calc<m,m/4>().calc(x,y,digits);
    return 0;
}

ということで，１６進の覆面算プログラムを作業の合間に組んでみた．８重の for ループでぶん回す．裏で別のプログラムが走ってるけど Pen4 2.6C で0.8秒ぐらい．ま，こんなもんか．

#include <iostream>
using namespace std;
int main(int argc, char *argv[])
{
    cout << hex;
    unsigned int digits[16];
    for(int i = 0; i < 16; i++){
        digits[i] = i;
    }
    unsigned int x = 0;
    unsigned int y = 0;
    unsigned int d0 = digits[0];
    for(int i0 = 0; i0 < 16; i0++){
        digits[0] = digits[i0];
        digits[i0] = d0;
        x = digits[0];
        unsigned int d1 = digits[1];
        for(int j0 = 1; j0 < 16; j0++){
            digits[1] = digits[j0];
            digits[j0] = d1;
            y = digits[1];
            if(y >= x) {            // omit symmetric pattern
                digits[j0] = digits[1];
                continue;
            }
            unsigned int z1 = (x * y) & 0xf;
            int k1 = 2;
            for(; k1 < 16; k1++) if (digits[k1]==z1) break;
            if(k1 >= 16) {
                digits[j0] = digits[1];
                continue;
            }
            digits[k1] = digits[2];
            digits[2] = z1;
            unsigned int d3 = digits[3];
            for(int i1 = 3; i1 < 16; i1++){
                digits[3] = digits[i1];
                digits[i1] = d3;
                x = (digits[3]<<4) | (x&0xf);
                unsigned int d4 = digits[4];
                for(int j1 = 4; j1 < 16; j1++){
                    digits[4] = digits[j1];
                    digits[j1] = d4;
                    y = (digits[4]<<4) | (y&0xf);
                    unsigned int z2 = ((x * y)>>4) & 0xf;
                    int k2 = 5;
                    for(; k2 < 16; k2++) if (digits[k2]==z2) break;
                    if(k2 >= 16) {
                        digits[j1] = digits[4];
                        continue;
                    }
                    digits[k2] = digits[5];
                    digits[5] = z2;
                    unsigned int d6 = digits[6];
                    for(int i2 = 6; i2 < 16; i2++){
                        digits[6] = digits[i2];
                        digits[i2] = d6;
                        x = (digits[6]<<8) | (x&0xff);
                        unsigned int d7 = digits[7];
                        for(int j2 = 7; j2 < 16; j2++){
                            digits[7] = digits[j2];
                            digits[j2] = d7;
                            y = (digits[7]<<8) | (y&0xff);
                            unsigned int z3 = ((x * y)>>8)&0xf;
                            int k3 = 8;
                            for(; k3 < 16; k3++) if (digits[k3]==z3) break;
                            if(k3 >= 16) {
                                digits[j2] = digits[7];
                                continue;
                            }
                            digits[k3] = digits[8];
                            digits[8] = z3;
                            unsigned int d9 = digits[9];
                            for(int i3 = 9; i3 < 16; i3++){
                                digits[9] = digits[i3];
                                digits[i3] = d9;
                                x = (digits[9]<<12) | (x&0xfff);
                                unsigned int d10 = digits[10];
                                for(int j3 = 10; j3 < 16; j3++){
                                    digits[10] = digits[j3];
                                    digits[j3] = d10;
                                    y = (digits[10]<<12) | (y&0xfff);
                                    unsigned int z4 = ((x * y) >>12)&0xfffff;
                                    unsigned int ds[5];
                                    for(int i = 0; i < 5; i++){
                                        ds[i] = digits[i+11];
                                    }
                                    int i = 0;
                                    for(; i < 5; i++){
                                        int k4 = i;
                                        unsigned int z5 = z4&0xf;
                                        z4>>=4;
                                        for(; k4 < 5; k4++) if (ds[k4]==z5) break;
                                        if(k4 >= 5) break;
                                        ds[k4] = ds[i];
                                    }
                                    if(i >= 5 && ds[4]!=0){
                                        cout << x << " * " << y << " = " << x*y << endl;
                                        cout << y << " * " << x << " = " << x*y << endl;
                                    }
                                    digits[j3] = digits[10];
                                }
                                digits[10] = d10;
                                digits[i3] = digits[9];
                            }
                            x &=0xfff;
                            digits[9] = d9;
                            digits[8] = digits[k3];
                            digits[k3] = z3;
                            digits[j2] = digits[7];
                        }
                        digits[7] = d7;
                        digits[i2] = digits[6];
                    }
                    x &=0xff;
                    digits[6] = d6;
                    digits[5] = digits[k2];
                    digits[k2] = z2;
                    digits[j1] = digits[4];
                }
                digits[4] = d4;
                digits[i1] = digits[3];
            }
            x &=0xf;
            digits[3] = d3;
            digits[2] = digits[k1];
            digits[k1] = z1;
            digits[j0] = digits[1];
        }
        digits[1] = d1;
        digits[i0] = digits[0];
    }
}

なげぇなぁ．

比較の仕方がデータによらないソートである Bitonic Sort を Haskell で書いてみた．要素数 2^n 限定．だけ．

import System.Random
randInt :: IO (Int)
randInt = getStdRandom (randomR (0,1000))
randInts n = randInts' n []
 where
 randInts' n xs = 
  if n == 0 then return xs
  else
   randInt >>= (\x -> randInts' (n-1) (x:xs))
  
bimerge xs ys 1 dir = unzip (zipWith (\x y -> if (x < y) == dir then (x,y) else (y,x)) xs ys)
bimerge xs ys n dir = (xxs++xys, yxs++yys)
	where
	(xs', ys') = unzip (zipWith (\x y -> if (x < y) == dir then (x,y) else (y,x)) xs ys)
	n2 = (n `div` 2)
	(xxs, xys) = bimerge (take n2 xs') (drop n2 xs') n2 dir
	(yxs, yys) = bimerge (take n2 ys') (drop n2 ys') n2 dir
  
bisort xs = bisort' xs (length xs) True
	where
	bisort' zs 1 dir = zs
	bisort' zs n dir = xs'' ++ ys''
		where
		n2 = (n `div` 2)
		xs' = bisort' (take n2 zs) n2 dir
		ys' = bisort' (drop n2 zs) n2 (not dir)
		(xs'', ys'') = bimerge xs' ys' n2 dir
  
 -- randInts 1024 >>= (\x -> print (bisort x))
 -- randInts 1024 >>= (\xs -> print (snd (foldl (\(x, f) y -> (y, and [f, (x <= y)])) (-1, True) (bisort xs))))

画像回転で最内のループを１６個展開してみたら１．２倍速くなった．block での DC の最内ループも 32個展開してとうとう毎秒 300回を超えた．いやあ，これだけ単純な計算だと展開したほうがやっぱ速いなぁ．ということで，これにてプログラム凍結．ソースとか

32x32 の変換テーブルを２段階用意して画像回転してみたら見事に１．５倍の速さに．もちろん，元画像は長さを1.5倍して if 分によるチェックをはずしてある．どうやらテーブルを引くことで演算がほぼいらなくなるので速くなるのと，アクセスが32x32のブロック単位になるため，画像がどの方向でもキャッシュのあたりはずれがあまり変わらないことによるようだ．ただ，変換テーブルのところで小数点以下を切り捨てているので今の方法では画面が汚い...　改良できるのだろうか？

ソースとか

画像回転で毎回アドレスをキャストやシフト演算で求めるよりも，テーブルを持っておいてぶん回したほうが速いらしい．ただし，全部のテーブルを持っているとメモリが死ねるので，32x32 の外側ブロックの回転と32x32の内側の回転にわけるらしい．しかも，バウンダリのチェックもはずされていたので...　速くなるのかなぁ？

画像回転でデフォルト画像としてパターン画像を生成していたのだが，友人に美しいパターンの生成方法を教わった．生成式は x 座標と y 座標の XOR で．これをやると見事なパターンがでる．

モノクロならこんな感じで．（パターンを描いたスクリーンショット）

    for(int i = 0 ; i < ImageHeight ; i++) {
      for(int j = 0 ; j < ImageWidth ; j++) {
        *p++ = (i^j)&0xff;
      }
    }

色つけるならこんな感じで．

    for(int i = 0 ; i < ImageHeight ; i++) {
      int r = (i * 0xFF) / ImageHeight;
      for(int j = 0 ; j < ImageWidth ; j++) {
        int y = (i^j)&0xff;
        int ry = r * y/0xff;
        int gy = (( j * 0xFF ) / ImageWidth) * y/0xff;
        int by = (~r&0xff) * y/0xff;
        *p++ = 0xff000000 | ry << 16 | gy << 8| by;
      }
    }

画像回転のメインルーチンで double による計算をしていたのだが，この部分を固定小数演算に置き換えたら劇的に速くなった．キャスト要らずだし，やっぱ整数演算は速いなぁと思っていたら，コンパイル時に -msse2 をつければ double でも SSE2 命令で速くなった．同様に，double を float に置き換えても -msse (SSE) で速くなる．いまさらながらだが，これらの拡張命令のすごさを体験してみたり．ちなみに fps はノートPC で30に達した．描画しなければ 90 だけど．

ということで，記念にソースと実行ファイル一式を置いておこう．馬鹿みたいに重たいスクリーンセーバーもどきか？

画像回転でOpenGLを使っていたのだが，どうやらビットマップよりテクスチャ貼り付けのほうが速いらしい．というわけで，描画を glDrawPixels からテクスチャに切り替え．

    glEnable( GL_TEXTURE_2D );
    glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
    glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
    glTexImage2D(GL_TEXTURE_2D , 0 , GL_RGBA , ImageWidth , ImageHeight, 0 , GL_RGBA, GL_UNSIGNED_BYTE , rotatedPixels);
    glBegin(GL_POLYGON);
        glTexCoord2f(0 , 0); glVertex2f(whiteSpace ,whiteSpace);
        glTexCoord2f(0 , 1); glVertex2f(whiteSpace , ImageHeight+whiteSpace);
        glTexCoord2f(1 , 1); glVertex2f(ImageWidth + whiteSpace, ImageHeight+whiteSpace);
        glTexCoord2f(1 , 0); glVertex2f(ImageWidth + whiteSpace, whiteSpace);
    glEnd();
    glDisable( GL_TEXTURE_2D );

1024x1024 の画像回転プログラム SDL+OpenGL 版（もちろん回転は自前）．Java でアルゴリズムを確かめた後にノートＰＣ上でC++移植していたため 10fps しか出てなかった．しかし，メインマシンに移したら 40fps で動くではないか．ＯｐｅｎＧＬ速いなぁ，そしてノートはやっぱ遅いんだなぁ．さーて任意の画像を白黒で回転できるようにしよう．

どうせ画像を回転させるなら png でも読み込ませようと考えて， libpng を使ってみることに．面倒なのでフルカラー png のみをひとつの大きな位置次元配列に展開する関数を作った．

// full color PNG を一列のバッファにして返す
int readPNG(const char *fname, unsigned char**buf, int *iHeight, int *iWidth, bool *bAlpha, int *iRowbytes)
{
    const int number = 8;
    unsigned char header[number];
    FILE *fp = fopen(fname, "rb");
    if (!fp) return -1;
    fread(header, 1, number, fp);
    bool is_png = !png_sig_cmp(header, 0, number);   // Signature check
    if (!is_png) return -1;
    png_structp png_ptr = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
    if (!png_ptr) return -1;
    png_infop info_ptr = png_create_info_struct(png_ptr);
    if (!info_ptr){
        png_destroy_read_struct(&png_ptr,
           (png_infopp)NULL, (png_infopp)NULL);
        return -1;
    }
    png_infop end_info = png_create_info_struct(png_ptr);
    if (!end_info){
        png_destroy_read_struct(&png_ptr, &info_ptr,
          (png_infopp)NULL);
        return -1;
    }    
    png_init_io(png_ptr, fp);                 // init
    png_set_sig_bytes(png_ptr, number);       // inform libpng read-bytes
    png_read_info(png_ptr, info_ptr);         // get info    
    png_uint_32 width, height;
    int bit_depth, color_type, interlace_type;
    int compression_type, filter_method;
    png_get_IHDR(png_ptr, info_ptr, &width, &height, &bit_depth, &color_type, &interlace_type, &compression_type, &filter_method);
    int channels = png_get_channels(png_ptr, info_ptr);
    int rowbytes = png_get_rowbytes(png_ptr, info_ptr);
    switch(color_type){
    case PNG_COLOR_TYPE_RGB:             // RGB only
    case PNG_COLOR_TYPE_RGB_ALPHA:
        break;
    case PNG_COLOR_TYPE_GRAY:
    case PNG_COLOR_TYPE_PALETTE:
    case PNG_COLOR_MASK_PALETTE:
        //    case PNG_COLOR_MASK_COLOR:
    case PNG_COLOR_MASK_ALPHA:
        return -1;
    }
    if(bit_depth!=8) return -1;
    *buf = new unsigned char[rowbytes*height];
    unsigned char **rows = new unsigned char*[height];
    for(int i = 0; i < height; i++){
        rows[i] = &(*buf)[i*rowbytes];
    }
    png_read_image(png_ptr, rows);
    png_read_end(png_ptr, end_info);
    png_destroy_read_struct(&png_ptr, &info_ptr, &end_info);
    *iHeight = height;
    *iWidth = width;
    *bAlpha = color_type == PNG_COLOR_TYPE_RGB_ALPHA;
    *iRowbytes = rowbytes;
    delete rows;
    return 0;
}

コンパイルは png.h をインクルードしてから -lpng12 -lz -lm をつければいいはず．

ちなみに，Cygwin には mingw の libpng がないのでソースから作る必要があるっぽい．cygwin の libpng を使ったらプログラムが止まらなくなったりしたし．作るのは普通どおりにただ cygwin 使うなと指定してやれば良いみたい．

CFLAGS="-Wall -O3 -funroll-loops -mno-cygwin" ./configure --libdir=/lib/mingw/ --includedir=/usr/include/mingw/
make
make install

Javaではいささか遅すぎるので， SDL + OpenGL で書き直し．

とりあえず，SDL日本語ドキュメントを参照しつつてきとうに組んでみた．画像表示部分には OpenGL を使うので，そっちもてきとうに調べるべし．ちなみに，SDL を入れて C++ でコンパイルしようとしたら SDL_audio.h の 97 行目で文句を言われたので，

 <- void (SDLCALL *filters[10])(struct SDL_AudioCVT *cvt, Uint16 format);
 -> void (SDLCALL *filters)(struct SDL_AudioCVT *cvt, Uint16 format);

こんな書き換えをしてやった．audio 使わない限り大丈夫でしょう．

んで，プログラムの感じは以下のとおりで．まず main 関数は SDL の初期化(VIDEOのみ)をして GL の属性を指定して，Window を作ったら他の初期化をしてからループしてメッセージ処理と描画処理をする．

int main(int argc,char *argv[])
{
     /* SDL初期化 */ 
    if(SDL_Init(SDL_INIT_VIDEO)<0){
        cerr << "SDL Initialization failed." << endl; 
        exit(-1);
    }
    SDL_GL_SetAttribute( SDL_GL_RED_SIZE, 8 );
    SDL_GL_SetAttribute( SDL_GL_GREEN_SIZE, 8 );
    SDL_GL_SetAttribute( SDL_GL_BLUE_SIZE, 8 );
    SDL_GL_SetAttribute( SDL_GL_DEPTH_SIZE, 1 );
    SDL_GL_SetAttribute( SDL_GL_DOUBLEBUFFER, 1 );
    /* 画面作成 */
    int flags = SDL_OPENGL;// | SDL_FULLSCREEN;
    screen = SDL_SetVideoMode(WindowSizeX,WindowSizeY,ColorDepth,flags);
    if(screen == NULL){    
        cerr << "Window Creation failed." << endl; 
        SDL_Quit();
        exit(-2);
    }
    init();
    /* メインループ */
    do{
        processEvent();
        updateScreen();
    } while(!done); 
    SDL_Quit();
    return 0;  
}

メッセージ処理はいつもの Windows のと変わるはずも無く．

int processEvent()
{
    SDL_Event e;
    if (SDL_PollEvent(&e) != 0){  // event がある？
        switch(e.type){
        case SDL_QUIT:
            done = true;
            break;
        case SDL_KEYDOWN:
            switch(e.key.keysym.sym){
            case SDLK_1:                 // アルゴリズム選択
                rotType = 1;
                break;
            case SDLK_2:
                rotType = 2;
                break;
            case SDLK_3:
                rotType = 3;
                break;
            case SDLK_4:
                rotType = 4;
                break;
            case SDLK_p:                 // 描画・非描画フリップ
                painting = !painting;
                break;
            case SDLK_u:                // 描画スキップ
                stepsPerDraw++;
                break;
            case SDLK_d:
                if(stepsPerDraw>1) stepsPerDraw--;
                break;
            case SDLK_f:                 // fps を描画
                fpsDrawFlag = true;
                break;
            case SDLK_q:
            case SDLK_ESCAPE:
                done = true;
                break;
            }
            break;
        }
    }
}

描画部分は fps を計算しつつ適当にと．ゲームなら fps の固定のための wait があるでしょう．

int updateScreen()
{
    static int steps = 0;
    static int time = SDL_GetTicks();
    if(steps%stepsPerDraw==0)
        drawImage();
    rotate();
    steps++;
    int t = SDL_GetTicks();
    if(t - time >= 1000 / freqFPSUpdate){
        totalSteps-=stepsHist[sp];
        totalMillis-=millisHist[sp];
        stepsHist[sp] = steps;
        millisHist[sp] = (int)(t - time);
        totalSteps+=stepsHist[sp];
        totalMillis+=millisHist[sp];
        sp = (sp + 1) % millisHistLength;
        steps = 0;
        fps = totalSteps * 100000 / totalMillis; // fps*100です
        time = t;
    }
}

OpenGL での描画はこんな感じで．glDrawPixels でフレームバッファをダイレクトに転送．その際，GL_LUMINANCE で輝度情報のみにしているが，GL_RGB とかにすると3byte / 1 ドット でカラーの転送とかできる．

画像転送後には文字列を描画している．

void drawImage()
{
    static char str[64]={0};
    static char *s1 = str;
    sprintf(s1, "%d.%d fps (skip %d) algo-%d", fps/100, fps%100, stepsPerDraw-1, rotType);
    if(painting){
        glClear(GL_COLOR_BUFFER_BIT);
        //glDrawPixels(ImageWidth , ImageHeight, GL_BLUE, GL_UNSIGNED_BYTE , rotatedPixels);
        glDrawPixels(ImageWidth , ImageHeight, GL_LUMINANCE, GL_UNSIGNED_BYTE , rotatedPixels);
        //glFlush();
        drawString(s1);
        SDL_GL_SwapBuffers();
    } else {
        if(fpsDrawFlag){
            glClear(GL_COLOR_BUFFER_BIT);
            fpsDrawFlag = false;
            drawString(s1);
            SDL_GL_SwapBuffers();
        }
    }
}

んで，文字列の描画には ASCII だけでいいなら glut のビットマップフォントが使えるのでそれを使う．font_init で各キャラクタの書き方を指定しておいて，あとでglCallListsで使う．font_init のループでは文字のASCIIコードに対してその文字の描画コードを対応付けているので，glCallLists でASCIIコード列であるところの文字列を渡すとうまい具合に文字が表示できる(リスト生成時にスクランブルしとくと楽しいかも)．ちなみに，日本語も使えるようにするにはここら辺が参考になるかも．

void font_init()
{
    base = glGenLists(128);
    for(int i=0;i<128;i++){
        glNewList(base+i, GL_COMPILE);
        glutBitmapCharacter(GLUT_BITMAP_HELVETICA_18, i);
        glEndList();
    }
    glListBase(base);
}
void drawString(const char *out, float r=1.0f, float g=1.0f, float b=1.0f)
{
     glPushAttrib(GL_ALL_ATTRIB_BITS);
     glColor3f(r, g, b);
     glRasterPos2i(whiteSpace/2, WindowSizeY-statusFieldY+5);
     glCallLists( strlen(out) , GL_BYTE, out);
     glPopAttrib();
}

あとはまあ，初期化部分ののこりをば．glOrtho のことろは 2D の描画の座標とスクリーン上の座標のスケールをあわせるためのもの．これをやっとかないと glRasterPos2i で予想以上に座標が動きすぎる．

void init()
{
     glViewport(0,0, WindowSizeX, WindowSizeY);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity(); 
    glOrtho(0.0, (GLdouble) WindowSizeX , 0.0, (GLdouble) WindowSizeY , -1.0, 1.0);
    glRasterPos2i(whiteSpace, whiteSpace);
    fps = 0;
    steps = 0;
    totalSteps = 0;
    totalMillis = 0;
    for(int i = 0; i < millisHistLength; i++){
        stepsHist[i] = 0;
        millisHist[i] = 0;
    }
    sp = 0;
    font_init();
    glDisable( GL_DEPTH_TEST );
    glDisable( GL_LIGHTING );
}

んで，最後にコンパイルはこんな感じで．オプションが長い．．．

g++ ImageRotator.cpp -o ImageRotator `sdl-config --cflags --libs` -lglut32 -lopengl32 -lglu32 -mno-cygwin -O3

とりあえず動くものをここにおいておこう．インクルードが <SDL/SDL.h> でなく <SDL.h> でないとだめかも．

少し前に画像の周りを黒く塗っておいて if を使わなければ速いだろうと予想したのだが，結果としてうまくいかん．各行の元画像からのコピーができる両端を先に計算しておくタイプのほうが速い．よくよく考えると，両端の計算は縁の長さ程度のオーダーですむ上に，その領域から外れている部分は黒く塗りつぶすことが確定しているので元画像へのアクセス不要．一方，縁を黒く塗りつぶして大きくした画像を使った場合は，縁の計算がいらない代わりに黒い部分も元画像にアクセスするためキャッシュが効かない状況では不利．

キャッシュをちゃんと考えないといけない問題で下手に番兵を使うと痛い目を見るいい例だな．以後気をつけよう．

今回のプログラムの課題は 1024x1024 の画像を1度ずつ回転させろとこのこと．速度重視で秒間に何回まわせるかを追求する．で，早速速度度外視の Java でアルゴリズムの検証をば．

まず，ピクセルは回転先画像の各ピクセルからから元の画像の対応するピクセルをコピーすれば穴はなくなると．んで，対応するピクセルが存在するかの判定を毎回するのは馬鹿なので，各行に対して対応するピクセルが存在する範囲をあらかじめ求めておいてジャンプを減らすと．これをやると 16fps から 19fps に速度アップ（でもあまり大きく上がらないなぁ．）

もうひとつ，普通に走らせていると90度ごとに fps が遅くなることに気づく．これは，90度回転あたりで row major に格納されているデータを column major でアクセスし，キャッシュミスが大量に発生することによる．ということで，キャッシュを当てやすくするように90度回転したデータをも保持しておいて，角度の近いほうからコピーするように改善．これで結構速度が安定してきた．

あとは，メモリを倍喰うけど元画像の周りに黒塗りの領域を作っておいて if の判定自体をなくしてしまうのと，これに伴って90度回転以外の中間画像も作っておいてキャッシュを当てやすくする．んで，効果が確認できたら速度重視のC++に移ると．

なんとなく Linux でも Windows でも動くGUI プログラムを作りたかったので wxWidgets をインストール．Windows でも cygwin(MinGW) から使うので，展開後に ./configure; make; make install .　まあ，結果として問題なく動くのだが，configure の途中でWindows の find が呼ばれてこけたのは少々痛い．パスの順番を cygwin 優先にして解決したけど，find とか sort とか windows につけるなよと思う．

とりあえず，フレームを出すだけのプログラムとコンパイルコマンドを書いておこう．

// g++ hello.cpp -o hello `wx-config --cxxflags --libs` -mwindows
#include "wx/wxprec.h"
#ifndef WX_PRECOMP
    #include "wx/wx.h"
#endif
 
class MyApp : public wxApp
{
public:
    virtual bool OnInit();
};
  
IMPLEMENT_APP(MyApp)
  
bool MyApp::OnInit()
{
    // create the main application window
	wxFrame *frame = new wxFrame(NULL,
                                 -1,
                                 wxT("Hello"),
                                 wxPoint(100, 100),
                                 wxSize(600, 480));
    frame->Show(true);
    return true;
}

研究室のメンバーが Haskell でリストをワンパスで２分割するプログラムを書いていたので便乗．

halfSplit l = let (len, ret) = halfSplit' (div len 2) l in ret
 where
  halfSplit' _ []     = (0, ([],[]))
  halfSplit' n (x:xs)= 
    let (len, ps) = halfSplit' (n-1) xs
    in (len + 1, if n > 0 then (x:fst ps, snd ps) else (fst ps, x:snd ps))
 
halfSplit2 l = let (len, ret) = halfSplit' (div len 2) l in ret
 where
  halfSplit' _ []     = (0, ([],[]))
  halfSplit' n (xxs@(x:xs))= 
    let (len, ps) = halfSplit' (n-1) xs
    in (len + 1, if n > 0 then (x:fst ps, snd ps) else ([], xxs))

halfSplit だとリストの後ろ半分も再構成しているが，halfSplit2 のようにすると後ろ半分の再構成がない分簡約ステップ数が減る．実際，Hugs で :set +s して簡約数とかを見てみると，

Main> halfSplit [1..100]
(4619 reductions, 8790 cells)
Main> halfSplit2 [1..100]
(3832 reductions, 7950 cells)   

のようにそれなりに差が出る．

Haskell でプログラムを動かしたときに，その実行とレースを取りたいことが時々？ある．でも，一般にHaskellで文字列を出力しようとするとモナドが出てきて面倒．

で，簡単にそれをやる方法があった．Debug.Trace.trace という関数でそれが簡単にできる．例えば mis の途中結果を知りたければ

import Debug.Trace
mis [] i = i
mis (x:xs) i = trace (show x ++" with "++show i) (mis xs (max i 0 + x))

とすれば途中結果が見れる．かなりありがたいかも．ただ，表示ために本来計算しない部分を計算したりすると動作が変わるので注意．

D3DXFont を使うと GDI を使わずに DirectX の世界だけで文字列を表示できる．で，早速それを使って描いてみたわけだが文字の縁がないので少々見にくい．そこで縁取りを入れてみることにした．

やり方は一番単純なやつで，縁の色で上下左右にずらした文字を先に描画，そのあと目的のを真ん中にと．これをやると５回描画が必要だけど，たかが文字列なので気にしない．他にうまいやり方ってあるのかなぁ？

DirectX でプログラムを組んで遊んでいるわけだが，fps の固定にも少々てこずったり．時間間隔を得るのに timeGetTime() を使ってたら分解能が 5ms でfps 固定できないし...　結局, QueryPerformanceCounter, QueryPerformanceFrequency を使うはめになった．こっちの精度はばか高いので，fps もかなり固定できるようになった．ただ，まだ 0.02 程度の揺らぎがあるけど...

DirectX のプログラムを書こうとしてクオータニオンが何なのか分からず調べることに．結果としてはただの四元数のことだったのだが，これが三次元上の任意軸の上の回転を現していることははじめて知った．ということで，検索に引っかかったよさげな解説書のリンクを張っておこう．阪大基礎工の方が書かれたらしい．

Haskell を使って「冪の下の桁から辞書順ソート」を組んでみた．メモリを節約するアイデアは Java のときと同様，必要な下の桁だけ計算すること．Java の実装では必要になった上位桁の再計算を明示的に書かなければならなかったが，Haskell だと遅延評価で勝手に必要なところだけを計算してソートしてくれる．おかげでコードもすっきり．いやぁ，遅延評価って本当に良いですねぇ．

> module Main (main) where
> import List
> import System
> import Char
> atoi :: String -> Int
> atoi = foldl (\x y-> x*10 + ord y - ord '0') 0
> main = getArgs >>= \lines -> let (n,b) = getOpt lines (10000, 2) in printans n b
> getOpt :: [String] -> (Int, Int) -> (Int, Int)
> getOpt [] p = p
> getOpt [x] p = p
> getOpt (x:y:zs) (n, b) = case x of 
>                               "-b" -> getOpt zs (n, atoi y)
>                               "-n" -> getOpt zs (atoi y, b)
>                               _    -> getOpt (y:zs) (n, b)
> printans n b = putStr (sortedliststring n b)
> sortedliststring n b = (foldr (\x y-> show b ++ "^"++(show x) ++"\n"++y) [] (sortedlist n b))
> sortedlist n b = map (\(x,y)->y) (List.sort(take n (powers b)))
> powers k = pows ([1], 0) 
>  where 
>  pows (an, n) = (map nrev p, n+1):pows (p, n+1)
>   where p = mmult an 0
>  base = 1000
>  mmult [] c   = if c == 0 then [] else (c `mod` base):mmult [] (c `div` base)
>  mmult (x:xs) c = let v = x * k + c in (v `mod` base):mmult xs (v `div` base)
>  nrev x = nrev' base 0 x
>   where 
>   nrev' 1 r n = r
>   nrev' b r n = nrev' (b `div` 10) (r * 10 + (n `mod` 10)) (n `div` 10) 

ちなみに，b = 2, n = 1000000 でうちのマシンだと 46秒 + Mem500M ていど．Java で 25秒 + Mem80M ていどなので，まあスピードはこんなもんでしょう．メモリは簡約を続けるために覚えとく分が多いからしょうがないし．

2^2000000 までで 2分+Mem1G. 限界...　ついでなのでソースを置いておこう．Haskellバージョン, Javaバージョン

Python Challenge の Lebel 9 で Tk を使おうと思ったのだがうまく動かない...　Tk の Window はでるけど Widget が描画されないし．何が悪いのやら？

階乗のソートをやろうとしている人間がいたので講義中に実装してみた．アイデアは最下位のゼロの数とそれに続くいくつかの桁数を保持してソート．計算中に 5 の倍数がきたときに末尾のゼロを捨てるので，その処理をしても比較に使う桁数を正確に保持していられるように工夫が必要．とりあえず 10000000 でも動くので良し．

ということで，横でゲームをしつつ冪を下位桁からの辞書順にソートするプログラムを組んでみた．冪をソートするときにはその桁数が大きすぎてメモリに載らないことが問題になる．そこで，基本的にメモリの節約をメインに考えて作成した．節約のための重要点は，冪をソートするといってもその全ての桁が必要ではなく，基本的には下位の桁の一部しか必要でないこと．そして，冪の下位桁は上位桁を計算しなくても正確に求められ，再計算も簡単であること．これらに基づき，最初は下位の一部を計算しておき，ソートのときに区別が付かなければより上位の桁まで再計算するようにした．

実際に計算してみると，２の冪の場合には冪が20000000になっても下位16桁だけでソートが完了．思いのほか下位桁がバラバラのようである．このときのメモリ使用量は 1G強．時間は...　分単位．それなりに動くプログラムが出来上がって満足．単位は去年とってしまったけど課題を提出しても良いかなぁ？

講義（単位取得済み）の課題で適当な基数の冪を下の桁からの辞書順でソートするプログラムを作れというものが出た．ついでに，階乗も同様にソートせよとの指示も出ている．

とりあえず，階乗を下の桁からソートするのは単純なのでプログラムする気力なし．やり方は以下のとおり．階乗の場合は 5 の倍数が現れるごとに 0 が少なくとも一つ増えるので，再開桁から連続する 0 の数が等しい５個程度の組をソートしておけば，あとはその組を 0 の数の多い順に並べるだけである．よって，このソートは入力サイズを n とすると O(n) で終了する．ちなみに，５個程度の組をソートするのに必要な有効桁数は１０桁もあれば足りるはずなので，

階乗を順に生成するのも n に関係ないオーダでできる．

ということで，ターゲットは冪のソートで，基数最下位桁が 0 でないものに．とりあえず案はできているのでコードを書きますか．