SHOEISHA iD

※旧SEメンバーシップ会員の方は、同じ登録情報(メールアドレス&パスワード)でログインいただけます

DeveloperZine(デベロッパージン)- エンジニアの意思決定を支える技術情報メディア ProductZine

CodeZine編集部では、現場で活躍するデベロッパーをスターにするためのカンファレンス「Developers Summit」や、エンジニアの生きざまをブーストするためのイベント「Developers Boost」など、さまざまなカンファレンスを企画・運営しています。

特集記事

<numeric>数値演算アルゴリズムひとめぐり

partial_sum(部分和)とadjacent_difference(隣接差)

 partial_sumはaccumulateのバリエーション、accumulateがX[N]の総和を求めるのに対し、partian_sumはY[i]=X[0]からX[i]までの和となるよう、Y[i]を埋めてくれます。

template<class InputIterator, OutputIterator>
OutputIterator partial_sum(InputIterator first, InputIterator last,
                           OutputIterator result);
template<class InputIterator, OutputIterator,
         class BinaryOperation>
OutputIterator partial_sum(InputIterator first, InputIterator last,
                           OutputIterator result
                           BinaryOperation binary_op);

 adjacent_differenceはpartial_sumの逆、Y[i] = X[i] - X[i-1]、つまり直前の値との差を求めます(ただしY[0] = X[0])。

 ではサンプル、有理数列の和1/2 + 1/4 + 1/8 + ……が次第に1に近づく様子を観察します。

list-9
#include <numeric>

#include <iostream>
#include <vector>
#include <iterator>

#include "sequence.h"
#include "rational.h"

template<typename T>
std::ostream& operator<<(std::ostream& stream, const rational<T>& r) 
  { return stream <<  right << r.first << '/' << left << r.second; }

using namespace std;

// 1/2 + 1/4 + 1/8 + 1/16 …… は 次第に 1 に近づく

int main() {
  const size_t N = 16;

  vector<double> c(N);
  iota(begin(c), end(c), sequence<double>(0.5, [](double x) { return x*0.5;}));

  // 第i項までの和を dout に求める
  vector<double> dout;
  partial_sum(begin(c), end(c), back_inserter(dout));
  for ( auto item : dout ) {
    cout << item << endl;
  }

  // 同じことを有理数について行う
  vector<rational<>> r(N, make_rational(0L));
  // r[] = 1/2, 1/4, 1/8 ……
  iota(begin(r), end(r),
       sequence<rational<>>(make_rational(1L,2L), 
       [](const rational<>& x) 
         { return make_rational(numerator(x),denominator(x)*2L); }));

  // 第i項までの和を rout に求める
  vector<rational<>> rout;
  partial_sum(begin(r), end(r), back_inserter(rout), 
              [](const rational<>& a, const rational<>& b) // 有理数の和 
                { return make_rational(numerator(a)*denominator(b) + numerator(b)*denominator(a),
                                       denominator(a)*denominator(b)); });
  for ( auto item : rout ) {
    cout << left << setw(10) << dout[i] << setw(10) << rout[i] << endl;
  }

}
fig-3
fig-3

 以前紹介した並列STLにも、accumulate,partial_sumに相当する並列アルゴリズムが導入されています。それが reduceとinclusive_scan、さらにinclusive_scanの姉妹品(?)exclusive_scan。inclusive_scanとexclusive_scanとはちょっとだけ違っていて、入力をX[],得られた部分和(の列)をY[]とすると、

 inclusive_scanだと、

Y[0] = X[0]
Y[1] = X[0] + X[1]
...
Y[i] = X[0] + X[1] + .... + X[i]

 対してexclusive_scanは、

Y[0] = 0
Y[1] = 0 + X[0]
Y[2] = 0 + X[0] + X[1]
...
Y[i] = 0 + X[0] + X[1] + .... + X[i]

 つまり、exclusive_scanのi番目の結果Y[i]はX[i]を含まず、X[i-1]までの和になってます。

 reduce, inclusive_scan(exclusive_scan)は複数のスレッドで並列に処理を行います。複数の演算が同時に行われるため、演算の順序が一意に定まりません。従ってreduce, inclusive_scan(exclusive_scan)内で呼び出されるoperator+, operator*および二項関数オブジェクト:binary_opが、

 交換則:

x + y = y + x
x * y = y * x
binary_op(x,y) = binary_op(y,x)

 結合則:

(x + y) + z = x + (y + z)
(x * y) * z = x * (y * z)
binary_op(binary_op(x,y),z) = binary_op(x,binary_op(y,z))

を満たさないと結果の一意性が保障されません。

 exclusive_scanを使ったちょっと面白いサンプル:移動平均(simple moving average)を紹介します。

 exclusive_scanの結果Y[i]はX[0]からX[i-1]までの和、Y[i]からwだけ離れたY[i+w]にはX[0]からX[i+w-1]までの和ですから、その差:Y[i+w]-Y[i]はX[i]からX[i+w-1]までの和となります。これをwで割るとX[i]からw個分の平均が求まるってスンポーです。

list-10 exclusive_scanによる移動平均
#include <experimental/algorithm>
#include <experimental/numeric>

#include <iostream>
#include <algorithm>
#include <random>
#include <vector>
#include <iterator>

using namespace std;
using namespace std::experimental;

int main() {
  
  mt19937 gen; // メルセンヌ・ツイスタ
  normal_distribution<float> dist; // (平均0,標準偏差1の)気温の揺れ
  float t = 25.0f;
  
  const size_t N = 31;
  vector<float> temp; // 日々の気温データ

  generate_n(back_inserter(temp), N, [&]() { return t += dist(gen);});

  vector<float> sum; // 部分和
  parallel::exclusive_scan(parallel::par, begin(temp), end(temp), back_inserter(sum), 0.0f);

  vector<float> sma; // simple-moving-average
  size_t w = 3; // window size
  parallel::transform(parallel::par, 
                      begin(sum)+w, end(sum), begin(sum), back_inserter(sma), 
                      [=](float x, float y) { return (x-y)/w; }); 

  for ( unsigned int i = 0U; i < sma.size(); ++i ) {
    cout <<  temp[i] << '\t' << sma[i] << endl;
  }

}

 ……<numeric>が提供してくれている関数はたったの5つながら、従来for-loopでくるくる回していた計算が関数呼び出しイッパツで済むのは快感です。サンプルで示したとおり、統計がらみの計算にはうってつけですね。

この記事は参考になりましたか?

連載通知を行うには会員登録(無料)が必要です。
既に会員の方はを行ってください。
特集記事連載記事一覧

もっと読む

この記事の著者

επιστημη(エピステーメー)

C++に首まで浸かったプログラマ。Microsoft MVP, Visual C++ (2004.01~2018.06) "だった"りわんくま同盟でたまにセッションスピーカやったり中国茶淹れてにわか茶...

※プロフィールは、執筆時点、または直近の記事の寄稿時点での内容です

この記事は参考になりましたか?

この記事をシェア

CodeZine(コードジン)
https://codezine.jp/article/detail/8778 2015/07/09 14:00

イベント

CodeZine編集部では、現場で活躍するデベロッパーをスターにするためのカンファレンス「Developers Summit」や、エンジニアの生きざまをブーストするためのイベント「Developers Boost」など、さまざまなカンファレンスを企画・運営しています。

新規会員登録無料のご案内

  • ・全ての過去記事が閲覧できます
  • ・会員限定メルマガを受信できます

メールバックナンバー