GCJ 2011 Round2 Spinning Blade復習

問題概要

R*Cの単位ブロックからなる長方形の鉄板がある。この鉄板から1辺がkの正方形を切り出す。さらに、4隅のブロックを切り落として、ブレードを作る。この鉄板は厚さが不均一なため、重心とブレードの中点が一致しないことがある。各ブロックの重さが与えられるので、重心とブレードの中点が一致するような最大のブレードのサイズを求める。

方針

正方形の左上のブロックを選ぶ。これを(x, y)とする。

　(x+k, y+k)が(r, c)を超えない間、kをインクリメント。

　　(x, y), (x+k, y+k)のブレードの重心および中点を計算する。一致すれば、答えに反映する。

これで、smallは通る。全体の処理は、O(N^5)。

largeを通すには、計算量を減らす必要がある。(x, y), (x+k, y+k)のブレードの重心計算は、(x, y), (x+k-1, y+k-1)のブレードの重心情報を利用すれば、O(N)で出来る。これは良くやる動的計画。中点は、((x+x+k)/2, (y+y+k)/2)なので、O(1)。これで全体の計算量はO(N^4)に落ちる。結構きわどいが、これでも通る。

Editorial解は、O(N^3)だった。ACrushとかrngさんとかもO(N^3)でやっていた。この方法は、非常に興味深い。

処理の概要は以下のとおり。

1. (0, 0), (x, y), x in {0,1,..,R-1}, y in {0,1,...C-1}のブレードの重心を予め計算しておく。このとき、(0,0),(x, y)の重心は、(0,0),(x-1,y-1)と、(0,0), (x-1,y)と、 (0,0),(x, y-1)の情報を利用すれば、O(1)で出来る。よって、1.の処理はO(N^2)
2. (x, y)(x+k, y+k)のブレードの重心は、(x, y) = (0, 0)の場合は、1.をそのまま使える。そうでない場合は、(0,0), (x+k, y+k)と、(0,0),(x+k, y-1)と、(0, 0),(x-1,y+k)と、(0,0),(x-1,y-1)の情報を使ってO(1)で求めることができる。よって2.の処理は、(x,y,k)のループのみとなり、O(N^3)

コードに落とすときの工夫だが、2.の場合分け（x=0とかy=0とか）をごちゃごちゃしなくていいように、上行と左列にdummyのバッファを入れておくとよい。

以下ソース。

int mass[512][512];
LL rec_m[512][512];
LL rec_gx[512][512];
LL rec_gy[512][512];

void gcjMain() {
    int t;

    cin >> t;
    REP(ii, t) {
        int r, c, d;

        cin >> r >> c >> d;
        vector<string> mtx(r);
        REP(i, r)
            cin >> mtx[i];

        REP(i, r) REP(j, c)
            mass[i][j] = mtx[i][j] - '0';

        memset(rec_m, 0, sizeof(rec_m));
        memset(rec_m, 0, sizeof(rec_gx));
        memset(rec_m, 0, sizeof(rec_gy));

        REP(i, r) REP(j, c) {
            rec_m[i+1][j+1] = rec_m[i][j+1] + rec_m[i+1][j] - rec_m[i][j] + mass[i][j];
            rec_gx[i+1][j+1] = rec_gx[i][j+1] + rec_gx[i+1][j] - rec_gx[i][j] + i * mass[i][j];
            rec_gy[i+1][j+1] = rec_gy[i][j+1] + rec_gy[i+1][j] - rec_gy[i][j] + j * mass[i][j];
        }

        int ret = -1;
        REP(i, r) REP(j, c) {
            for (int k = 2; i + k < r && j + k < c; k++) {
                LL ms;
                LL gx;
                LL gy;

                ms = rec_m[i+1+k][j+1+k] + rec_m[i][j] - rec_m[i+1+k][j] - rec_m[i][j+1+k];
                gx = rec_gx[i+1+k][j+1+k] + rec_gx[i][j] - rec_gx[i+1+k][j] - rec_gx[i][j+1+k];
                gy = rec_gy[i+1+k][j+1+k] + rec_gy[i][j] - rec_gy[i+1+k][j] - rec_gy[i][j+1+k];

                ms -= mass[i][j] + mass[i][j+k] + mass[i+k][j] + mass[i+k][j+k];
                gx -= i*mass[i][j] + i*mass[i][j+k] + (i+k)*mass[i+k][j] + (i+k)*mass[i+k][j+k];
                gy -= j*mass[i][j] + (j+k)*mass[i][j+k] + j*mass[i+k][j] + (j+k)*mass[i+k][j+k];

                if (ms * (i + i + k) == 2 * gx && ms * (j + j + k) == 2 * gy)
                    ret = max(ret, k+1);
            }
        }

        if (ret == -1)
            printf("Case #%d: IMPOSSIBLE\n", ii+1);
        else
            printf("Case #%d: %d\n", ii+1, ret);
    }
}

