Compiler Explorer

Source code

//来自github的爱！
//https://github.com/jbapple/128-bit-literals
namespace suffix128 {

// Returns the value of the hexadecimal character
constexpr unsigned __int128 CharValue(char c) {
return (c >= '0' && c <= '9')
			? (c - '0')
			: ((c >= 'a' && c <= 'f') ? (10 + (c - 'a')) : (10 + (c - 'A')));
}

static constexpr unsigned __int128 MAX128 = ~static_cast<unsigned __int128>(0);

// ValidateU128Helper<BASE, c1, c2, ... , cn>(v) returns true iff cn + ... + c2
// * BASE^(n-2) + c1 * BASE^(n-1) + v * BASE^n is a valid 128-bit unsigned
// number when interpreted in base BASE.
template <int BASE>
constexpr bool ValidateU128Helper(unsigned __int128) {
return true;
}

template <int BASE, char C, char... CS>
constexpr bool ValidateU128Helper(unsigned __int128 accumulate) {
return (C == '\'') ? ValidateU128Helper<BASE, CS...>(accumulate)
					: ((accumulate <= MAX128 / BASE) &&
						(BASE * accumulate <= MAX128 - CharValue(C)) &&
						ValidateU128Helper<BASE, CS...>(accumulate * BASE +
														CharValue(C)));
}

// ValidateU128<BASE, c1, c2, ... , cn>(v) returns true iff cn + ... + c2 *
// BASE^(n-2) + c1 * BASE^(n-1) is a valid 128-bit unsigned number when
// interpreted in base BASE.
template <int BASE, char... CS>
constexpr bool ValidateU128() {
return ValidateU128Helper<BASE, CS...>(0);
}

// MakeU128Helper<BASE, c1, c2, ... , cn>(v) returns cn + ... + c2 *
// BASE^(n-2) + c1 * BASE^(n-1) + result * BASE^n.
template <int BASE>
constexpr unsigned __int128 MakeU128Helper(unsigned __int128 result) {
return result;
}

template <int BASE, char C, char... CS>
constexpr unsigned __int128 MakeU128Helper(unsigned __int128 result) {
return MakeU128Helper<BASE, CS...>(
	(C == '\'') ? result : (result * BASE + CharValue(C)));
}

// MakeU128<BASE, c1, c2, ... , cn>(v) returns cn + ... + c2 * BASE^(n-2) + c1 *
// BASE^(n-1).
template <int BASE, char... CS>
constexpr unsigned __int128 MakeU128() {
return MakeU128Helper<BASE, CS...>(0);
}

template <char... CS>
struct StaticU128 {
static constexpr bool IS_VALID = ValidateU128<10, CS...>();
static constexpr unsigned __int128 PAYLOAD = MakeU128<10, CS...>();
};

template <char... CS>
struct StaticU128<'0', 'x', CS...> {
static constexpr bool IS_VALID = ValidateU128<16, CS...>();
static constexpr unsigned __int128 PAYLOAD = MakeU128<16, CS...>();
};

template <char... CS>
struct StaticU128<'0', 'X', CS...> {
static constexpr bool IS_VALID = ValidateU128<16, CS...>();
static constexpr unsigned __int128 PAYLOAD = MakeU128<16, CS...>();
};

template <char... CS>
struct StaticU128<'0', 'b', CS...> {
static constexpr bool IS_VALID = ValidateU128<2, CS...>();
static constexpr unsigned __int128 PAYLOAD = MakeU128<2, CS...>();
};

template <char... CS>
struct StaticU128<'0', 'B', CS...> {
static constexpr bool IS_VALID = ValidateU128<2, CS...>();
static constexpr unsigned __int128 PAYLOAD = MakeU128<2, CS...>();
};

template <char... CS>
struct StaticU128<'0', CS...> {
static constexpr bool IS_VALID = ValidateU128<8, CS...>();
static constexpr unsigned __int128 PAYLOAD = MakeU128<8, CS...>();
};

static constexpr unsigned __int128 SIGNED_LIMIT =
	(static_cast<unsigned __int128>(1) << 127) - 1;

template <char... CS>
struct Static128 {
static constexpr bool IS_VALID =
	ValidateU128<10, CS...>() && (MakeU128<10, CS...>() <= SIGNED_LIMIT);
static constexpr __int128 PAYLOAD =
	static_cast<__int128>(MakeU128<10, CS...>());
};

template <char... CS>
struct Static128<'0', 'x', CS...> {
static constexpr bool IS_VALID =
	ValidateU128<16, CS...>() && (MakeU128<16, CS...>() <= SIGNED_LIMIT);
static constexpr __int128 PAYLOAD =
	static_cast<__int128>(MakeU128<16, CS...>());
};

template <char... CS>
struct Static128<'0', 'X', CS...> {
static constexpr bool IS_VALID =
	ValidateU128<16, CS...>() && (MakeU128<16, CS...>() <= SIGNED_LIMIT);
static constexpr __int128 PAYLOAD =
	static_cast<__int128>(MakeU128<16, CS...>());
};

template <char... CS>
struct Static128<'0', 'b', CS...> {
static constexpr bool IS_VALID =
	ValidateU128<2, CS...>() && (MakeU128<2, CS...>() <= SIGNED_LIMIT);
static constexpr __int128 PAYLOAD =
	static_cast<__int128>(MakeU128<2, CS...>());
};

template <char... CS>
struct Static128<'0', 'B', CS...> {
static constexpr bool IS_VALID =
	ValidateU128<2, CS...>() && (MakeU128<2, CS...>() <= SIGNED_LIMIT);
static constexpr __int128 PAYLOAD =
	static_cast<__int128>(MakeU128<2, CS...>());
};

template <char... CS>
struct Static128<'0', CS...> {
static constexpr bool IS_VALID =
	ValidateU128<8, CS...>() && (MakeU128<8, CS...>() <= SIGNED_LIMIT);
static constexpr __int128 PAYLOAD =
	static_cast<__int128>(MakeU128<8, CS...>());
};

}  // namespace suffix128

