Wyhash & Random Number Generator

rng-example.cpp

#include <iostream>
#include <map>
using namespace std;

int main() {
	u64 const REQS = 100000000;
	u64 const SERVERS = 250;
	f64 const MEAN_COST = 10;
	f64 const RED = 35 / MEAN_COST / (REQS / SERVERS);

	map<i64, i64> hists;
	map<i64, f64> servers, costs;
	for (u64 n = 0; n < REQS; n += 1) {
		i64 uid = i64(rng.gau(10000.0, 2500.0));
		i64 rid = i64(rng.gau(10000.0, 2500.0));
		f64 cost = costs.count(rid) ? costs[rid] : (costs[rid] = rng.gau(MEAN_COST, 2.5));
		servers[Wy::poorrange(Wy::mix(rid, 0), SERVERS)] += cost;
	}
	for (auto [k, c] : servers) hists[i64(c * RED)] += 1;
	// for (auto [k, c] : servers) cout << k << ": " << string(i64(c * RED), '=') << endl;
	for (auto [k, n] : hists) cout << k << ": " << string(n, '=') << endl;
}

wy.hpp

#include <cstring>
#include <cstddef>
#include <cstdint>
#include <cmath>
#include <type_traits>
#include <random>
#include <algorithm>
#include <array>
#include <memory>

using u8 = uint8_t;
using i8 = int8_t;
using u16 = uint16_t;
using i16 = int16_t;
using u32 = uint32_t;
using i32 = int32_t;
using u64 = unsigned long long;
using i64 = signed long long;
using u128 = __uint128_t;
using i128 = __int128_t;
using usize = size_t;
using isize = typename std::make_signed_t<usize>;
using f32 = float;
using f64 = double;
using f80 = long double;


template<class Rng, class Iter>
static inline void shuffle(Rng rng, Iter const a, Iter const b) {
	auto const n = std::distance(a, b);
	for (auto i = n; i >= 2; i -= 1)
		std::swap(a[i - 1], a[rng(i)]);
}

#ifndef WYHASH_LITTLE_ENDIAN
	#if defined(_WIN32) || defined(__LITTLE_ENDIAN__) || (defined(__BYTE_ORDER__) && __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__)
		#define WYHASH_LITTLE_ENDIAN 1
	#elif defined(__BIG_ENDIAN__) || (defined(__BYTE_ORDER__) && __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__)
		#define WYHASH_LITTLE_ENDIAN 0
	#else
		#warning could not determine endianness! Falling back to little endian.
		#define WYHASH_LITTLE_ENDIAN 1
	#endif
#endif

namespace WyTraits {
	template<class T, class U> concept StrictType =
		std::same_as<std::remove_cvref_t<T>, std::remove_cvref_t<U>>;

	template<typename T> concept NotPointer = 
		std::is_pointer_v<T> == false;

	template<typename T> concept Container = requires(T v) {
		std::distance((u8 const*)std::data(v), (u8 const*)std::next(std::data(v), std::size(v)));
	};

	template<typename T> concept Trivial = 
		NotPointer<T> && std::is_trivial_v<T> && not Container<T>;
};

struct Wy {
	static constexpr std::array<u64, 4> secrets = {
		0xA0761D6478BD642FULL,
		0xE7037ED1A0B428DBULL,
		0x8EBC6AF09C88C6E3ULL,
		0x589965CC75374CC3ULL,
	};
	static constexpr f64 gau_norm = 1.0 / (u64(1) << 20);
	static constexpr f64 uni_norm = 1.0 / (u64(1) << 52);
	static constexpr u64 CL = 0xA0761D6478BD642FULL;
	static constexpr u64 CH = 0xE7037ED1A0B428DBULL;
	u64 g;

	Wy(): g(std::random_device{}()) {}
	Wy(u64 const seed): g(seed) {}

	static inline u64 min() {
		return std::numeric_limits<u64>::min();
	}
	static inline u64 max() {
		return std::numeric_limits<u64>::max();
	}