サルでも分かるFFT（4）

前回「サルでも分かるFFT(3)」で、FFTのC++による実装を書きました。この実装は、O(n log n)で動くのですが、まだまだ遅いです。

• 入力用と出力用バッファのデータ移動
• 出力用バッファは毎回newしている
• ビット逆転で毎回Nビットなめている
• 同じWを何回も計算している

この辺の処理を工夫して、定数項を小さくします。

　まず、一つ目と二つ目は、in-placeに変換を行うことで解決できます。つまり、入力のデータを直接いじって入力用バッファに結果を書き込みます。三つ目は、ループを一つ増やしてinnermostなloopの外で計算することで解決します。（余談ですが、exp()とかhypot()は相当重い処理です。）

四つ目は、ビット逆転では、LSBとMSBが逆転しているということを利用すれば、効率的に書くことができます。（spaghetti sourceのを参考にしました。）

#include <complex>

typedef complex<double> comp;

/*
 * sgn -1 : time to freq
 *      1 : freq to time
 */
void FFT(comp x[], int n, int sgn) {
    comp theta = comp(0, sgn*2*PI/n);

    for (int m = n; m >= 2; m>>=1) {
        int mh = m >> 1;

        for (int i = 0; i < mh; i++) {
            comp W = exp(1.0*i*theta);

            for (int j = i; j < n; j += m) {
                int k = j + mh;
                comp x_t = x[j] - x[k];
                x[j] += x[k];
                x[k] = W * x_t;
            }
        }
        theta *= 2;
    }

    int i = 0;
    for (int j = 1; j < n - 1; j++) {
        for (int k = n >> 1; k > (i ^= k); k >>= 1);
        if (j < i) swap(x[i], x[j]);
    }
}

前回は、時間間引き（時間成分を偶数、奇数に分ける）のFFTを書きましたが、上のソースは周波数間引き（周波数成分を偶数、奇数に分ける）で書いています。どちらも同じように書けますが、バタフライ演算の向きが逆になります。

ビット逆転のところが興味深いですが、ビット演算を使わない書くと処理の意味が分かります。LSBとMSBを逆に考えて1をインクリメントしています。

void bit_reverse(int x[], int N) {
    int r = 0;

    for (int i = 1; i < N-1; i++) {
        int k = N/2;

        while (r >= k) {
            r -= k;
            k /= 2;
        }
        r += k;

        if (r > i)
            swap(x[i], x[r]);
    }
}

さて、最終テーマであるFFTを用いた多倍長整数乗算の高速化です。

普通に多倍長整数の乗算をやると、O(n^2)の時間が必要です（nは桁数）。100万桁同士の掛け算をすると、30分~1時間程度かかるのではないかと思います。FFTを使うと、O(n log n)で乗算ができます。100万桁同士の掛け算でも数秒で終わる計算量です。

数学的な詳しい解説は、ここのサイトを参照。ここでは、大まかな処理だけをまとめます。

「文字列xとyが与えられる。z = x * yを求める。」とします。