template <char... CS>
constexpr unsigned __int128 operator"" _u128() {
static_assert(suffix128::StaticU128<CS...>::IS_VALID,
				"Invalid characters or number too large");
return suffix128::StaticU128<CS...>::PAYLOAD;
}

template <char... CS>
constexpr __int128 operator"" _128() {
static_assert(suffix128::Static128<CS...>::IS_VALID,
				"Invalid characters or number too large");
return suffix128::Static128<CS...>::PAYLOAD;
}

//代码用于寻找合适的快速逆开方魔数用于记录在elc的float_info.hpp
//或许我该把这个蠢东西放在elc的github action中不定时自动运行并更新文件内容
//剩下的一切交给虚空和他妈的庞加莱

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

typedef __int128 int128_t;
typedef unsigned __int128 uint128_t;
typedef long double float128_t;
static_assert(sizeof(float128_t)*8==128);
decltype(cout)& operator<< (decltype(cout)&stm,int128_t num){
	if(num<0){stm << "-0x";num=-num;}
	if(num>=16){stm << (int128_t)(num/16);num%=16;}
	stm << (char)(num+(num<10?'0':('A'-10)));
	return stm;
}

template<class T,class U>
T union_cast(U data){//酷酷的小语法糖，和它在elc中的版本一样酷。
	union{U a;T b;}t{data};
	return t.b;
}
float128_t quick_invsqrt_with_magic_num(int128_t mnum,float128_t X){
	auto x= union_cast<float128_t>(mnum - (union_cast<int128_t>(X) >> 1));
	//x=x*(1.5-x*x*X/2);//有时候这会有用。
	return x;
}

//从汉堡反推出其原料：汉堡矿！
int128_t rewind_get_magic_number(float128_t X){
	float128_t x=1/sqrt((double)X);
	//反向合成
	int128_t mnum=union_cast<int128_t>(X)>>1;
	mnum+=union_cast<int128_t>(x);
	return mnum;
}

