libsimple_support/unit_tests/math.cpp

// TODO: need a test file per header and include the header first, to detect missing includes
#define SIMPLE_SUPPORT_ALGORITHM_NUMERIC_CLEVER_REPEAT

#include "simple/support/math.hpp"
#include "simple/support/random.hpp"

#include <cassert>
#include <random>
#include <iostream>
#include <iomanip>

#include <chrono>

using namespace simple::support;
using namespace std::chrono;

void FloatEquality()
{
	// ummm how exactly do i do this?

	auto seed = std::random_device{}();
	std::cout << "Float equality test seed: " << std::hex << std::showbase << seed << std::endl;
	random::engine::tiny<unsigned long long> random{seed};
	std::uniform_real_distribution<double> number{};

	for(int i = 0; i < 10'000'000; ++i)
	{
		double a = number(random);
		double b = number(random);

		// someting like this??
		assert( almost_equal(a/b * b, a) );
		assert( almost_equal(b/a * a, b) );
		assert( almost_equal(std::sqrt(a)*std::sqrt(a), a) );
		assert( almost_equal(std::sqrt(b)*std::sqrt(b), b) );
		assert(	!almost_equal(a,b) ); // ummmm... well...  i mean... if this fails - you win! =)
	};

}

void FloatTruncation()
{
	auto seed = std::random_device{}();
	std::cout << "Float trunc test seed: " << std::hex << std::showbase << seed << std::endl;
	random::engine::tiny<unsigned long long> random{seed};
	std::uniform_real_distribution<long double> sign{-1,1};

	int trunced = 0;
	std::uniform_real_distribution<long double> number{
		std::numeric_limits<long double>::min(),
		(long double)std::numeric_limits<std::intmax_t>::max(),
	};
	for(int i = 0; i < 10'000'000; ++i)
	{
		long double a = std::copysign(number(random), sign(random));
		assert( trunc_up_to_intmax(a) == std::trunc(a) );
		if(trunc_up_to_intmax(a) != a)
			++trunced;
	};
	assert(trunced > 0);

	trunced = 0;
	std::uniform_real_distribution<long double> bigger_number{
		(long double)std::numeric_limits<std::intmax_t>::max(),
		std::numeric_limits<long double>::max(),
	};
	for(int i = 0; i < 10'000'000; ++i)
	{
		long double a = std::copysign(bigger_number(random), sign(random));
		assert( trunc_up_to_intmax(a) == std::trunc(a) );
		if(trunc_up_to_intmax(a) != a)
			++trunced;
	};
	assert(0 == trunced);

}

void BabelonianSquareRoot()
{
	auto seed = std::random_device{}();
	std::cout << "Babelonian square root test seed:     " << std::hex << std::showbase << seed << std::dec;
	random::engine::tiny<unsigned long long> random{seed};
	std::uniform_real_distribution<double> number{};

	double good = 0;
	double bad = 0;
	int awful = 0;
	auto then = high_resolution_clock::now();
	for(int i = 0; i < 10'000'000; ++i)
	{
		double n = number(random);
		if((std::abs(babelonian_root2_f(n) - std::sqrt(n)) <
			(std::sqrt(n) * std::numeric_limits<double>::epsilon())))
			++good;
		else
			++bad;

		if(!almost_equal(babelonian_root2_f(n) , std::sqrt(n)))
			++awful;
	}
	std::cout << " time: " << (high_resolution_clock::now() - then).count() << std::endl;
	assert( (good / (good+bad)) > 0.99 );
	assert( awful == 0 );

	// edge cases
	assert( std::isnan(babelonian_root2_f(-1.f)) );
	assert( std::isnan(babelonian_root2_f(0.f/0.f)) );

	float x = 0;
	assert( babelonian_root2_f(x) == x );
	x = std::numeric_limits<float>::min()/2;
	assert( babelonian_root2_f(x) == x );
	x = std::numeric_limits<float>::max()*2;
	assert( babelonian_root2_f(x) == x );
}