1. 上記のFFTが使用できるように信号長を2のべき乗とします。N >= max(len(x), len(y))となるような最小の2のべき乗Nを求めます。
2. 文字列を信号列に変換します。x_t[]、y_t[]にそれぞれx, yのdigitsをLSBから詰めていきます。x_t[i], i=len(x), len(x)+1, ..., 2*N-1、および、y_t[i], i=len(y), len(y)+1, ...,2*N-1には0を入れます。
3. x_t[]、y_t[]を長さ2*Nの信号としてフーリエ変換します。
4. z_t[i] = x_t[i] * y_t[i], i = 0,1,..., 2*N-1とします。
5. z_t[]を逆フーリエ変換します。
6. z_t[]はx*yと一致しています。ただし、繰り上がりの処理を忘れずに。

以下ソース。SPOJのNot So Fast Multiplicationを解きました。

char num1[1<<16], num2[1<<16];
comp xt[1<<16], yt[1<<16], zt[1<<16];
int carried[1<<16];
char ans[1<<16];

/*
 * conversion from a string to a signal sequence
 */
void convert(const char x[], comp *ret, int len) {
    for (int i = 0; i < len; i++)
        ret[i] = x[i]-'0';

    for (int i = 0; i < len; i++)
        ret[len+i] = 0;
}

/*
 * conversion from a signal sequence to an integer
 */
void toInteger(int N) {
    for (int i = 0; i < N; i++)
        carried[i] = (int)(zt[i].real()/(N) + 0.5);

    for (int i = 0; i < N; i++) {
        if (carried[i] >= 10) {
            carried[i+1] += carried[i] / 10;
            carried[i] %= 10;
        }
    }

    int last = N-1;
    for (; last >= 0 && !carried[last]; last--)
        ;

    for (int i = 0; i <= last; i++)
        ans[last-i] = (char)(carried[i]+'0');
    if (last == -1)
        ans[0] = '0';
}

/*
 * O(n log n) multiplication
 */
void multiply(char x[], char y[]) {
    int len1 = strlen(x);
    int len2 = strlen(y);

    int N = max(len1, len2);

    int b = 1;
    while (b < N)
        b *= 2;
    N = b;

    reverse(x, x+len1);
    reverse(y, y+len2);

    for (int i = len1; i < N; i++)
        x[i] = '0';
    for (int i = len2; i < N; i++)
        y[i] = '0';

    convert(x, xt, N);
    convert(y, yt, N);

    FFT(xt, 2*N, -1);
    FFT(yt, 2*N, -1);

    for (int i = 0; i < 2*N; i++)
        zt[i] = xt[i] * yt[i];

    FFT(zt, 2*N, 1);
    toInteger(2*N);
}

int main() {
    int n;
    scanf("%d", &n);

    while (n--) {
        scanf("%s %s", num1, num2);

        memset(ans, 0, sizeof(ans));
        multiply(num1, num2);
        printf("%s\n", ans);
    }
    return 0;

}

サルでも分かるFFT（3）

ようやくFFTが分かった。結構長い道のりだった。

まずは、普通のO(n^2)のDFTを行列表現で書いてみます。このとき、eのべき乗のところを式で書くとイメージが大変なので、極座標系の矢印で書いてあげると分かりやすいです。

この式をじっくり見ると、「サルでも分かるFFT(2)」で紹介したk成分と(N-k)成分の対称性が見えます。Fourier visualizerがはき出した極座標系の周波数スペクトラムとの対比がよく分かると思います。

さて、FFTですが、まず偶数番と奇数番の周波数成分に分けます。偶数番成分に注目すると、

と書けますが、矢印の部分の4*8の行列は真ん中を中心に対象であることが分かります。よってこれは、

と書けます。この時点で演算数が半分になりました。これを繰り返せば分割統治の容量でO(N^2)からO(N log N)まで落とせそうです。

　同様に、奇数番のものは以下のように表すことができます。