	static inline u64 _mix(u64 const l, u64 const h) {
		u128 const v = l * u128(h);
		return u64(v) ^ u64(v >> 64);
	}
	static inline u64 mix(u64 const l, u64 const h) {
		return _mix(l ^ CL, h ^ CH);
	}
	static inline u64 mix(u64 const x) {
		return _mix(x, 0);
	}
	static inline u64 mix2(u64 const l, u64 const h) {
		u128 const v = (l ^ CL) * u128(h ^ CH);
		return mix(u64(v), u64(v >> 64));
	}
	static inline u64 mix2(u64 const x) {
		return mix2(x, 0);
	}

	static inline u64 poorrange(u64 const x, u64 const n) {
		return x * u128(n) >> 64;
	}

	#if (WYHASH_LITTLE_ENDIAN)
		static inline u64 rot8(u8 const* const p) { u64 v; memcpy(&v, p, 8); return v; }
		static inline u64 rot4(u8 const* const p) { u32 v; memcpy(&v, p, 4); return v; }
	#elif defined(__GNUC__) || defined(__INTEL_COMPILER) || defined(__clang__)
		static inline u64 rot8(u8 const* const p) { u64 v; memcpy(&v, p, 8); return __builtin_bswap64(v); }
		static inline u64 rot4(u8 const* const p) { u32 v; memcpy(&v, p, 4); return __builtin_bswap32(v); }
	#elif defined(_MSC_VER)
		static inline u64 rot8(u8 const* const p) { u64 v; memcpy(&v, p, 8); return _byteswap_uint64(v); }
		static inline u64 rot4(u8 const* const p) { u32 v; memcpy(&v, p, 4); return _byteswap_ulong(v); }
	#endif
	static inline u64 rot3(u8 const* const p, usize const k) { return (u64(p[0]) << 16) | (u64(p[k >> 1]) << 8) | p[k - 1]; }

	#define ADD_HASH_WRAPPERS(hash) \
		template<class Iter> \
		static inline u64 hash(Iter const begin, Iter const end, u64 seed = 0) { \
			u8 const* p = (u8 const*)std::addressof(*begin); \
			usize const n = std::distance(p, (u8 const*)std::next(std::addressof(*begin), std::distance(begin, end))); \
			return hash(p, n, seed); \
		} \
		template<class T> \
		static inline u64 hash(T const& cont, u64 seed = 0) requires WyTraits::Container<T> { \
			u8 const* p = (u8 const*)std::data(cont); \
			usize const n = std::distance(p, (u8 const*)std::next(std::data(cont), std::size(cont))); \
			return hash(p, n, seed); \
		} \
		template<class T> \
		static inline u64 hash(T data, u64 seed = 0) requires WyTraits::Trivial<T> { \
			return hash(std::addressof(data), sizeof(data), seed); \
		}

	static inline u64 hash(void const* rawptr, usize const n, u64 seed = 0) {
		u8 const* p = (u8 const*)rawptr;
		seed ^= _mix(seed ^ secrets[0], secrets[1]);
		u64 l, h;
		if (n <= 16) [[likely]] {
			if (n >= 4) {
				h = (rot4(p + n - 4) << 32) | rot4(p + n - 4 - ((n >> 3) << 2));
				l = (rot4(p) << 32) | rot4(p + ((n >> 3) << 2));
			} else {
				h = 0;
				if (n > 0) [[likely]] {
					l = rot3(p, n);
				} else {
					l = 0;
				}
			}
		} else {
			usize m = n; 
			if(m > 48) [[unlikely]] {
				u64 s1 = seed, s2 = seed;
				do [[likely]] {
					seed = _mix(rot8(p) ^ secrets[1], rot8(p + 8) ^ seed);
					s1 = _mix(rot8(p + 16) ^ secrets[2], rot8(p + 24) ^ s1);
					s2 = _mix(rot8(p + 32) ^ secrets[3], rot8(p + 40) ^ s2);
					p += 48;
					m -= 48;
				} while(m > 48);
				seed ^= s1 ^ s2;
			}
			if (m > 16) [[unlikely]] {
				seed = _mix(rot8(p) ^ secrets[1], rot8(p + 8) ^ seed);
				m -= 16;
				p += 16;
				if (m > 16) [[unlikely]] {
					seed = _mix(rot8(p) ^ secrets[1], rot8(p + 8) ^ seed);
					m -= 16;
					p += 16;
				}
			}
			l = rot8(p + m - 16);
			h = rot8(p + m - 8);
		}
		u128 const v = (l ^ secrets[1]) * u128(h ^ seed);
		return _mix(u64(v) ^ secrets[0] ^ u64(n), u64(v >> 64) ^ secrets[1]);
	}
	ADD_HASH_WRAPPERS(hash);

