#ifndef MTGP64EX_H
#define MTGP64EX_H

#include <stdint.h>
#include <string>
#include "mtgp_param.h"

namespace mtgp {
    using namespace std;

    /**
     * @class mtgp64ex
     * @brief Mersenne Twister for Graphic Processors 64-bit executor
     */
    class mtgp64ex {
    public:
	mtgp64ex(const mtgp_param<uint64_t> & param, uint64_t seed);
	mtgp64ex(const mtgp_param<uint64_t> & param,
		 const uint64_t *seed, int size);
	mtgp64ex(const mtgp_param<uint64_t> & param,
		 const string& seed);
	mtgp64ex(const mtgp_param<uint64_t> & param);
	~mtgp64ex();
	uint64_t next() {
	    idx = (idx + 1) & array_size_mask;
	    uint64_t r = recur(array[(array_size - status_size + idx)
				   & array_size_mask],
			       array[(array_size - status_size + idx + 1)
				     & array_size_mask],
			       array[(array_size - status_size + idx + pos)
				     & array_size_mask]);
	    array[idx] = r;
	    r = temper(r, array[(array_size - status_size + idx + pos - 1)
				& array_size_mask]);
	    return r;
	};
	double next_float() {
	    return next_float_close1_open2() - 1.0;
	};
	double next_float_close_open() {
	    return next_float();
	};
	double next_float_open_close() {
	    return 2.0 - next_float_close1_open2();
	};
	double next_float_open_open() {
	    idx = (idx + 1) & array_size_mask;
	    uint64_t r = recur(array[(array_size - status_size + idx)
				   & array_size_mask],
			       array[(array_size - status_size + idx + 1)
				     & array_size_mask],
			       array[(array_size - status_size + idx + pos)
				     & array_size_mask]);
	    array[idx] = r;
	    double f
		= float_temper_open_open(
		    r,
		    array[(array_size - status_size + idx + pos - 1)
			  & array_size_mask]);
	    return f;
	};
	double next_float_close1_open2() {
	    idx = (idx + 1) & array_size_mask;
	    uint64_t r = recur(array[(array_size - status_size + idx)
				   & array_size_mask],
			       array[(array_size - status_size + idx + 1)
				     & array_size_mask],
			       array[(array_size - status_size + idx + pos)
				     & array_size_mask]);
	    array[idx] = r;
	    double f = float_temper(
		r,
		array[(array_size - status_size + idx + pos - 1)
		      & array_size_mask]);
	    return f;
	};
	void reseed(uint64_t seed);
	void reseed(const std::string& seed);
	void reseed(const uint64_t *seed, int size);
	int get_mexp() const;
	const std::string& get_sha1() const;
	void print(ostream& os);
    private:
	uint64_t recur(uint64_t p1, uint64_t p2, uint64_t p3) {
	    uint64_t x = (p1 & mask) ^ p2;
	    uint64_t y = p3;
	    uint64_t r;
	    uint32_t xh = x >> 32;
	    uint32_t xl = x & UINT32_C(0xffffffff);
	    uint32_t yh = y >> 32;
	    uint32_t yl = y & UINT32_C(0xffffffff);
	    xh ^= xh << sh1;
	    xl ^= xl << sh1;
	    yh = xl ^ (yh >> sh2);
	    yl = xh ^ (yl >> sh2);
	    r = ((uint64_t)yh << 32) | yl;
	    r ^= rec_tbl[yl & 0x0f];
	    return r;
	}
	uint64_t temper(uint64_t v, uint64_t h) {
	    h = h ^ (h >> 16);
	    h = h ^ (h >> 8);
	    return v ^ tmp_tbl[h & 0x0f];
	};
	double float_temper(uint64_t v, uint64_t h) {
	    union {
		uint64_t u;
		double f;
	    } w;
	    h = h ^ (h >> 16);
	    h = h ^ (h >> 8);
	    w.u = (v >> 12) ^ flt_tbl[h & 0x0f];
	    return w.f;
	};
	double float_temper_open_open(uint64_t v, uint64_t h) {
	    union {
		uint64_t u;
		double f;
	    } w;
	    h = h ^ (h >> 16);
	    h = h ^ (h >> 8);
	    w.u = ((v >> 12) ^ flt_tbl[h & 0x0f]) | 1;
	    return w.f - 1.0;
	};
	void setup(const mtgp_param<uint64_t>& param);
	uint64_t *array;
	uint64_t rec_tbl[16];
	uint64_t tmp_tbl[16];
	uint64_t flt_tbl[16];
	uint64_t mask;
	int mexp;
	int sh1;
	int sh2;
	int pos;
	int status_size;
	int array_size;
	int array_size_mask;
	int idx;
	std::string sha1;
    };
}

#endif
