物理の駅 Physics station by 現役研究者

テクノロジーは共有されてこそ栄える

2個の数列(std::vector<int>とか)から重複とかを探す

2個の数列から重複等を探す方法を紹介する。Linuxのjoinコマンドを高速化する場合などに使える。

続きは 2個の自作クラスの配列(std::vector<MyClass>)から重複とかを探す - 物理の駅 Physics station by 現役研究者

#include <algorithm>に便利な関数が用意されている。 ソートして、重複を削除した後、set_intersection, set_union, set_differenceを使ってみよう。

参照

Tech Tips: 知ってると便利なSTL(10) set_intersection, set_union, set_difference

乱数発生器についてのコードを修正しました(2016/11/30)

list1 = {5,7,4,2,9}
list2 = {4,3,8,3,2}
↓ソート
list1 = {2,4,5,7,9}
list2 = {2,3,3,4,8}
↓重複を削除
list1 = {2,4,5,7,9}
list2 = {2,3,4,8}
↓比較
積集合 = {2,4}
和集合 = {2,3,4,5,7,8,9}
差集合1 = {5,7,9}
差集合2 = {3,8}
#include <random>
#include <vector>
#include <algorithm>
#include <iostream>
#include <iterator>

void main()
{
    const int Ntrk = 10000000;

    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_int_distribution<> dis(0, Ntrk - 1);
    std::vector<unsigned int> list1, list2;

    for (int i = 0; i < Ntrk; ++i) {
        list1.emplace_back(dis(gen)); //乱数を作成
        list2.emplace_back(dis(gen)); //乱数を作成
    }

    std::sort(list1.begin(), list1.end()); //ソート
    std::sort(list2.begin(), list2.end()); //ソート

    std::cout << list1.size() << " -> ";
    //重複を削除
    list1.erase(std::unique(list1.begin(), list1.end()), list1.end());
    std::cout << list1.size() << std::endl;

    std::cout << list2.size() << " -> ";
    //重複を削除
    list2.erase(std::unique(list2.begin(), list2.end()), list2.end()); 
    std::cout << list2.size() << std::endl;

    std::vector<unsigned int> list_inter; 
    //list1とlist2の両方にある値 (積集合)
    std::set_intersection(list1.begin(), list1.end(), list2.begin(), list2.end(), std::back_inserter(list_inter));

    std::vector<unsigned int> list_union; 
    //list1かlist2のどちらかに存在する値(和集合)
    std::set_union(list1.begin(), list1.end(), list2.begin(), list2.end(), std::back_inserter(list_union));

    std::vector<unsigned int> list_dif1; 
    //list1にあってlist2にない値(差集合)
    std::set_difference(list1.begin(), list1.end(), list2.begin(), list2.end(), std::back_inserter(list_dif1));

    std::vector<unsigned int> list_dif2; 
    //list2にあってlist1にない値(差集合)
    std::set_difference(list2.begin(), list2.end(), list1.begin(), list1.end(), std::back_inserter(list_dif2));

    //結果の確認
    std::cout << "dif1 + dif2 + inter = union" << std::endl;
    printf_s("%u + %u + %u = %u\n", list_dif1.size(), list_dif2.size(), list_inter.size(), list_union.size());
}

δtanθ - tanθグラフをδθ - tanθグラフに変換する

より正しい記述へ

phst.hateblo.jp

tanθを0.1刻みで、δtanθの値が一律で0.1のときにδθを計算するコードの例。

#include <vector>

void main()
{
    std::vector<double> x, y; //tanθ δtanθ
    for (int i = 0; i < 20; i++) {
        x.emplace_back(i * 0.1 + 0.05);
        y.emplace_back(0.1);
    }

    for (int i = 0; i < 20; i++) {
        double dtheta = atan(x.at(i) + y.at(i) / 2) - atan(x.at(i) - y.at(i) / 2); //δθを計算
        printf("%8.5f %8.5f %8.5f\n", x.at(i), y.at(i), dtheta);
    }
}

f:id:onsanai:20160918044320p:plain

OpenMPでfor文を高速化する

