diff --git a/lean/Makefile b/lean/Makefile
index 61a42e0..18b51bb 100644
--- a/lean/Makefile
+++ b/lean/Makefile
@@ -1,5 +1,7 @@
 SRC_C   = $(wildcard src/*.c)
 SRC_CXX = $(wildcard src/*.cpp)
+SRC_CXX+= $(wildcard src/rng/*.cpp)
+
 INCLUDE = include
 
 PROJECT=shuffle3
@@ -41,6 +43,7 @@ debug: | dirs $(PROJECT)-debug
 
 dirs:
 	@mkdir -p obj/c{,xx}/src
+	@mkdir -p obj/cxx/src/rng
 
 obj/c/%.o: %.c
 	$(CC) -c $< $(CFLAGS) -o $@ $(LDFLAGS)
diff --git a/lean/include/rng.h b/lean/include/rng.h
new file mode 100644
index 0000000..8334d6c
--- /dev/null
+++ b/lean/include/rng.h
@@ -0,0 +1,6 @@
+#ifndef _RNG_H
+#define _RNG_H
+
+// C interface to C++ RNG impls
+
+#endif /* _RNG_H */
diff --git a/lean/include/rng/frng.hpp b/lean/include/rng/frng.hpp
new file mode 100644
index 0000000..7a634a6
--- /dev/null
+++ b/lean/include/rng/frng.hpp
@@ -0,0 +1,17 @@
+
+#include "impl.hpp"
+
+namespace rng
+{
+	struct frng : public RNG 
+	{
+		inline frng(double s1, double s2) : state({s1, s2}){}
+		inline frng(const std::array<double, 2>& ar) : state(ar){}
+		inline frng(std::array<double, 2>&& ar) : state(ar){}
+		inline frng(const double (&ar)[2]) : state({ar[0], ar[1]}) {}
+	protected: 
+		double sample();
+	private:
+		std::array<double, 2> state;
+	}; 
+}
diff --git a/lean/include/rng/impl.hpp b/lean/include/rng/impl.hpp
new file mode 100644
index 0000000..ca10bd8
--- /dev/null
+++ b/lean/include/rng/impl.hpp
@@ -0,0 +1,34 @@
+#pragma once
+
+#include <array>
+#include <vector>
+#include <cstdint>
+
+/// Base class for RNG impls
+struct RNG {
+	virtual unsigned char byte();
+	virtual void bytes(unsigned char* ptr, std::size_t len);
+
+	template<std::size_t N>
+	inline void bytes(unsigned char (&arr)[N]) { return bytes(arr, N); }
+	template<std::size_t N>
+	inline void bytes(std::array<unsigned char, N>& array) { return bytes(&array[0], N); }
+	inline void bytes(std::vector<unsigned char>& vec) { return bytes(&vec[0], vec.size()); }
+
+	bool chance();
+	virtual bool chance(double chance);
+
+	virtual std::int32_t next_int();
+	inline std::int32_t next_int(std::int32_t max) { return next_int(0, max); }
+	std::int32_t next_int(std::int32_t min, std::int32_t max);
+
+	virtual std::int64_t next_long();
+
+       	inline std::int64_t next_long(std::int64_t max) { return next_long(0, max); }
+	std::int64_t next_long(std::int64_t min, std::int64_t max);
+
+	inline float  next_float() { return (float)sample(); }
+	inline double next_double() { return sample(); }
+protected:
+	virtual double sample() = 0;
+};
diff --git a/lean/src/rng.cpp b/lean/src/rng.cpp
new file mode 100644
index 0000000..1fd804b
--- /dev/null
+++ b/lean/src/rng.cpp
@@ -0,0 +1,47 @@
+
+#include <rng/impl.hpp>
+#include <cmath>
+
+inline unsigned char RNG::byte()
+{
+	return (unsigned char)next_int(255);
+}
+
+void RNG::bytes(unsigned char* ptr, std::size_t len)
+{
+	for(std::size_t i=0;i<len;i++)
+	{
+		ptr[i] = byte();
+	}
+}
+
+bool RNG::chance()
+{
+	return next_long() <= 0;
+}
+
+bool RNG::chance(double d)
+{
+	if (d<=0) return false;
+	return sample() <= d;
+}
+
+std::int32_t RNG::next_int(std::int32_t min, std::int32_t max)
+{
+	return (std::int32_t)floor(sample() * (double)(max-min))+min;
+}
+
+std::int32_t RNG::next_int()
+{
+	return next_int(INT32_MIN, INT32_MAX);
+}
+
+std::int64_t RNG::next_long(std::int64_t min, std::int64_t max)
+{
+	return (std::int64_t)floor(sample() * (double)(max-min))+min;
+}
+
+std::int64_t RNG::next_long()
+{
+	return next_long(INT64_MIN, INT32_MAX);
+}
diff --git a/lean/src/rng/frng.cpp b/lean/src/rng/frng.cpp
new file mode 100644
index 0000000..f20cb1f
--- /dev/null
+++ b/lean/src/rng/frng.cpp
@@ -0,0 +1,48 @@
+#include <rng/impl.hpp>
+#include <rng/frng.hpp>
+#include <cmath>
+
+template<std::size_t N>
+constexpr inline double dot(const std::array<double, N>& v, const std::array<double, N>& u)
+{
+	double res=0;
+	for(std::size_t i=0;i<N;i++)
+	{
+		res += v[i] * u[i];
+	}
+	return res;
+}
+
+inline constexpr double fract(double x)
+{
+	return x - floor(x);
+}
+
+namespace rng {
+	inline constexpr double sample_double(const std::array<double, 2>& state)
+	{
+		const constexpr std::array<double, 2> vec2 = {  12.9898, 78.223  };
+		return fract(sin(dot(state, vec2)) * 43758.5453);
+	}
+	inline void update_state(std::array<double, 2>& state, double r)
+	{
+		float v1 = (float)state[0];
+		float v2 = (float)state[1];
+
+		std::array<double, 2> nvec = {
+			r,
+			(double)v2,
+		};
+
+		state[0] = sample_double(nvec);
+
+		nvec[1] = (double)v1;
+		state[1] = sample_double(nvec);
+	}
+	double frng::sample()
+	{
+		double res = sample_double(state);
+		update_state(state, res);
+		return res;
+	}
+}