AlgorithmEquidistribution.hpp

[詳解]

#ifndef MTTOOLBOX_ALGORITHM_EQUIDISTRIBUTION_HPP
#define MTTOOLBOX_ALGORITHM_EQUIDISTRIBUTION_HPP

#if __cplusplus >= 201103L
#include <memory>
#else
#if defined(__clang__)
#if __has_include(<tr1/memory>)
#define MTTOOLBOX_USE_TR1
#include <tr1/memory>
#else
#include <memory>
#endif // __has_indlude
#else  // not clang
#define MTTOOLBOX_USE_TR1
#include <tr1/memory>
#endif // clang
#endif // cplusplus version
#include <stdexcept>
#include <MTToolBox/EquidistributionCalculatable.hpp>
#include <MTToolBox/period.hpp>
#include <MTToolBox/util.hpp>

namespace MTToolBox {
#if defined(MTTOOLBOX_USE_TR1)
    using std::tr1::shared_ptr;
#else
    using std::shared_ptr;
#endif

    template<typename U, typename V = U>
    class linear_generator_vector {
    public:

        typedef EquidistributionCalculatable<U, V> ECGenerator;

        linear_generator_vector<U, V>(const ECGenerator& generator) {
            shared_ptr<ECGenerator> r(generator.clone());
            rand = r;
            //rand->seed(1);
            count = 0;
            zero = false;
            setZero(next);
        }

        linear_generator_vector<U, V>(const ECGenerator& generator,
                                   int bit_pos) {
            shared_ptr<ECGenerator> r(generator.clone());
            rand = r;
            rand->setZero();
            count = 0;
            zero = false;
            next = getOne<U>() << (bit_size<U>() - bit_pos - 1);
#if defined(DEBUG)
            cout << "DEBUG:" << dec << bit_pos << ":"
                 << hex << next << endl;
#endif
        }

        void add(const linear_generator_vector<U, V>& src);
        void next_state(int bit_len);
        void debug_print();

        shared_ptr<ECGenerator> rand;
        int count;

        bool zero;

        U next;
    };

    template<typename U, typename V = U>
    class AlgorithmEquidistribution {

        typedef linear_generator_vector<U, V> linear_vec;

        /*
         *\japanese
         * 均等分布次元計算可能な疑似乱数生成器
         *\endjapanese
         *\english
         * Pseudo random number generator which can calculate dimension
         * of equi-distribution.
         *\endenglish
         */
        typedef EquidistributionCalculatable<U, V> ECGenerator;

    public:

        AlgorithmEquidistribution(const ECGenerator& rand, int bit_length,
                                  int mexp = 0) {
            bit_len = bit_length;
            size = bit_len + 1;
            basis = new linear_vec * [size];
            if (mexp == 0) {
                stateBitSize = rand.bitSize();
            } else {
                stateBitSize = mexp;
            }
            for (int i = 0; i < bit_len; i++) {
                basis[i] = new linear_vec(rand, i);
            }
            basis[bit_len] = new linear_vec(rand);
            basis[bit_len]->next_state(bit_len);
        }

        ~AlgorithmEquidistribution() {
            for (int i = 0; i < size; i++) {
                delete basis[i];
            }
            delete[] basis;
        }

        int get_all_equidist(int veq[]);
        int get_equidist(int *sum_equidist);
    private:
        int get_equidist_main(int bit_len);
        void adjust(int new_len);

        linear_vec **basis;

        int bit_len;

        int stateBitSize;

        int size;
    };

    template<typename U, typename V>
    void AlgorithmEquidistribution<U, V>::adjust(int new_len) {
        using namespace std;
        U tmp;
        setZero(tmp);
        U mask = ~tmp << (bit_size<U>() - new_len);
        for (int i = 0; i < size; i++) {
            basis[i]->next = basis[i]->next & mask;
            if (isZero(basis[i]->next)) {
                basis[i]->next_state(new_len);
            }
        }
    }

#if defined(DEBUG)