openmpでfor文を高速化してみよう。OpenMPを有効にするには、Visual Studioのプロジェクトのプロパティページを開いて、C/C++の言語のOpenMPのサポートをはい(/openmp)にする必要がある。

例として1から50000までの数が素数(prime)かどうかを調べてみよう。素数の調べ方は2以上候補値未満の数で割った余りを調べる愚直な方法で実装した。

#include <omp.h>
#include <vector>
#include <iostream>
#include <chrono>
#include <algorithm>

int main()
{
    auto *mylock = new omp_lock_t;
    omp_init_lock(mylock); //lockを初期化

    std::vector<int> result; //素数の結果

    auto start = std::chrono::system_clock::now(); //開始時間

#pragma omp parallel for num_threads(12)
    for (int i = 0; i < 100; i++) {
        std::vector<int> list; //素数の中間ファイル
        
        for (int j = 0; j < 5000; j++) {
            int candidate = i + j * 100 + 1;
            if (candidate < 2) { continue; }

            bool prime_flag = true;
            for (int k = 2; k < candidate; k++) {
                if ((candidate%k) == 0) { //2以上候補値未満の数で割った余り
                    prime_flag = false; //素数ではない
                    break;
                }
            }
            if (prime_flag) {
                list.emplace_back(candidate);
            }
        }
        omp_set_lock(mylock); //std::vector<T>の変更操作をするためlock
        for (auto p = list.begin(); p != list.end(); p++) {
            result.emplace_back(*p);
        }
        omp_unset_lock(mylock); //lockを解除
    }
    auto end = std::chrono::system_clock::now(); //終了時間

    omp_destroy_lock(mylock); //lock破棄

    std::cout << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() << " [msec]" << std::endl;

    std::cout << "max prime is " << *std::max_element(result.begin(), result.end()) << std::endl; //調査した中で最大の素数を出力
    ::system("pause");
}

当環境で 37秒かかったものが12スレッド並列で行うと6.5秒に短縮された。オーバーヘッドがあるため、単純に1/12にはならない。 std::vector<T>の変更操作はスレッドセーフではないため、omp_lock_tを使ってlockとlock解除を行っている。

素数の調べ方は最適化しておらず、世の中にはもっと高速な方法が沢山ある。一般的に、openmpで高速化するのはアルゴリズムを改善した後のほうが良いだろう。

スレッド数を指定しない場合は以下のように記述することも可能。この場合は全スレッドを使って計算が行われる。

#pragma omp parallel for

Visual Studio 2013でstd::min std::maxが使えなくなった時の対処

std::max - cppreference.com

Visual Studio 2013以降、Visual Studio 2015も完全にalgorithmに移動したらしいので以下の記述が必要。

#include <algorithm> //必要
int main()
{
    int i = std::max(1,2);
    int j = std::min(3,4);
}

ちなみにVisual Studio 2012以下は以下のコードで動く。

#include <xutility> //非推奨
int main()
{
    int i = std::max(1, 2);
    int j = std::min(3, 4);
}

VC10, VC11ではalgorithm 経由で xutilityに記述、VC12, VC14ではalgorithm に直接記述されるようになったようだ。VC12, VC14では、別のincludeファイルからxutilityが参照されている場合は、明示的にalgorithmを追加してやる必要がある。

メモリリークしているソースコードの行を特定する

#include <cstdlib>
#include <new>
#include <memory>

#include <crtdbg.h>

#define _CRTDBG_MAP_ALLOC

#define new ::new(_NORMAL_BLOCK, __FILE__, __LINE__)

void main()
{
    int *i = new int;
    _CrtDumpMemoryLeaks();
    return;
}

これをデバックモードで実行すると

Detected memory leaks!
Dumping objects ->
c:\users\documents\source.cpp(13) : {170} normal block at 0x00000200AD5747B0, 4 bytes long.
 Data: <    > CD CD CD CD 
Object dump complete.

ちゃんとint *i = new int; である13行目でエラーが出ていることが分かる。

参照

プログラマーの友 第八報:メモリリークと crtdbg.h