#pragma once #include "gp/algorithm/repeat.hpp" #include #include #include namespace gp{ template constexpr T pi; template<> constexpr float pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062; template<> constexpr double pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062L; template T abs(T); template<> float abs(float value) { static_assert(sizeof(float) == 4, "bad float size"); union { float fp; uint32_t ab; } p; p.fp = value; p.ab &= 0x7fFFffFF; return p.fp; } template<> double abs(double value) { static_assert(sizeof(double) == 8, "bad double size"); union { double fp; uint64_t ab; } p; p.fp = value; p.ab &= 0x7fFFffFFffFFffFF; return p.fp; } template T floor(T); template<> float floor(float value) { static_assert(sizeof(float) == 4, "bad float size"); if( value >= 16777216 || value <= std::numeric_limits::min() || value != value ) { return value; } int32_t ret = value; float ret_d = ret; if(value == ret_d || value >= 0) { return ret; } else { return ret-1; } } template<> double floor(double value) { static_assert(sizeof(double) == 8, "bad double size"); if( value >= 9007199254740992 || value <= std::numeric_limits::min() || value != value ) { return value; } int64_t ret = value; double ret_d = ret; if(value == ret_d || value >= 0) { return ret; } else { return ret-1; } } template T sign(T); template<> float sign(float value) { static_assert(sizeof(float) == 4, "bad float size"); if(!value) return 0; union { float fp; uint32_t ab; } p; p.fp = value; p.ab &= 0x7fFFffFF; return value/p.fp; } template<> double sign(double value) { static_assert(sizeof(double) == 8, "bad double size"); if(!value) return 0; union { double fp; uint64_t ab; } p; p.fp = value; p.ab &= 0x7fFFffFFffFFffFF; return value/p.fp; } template T sin_taylor(T value) { const T acc = T{1}/T{accuracy}; T B = value; T C = 1; T ret = B/C; for(size_t i = 1; (i < steps) && (abs<>(B/C) > acc); ++i) { B *= -1*value*value; C *= 2*i*(2*i+1); ret += B/C; } return ret; } float sin(float v) { v += pi; v = v - 2*pi*floor(v/(2*pi)); v -= pi; float s = sign(v); v *= s; return sin_taylor<10>(v)*s; } double sin(double v) { v += pi; v = v - 2*pi*floor(v/(2*pi)); v -= pi; float s = sign(v); v *= s; return sin_taylor<10>(v)*s; } // TODO: replace with an actual implementation float cos(float v) { return sin(v+pi/2); } // TODO: replace with an actual implementation double cos(double v) { return sin(v+pi/2); } // TODO: replace with an actual implementation float tan(float v) { return sin(v)/cos(v); } // TODO: replace with an actual implementation double tan(double v) { return sin(v)/cos(v); } template float isqrt(float v) { int32_t i; float x2, y; constexpr float threehalfs = 1.5F; x2 = v * 0.5F; y = v; i = * ( int32_t * ) &y; i = 0x5F375A86 - ( i >> 1 ); y = * ( float * ) &i; gp::repeat(cycles, [&](){ y = y * ( threehalfs - ( x2 * y * y ) ); }); return y; } template double isqrt(double v) { int64_t i; double x2, y; constexpr double threehalfs = 1.5F; x2 = v * 0.5F; y = v; i = * ( int64_t * ) &y; i = 0x5FE6EB50C7B537A9 - ( i >> 1 ); y = * ( double * ) &i; gp::repeat(cycles, [&](){ y = y * ( threehalfs - ( x2 * y * y ) ); }); return y; } float fast_isqrt(float v) {return isqrt<1>(v);} double fast_isqrt(double v) {return isqrt<1>(v);} }