これをゴリゴリやり続けると、最終的に以下の要素に辿り着きます。

ここまでくれば、あとは近い塊から順に例の「バタフライダイアグラム」(参照元：wikipedia、作成者: たまなるたみ氏)に当てはめればいいです。

左側のx[i]の並びですが、ビット逆転（ビット列に直して、reverseをする）が使われています。これは何故？と思いましたが、偶数と奇数の成分に逐次分けていったことを思い出せば理解できます。

まず、上半分は、偶数です。さらにその中の上半分の2個は、4で割れる数（下二桁のビットが00のもの）です。逆に下半分の2個は、ビットが10のものです。

これを再帰的に考えていけば、LSBとMSBの関係が逆転していることが分かります。それでビット逆転を使うと都合がいいわけです。

ほんで、C++のソース。最初pythonで書こうとしていましたが、C++のcomplexライブラリの方が使いやすかったので、C++で書きました。FFTとInverse FFTは、非常に似た処理を行うので、共通部分はgenFFT()にまとめています。ソースのコンパクトさよりも分かりやすさを追求したので、比較的読みやすい（はず？）です。

typedef double Amp;               // amplitude at each sampling time
typedef complex<double> Freq;     // constituent frequencies

int bit_reverse(int x, int n) {
    int rev = 0;

    for (--n; n > 0; n >>= 1, x >>= 1)
        rev = (rev << 1) | (x & 1);

    return rev;
}

vector<Freq> genFFT(const vector<Freq> &x, int sgn) {
    int N = x.size();
    vector<Freq> X(N);

    // set init values with bit reversed
    for (int i = 0; i < N; i++) {
        int rev = bit_reverse(i, N);
        X[i] = x[rev];
    }

    // butterfly diagrams
    for (int b = 2; b <= N; b <<= 1) {
        for (int i = 0; i < N; i++) {
            if (i % b >= b/2)
                continue;

            int k = i + b/2;
            Freq Wi = exp(Freq(0, sgn*2*PI*(i%b)/b));
            Freq Wk = exp(Freq(0, sgn*2*PI*(k%b)/b));

            X[i] += Wi*X[k];
            X[k] = Wk*X[k] + (X[i] - Wi*X[k]);
        }
    }

    return X;
}

vector<Freq> FFT(const vector<Amp> &x) {
    int N = x.size();

    // convert real numbers to complex numbers
    vector<Freq> x_t(N);
    for (int i = 0; i < N; i++)
        x_t[i] = Freq(x[i], 0);

    return genFFT(x_t, -1);
}

vector<Amp> IFFT(const vector<Freq> &X) {
    int N = X.size();

    vector<Freq> x_t = genFFT(X, 1);

    // convert complex numbers to real numbers
    // also normalize the amplituds (devide by N)
    vector<Amp>x(N);
    for (int i = 0; i < N; i++)
        x[i] = x_t[i].real()/N;

    return x;
}

FFTが出来たので、次回「多倍長整数の除算をFFTで高速化する」でFFTシリーズをシメたいと思います。

サルでも分かるFFT（2）

　ふーむ、前回の「サルでも分かるFFT(1)」から時間が空きましたが、第２弾。（※サルは無知な自分のこと。）

　結局のところ、DFTでやるのは、単なる一次変換なので、実装するのはそんなに難しくないです。行列の積と、複素数の扱いくらいが書ければ、簡単にできます。（難しいのは、えられた結果の物理的か解釈だと思います。）また、FFTは、DFTを高速に実施するための一手法であると、前回書きました。

　前回の疑問点が解決したので、それを書きます。

1. 連続信号を離散でとっているので、周波数成分分解のパターンはいくつか存在するはず。DFTでは、解が一意に求まるが、どんな解が選ばれてるの？

これは、実際にvisualizingするツールを作っているときにいろいろ気付きました。

DFTによってえられる周波数成分は、以下の2式の条件を満たしています。

1つ目の式は、周波数スペクトラムが、ある周波数を中心に左右対称になっていることを表しています。