	static inline u64 poorhash(u8 const* p, usize n, u64 seed = 0) {
		if(n < 4) [[unlikely]] {
			if (!n) [[unlikely]] return 0;
			return ((u32(p[0]) << 16) | (u32(p[n >> 1]) << 8) | p[n-1]) * CL;
		}
		return u64(rot4(p) + rot4(p + (n >> 1) - 2)) * rot4(p + n - 4);
	}
	ADD_HASH_WRAPPERS(poorhash);

	#undef ADD_HASH_WRAPPERS



	inline u64 operator()() {
		g += CL;
		return _mix(g, g ^ CH);
	}

	template<class T> requires WyTraits::StrictType<T, u64>
	inline T operator()(T const n) {
		u128 m = u64(operator()()) * u128(n);
		if (u64(m) < n) {
			u64 t = -n;
			if (t >= n) {
				t -= n;
				if (t >= n) t %= n;
			}
			while (u64(m) < t) m = u64(operator()()) * u128(n);
		}
		return m >> 64;
	}

	template<class T> requires WyTraits::StrictType<T, u32>
	inline T operator()(T const n) {
		u64 m = u32(operator()()) * u64(n);
		if (u32(m) < n) {
			u32 t = -n;
			if (t >= n) {
				t -= n;
				if (t >= n) t %= n;
			}
			while (u32(m) < t) m = u32(operator()()) * u64(n);
		}
		return m >> 32;
	}

	template<class T> requires WyTraits::StrictType<T, i64>
	inline i64 operator()(T const n) { return i64(operator()(u64(n))); }
	template<class T> requires WyTraits::StrictType<T, i32>
	inline i32 operator()(T const n) { return i32(operator()(u32(n))); }

	template<class T> requires WyTraits::StrictType<T, f64>
	inline T operator()(T const n) {
		u64 const r = operator()();
		return n * (uni_norm * (r >> 12));
	}

	template<class T>
	inline T operator()(T const start, T const end) {
		return operator()(end - start) + start;
	}