    template<typename U, typename V>
    void linear_generator_vector<U, V>::debug_print() {
        using namespace std;

        cout << "debug ====" << endl;
        cout << "count = " << dec << count << endl;
        cout << "zero = " << zero << endl;
        cout << "next = " << hex << next << endl;
        cout << "debug ====" << endl;
    }
#else
    template<typename U, typename V>
    void linear_generator_vector<U, V>::debug_print() {
    }
#endif

    template<typename U, typename V>
    int AlgorithmEquidistribution<U, V>::get_all_equidist(int veq[]) {
        using namespace std;

        int sum = 0;

        veq[bit_len - 1] = get_equidist_main(bit_len);
#if defined(DEBUG)
        for (int i = 0; i < size; i++) {
            basis[i]->debug_print();
        }
#endif
        sum += stateBitSize / bit_len - veq[bit_len - 1];
        bit_len--;
        for (; bit_len >= 1; bit_len--) {
            adjust(bit_len);
            veq[bit_len - 1] = get_equidist_main(bit_len);
            sum += stateBitSize / bit_len - veq[bit_len - 1];
        }
        return sum;
    }

    template<typename U, typename V>
    int AlgorithmEquidistribution<U, V>::get_equidist(int *sum_equidist) {
        using namespace std;

        int veq = get_equidist_main(bit_len);
        int sum = 0;
        bit_len--;
        for (; bit_len >= 1; bit_len--) {
            adjust(bit_len);
            sum += stateBitSize / bit_len - get_equidist_main(bit_len);
        }
        *sum_equidist = sum;
        return veq;
    }

    template<typename U, typename V>
    void linear_generator_vector<U, V>::add(
        const linear_generator_vector<U, V>& src) {
        using namespace std;

        rand->add(*src.rand);
        next ^= src.next;
    }

    template<typename U, typename V>
    void linear_generator_vector<U, V>::next_state(int bit_len) {
        using namespace std;

        if (zero) {
            return;
        }
        int zero_count = 0;
        next = rand->generate(bit_len);
        count++;
        while (isZero(next)) {
            zero_count++;
            if (zero_count > rand->bitSize() * 2) {
                zero = true;
                if (rand->isZero()) {
                    zero = true;
                }
                break;
            }
            next = rand->generate(bit_len);
            count++;
        }
    }

