累積和の使い方を学ぶ

ちょっととあるところで出くわしたので。
アルゴリズム苦手勢として、やったことはまとめておきたい。

あるデータ列Xがあって、Xのk個ずつの区間でのそれぞれの和を必要とする問題。(k≦|X|=N)

安直にやったら

t_element x[N], r[N-K];
t_counter i, j;

for (i = 0; i < N - K; i++) {
    r[i] = 0;
    for (j = 0; j < K; j++) {
        r[i] += x[i + j];
    }
}

みたいな感じな気がする(コンパイルしてない)。でも超絶遅い。

でぐぐってみたら http://togetter.com/li/617816 の記事というかまとめを発見。累積和なるものを知る。

要は、

こんなことすると速くなるで、ってこと。

残念ながらそれまでの自分の頭のなかは画像の一番最後辺りの

# zsh
echo $(($(echo "${${x[$p,$q]}// /+}")))

だった。

---

Matlab系だと、(下はOctave 3.8.1 on OSX 10.10で実行)

% Configure
N = 10;
K = 3;

% Initialize data
x = 1 : N;

% Initialize result container
m = zeros(N - K + 1, 1);

% Calculate
for i = 1 : N - K + 1;
    m(i) = sum(x(i : i + K - 1));
end

が普通で、探したらcumsum関数があったので、それを使って書くと、

% Configure
N = 10;
K = 3;

% Initialize data
x = 1 : N;

% Prepare cumulative sum
y = [0, cumsum(x)];

% Initialize result container
m = zeros(N - K + 1, 1);

% Calculate
for i = 1 : N - K + 1;
    m(i) = y(i + K) - y(i);
end

こんな感じ。

以下のコードでベンチマークとってみた。

if 0 == 1;
    printf("This method is just trick\n");
end

function m = test1(N, K, x)
    % Initialize data
    x = 1 : N;

    % Initialize result container
    m = zeros(N - K + 1, 1);

    % Calculate
    for i = 1 : N - K + 1;
        m(i) = sum(x(i : i + K - 1));
    end
endfunction

function m = test2(N, K, x)
    % Prepare cumulative sum
    y = [0, cumsum(x)];

    % Initialize result container
    m = zeros(N - K + 1, 1);

    % Calculate
    for i = 1 : N - K + 1;
        m(i) = y(i + K) - y(i);
    end
endfunction

function y = benchmark(N, K)
    % Initialize
    x = 1 : N;
    t = [0, 0];

    s_t = time();
    test1(N, K, x);
    t(1) = time() - s_t;

    s_t = time();
    test2(N, K, x);
    t(2) = time() - s_t;

    printf("N = %d, K = %d\n    test1 : %10.5fs\n    test2 : %10.5fs\n", N, K, t(1), t(2));

    y = t(1) / t(2);
endfunction

function run_bench()
    U   = 16 + 1;
    N   = 2 ** U;
    res = zeros(U, 2);

    for i = 0 : U - 1;
        K = 2 ** i;
        r = benchmark(N, K);
        res(i + 1, :) = [K, r];
    end

    plot(res(:, 1), res(:, 2));
endfunction

で、結果が

N = 131072, K = 1
    test1 :    3.19987s
    test2 :    2.36953s
N = 131072, K = 2
    test1 :    3.25742s
    test2 :    2.04249s
N = 131072, K = 4
    test1 :    3.29297s
    test2 :    1.96278s
N = 131072, K = 8
    test1 :    3.28375s
    test2 :    2.03915s
N = 131072, K = 16
    test1 :    3.28769s
    test2 :    1.98568s
N = 131072, K = 32
    test1 :    3.50897s
    test2 :    2.15263s
N = 131072, K = 64
    test1 :    3.68561s
    test2 :    2.08608s
N = 131072, K = 128
    test1 :    3.96923s
    test2 :    2.23401s
N = 131072, K = 256
    test1 :    3.99468s
    test2 :    2.22413s
N = 131072, K = 512
    test1 :    4.07600s
    test2 :    1.94705s
N = 131072, K = 1024
    test1 :    4.30081s
    test2 :    2.03379s
N = 131072, K = 2048
    test1 :    5.23567s
    test2 :    1.99727s
N = 131072, K = 4096
    test1 :    6.15645s
    test2 :    1.98333s
N = 131072, K = 8192
    test1 :    9.03844s
    test2 :    1.84308s
N = 131072, K = 16384
    test1 :   14.21577s
    test2 :    1.75381s
N = 131072, K = 32768
    test1 :   23.12840s
    test2 :    1.48494s
N = 131072, K = 65536
    test1 :   28.41491s
    test2 :    1.22984s

なぜかOctaveからプロットの軸のタイトルが設定できなかったんで書くと、

• x軸 : K (2^{i})
• y軸 : 高速化割合(%)

Kが小さいとあまりありがたみないけど、大きくなるにつれて普通のと比べて速いことがわかる。
というか、

N = 131072, K = 2048
    test1 :    5.23567s
    test2 :    1.99727s
N = 131072, K = 4096
    test1 :    6.15645s
    test2 :    1.98333s
N = 131072, K = 8192
    test1 :    9.03844s
    test2 :    1.84308s
N = 131072, K = 16384
    test1 :   14.21577s
    test2 :    1.75381s
N = 131072, K = 32768
    test1 :   23.12840s
    test2 :    1.48494s
N = 131072, K = 65536
    test1 :   28.41491s
    test2 :    1.22984s

あたりで(自分の環境で)どれも2秒切ってるの見てすげぇって思った。
記事にもあったとおり、累積和を求めるのにO(n)、あとは各要素はO(1)で求められるのでこんな結果になるはず。

上のプログラムはここでも公開中。