	#define ZIGGURAT_BLOCK_COUNT  128
	#define ZIGGURAT_R            3.442619855899
	#define ZIGGURAT_AREA         0.00991256303526217
	f64 gau() {
		static constexpr f64 xtbl[ZIGGURAT_BLOCK_COUNT + 1] = {
			3.7130862467425505002,  3.4426198558990002141,  3.2230849845811415655,
			3.0832288582168683178,  2.9786962526477802626,  2.8943440070215289417,
			2.8231253505489104505,  2.7611693723871768569,  2.7061135731218195488,
			2.6564064112613596791,  2.6109722484318473867,  2.5690336259249377804,
			2.5300096723888274575,  2.4934545220953721056,  2.4590181774118304858,
			2.4264206455337498092,  2.3954342780110624567,  2.3658713701176385946,
			2.3375752413392367757,  2.3104136836987629877,  2.2842740596774717687,
			2.2590595738691980898,  2.2346863955909785915,  2.2110814088787025256,
			2.1881804320760482874,  2.1659267937489210532,  2.1442701823603944611,
			2.1231657086739761375,  2.1025731351892376075,  2.0824562379920159572,
			2.0627822745083079781,  2.0435215366550671945,  2.02464697337738464,
			2.0061338699634712057,  1.9879595741276190335,  1.970103260854325633,
			1.9525457295535557645,  1.9352692282966217352,  1.9182573008645087409,
			1.901494653105150201,   1.8849670357077579208,  1.8686611409944875462,
			1.85256451172808978,    1.836665460258444682,   1.820952996596124418,
			1.8054167642192271437,  1.7900469825998572837,  1.7748343955860681476,
			1.759770224899592117,   1.744846128113799022,   1.730054160563729182,
			1.7153867407136660361,  1.7008366185699153039,  1.6863968467791665695,
			1.6720607540975995775,  1.6578219209540228096,  1.6436741568628672194,
			1.6296114794706331175,  1.6156280950431594068,  1.6017183802213763588,
			1.5878768648905743355,  1.5740982160229990416,  1.5603772223661671603,
			1.5467087798599086224,  1.5330878776740417546,  1.5195095847659385591,
			1.5059690368632017154,  1.492461423781352492,   1.4789819769899226198,
			1.4655259573427090736,  1.4520886428892227915,  1.4386653166845615459,
			1.4252512545140580968,  1.4118417124470556967,  1.3984319141310033174,
			1.3850170377326498361,  1.3715922024273405899,  1.3581524543301413122,
			1.3446927517535449681,  1.3312079496656250566,  1.3176927832094120774,
			1.3041418501286148324,  1.2905495919261944504,  1.2769102735601534082,
			1.2632179614546188429,  1.2494664995730662138,  1.2356494832633604375,
			1.2217602305399941631,  1.2077917504159472184,  1.1937367078331260206,
			1.1795873846639857163,  1.1653356361647499995,  1.1509728421488649719,
			1.1364898520131583304,  1.1218769225825397928,  1.1071236475340335836,
			1.092218876907274927,   1.0771506248928932603,  1.0619059636948215974,
			1.0464709007640424776,  1.0308302360681926846,  1.0149673952513273978,
			0.99886423349298048002, 0.98250080351542590229, 0.96585507940114656567,
			0.94890262551130311053, 0.93161619661514749602, 0.91396525102302894616,
			0.89591535258093435434, 0.87742742911291993213, 0.85845684319380943794,
			0.83895221429757360632, 0.81885390670035329563, 0.79809206064405291414,
			0.7765839878947558006,  0.7542306644540515137,  0.73091191064248450804,
			0.70647961133543202283, 0.68074791866914985405, 0.6534786387399702523,
			0.62435859733604526234, 0.59296294247144254452, 0.55869217840817897436,
			0.52065603876205379663, 0.47743783729668198834, 0.42654798635541407714,
			0.36287143109701985866, 0.27232086481394562893, 0,
		};

		f64 x, y, f0, f1;
		while (1) {
			const u64 u = operator()();
			const u64 idx = u & (ZIGGURAT_BLOCK_COUNT - 1);
			const f64 uf = (2.0 * (uni_norm * (u >> 12)) - 1.0) * xtbl[idx];

			if (std::fabs(uf) < xtbl[idx + 1]) return uf;

			if (idx == 0) {
				do {
					x = std::log(operator()(1.0 / ZIGGURAT_R));
					y = std::log(operator()(1.0));
				} while (-(y + y) < x * x);
				return uf < 0 ? x - ZIGGURAT_R : ZIGGURAT_R - x;
			}

			y = uf * uf;
			f0 = std::exp(-0.5 * (xtbl[idx]     * xtbl[idx]     - y));
			f1 = std::exp(-0.5 * (xtbl[idx + 1] * xtbl[idx + 1] - y));
			if (operator()(f1, f0) < 1.0) return uf;
		}
	#undef RNG2F
	}
	inline f64 gau(f64 const mean, f64 const std) {
		return gau() * std + mean;
	}


	inline u64 poorrange(u64 const n) {
		return poorrange(operator()(), n);
	}
	inline f64 clt() {
		u64 const r = operator()();
		return ((r & 0x1FFFFF) + ((r >> 21) & 0x1FFFFF) + ((r >> 42) & 0x1FFFFF)) * gau_norm - 3.0;
	}
	inline f64 clt(f64 const mean, f64 const std) {
		return clt() * std + mean;
	}
};

thread_local static inline Wy rng;