debias_bits.h

#pragma once
#include <cstdint>

namespace rng::utils {
	struct debiased_u64 {
		uint64_t bits = {};
		uint32_t bit_count = {};
	};

	struct debiased_bool {
		bool bit = {};
		bool bit_count = {};
	};

	debiased_bool debias_bit(bool bit_0, bool bit_1) {
		debiased_bool result = {};
		result.bit_count = bit_0 != bit_1;
		result.bit = result.bit_count && bit_0;
		return result;
	}

	// vonn Neumann's debiasing coin method
	// Consider a fair coin toss eg:
	// where we have 50% probability of flipping heads and
	// where we have 50% probability of flipping tails
	// An unfair (weighted) coin toss being
	// where the probability of tossing heads or tails is not 50%
	// We would like a way of "debiasing" a biased coin assuming we're handed
	// a coin where we do not know whether it is biased or not
	// 
	// Given the probability of heads (p) then the probability of tails is (1-p)
	// Note the probabilities of results when tossing the coin twice
	// 
	// Heads Heads -> p * p
	// Heads Tails -> p * (1-p)
	// Tails Heads -> (1-p) * p
	// Tails Tails -> (1-p) * (1-p)
	// 
	// We can see that Heads Tails and Tails Heads have equal probability of occuring while
	// Heads Heads and Tails Tails are both different
	// 
	// Method: toss the (potentially) weighted coin twice
	// if we see Heads Heads or Tails Tails, toss the coin twice again
	// Map Heads Tails result to Heads*
	// Map Tails Heads result to Tails*
	// 
	// Furthur reading: https://www.eecs.harvard.edu/~michaelm/coinflipext.pdf
	
	// Apply von Neuman's trick assuming the unsigned integer represents a series of coin tosses
	// debiased_u64.bits -> the resulting "coin tosses"
	// debiased_u64.bit_count -> the # of "debiased coin tosses" detected
	constexpr debiased_u64 debias_bit_positions(uint64_t data_0, uint64_t data_1) noexcept {
		debiased_u64 result = {};
		
		uint64_t xor_result = data_0 ^ data_1;

		for (uint32_t i = 0; i < (8 * sizeof(uint64_t)); i++) {
			uint64_t flag = (1ull << i) & xor_result;
			uint64_t bit_flag = (1ull << i) & data_0;

			uint64_t keep_bit = flag > 0;
			uint64_t bit = bit_flag > 0;

			result.bits = result.bits | ((bit << result.bit_count) * keep_bit);
			result.bit_count += keep_bit;
		}

		return result;
	}

	// used in debias_bit_positions_interleaved, skips 0'd bits from bitwise & operation
	constexpr debiased_u64 debias_bit_positions_skip(uint64_t data_0, uint64_t data_1, uint32_t start) noexcept {
		debiased_u64 result = {};

		uint64_t xor_result = data_0 ^ data_1;

		for (uint32_t i = start; i < (8 * sizeof(uint64_t)); i+=2) {
			uint64_t flag = (1ull << i) & xor_result;
			uint64_t bit_flag = (1ull << i) & data_0;

			uint64_t keep_bit = flag > 0;
			uint64_t bit = bit_flag > 0;

			result.bits = result.bits | ((bit << result.bit_count) * keep_bit);
			result.bit_count += keep_bit;
		}

		return result;
	}

	constexpr debiased_u64 debias_bit_positions_interleaved(uint64_t data_interleaved) noexcept {
		// deinterleave even and odd bits
		uint64_t a = (data_interleaved & 0xAAAAAAAAAAAAAAAAULL) >> 1; //odd bits
		uint64_t b = (data_interleaved & 0x5555555555555555ULL); //even bits
		return debias_bit_positions_skip(a, b, 0u);
	}
}