2つ目の式は、中心周波数から同じ距離にある対象な周波数成分は、反対の位相を持っているということを意味します。

以下は、私が作ったvisualizer(※Google Chromeでの閲覧推奨)が吐き出した周波数スペクトラムです。極座標系にプロットしているので、これを見ると、式の意味が一目瞭然です。（上：時間軸波形、下：対応する周波数スペクトラム）

2. F[i]のiって何？F[1]は周波数が1っていうのは分かるけど、これは1[Hz]?それともサンプリング間隔に対して、1倍の周波数を持っているという意味?

これは、サンプリング間隔に対する周波数です。実際の物理的な時間単位とは違います。

すべてのサンプリングをとるのに要した時間を1[s]と考えたときに、F[i]は、i[Hz]の周波数成分となっています。

3. その他

DFTしてえられた周波数成分（F[i], i=0,1, ... N.）を実際の物理的な波に変換するときどうすればいいか？

F[i]に対応する正弦波は以下の式で表されます。

さてさて、記事のタイトルが「サルでも分かるFFT」なので、FFTをやらないといけません。FFTの効果を実感したいなら、サンプル数が10^5とか10^6とかないといけません。(O(N^2)とO(N*log N)の対比なので。)

音楽CDのデータをとってきてそれにDFTかけてみる。

にチャレンジしようと思いましたが、データの抽出に時間がかかりそうなので、「FFTを利用した多倍長整数乗算の高速化」にトライしようかと思います。

ビットを操るgccビルトイン関数たち

gccのビルトイン関数の中で、ビットに関するもののまとめ。

function	feature
__builtin_clz(unsigned int)	count leading zeros
__builtin_ctz(unsigned int)	count trailing zeros
__builtin_ffs(unsigned int)	find first set bit (1-based)
__builtin_popcount(unsigned int)	population count
__builtin_parity(unsigned int)	parity of bits set

■使い方とか

1. _builtin_clz(x)は、xを二進数にしたときの頭の0の数を返します。

10は、2進数だと、

00000000000000000000000000001010
なので、28が返ってきます。

ぱっと思いつく応用は、2^i <= x < 2^{i+1}となるiを求めるとか。

2. __builtin_ctz(x)は、xを二進数にしたときのお尻の0の数を返します。

00000000000000000000000000001010

なので、10の場合は、1です。

これを使うと、一発で、lsbの位置(o-based)が分かります。ただし、lsbを使いたい場合は、1<< lsbという使い方が多いと思うので、そういう場合は、x & -xの方がやすい。

xは、2で何回割り切れる?かを調べたいとき、これを使えば一発。

3. __builtin_ffs(x)は、lsbの位置を1-basedの数字で返します。

__builtin_ffs(x) = 1 + __builtin_ctz(x)

の関係がなりたちます。特にこれと言って便利なシチュエーションは浮かびません。すいません。

4. __builtin_popcount(x)は、1が立っているビット数を返します。

10の場合は、2が返ってきます。この関数はビットで集合を表す場合に、重宝します。

5. __builtin_parity(x)は、1が立ってるビットのパリティなので、以下の関係が成り立ちます。

__builtin_parity(unsigned int) = __builtin_popcount(x) % 2

これは、包除原理の足すか引くかの判定のときに、使えます。

知ってると便利なSTL(9) remove, erase

remove()とerase()は紛らわしい。

この２つの関数がやることは全然違う。

remove()は、"取り除く"という意味なので、指定した値を（意味のある）区間から取り除く。

erase()は、"消去する"という意味なので、指定した範囲のメモリを消去する。

　まず、remove()について詳しくみていく。以下の例を考える。

int x[] = {1,2,3,4,3,6};

int main() {
  int sz = SIZE(x);

  remove(x, x+sz, 3);

  for (int i = 0; i < sz; i++)
    cout << x[i] << endl;
}