    template<typename U, typename V>
    int AlgorithmEquidistribution<U, V>::get_equidist_main(int v) {
        using namespace std;
        using namespace NTL;
        int bit_len = v;
        int pivot_index;
        int old_pivot = 0;

        pivot_index = calc_1pos(basis[bit_len]->next);
        while (!basis[bit_len]->zero) {
#if defined(DEBUG)
            if (pivot_index == -1) {
                cout << "pivot_index = " << dec << pivot_index << endl;
                cout << "zero = " << basis[bit_len]->zero << endl;
                cout << "next = " << hex << basis[bit_len]->next << endl;
                throw new std::logic_error("pivot error 0");
            }
            if (pivot_index >= bit_len) {
                cout << "pivot_index = " << dec << pivot_index << endl;
                cout << "bit_len = " << bit_len << endl;
                cout << "next = " << hex << basis[bit_len]->next << endl;
                throw new std::logic_error("pivot error 0.1");
            }
            if (pivot_index != calc_1pos(basis[pivot_index]->next)) {
                cerr << "pivot error 1" << endl;
                cerr << "pivot_index:" << dec << pivot_index << endl;
                cerr << "calc_1pos:" << dec
                     << calc_1pos(basis[pivot_index]->next) << endl;
                cerr << "next:" << hex << basis[pivot_index]->next << endl;
                throw new std::logic_error("pivot error 1");
            }
#endif
            // アルゴリズムとして、全部のcount を平均的に大きくしたい。
            // 従って count の小さい方を変化させたい
            if (basis[bit_len]->count > basis[pivot_index]->count) {
                swap(basis[bit_len], basis[pivot_index]);
            }
#if defined(DEBUG)
            cout << "before add bit_len next = " << hex
                 << basis[bit_len]->next << "  count = " << dec
                 << basis[bit_len]->count << endl;
            cout << "before add pivot   next = " << hex
                 << basis[pivot_index]->next << "  count = " << dec
                 << basis[pivot_index]->count << endl;
#endif
            basis[bit_len]->add(*basis[pivot_index]);
#if defined(DEBUG)
            cout << "after  add bit_len next = " << hex
                 << basis[bit_len]->next << "  count = " << dec
                 << basis[bit_len]->count << endl;
#endif
            // add の結果 next の最後の1 は必ず 0 になる。
            // 全部0なら次の状態に進める。（内部でcount が大きくなる）
            if (isZero(basis[bit_len]->next)) {
                basis[bit_len]->next_state(bit_len);
                pivot_index = calc_1pos(basis[bit_len]->next);
#if defined(DEBUG)
                cout << "zero" << endl;
                cout << "pivot_index = " << dec << pivot_index << endl;
                if (pivot_index >= bit_len) {
                    cout << "pivot_index = " << dec << pivot_index << endl;
                    cout << "bit_len = " << bit_len << endl;
                    cout << "next = " << hex << basis[bit_len]->next << endl;
                    throw new std::logic_error("pivot error 1.1");
                }
                if (basis[bit_len]->zero) {
                    cout << "loop exit condition" << endl;
                } else if (pivot_index == -1) {
                    cerr << "pivot error 1.1" << endl;
                    throw new std::logic_error("pivot error 1.2");
                }
#endif
            // 全部0でなければ、pivot_index は小さくなる。
            // pivot_index が 0 になれば最上位bit のみ1なので
            // 次の add で全部0になる。
            } else {
                old_pivot = pivot_index;
                pivot_index = calc_1pos(basis[bit_len]->next);
                if (pivot_index >= bit_len) {
                    cout << "pivot_index = " << dec << pivot_index << endl;
                    cout << "bit_len = " << bit_len << endl;
                    cout << "next = " << hex << basis[bit_len]->next << endl;
                    throw new std::logic_error("pivot error 2");
                }
                if (old_pivot <= pivot_index) {
                    cerr << "pivot error 2" << endl;
                    cerr << "old_pivot = " << dec << old_pivot << endl;
                    cerr << "pivot_index = " << dec << pivot_index << endl;
                    throw new std::logic_error("pivot error 2.1");
                }
            }
        }

        // 計算終了したので最長のベクトルを求める。（長いとはcountが少ないこと）
        int min_count = basis[0]->count;
        for (int i = 1; i < bit_len; i++) {
            if (min_count > basis[i]->count) {
                min_count = basis[i]->count;
            }
        }
#if defined(DEBUG)
        cout << "min_count = " << min_count << endl;
        cout << "stateBitSize = " << stateBitSize << endl;
        cout << "bit_len = " << bit_len << endl;
        cout << "theoretical bound " << (stateBitSize / bit_len) << endl;
#endif
        if (min_count > stateBitSize / bit_len) {
            cerr << basis[0]->rand->getParamString() << endl;
            for(int i = 0; i < size; i++) {
                basis[i]->debug_print();
            }
            cerr << "min_count = " << min_count << endl;
            cerr << "stateBitSize = " << stateBitSize << endl;
            cerr << "bit_len = " << bit_len << endl;
            cerr << "over theoretical bound " << (stateBitSize / bit_len)
                 << endl;
            throw new std::logic_error("over theoretical bound");
        }
        return min_count;
    }
}
#if defined(MTTOOLBOX_USE_TR1)
#undef MTTOOLBOX_USE_TR1
#endif
#endif // MTTOOLBOX_ALGORITHM_EQUIDISTRIBUTION_HPP