|
|
- #pragma once
- #include <stddef.h>
- #include <stdint.h>
- #include <limits>
- #include "gp/algorithm/repeat.hpp"
-
- namespace gp{
-
- template<typename T>
- constexpr T pi;
-
- template<>
- constexpr float pi<float> = 3.1415926535897932384626433832795028841971693993751058209749445923078164062;
-
- template<>
- constexpr double pi<double> = 3.1415926535897932384626433832795028841971693993751058209749445923078164062;
-
- template<typename T>
- T abs(T);
-
- template<>
- float abs<float>(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>(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<typename T>
- T floor(T);
-
- template<>
- float floor<float>(float value) {
- static_assert(sizeof(float) == 4, "bad float size");
- if(
- value >= std::numeric_limits<int32_t>::max()
- || value <= std::numeric_limits<int32_t>::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>(double value) {
- static_assert(sizeof(double) == 8, "bad double size");
- if(
- value >= std::numeric_limits<int64_t>::max()
- || value <= std::numeric_limits<int64_t>::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<typename T>
- T sign(T);
-
- template<>
- float sign<float>(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>(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<size_t steps, typename T, size_t accuracy = 1000000>
- 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<float>;
- v = v - 2*pi<float>*floor(v/(2*pi<float>));
- v -= pi<float>;
- float s = sign(v);
- v *= s;
- return sin_taylor<10>(v)*s;
- }
-
- double sin(double v) {
- v += pi<double>;
- v = v - 2*pi<double>*floor(v/(2*pi<double>));
- v -= pi<double>;
- 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<float>/2);
- }
-
- // TODO: replace with an actual implementation
- double cos(double v) {
- return sin(v+pi<float>/2);
- }
-
- template<size_t cycles = 5>
- 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<size_t cycles = 5>
- 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);}
- }
|