この例では、{1, 2, 3, 4, 3, 6}から、3を取り除く。実際にやるのは、x[0]から、x + szまでを走査して、3じゃない値だけをx[0], x[1], ..にコピーし直している。

よって、表示されるのは、{1, 2, 4, 6, 3, 6}となる。

参考までに、この部分のソース(gccの実装)を抜粋すると、

 template<typename _InputIterator, typename _OutputIterator, typename _Tp>
   _OutputIterator
   remove_copy(_InputIterator __first, _InputIterator __last,
        _OutputIterator __result, const _Tp& __value)
   {
     for ( ; __first != __last; ++__first)
    if (!(*__first == __value))
      {
        *__result = *__first;
        ++__result;
      }
     return __result;
   }

という感じ。

remove()は、不要な値を除いてコピーをし終えた位置のイテレーターを返すので、上記の例で、{1,2,4,6}までを表示させたければ、

int x[] = {1,2,3,4,3,6};

int main() {
  int sz = SIZE(x);

  int *p = x;
  int *q = remove(x, x+sz, 3);

  while (p != q)
    cout << *p++ << endl;
}

と書けばよい。vectorについても同様なので、

vector<int> x;

int main() {
  x.push_back(1);
  x.push_back(2);
  x.push_back(3);
  x.push_back(4);
  x.push_back(3);
  x.push_back(6);

  vector<int>::iterator p = x.begin();
  vector<int>::iterator q = remove(x.begin(), x.end(), 3);

  while (p != q)
    cout << *p++ << endl;
}

と書くと、remove後のコンテナ内の値を表示させることができる。この時点では、コンテナのサイズや終端イテレーターの位置は更新されていないことに注意。x.end()は、x[5]の次の位置を指している。

　実際に、メモリ上から値を削除して、コンテナのサイズや終端イテレーターの位置を更新したい場合は、erase()を使う。

erase()には、削除したいイテレーターの位置、またはその範囲を指定することができる。

vector<int> x;
int main() {
  // 1 2 3 4 3 6
  x.push_back(1);
  x.push_back(2);
  x.push_back(3);
  x.push_back(4);
  x.push_back(3);
  x.push_back(6);

  x.erase(x.begin()+2);    // erase x[2]
  for (int i = 0; i < (int)x.size(); i++)
    cout << x[i] << endl;     // 1 2 4 3 6

  x.erase(x.begin()+1, x.begin()+4); // erase the elements [x[1], x[4])
  for (int i = 0; i < (int)x.size(); i++)
    cout << x[i] << endl;     // 1 6
}

という具合。

以下のように、remove()とerase()を併用すると、値を取り除いて範囲外に出すのと、メモリからの消去がまとめて出来る。

int main() {
  vector<int> x;

  x.push_back(1);
  x.push_back(2);
  x.push_back(3);
  x.push_back(4);
  x.push_back(3);
  x.push_back(6);

  x.erase(remove(x.begin(), x.end(), 3), x.end());
  // 1 2 4 6
  for (int i = 0; i < (int)x.size(); i++)
    cout << x[i] << endl;

  return 0;
}

また、文字列から指定した文字を削除する関数は、以下のように書くことができる。

void removeAll(string &s, char t) {
  s.erase(remove(s.begin(), s.end(), t), s.end());
}

サーバー奮闘記（14）schemeインストール

Lispの方言のschemeを入れた。

MIT schemeという処理系をインストールした。

tar xzf mit-scheme-9.0-i386.tar.gz

cd mit-scheme-9.0.1/src

./configure

make compile-microcode

make install

でエラーなくすんなりと入った。

（いつもならコマンドを叩くたびにエラーが出て大変なんだが。）

configureの実体はシェルスクリプトだということをはじめて知った。

makefileの後のcompile-microcodeとかinstallとかは、phony targetっていうやつですね？？

最近makefile書いたから、分かるのさ。

まー何にせよ、ガイドに従うだけじゃなくて、

中身がどうなってるのに興味を持つことは大切だと思います。

schemeを導入したのは、SICPを読むため。