#ifndef SEARCH_TEMPER_HPP
#define SEARCH_TEMPER_HPP
/**
 * @file search_temper.hpp
 *
 * @brief search tempering parameters so that the random number
 * generator has a good equidistribution property.
 *
 * @author Mutsuo Saito (Hiroshima University)
 * @author Makoto Matsumoto (Hiroshima University)
 *
 * Copyright (C) 2010 Mutsuo Saito, Makoto Matsumoto and
 * Hiroshima University. All rights reserved.
 *
 * The new BSD License is applied to this software, see LICENSE.txt
 */
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <unistd.h>
#include <tr1/memory>
#include "linear_generator.hpp"
#include "shortest_basis.hpp"

namespace MTToolBox {
    /**
     * @class search_temper
     *
     * @brief search tempering parameters so that the random number
     * generator has a good equidistribution property.
     *
     * This class is expected to be a strategy pattern.
     *
     * @tparam G generator class.
     * \b G must have get_mexp() method,
     * \b G must have set_temper_bit(int pos, int param_pos) and
     * reset_temper_bit(int pos, int param_pos).
     * \b G must have setup_temper().
     * and \b G must be acceptable by shortest_basis class.
     * @tparam T type of the output of \b G.
     * @tparam bit_len bit length of a tempering parameter.
     * @tparam param_num number of tempering parameters.
     * @tparam try_bit_len bit length of the output of \b G to
     * be tempered. Only \b try_bit_len MSBs of the output are tempered.
     * @tparam step Equidistribution property of v-bit accuracy is calculated
     * each \b step bits pattern of tempering parameter. Large \b step will
     * cause low performance.
     */
    template<typename G, typename T,
	     int bit_len, int param_num, int try_bit_len, int step>
    class search_temper {
    public:
	/**
	 * search tempering parameters.
	 * @returns 0
	 */
	int operator()(G& rand, bool verbose = false) {
	    using namespace std;
	    if (verbose) {
		cout << "searching..." << endl;
	    }
	    int delta = 0;
	    for (int p = 0; p < try_bit_len; p += step) {
		int max_depth = p + step;
		if (max_depth > try_bit_len) {
		    max_depth = try_bit_len;
		}
		for (int i = 0; i < param_num; i++) {
		    delta = search_best_temper(rand, p, i, max_depth, verbose);
		}
	    }
	    if (verbose) {
		cout << "delta = " << delta << endl;
	    }
	    return 0;
	};
    private:
	/**
	 * search for one tempering parameter. generate all bit pattern
	 * for v-bit to max_v_bit of the tempering parameter, and calculate
	 * equidistribution property, then select the best bit pattern.
	 * This function calls itself recursively.
	 *
	 * @param rand linear generator
	 * @param v_bit the bit of tempering parameter is searched
	 * @param param_pos position of current tempering parameter.
	 * @param max_v_bit search stops at this bit.
	 * @param verbose if true output searching process
	 * @return delta an integer which shows equidistribution property.
	 */
	int search_best_temper(G& rand, int v_bit,
			       int param_pos, int max_v_bit, bool verbose) {
	    using namespace std;
	    int mexp = rand.get_mexp();
	    int delta1;
	    int delta2;
	    int eq;
	    int sum = 0;
	    rand.reset_temper_bit(bit_len - v_bit, param_pos);
	    if (v_bit >= max_v_bit) {
		rand.setup_temper();
		eq = get_equidist(rand, v_bit, &sum);
		delta1 = sum + mexp / v_bit - eq;
	    } else {
		delta1 = search_best_temper(rand, v_bit + 1, param_pos,
					    max_v_bit, verbose);
	    }
	    rand.set_temper_bit(bit_len - v_bit, param_pos);
	    if (v_bit >= max_v_bit) {
		rand.setup_temper();
		eq = get_equidist(rand, v_bit, &sum);
		delta2 = sum + mexp / v_bit - eq;
	    } else {
		delta2 = search_best_temper(rand, v_bit + 1, param_pos,
					    max_v_bit, verbose);
	    }
	    if (verbose) {
		cout << param_pos << ":" << v_bit << ":"
		     << delta1 << "," << delta2 << endl;
	    }
	    if (delta1 > delta2) {
		return delta2;
	    } else {
		rand.reset_temper_bit(bit_len - v_bit, param_pos);
		return delta1;
	    }
	};

	/**
	 * wrapper of shortest_basis#get_equidist()
	 *
	 * @param rand linear generator
	 * @param bit_len_ \b bit_len_ MSBs of equidistribution
	 * property of the generator is calculated.
	 * @param sum summation of equidistribution property
	 * from 0 to \b bit_len_ -1 MSBs.
	 * @returns equidistribution property at \b bit_len_-th
	 * MSB.
	 */
	int get_equidist(G& rand, int bit_len_, int *sum) {
	    G r(rand);
	    shortest_basis<G, T> sb(r, bit_len_);
	    return sb.get_equidist(sum);
	};
    };
}

#endif