//我不打算启用最下面两行的测试用例，现有的数足够评价你获得的魔数好坏了
//对于极小的测试数据 它们会使得魔数的品质因素p比例失衡 脱离我原有的构想。
float128_t test_nums[]={
	2,7,17,13,72,73,2003,2011,2017,2027,1989,
	72893,72901,72907,72911,72923,72931,72937,72949,72953,
	9.99e+37,72e72,17e17,114514e191,980,
	1.0/3,1.0/7,1.0/17,1.0/13,1.0/72,1.0/73,1.0/2003,1.0/2011,1.0/2017,1.0/2027,1.0/1989,
	//1.0/72893,1.0/72901,1.0/72907,1.0/72911,1.0/72923,1.0/72931,1.0/72937,1.0/72949,1.0/72953,
	//1.0/9.99e+37,1.0/72e72,1.0/17e17,1.0/114514e191,1.0/980,
};

float128_t get_q_of_magic_num(int128_t mnum){
	float128_t result=0;
	for(auto&tnum:test_nums){
		auto tres=1/sqrt(tnum);
		auto res=quick_invsqrt_with_magic_num(mnum,tnum);
		if(isnan(res))
			result+=13;
		else{
			auto squrt=abs(tres-res)/tnum;
			result+=squrt;
		}
	}
	return result;
}

#include <random>
mt19937_64 rand_engene(72);//mt19937不是完美的随机数引擎（众所周知 它没有通过TestU01）
//但是elc在gcc下无法编译，所以让我们和elc::rand挥手说再见

int128_t rand_128(){//我知道的 不可以泛型 真蠢
	uint128_t first=rand_engene();
	return int128_t((first<< 64)+rand_engene());
}

int128_t average(int128_t a,int128_t b){//我们需要防止溢出
	return (a>>1)+(b>>1)+(a&b&1);
}

void log_magic_num(int128_t mnum){
	size_t bad_num=0;
	const auto q_of_i=get_q_of_magic_num(mnum);
	cout << mnum << " : " << q_of_i << endl;
	for(auto&tnum:test_nums){
		auto res=quick_invsqrt_with_magic_num(mnum,tnum);
		if(isnan(res)){
			if(!bad_num)
				cout << "not good enough on";
			cout << ' ' << tnum;
			bad_num++;
		}
	}
	cout << endl;
	constexpr auto test_size=sizeof(test_nums)/sizeof(int128_t);
	cout << "bad rate " << bad_num << "/" << test_size
		 << " (" << float128_t(bad_num)/test_size*100 << "%)" << endl;
	
}

constexpr auto r0=0.432744889959443195468521587014f128;

#include<chrono>

//一个魔鬼的想法是当你将这一切交给随机数和虚空后
//我不认为float128有完美的快速逆开方魔数存在使其对于你能想到的任何测试值都不产生NaN
//如果你找到了 请联系我。

int main() {
	cout << (190_128 >> 1);
	cout << int128_t(r0*(1_u128 << 23)+0.5);
	cout << endl;
	cout << (24574_128 >> 1);
	cout << int128_t(r0*(1_u128 << 112)+0.5);
	cout << endl;//理论最佳值？（来自FAST INVERSE SQUARE ROOT form CHRIS LOMONT）
	auto best_mnum=0x2FFF6EC85E7DE30DAABC6027118577EF_128;
	auto best_q=get_q_of_magic_num(best_mnum);
log:
	log_magic_num(best_mnum);
sleep:
	//_sleep(1000);//休息一下！
reset_seed:
	rand_engene.seed(chrono::system_clock::now().time_since_epoch().count());
loop:
	for(unsigned i=0;i<1000;i++){
		auto mnum=rand_128();
		auto q_of_i=get_q_of_magic_num(mnum);
		if(q_of_i<best_q){
			best_q=q_of_i;
			best_mnum=mnum;
			//向数值两侧扩展，寻找更好的魔数
			const auto step=0xbeadbeef;
			auto tmp_mnum=best_mnum;
			for(int i=0;i<100;i++){
				tmp_mnum+=step;
				auto q_of_i=get_q_of_magic_num(tmp_mnum);
				if(q_of_i<best_q){
					best_q=q_of_i;
					best_mnum=tmp_mnum;
				}
			}
			tmp_mnum=best_mnum;
			for(int i=0;i<100;i++){
				tmp_mnum-=step;
				auto q_of_i=get_q_of_magic_num(tmp_mnum);
				if(q_of_i<best_q){
					best_q=q_of_i;
					best_mnum=tmp_mnum;
				}
			}
			goto log;
		}
	}
	//goto sleep;
	//*/
}