void Babelonian4thRoot()
{
	auto seed = std::random_device{}();
	std::cout << "Babelonian 4th root test seed:        " << std::hex << std::showbase << seed << std::dec;
	random::engine::tiny<unsigned long long> random{seed};
	std::uniform_real_distribution<double> number{};

	double good = 0;
	double bad = 0;
	int awful = 0;
	auto then = high_resolution_clock::now();
	for(int i = 0; i < 10'000'000; ++i)
	{
		double n = number(random);
		if((std::abs(babelonian_root2_f(babelonian_root2_f(n)) - std::pow(n,1.0/4)) <
			(std::pow(n,1.0/4) * std::numeric_limits<double>::epsilon())))
			++good;
		else
			++bad;

		if(!almost_equal(babelonian_root2_f(babelonian_root2_f(n)) , std::pow(n,1.0/4)))
			++awful;
	}
	std::cout << " time: " << (high_resolution_clock::now() - then).count() << std::endl;
	assert( (good / (good+bad)) > 0.99 );
	assert( awful == 0 );

	// edge cases
	assert( std::isnan(babelonian_root2_f(-1.f)) );
	assert( std::isnan(babelonian_root2_f(0.f/0.f)) );

	float x = 0;
	assert( babelonian_root2_f(x) == x );
	x = std::numeric_limits<float>::min()/2;
	assert( babelonian_root2_f(x) == x );
	x = std::numeric_limits<float>::max()*2;
	assert( babelonian_root2_f(x) == x );
}

template <unsigned Degree>
void NewtonianRoot()
{
	double power = 1.0/Degree;
	auto seed = std::random_device{}();
	std::cout << "Newtonian root of degree " << std::dec << Degree <<  " test seed: " << std::hex << std::showbase << seed << std::dec;
	random::engine::tiny<unsigned long long> random{seed};
	std::uniform_real_distribution<double> number{};

	// with all the powers involved newton sometimes won't even terminate with default definition check, so instead have to use a convergence check
	// TODO: figure out what's wrong with the default, or make this default
	auto cnvrg = [prev = 0.0](auto guess) mutable
	{
		auto ret = almost_equal(prev, guess);
		prev = guess;
		return ret;
	};

	double good = 0;
	double bad = 0;
	int awful = 0;
	constexpr int count = 10'000'000;
	auto then = high_resolution_clock::now();
	for(int i = 0; i < count; ++i)
	{
		double n = number(random);
		if((std::abs(newtonian_root_f<Degree>(n,cnvrg) - std::pow(n,power)) <
			(std::pow(n,power) * std::numeric_limits<double>::epsilon())))
			++good;
		else
			++bad;

		if(!almost_equal(newtonian_root_f<Degree>(n,cnvrg) , std::pow(n,power)))
			++awful;
	}
	std::cout << " time: " << (high_resolution_clock::now() - then).count() << std::endl;
	assert( (good / (good+bad)) > 0.99 );
	assert( awful == 0 );

	// edge cases
	assert( (Degree % 2 == 1) || std::isnan(newtonian_root_f<Degree>(-1.f)) );
	assert( std::isnan(newtonian_root_f<Degree>(0.f/0.f)) );

	float x = 0;
	assert( newtonian_root_f<Degree>(x) == x );
	x = std::numeric_limits<float>::min()/2;
	assert( newtonian_root_f<Degree>(x) == x );
	x = std::numeric_limits<float>::max()*2;
	assert( newtonian_root_f<Degree>(x) == x );
}

constexpr bool Constexpr()
{
	void(root2(2.0));
	return true;
}

int main()
{
	FloatEquality();
	FloatTruncation();
	BabelonianSquareRoot();
	Babelonian4thRoot();
	NewtonianRoot<1>();
	NewtonianRoot<2>();
	NewtonianRoot<3>();
	NewtonianRoot<4>();
	NewtonianRoot<5>();
	static_assert(Constexpr(),"");
	return 0;
}