From 89eb1775dbe8cd43f6971020a0e33194380f54a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Branimir=20Karad=C5=BEi=C4=87?= Date: Thu, 8 Oct 2026 19:16:41 -0700 Subject: [PATCH] Added a constexpr double math surface, and made fromString exact. (#451) --- include/bx/constants.h | 90 ++ include/bx/inline/math.inl | 2503 ++++++++++++++++++++++++----- include/bx/inline/simd128_ref.inl | 119 +- include/bx/inline/simd256_ref.inl | 58 +- include/bx/inline/simd_ni.inl | 264 ++- include/bx/math.h | 618 ++++++- src/crtnone.cpp | 54 +- src/dtoa.cpp | 907 +++++------ src/hash.cpp | 6 +- src/math.cpp | 18 - src/string.cpp | 10 +- tests/hash_test.cpp | 6 +- tests/math_test.cpp | 620 ++++++- tests/simd_test.cpp | 126 +- tests/string_test.cpp | 150 ++ 15 files changed, 4382 insertions(+), 1167 deletions(-) diff --git a/include/bx/constants.h b/include/bx/constants.h index 5515f5c..b08330d 100644 --- a/include/bx/constants.h +++ b/include/bx/constants.h @@ -83,6 +83,51 @@ namespace bx /// Floating-point infinity. // constexpr float kFloatInfinity; + /// The reciprocal of ln(2). 1/ln(2) + constexpr float kFloatInvLn2 = kInvLogNat2; + + /// The reciprocal of ln(10). 1/ln(10) + constexpr float kFloatInvLn10 = 1.0f/kLogNat10; + + /// Leading part of ln(2), with the trailing mantissa bits cleared. + constexpr float kFloatLn2Hi = 6.9314575195e-01f; + + /// What is left of ln(2) after kFloatLn2Hi. + constexpr float kFloatLn2Lo = 1.4286067653e-06f; + + /// Largest argument exp can answer with a finite number. + constexpr float kFloatExpOverflow = 88.7228394f; + + /// Smallest argument exp can answer with a non-zero number. + constexpr float kFloatExpUnderflow = -103.972084f; + + /// Exponent of the square root of the scale expHalf splits its result by. + constexpr int32_t kFloatExpHalfShift = 235; + + /// Argument reduction expHalf pairs with kFloatExpHalfShift. ln(2)*kFloatExpHalfShift/2 + constexpr float kFloatExpHalfLn2 = 1.62889587431602213e+02f; + + /// Below this the hyperbolic functions take their series instead of exp. + constexpr float kFloatHypSmall = 0.0625f; + + /// Above this exp(-_a) no longer contributes to sinh and cosh. + constexpr float kFloatHypBig = 9.0f; + + /// Above this half of exp(_a) overflows, and expHalf takes over. + constexpr float kFloatHypNoFit = 88.0f; + + /// Above this tanh is 1. + constexpr float kFloatTanhLarge = 10.0f; + + /// Below this the arc hyperbolic functions answer with the argument. + constexpr float kFloatArcHypTiny = 0.000244140625f; + + /// Above this asinh and acosh are log(_a) + ln(2). + constexpr float kFloatArcHypHuge = 4096.0f; + + /// Seed for the reciprocal square root iteration. + constexpr uint32_t kFloatRsqrtSeed = UINT32_C(0x5f3759df); + /// constexpr uint8_t kDoubleSignNumBits = 1; constexpr uint8_t kDoubleSignBitShift = 63; @@ -104,6 +149,51 @@ namespace bx // Double-precision floating-point infinity. // constexpr double kDoubleInfinity; + /// The reciprocal of ln(2). 1/ln(2) + constexpr double kDoubleInvLn2 = 1.44269504088896338700; + + /// The reciprocal of ln(10). 1/ln(10) + constexpr double kDoubleInvLn10 = 4.34294481903251827651e-01; + + /// Leading part of ln(2), with the trailing mantissa bits cleared. + constexpr double kDoubleLn2Hi = 6.93147180369123816490e-01; + + /// What is left of ln(2) after kDoubleLn2Hi. + constexpr double kDoubleLn2Lo = 1.90821492927058770002e-10; + + /// Largest argument exp can answer with a finite number. + constexpr double kDoubleExpOverflow = 7.09782712893383973096e+02; + + /// Smallest argument exp can answer with a non-zero number. + constexpr double kDoubleExpUnderflow = -7.45133219101941108420e+02; + + /// Exponent of the square root of the scale expHalf splits its result by. + constexpr int32_t kDoubleExpHalfShift = 2043; + + /// Argument reduction expHalf pairs with kDoubleExpHalfShift. ln(2)*kDoubleExpHalfShift/2 + constexpr double kDoubleExpHalfLn2 = 1.41609968988396828e+03; + + /// Below this the hyperbolic functions take their series instead of exp. + constexpr double kDoubleHypSmall = 0.03125; + + /// Above this exp(-_a) no longer contributes to sinh and cosh. + constexpr double kDoubleHypBig = 22.0; + + /// Above this half of exp(_a) overflows, and expHalf takes over. + constexpr double kDoubleHypNoFit = 709.0; + + /// Above this tanh is 1. + constexpr double kDoubleTanhLarge = 20.0; + + /// Below this the arc hyperbolic functions answer with the argument. + constexpr double kDoubleArcHypTiny = 3.72529029846191406250e-09; + + /// Above this asinh and acosh are log(_a) + ln(2). + constexpr double kDoubleArcHypHuge = 2.68435456e+08; + + /// Seed for the reciprocal square root iteration. + constexpr uint64_t kDoubleRsqrtSeed = UINT64_C(0x5fe6eb50c7b537a9); + } // namespace bx #endif // BX_CONSTANTS_H_HEADER_GUARD diff --git a/include/bx/inline/math.inl b/include/bx/inline/math.inl index 4aa7ed4..7748138 100644 --- a/include/bx/inline/math.inl +++ b/include/bx/inline/math.inl @@ -109,48 +109,127 @@ namespace bx return tmp == kDoubleExponentMask; } - inline BX_CONSTEXPR_FUNC float truncRef(float _a) + /// True for an argument whose magnitude has no fractional part left. + template + inline BX_CONSTEXPR_FUNC bool isIntegralMagnitude(Ty _a) { - const uint32_t bits = floatToBits(_a); - const uint32_t sign = bits & kFloatSignMask; - const uint32_t magnitude = bits & ~kFloatSignMask; + using Traits = FloatT; + using Bits = typename Traits::Bits; - if (0x4b000000u <= magnitude) + constexpr Bits kIntegralExp = Bits(Traits::kExponentBias + Traits::kMantissaNumBits); + constexpr Bits kIntegral = kIntegralExp << Traits::kExponentBitShift; + + const Bits bits = Traits::toBits(_a); + const Bits magnitude = bits & ~Traits::kSignMask; + + return kIntegral <= magnitude; + } + + template + inline BX_CONSTEXPR_FUNC Ty truncT(Ty _a) + { + using Traits = FloatT; + using Bits = typename Traits::Bits; + using Int = typename Traits::Int; + + if (isIntegralMagnitude(_a) ) { return _a; } - const float tr = float(int(_a) ); - const uint32_t trBits = floatToBits(tr) | sign; - const float result = bitsToFloat(trBits); + const Bits sign = Traits::toBits(_a) & Traits::kSignMask; + const Int integral = Int(_a); + const Ty tr = Ty(integral); + const Bits trBits = Traits::toBits(tr) | sign; + const Ty result = Traits::fromBits(trBits); + return result; } + template + inline BX_CONSTEXPR_FUNC Ty floorT(Ty _a) + { + using Traits = FloatT; + using Bits = typename Traits::Bits; + using Int = typename Traits::Int; + + if (isIntegralMagnitude(_a) ) + { + return _a; + } + + const Bits sign = Traits::toBits(_a) & Traits::kSignMask; + const Int integral = Int(_a); + const Ty tr = Ty(integral); + const Ty fl = tr > _a ? tr - Ty(1.0) : tr; + const Bits flBits = Traits::toBits(fl) | sign; + const Ty result = Traits::fromBits(flBits); + + return result; + } + + template + inline BX_CONSTEXPR_FUNC Ty ceilT(Ty _a) + { + const Ty na = -_a; + const Ty fl = floorT(na); + const Ty result = -fl; + + return result; + } + + template + inline BX_CONSTEXPR_FUNC Ty roundT(Ty _a) + { + using Traits = FloatT; + using Bits = typename Traits::Bits; + + if (isIntegralMagnitude(_a) ) + { + return _a; + } + + constexpr Ty kMagic = Ty(Bits(1) << Traits::kMantissaNumBits); + + const Bits bits = Traits::toBits(_a); + const Bits sign = bits & Traits::kSignMask; + const Ty absA = Traits::fromBits(bits & ~Traits::kSignMask); + const Ty raised = absA + kMagic; + const Ty rd = raised - kMagic; + const Bits rdBits = Traits::toBits(rd) | sign; + const Ty result = Traits::fromBits(rdBits); + + return result; + } + + inline BX_CONSTEXPR_FUNC float truncRef(float _a) + { + return truncT(_a); + } + + inline BX_CONSTEXPR_FUNC double truncRef(double _a) + { + return truncT(_a); + } + inline BX_CONSTEXPR_FUNC float floorRef(float _a) { - const uint32_t bits = floatToBits(_a); - const uint32_t sign = bits & kFloatSignMask; - const uint32_t magnitude = bits & ~kFloatSignMask; + return floorT(_a); + } - if (0x4b000000u <= magnitude) - { - return _a; - } - - const float tr = float(int(_a) ); - const float fl = tr > _a ? tr - 1.0f : tr; - const uint32_t flBits = floatToBits(fl) | sign; - const float result = bitsToFloat(flBits); - - return result; + inline BX_CONSTEXPR_FUNC double floorRef(double _a) + { + return floorT(_a); } inline BX_CONSTEXPR_FUNC float ceilRef(float _a) { - const float na = -_a; - const float fl = floorRef(na); - const float result = -fl; - return result; + return ceilT(_a); + } + + inline BX_CONSTEXPR_FUNC double ceilRef(double _a) + { + return ceilT(_a); } #if BX_SIMD_SUPPORTED @@ -176,6 +255,17 @@ namespace bx return out; } + inline BX_CONST_FUNC double floorSimd(double _a) + { + const simd128_t aa = simd_splat(_a); + const simd128_t result = simd_f64_floor(aa); + + alignas(16) double out[2] = { 0.0, 0.0 }; + simd_st(out, result); + + return out[0]; + } + inline BX_CONST_FUNC float ceilSimd(float _a) { const simd128_t aa = simd_splat(_a); @@ -186,55 +276,27 @@ namespace bx return out; } -#endif // BX_SIMD_SUPPORTED - inline BX_CONSTEXPR_FUNC float floor(float _a) + inline BX_CONST_FUNC double ceilSimd(double _a) { -#if BX_SIMD_SUPPORTED - if (isConstantEvaluated() ) - { - return floorRef(_a); - } + const simd128_t aa = simd_splat(_a); + const simd128_t result = simd_f64_ceil(aa); - return floorSimd(_a); -#else - return floorRef(_a); -#endif // BX_SIMD_SUPPORTED + alignas(16) double out[2] = { 0.0, 0.0 }; + simd_st(out, result); + + return out[0]; } - - inline BX_CONSTEXPR_FUNC float ceil(float _a) - { -#if BX_SIMD_SUPPORTED - if (isConstantEvaluated() ) - { - return ceilRef(_a); - } - - return ceilSimd(_a); -#else - return ceilRef(_a); #endif // BX_SIMD_SUPPORTED - } inline BX_CONSTEXPR_FUNC float roundRef(float _a) { - const uint32_t bits = floatToBits(_a); - const uint32_t sign = bits & kFloatSignMask; - const uint32_t magnitude = bits & ~kFloatSignMask; + return roundT(_a); + } - if (0x4b000000u <= magnitude) - { - return _a; - } - - const float fl = floorRef(_a); - const float fr = _a - fl; - const bool odd = 0 != (int32_t(fl) & 1); - const bool up = fr > 0.5f || (fr == 0.5f && odd); - const float rd = up ? fl + 1.0f : fl; - const uint32_t rdBits = floatToBits(rd) | sign; - const float result = bitsToFloat(rdBits); - return result; + inline BX_CONSTEXPR_FUNC double roundRef(double _a) + { + return roundT(_a); } #if BX_SIMD_SUPPORTED @@ -248,21 +310,18 @@ namespace bx return out; } -#endif // BX_SIMD_SUPPORTED - inline BX_CONSTEXPR_FUNC float round(float _a) + inline BX_CONST_FUNC double roundSimd(double _a) { -#if BX_SIMD_SUPPORTED - if (isConstantEvaluated() ) - { - return roundRef(_a); - } + const simd128_t aa = simd_splat(_a); + const simd128_t result = simd_f64_round(aa); - return roundSimd(_a); -#else - return roundRef(_a); -#endif // BX_SIMD_SUPPORTED + alignas(16) double out[2] = { 0.0, 0.0 }; + simd_st(out, result); + + return out[0]; } +#endif // BX_SIMD_SUPPORTED inline BX_CONSTEXPR_FUNC float lerp(float _a, float _b, float _t) { @@ -278,63 +337,53 @@ namespace bx return (_value - _a) / (_b - _a); } - inline BX_CONSTEXPR_FUNC float sign(float _a) - { - return float( (0.0f < _a) - (0.0f > _a) ); - } - inline BX_CONSTEXPR_FUNC bool signBit(float _a) { const uint32_t bits = floatToBits(_a); return 0 != (bits & kFloatSignMask); } - inline BX_CONSTEXPR_FUNC float copySign(float _value, float _sign) + template + inline BX_CONSTEXPR_FUNC Ty fractT(Ty _a) { -#if BX_COMPILER_MSVC - const uint32_t magnitude = floatToBits(_value) & ~kFloatSignMask; - const uint32_t sign = floatToBits(_sign) & kFloatSignMask; - const uint32_t bits = magnitude | sign; - const float result = bitsToFloat(bits); + const Ty tr = trunc(_a); + const Ty result = _a - tr; + return result; -#else - return __builtin_copysign(_value, _sign); -#endif // BX_COMPILER_MSVC } - inline BX_CONSTEXPR_FUNC float abs(float _a) + template + inline BX_CONSTEXPR_FUNC Ty fmodT(Ty _a, Ty _b) { - return _a < 0.0f ? -_a : _a; + const Ty quotient = _a / _b; + const Ty whole = trunc(quotient); + const Ty result = nms(_b, whole, _a); + + return result; } - inline BX_CONSTEXPR_FUNC float square(float _a) + template + inline BX_CONSTEXPR_FUNC Ty modT(Ty _a, Ty _b) { - return _a * _a; + const Ty quotient = _a / _b; + const Ty whole = floor(quotient); + const Ty result = nms(_b, whole, _a); + + return result; } - inline BX_CONSTEXPR_FUNC float trunc(float _a) + template + inline BX_CONSTEXPR_FUNC Ty exp2T(Ty _a) { -#if BX_SIMD_SUPPORTED - if (isConstantEvaluated() ) - { - return truncRef(_a); - } - - return truncSimd(_a); -#else - return truncRef(_a); -#endif // BX_SIMD_SUPPORTED + return pow(Ty(2.0), _a); } - inline BX_CONSTEXPR_FUNC float fract(float _a) + template + inline BX_CONSTEXPR_FUNC Ty tanT(Ty _a) { - return _a - trunc(_a); - } - - inline BX_CONSTEXPR_FUNC float nms(float _a, float _b, float _c) - { - const float na = -_a; - const float result = mad(na, _b, _c); + const Ty sn = sin(_a); + const Ty cs = cos(_a); + const Ty result = sn / cs; return result; } @@ -399,16 +448,6 @@ namespace bx return _a * _b; } - inline BX_CONSTEXPR_FUNC float mad(float _a, float _b, float _c) - { - const simd32_t aa = simd32_ld(_a); - const simd32_t bb = simd32_ld(_b); - const simd32_t cc = simd32_ld(_c); - const simd32_t result = simd32_f32_madd(aa, bb, cc); - - return bitCast(result); - } - inline BX_CONSTEXPR_FUNC float rcp(float _a) { return 1.0f / _a; @@ -431,87 +470,6 @@ namespace bx BX_FP_PRECISE_BEGIN() - inline BX_CONSTEXPR_FUNC float mod(float _a, float _b) - { - const float quotient = _a / _b; - const float whole = floor(quotient); - const float result = nms(_b, whole, _a); - return result; - } - - inline BX_CONSTEXPR_FUNC float cos(float _a) - { - const float scaled = _a * (2.0f*kInvPi); - const float real = floor(scaled); - const float xx = _a - real * kPiHalf; - const int32_t bits = isFinite(real) - ? int32_t(mod(real, 4.0f) ) & 3 - : 0 - ; - - constexpr float kSinC2 = -0.16666667163372039794921875f; - constexpr float kSinC4 = 8.333347737789154052734375e-3f; - constexpr float kSinC6 = -1.9842604524455964565277099609375e-4f; - constexpr float kSinC8 = 2.760012648650445044040679931640625e-6f; - constexpr float kSinC10 = -2.50293279435709337121807038784027099609375e-8f; - - constexpr float kCosC2 = -0.5f; - constexpr float kCosC4 = 4.166664183139801025390625e-2f; - constexpr float kCosC6 = -1.388833043165504932403564453125e-3f; - constexpr float kCosC8 = 2.47562347794882953166961669921875e-5f; - constexpr float kCosC10 = -2.59630184018533327616751194000244140625e-7f; - - float c0 = xx; - float c2 = kSinC2; - float c4 = kSinC4; - float c6 = kSinC6; - float c8 = kSinC8; - float c10 = kSinC10; - - if (bits == 0 - || bits == 2) - { - c0 = 1.0f; - c2 = kCosC2; - c4 = kCosC4; - c6 = kCosC6; - c8 = kCosC8; - c10 = kCosC10; - } - - const float xsq = square(xx); - const float tmp0 = mad(c10, xsq, c8 ); - const float tmp1 = mad(tmp0, xsq, c6 ); - const float tmp2 = mad(tmp1, xsq, c4 ); - const float tmp3 = mad(tmp2, xsq, c2 ); - const float tmp4 = mad(tmp3, xsq, 1.0); - const float result = tmp4 * c0; - - return bits == 1 || bits == 2 - ? -result - : result - ; - } - - inline BX_CONSTEXPR_FUNC float acos(float _a) - { - constexpr float kAcosC0 = 1.5707288f; - constexpr float kAcosC1 = -0.2121144f; - constexpr float kAcosC2 = 0.0742610f; - constexpr float kAcosC3 = -0.0187293f; - - const float absa = abs(_a); - const float tmp0 = mad(kAcosC3, absa, kAcosC2); - const float tmp1 = mad(tmp0, absa, kAcosC1); - const float tmp2 = mad(tmp1, absa, kAcosC0); - const float tmp3 = tmp2 * sqrt(1.0f - absa); - const float negate = float(_a < 0.0f); - const float tmp4 = tmp3 - 2.0f*negate*tmp3; - const float result = negate*kPi + tmp4; - - return result; - } - inline void sinCosApprox(float& _outSinApprox, float& _outCos, float _a) { const float aa = _a - floor(_a*kInvPi2)*kPi2; @@ -527,243 +485,1956 @@ BX_FP_PRECISE_BEGIN() _outCos = cosA; } - inline BX_CONSTEXPR_FUNC float sin(float _a) + template + inline BX_CONSTEXPR_FUNC Ty expHalfT(Ty _a) { - return cos(_a - kPiHalf); - } + using Traits = FloatT; + using Bits = typename Traits::Bits; - inline BX_CONSTEXPR_FUNC float sinh(float _a) - { - return 0.5f*(exp(_a) - exp(-_a) ); - } + constexpr Bits kBits = Bits(Traits::kExponentBias + Traits::kExpHalfShift/2) << Traits::kExponentBitShift; + constexpr Ty kScale = Traits::fromBits(kBits); - inline BX_CONSTEXPR_FUNC float asin(float _a) - { - return kPiHalf - acos(_a); - } - - inline BX_CONSTEXPR_FUNC float cosh(float _a) - { - return 0.5f*(exp(_a) + exp(-_a) ); - } - - inline BX_CONSTEXPR_FUNC float tan(float _a) - { - return sin(_a) / cos(_a); - } - - inline BX_CONSTEXPR_FUNC float tanh(float _a) - { - const float tmp0 = exp(2.0f*_a); - const float tmp1 = tmp0 - 1.0f; - const float tmp2 = tmp0 + 1.0f; - const float result = tmp1 / tmp2; + const Ty reduced = _a - Traits::kExpHalfLn2; + const Ty ee = exp(reduced); + const Ty once = ee * kScale; + const Ty result = once * kScale; return result; } - inline BX_CONSTEXPR_FUNC float atan(float _a) + inline BX_CONSTEXPR_FUNC float expHalf(float _a) { - return atan2(_a, 1.0f); + return expHalfT(_a); } - inline BX_CONSTEXPR_FUNC float atan2(float _y, float _x) + inline BX_CONSTEXPR_FUNC double expHalf(double _a) { - const float ax = abs(_x); - const float ay = abs(_y); - const float maxaxy = max(ax, ay); - const float minaxy = min(ax, ay); - - const uint32_t ysign = floatToBits(_y) & kFloatSignMask; - - if (maxaxy == 0.0f) - { - return bitsToFloat(ysign); - } - - constexpr float kAtan2C0 = -0.013480470f; - constexpr float kAtan2C1 = 0.057477314f; - constexpr float kAtan2C2 = -0.121239071f; - constexpr float kAtan2C3 = 0.195635925f; - constexpr float kAtan2C4 = -0.332994597f; - constexpr float kAtan2C5 = 0.999995630f; - - const float mxy = minaxy / maxaxy; - const float mxysq = square(mxy); - const float tmp0 = mad(kAtan2C0, mxysq, kAtan2C1); - const float tmp1 = mad(tmp0, mxysq, kAtan2C2); - const float tmp2 = mad(tmp1, mxysq, kAtan2C3); - const float tmp3 = mad(tmp2, mxysq, kAtan2C4); - const float tmp4 = mad(tmp3, mxysq, kAtan2C5); - const float tmp5 = tmp4 * mxy; - const float tmp6 = ay > ax ? kPiHalf - tmp5 : tmp5; - const float tmp7 = _x < 0.0f ? kPi - tmp6 : tmp6; - - const uint32_t bits = floatToBits(tmp7) | ysign; - const float result = bitsToFloat(bits); - - return result; + return expHalfT(_a); } - inline BX_CONSTEXPR_FUNC float ldexp(float _a, int32_t _b) + template + inline BX_CONSTEXPR_FUNC Ty sinhT(Ty _a) { - const simd32_t ftob = simd32_splat(floatToBits(_a)); - const simd32_t signexpmask = simd32_splat(kFloatSignMask | kFloatExponentMask); - const simd32_t mantmask = simd32_splat(kFloatMantissaMask); - const simd32_t b = simd32_splat(_b); - const simd32_t masked = simd32_and(ftob, signexpmask); - const simd32_t expsign0 = simd32_x32_sra(masked, kFloatExponentBitShift); - const simd32_t tmp = simd32_i32_add(expsign0, b); - const simd32_t expsign1 = simd32_x32_sll(tmp, kFloatExponentBitShift); - const simd32_t mantissa = simd32_and(ftob, mantmask); - const simd32_t bits = simd32_or(mantissa, expsign1); + using Traits = FloatT; - return bitsToFloat(bits.u32); + constexpr Ty kSinhC3 = Ty(1.0)/Ty(6.0); + constexpr Ty kSinhC5 = Ty(1.0)/Ty(120.0); + constexpr Ty kSinhC7 = Ty(1.0)/Ty(5040.0); + + const Ty absA = abs(_a); + + Ty magnitude = Ty(0.0); + + if (absA < Traits::kHypSmall) + { + const Ty zz = square(absA); + const Ty inn0 = mad(zz, kSinhC7, kSinhC5); + const Ty inn1 = mad(zz, inn0, kSinhC3); + const Ty poly = mad(zz, inn1, Ty(1.0) ); + + magnitude = absA * poly; + } + else if (absA < Traits::kHypBig) + { + const Ty ep = exp(absA); + const Ty en = exp(-absA); + const Ty diff = ep - en; + + magnitude = Ty(0.5) * diff; + } + else if (absA < Traits::kHypNoFit) + { + const Ty ep = exp(absA); + + magnitude = Ty(0.5) * ep; + } + else + { + magnitude = expHalfT(absA); + } + + return signBit(_a) ? -magnitude : magnitude; } - inline BX_CONSTEXPR_FUNC float log(float _a) + template + inline BX_CONSTEXPR_FUNC Ty coshT(Ty _a) { - if (_a < 0.0f) + using Traits = FloatT; + + const Ty absA = abs(_a); + + if (absA < Traits::kHypBig) { - return bitsToFloat(kFloatSignMask | kFloatExponentMask | kFloatMantissaMask); + const Ty ep = exp(absA); + const Ty en = exp(-absA); + const Ty sum = ep + en; + + return Ty(0.5) * sum; } - if (_a == 0.0f) + if (absA < Traits::kHypNoFit) { - return -kFloatInfinity; + const Ty ep = exp(absA); + + return Ty(0.5) * ep; } - const simd32_t ftob = simd32_splat(floatToBits(_a)); - const simd32_t expmask = simd32_splat(kFloatExponentMask); - const simd32_t signmantmask = simd32_splat(kFloatSignMask | kFloatMantissaMask); - const simd32_t half = simd32_splat(UINT32_C(0x3f000000)); - - const simd32_t masked0 = simd32_and(ftob, expmask); - const simd32_t exp0 = simd32_x32_srl(masked0, kFloatExponentBitShift); - - int32_t exp = int32_t(exp0.u32) - 0x7e; - - const simd32_t masked1 = simd32_and(ftob, signmantmask); - const simd32_t bits = simd32_or(masked1, half); - float ff = bitsToFloat(bits.u32); - - if (ff < kSqrt2*0.5f) - { - ff *= 2.0f; - --exp; - } - - constexpr float kLogC0 = 6.666666666666735130e-01f; - constexpr float kLogC1 = 3.999999999940941908e-01f; - constexpr float kLogC2 = 2.857142874366239149e-01f; - constexpr float kLogC3 = 2.222219843214978396e-01f; - constexpr float kLogC4 = 1.818357216161805012e-01f; - constexpr float kLogC5 = 1.531383769920937332e-01f; - constexpr float kLogC6 = 1.479819860511658591e-01f; - constexpr float kLogNat2Lo = 1.90821492927058770002e-10f; - - ff -= 1.0f; - const float kk = float(exp); - const float hi = kk*kLogNat2; - const float lo = kk*kLogNat2Lo; - const float ss = ff / (2.0f + ff); - const float s2 = square(ss); - const float s4 = square(s2); - - const float tmp0 = mad(kLogC6, s4, kLogC4); - const float tmp1 = mad(tmp0, s4, kLogC2); - const float tmp2 = mad(tmp1, s4, kLogC0); - const float t1 = s2*tmp2; - - const float tmp3 = mad(kLogC5, s4, kLogC3); - const float tmp4 = mad(tmp3, s4, kLogC1); - const float t2 = s4*tmp4; - - const float t12 = t1 + t2; - const float hfsq = 0.5f*square(ff); - const float result = hi - ( (hfsq - (ss*(hfsq+t12) + lo) ) - ff); - - return result; + return expHalfT(absA); } - inline BX_CONSTEXPR_FUNC float exp(float _a) + template + inline BX_CONSTEXPR_FUNC Ty tanhT(Ty _a) { - if (abs(_a) <= kNearZero) - { - return _a + 1.0f; - } - - constexpr float expMin = log(kFloatSmallest); - - if (_a <= expMin) - { - return 0.0f; - } - - if (isInfinite(_a) && _a < 0.0f) - { - return 0.0f; - } - - constexpr float expMax = log(kFloatLargest); - - if (_a >= expMax) - { - return kFloatInfinity; - } + using Traits = FloatT; if (isNan(_a) ) { return _a; } - constexpr float kExpC0 = 1.66666666666666019037e-01f; - constexpr float kExpC1 = -2.77777777770155933842e-03f; - constexpr float kExpC2 = 6.61375632143793436117e-05f; - constexpr float kExpC3 = -1.65339022054652515390e-06f; - constexpr float kExpC4 = 4.13813679705723846039e-08f; - constexpr float kLogNat2Lo = 1.90821492927058770002e-10f; + const Ty absA = abs(_a); - const float kk = round(_a*kInvLogNat2); - const float hi = _a - kk*kLogNat2; - const float lo = kk*kLogNat2Lo; - const float hml = hi - lo; - const float hmlsq = square(hml); - const float tmp0 = mad(kExpC4, hmlsq, kExpC3); - const float tmp1 = mad(tmp0, hmlsq, kExpC2); - const float tmp2 = mad(tmp1, hmlsq, kExpC1); - const float tmp3 = mad(tmp2, hmlsq, kExpC0); - const float tmp4 = hml - hmlsq * tmp3; - const float tmp5 = hml*tmp4/(2.0f-tmp4); - const float tmp6 = 1.0f - ( (lo - tmp5) - hi); - const float result = ldexp(tmp6, int32_t(kk) ); + Ty magnitude = Ty(1.0); + + if (absA < Traits::kHypSmall) + { + const Ty zz = square(absA); + + Ty acc = Traits::kTanhPoly[0]; + + for (uint32_t ii = 1, num = uint32_t(BX_COUNTOF(Traits::kTanhPoly) ); ii < num; ++ii) + { + acc = mad(zz, acc, Traits::kTanhPoly[ii]); + } + + const Ty poly = mad(zz, acc, Ty(1.0) ); + + magnitude = absA * poly; + } + else if (absA < Traits::kTanhLarge) + { + const Ty ee = exp(Ty(2.0)*absA); + const Ty num = ee - Ty(1.0); + const Ty den = ee + Ty(1.0); + + magnitude = num / den; + } + + return signBit(_a) ? -magnitude : magnitude; + } + + inline BX_CONSTEXPR_FUNC float sinh(float _a) + { + return sinhT(_a); + } + + inline BX_CONSTEXPR_FUNC double sinh(double _a) + { + return sinhT(_a); + } + + inline BX_CONSTEXPR_FUNC float cosh(float _a) + { + return coshT(_a); + } + + inline BX_CONSTEXPR_FUNC double cosh(double _a) + { + return coshT(_a); + } + + inline BX_CONSTEXPR_FUNC float tanh(float _a) + { + return tanhT(_a); + } + + inline BX_CONSTEXPR_FUNC double tanh(double _a) + { + return tanhT(_a); + } + + template + inline BX_CONSTEXPR_FUNC Ty ldexpT(Ty _a, int32_t _b) + { + using Traits = FloatT; + using Bits = typename Traits::Bits; + + constexpr int32_t kMinStep = 1 - Traits::kExponentBias; + constexpr int32_t kMaxStep = Traits::kExponentBias; + constexpr int32_t kLimit = 2*Traits::kExponentBias + Traits::kMantissaNumBits + 2; + + const int32_t total = clamp(_b, -kLimit, kLimit); + + if (Ty(0.0) == _a + || isNan(_a) + || isInfinite(_a) ) + { + return _a; + } + + int32_t expA = 0; + frexp(_a, &expA); + + if (expA + total > Traits::kExponentBias + 1) + { + return signBit(_a) ? -Traits::kInfinity : Traits::kInfinity; + } + + const int32_t step0 = clamp(total, kMinStep, kMaxStep); + const int32_t rest0 = total - step0; + const int32_t step1 = clamp(rest0, kMinStep, kMaxStep); + const int32_t rest1 = rest0 - step1; + const int32_t step2 = clamp(rest1, kMinStep, kMaxStep); + + const Bits exp0 = Bits(step0 + Traits::kExponentBias); + const Bits exp1 = Bits(step1 + Traits::kExponentBias); + const Bits exp2 = Bits(step2 + Traits::kExponentBias); + + const Bits bits0 = exp0 << Traits::kExponentBitShift; + const Bits bits1 = exp1 << Traits::kExponentBitShift; + const Bits bits2 = exp2 << Traits::kExponentBitShift; + + const Ty scale0 = Traits::fromBits(bits0); + const Ty scale1 = Traits::fromBits(bits1); + const Ty scale2 = Traits::fromBits(bits2); + + const Ty tmp0 = _a * scale0; + const Ty tmp1 = tmp0 * scale1; + const Ty result = tmp1 * scale2; return result; } - inline BX_CONSTEXPR_FUNC float pow(float _a, float _b) + template + inline BX_CONSTEXPR_FUNC Ty logT(Ty _a) { - if (abs(_b) < kFloatSmallest) + using Traits = FloatT; + using Bits = typename Traits::Bits; + + if (isNan(_a) ) { - return 1.0f; + return _a; } - if (abs(_a) < kFloatSmallest) + if (_a < Ty(0.0) ) { - return 0.0f; + constexpr Bits kNanBits = Traits::kSignMask | Traits::kExponentMask | Traits::kMantissaMask; + + return Traits::fromBits(kNanBits); } - return copySign(exp(_b * log(abs(_a) ) ), _a); + if (Ty(0.0) == _a) + { + return -Traits::kInfinity; + } + + if (isInfinite(_a) ) + { + return _a; + } + + constexpr Ty kSqrt2Half = Ty(7.07106781186547524401e-01); + constexpr Ty kLogC0 = Ty(6.666666666666735130e-01); + constexpr Ty kLogC1 = Ty(3.999999999940941908e-01); + constexpr Ty kLogC2 = Ty(2.857142874366239149e-01); + constexpr Ty kLogC3 = Ty(2.222219843214978396e-01); + constexpr Ty kLogC4 = Ty(1.818357216161805012e-01); + constexpr Ty kLogC5 = Ty(1.531383769920937332e-01); + constexpr Ty kLogC6 = Ty(1.479819860511658591e-01); + + int32_t exp = 0; + Ty ff = frexp(_a, &exp); + + if (ff < kSqrt2Half) + { + ff += ff; + --exp; + } + + ff -= Ty(1.0); + + const Ty kk = Ty(exp); + const Ty hi = kk * Traits::kLn2Hi; + const Ty lo = kk * Traits::kLn2Lo; + const Ty ss = ff / (Ty(2.0) + ff); + const Ty s2 = square(ss); + const Ty s4 = square(s2); + + const Ty tmp0 = mad(kLogC6, s4, kLogC4); + const Ty tmp1 = mad(tmp0, s4, kLogC2); + const Ty tmp2 = mad(tmp1, s4, kLogC0); + const Ty t1 = s2*tmp2; + + const Ty tmp3 = mad(kLogC5, s4, kLogC3); + const Ty tmp4 = mad(tmp3, s4, kLogC1); + const Ty t2 = s4*tmp4; + + const Ty t12 = t1 + t2; + const Ty fsq = square(ff); + const Ty hfsq = Ty(0.5)*fsq; + const Ty result = hi - ( (hfsq - (ss*(hfsq+t12) + lo) ) - ff); + + return result; } - inline BX_CONSTEXPR_FUNC float exp2(float _a) + template + inline BX_CONSTEXPR_FUNC Ty log2T(Ty _a) { - return pow(2.0f, _a); + using Traits = FloatT; + + const Ty ln = log(_a); + const Ty result = ln * Traits::kInvLn2; + + return result; + } + + template + inline BX_CONSTEXPR_FUNC Ty log10T(Ty _a) + { + using Traits = FloatT; + + const Ty ln = log(_a); + const Ty result = ln * Traits::kInvLn10; + + return result; + } + + template + inline BX_CONSTEXPR_FUNC Ty log1pT(Ty _a) + { + using Traits = FloatT; + using Bits = typename Traits::Bits; + + if (Ty(0.0) == _a) + { + return _a; + } + + if (Ty(-1.0) > _a) + { + constexpr Bits kNanBits = Traits::kSignMask | Traits::kExponentMask | Traits::kMantissaMask; + + return Traits::fromBits(kNanBits); + } + + if (Ty(-1.0) == _a) + { + return -Traits::kInfinity; + } + + const Ty yy = Ty(1.0) + _a; + const Ty zz = yy - Ty(1.0); + const Ty ln = log(yy); + const Ty err = (zz - _a) / yy; + const Ty result = ln - err; + + return result; + } + + template + inline BX_CONSTEXPR_FUNC Ty expT(Ty _a) + { + using Traits = FloatT; + + if (isNan(_a) ) + { + return _a; + } + + if (_a > Traits::kExpOverflow) + { + return Traits::kInfinity; + } + + if (_a < Traits::kExpUnderflow) + { + return Ty(0.0); + } + + const Ty absA = abs(_a); + + if (absA <= Ty(kNearZero) ) + { + return _a + Ty(1.0); + } + + constexpr Ty kExpC0 = Ty( 1.66666666666666019037e-01); + constexpr Ty kExpC1 = Ty(-2.77777777770155933842e-03); + constexpr Ty kExpC2 = Ty( 6.61375632143793436117e-05); + constexpr Ty kExpC3 = Ty(-1.65339022054652515390e-06); + constexpr Ty kExpC4 = Ty( 4.13813679705723846039e-08); + + const Ty scaled = _a * Traits::kInvLn2; + const Ty kk = round(scaled); + const Ty hi = _a - kk*Traits::kLn2Hi; + const Ty lo = kk*Traits::kLn2Lo; + const Ty hml = hi - lo; + const Ty hmlsq = square(hml); + const Ty tmp0 = mad(kExpC4, hmlsq, kExpC3); + const Ty tmp1 = mad(tmp0, hmlsq, kExpC2); + const Ty tmp2 = mad(tmp1, hmlsq, kExpC1); + const Ty tmp3 = mad(tmp2, hmlsq, kExpC0); + const Ty tmp4 = hml - hmlsq * tmp3; + const Ty tmp5 = hml*tmp4/(Ty(2.0) - tmp4); + const Ty tmp6 = Ty(1.0) - ( (lo - tmp5) - hi); + const Ty result = ldexpT(tmp6, int32_t(kk) ); + + return result; + } + + template + inline BX_CONSTEXPR_FUNC Ty powT(Ty _a, Ty _b) + { + using Traits = FloatT; + using Bits = typename Traits::Bits; + + constexpr Bits kNanBits = Traits::kExponentMask | Traits::kMantissaMask; + + if (Ty(0.0) == _b) + { + return Ty(1.0); + } + + if (isNan(_a) + || isNan(_b) ) + { + return Traits::fromBits(kNanBits); + } + + const Ty magnitude = abs(_a); + const Ty halfB = Ty(0.5) * _b; + const Ty truncB = trunc(_b); + const Ty truncHalf = trunc(halfB); + const bool integral = _b == truncB; + const bool odd = halfB != truncHalf; + + if (_a < Ty(0.0) + && !integral) + { + return Traits::fromBits(kNanBits); + } + + if (Ty(0.0) == magnitude) + { + return _b < Ty(0.0) ? Traits::kInfinity : Ty(0.0); + } + + constexpr Ty kMaxIntegral = Ty(2147483648.0); + + const Ty absB = abs(_b); + + Ty mag = Ty(0.0); + + if (integral + && absB < kMaxIntegral) + { + Ty acc = Ty(1.0); + Ty base = magnitude; + int32_t count = int32_t(absB); + + while (count > 0) + { + if (0 != (count & 1) ) + { + acc = acc * base; + } + + base = base * base; + count >>= 1; + } + + mag = _b < Ty(0.0) ? Ty(1.0) / acc : acc; + } + else + { + const Ty ln = log(magnitude); + const Ty scaled = _b * ln; + + mag = exp(scaled); + } + + if (_a >= Ty(0.0) ) + { + return mag; + } + + return odd ? -mag : mag; + } + + inline BX_CONSTEXPR_FUNC float mad(float _a, float _b, float _c) + { + const simd32_t aa = simd32_ld(_a); + const simd32_t bb = simd32_ld(_b); + const simd32_t cc = simd32_ld(_c); + const simd32_t result = simd32_f32_madd(aa, bb, cc); + + return bitCast(result); + } + + inline BX_CONSTEXPR_FUNC double mad(double _a, double _b, double _c) + { + return _a*_b + _c; + } + + inline BX_CONSTEXPR_FUNC float nms(float _a, float _b, float _c) + { + const float na = -_a; + const float result = mad(na, _b, _c); + + return result; + } + + inline BX_CONSTEXPR_FUNC double nms(double _a, double _b, double _c) + { + const double na = -_a; + const double result = mad(na, _b, _c); + + return result; + } + + inline BX_CONSTEXPR_FUNC float abs(float _a) + { + const uint32_t bits = floatToBits(_a); + const uint32_t magnitude = bits & ~kFloatSignMask; + const float result = bitsToFloat(magnitude); + + return result; + } + + inline BX_CONSTEXPR_FUNC double abs(double _a) + { + const uint64_t bits = doubleToBits(_a); + const uint64_t magnitude = bits & ~kDoubleSignMask; + const double result = bitsToDouble(magnitude); + + return result; + } + + template + requires (isInteger() ) + inline constexpr Ty abs(Ty _a) + { + if constexpr (isSigned() ) + { + const Ty negated = Ty(-_a); + const Ty result = _a < Ty(0) ? negated : _a; + + return result; + } + else + { + return _a; + } + } + + inline BX_CONSTEXPR_FUNC float square(float _a) + { + return _a * _a; + } + + inline BX_CONSTEXPR_FUNC double square(double _a) + { + return _a * _a; + } + + inline BX_CONSTEXPR_FUNC float sign(float _a) + { + return float( (0.0f < _a) - (0.0f > _a) ); + } + + inline BX_CONSTEXPR_FUNC double sign(double _a) + { + const int32_t positive = 0.0 < _a; + const int32_t negative = 0.0 > _a; + const int32_t signum = positive - negative; + const double result = double(signum); + + return result; + } + + inline BX_CONSTEXPR_FUNC bool signBit(double _a) + { + const uint64_t bits = doubleToBits(_a); + + return 0 != (bits & kDoubleSignMask); + } + + inline BX_CONSTEXPR_FUNC float copySign(float _value, float _sign) + { +#if BX_COMPILER_MSVC + const uint32_t magnitude = floatToBits(_value) & ~kFloatSignMask; + const uint32_t sign = floatToBits(_sign) & kFloatSignMask; + const uint32_t bits = magnitude | sign; + const float result = bitsToFloat(bits); + return result; +#else + return __builtin_copysign(_value, _sign); +#endif // BX_COMPILER_MSVC + } + + inline BX_CONSTEXPR_FUNC double copySign(double _value, double _sign) + { + const uint64_t magnitude = doubleToBits(_value) & ~kDoubleSignMask; + const uint64_t sign = doubleToBits(_sign) & kDoubleSignMask; + const uint64_t bits = magnitude | sign; + const double result = bitsToDouble(bits); + + return result; + } + + inline BX_CONSTEXPR_FUNC float trunc(float _a) + { +#if BX_SIMD_SUPPORTED + if (isConstantEvaluated() ) + { + return truncRef(_a); + } + + return truncSimd(_a); +#else + return truncRef(_a); +#endif // BX_SIMD_SUPPORTED + } + + inline BX_CONSTEXPR_FUNC double trunc(double _a) + { + return truncRef(_a); + } + + inline BX_CONSTEXPR_FUNC float floor(float _a) + { +#if BX_SIMD_SUPPORTED + if (isConstantEvaluated() ) + { + return floorRef(_a); + } + + return floorSimd(_a); +#else + return floorRef(_a); +#endif // BX_SIMD_SUPPORTED + } + + inline BX_CONSTEXPR_FUNC double floor(double _a) + { +#if BX_SIMD_SUPPORTED + if (isConstantEvaluated() ) + { + return floorRef(_a); + } + + return floorSimd(_a); +#else + return floorRef(_a); +#endif // BX_SIMD_SUPPORTED + } + + inline BX_CONSTEXPR_FUNC float ceil(float _a) + { +#if BX_SIMD_SUPPORTED + if (isConstantEvaluated() ) + { + return ceilRef(_a); + } + + return ceilSimd(_a); +#else + return ceilRef(_a); +#endif // BX_SIMD_SUPPORTED + } + + inline BX_CONSTEXPR_FUNC double ceil(double _a) + { +#if BX_SIMD_SUPPORTED + if (isConstantEvaluated() ) + { + return ceilRef(_a); + } + + return ceilSimd(_a); +#else + return ceilRef(_a); +#endif // BX_SIMD_SUPPORTED + } + + inline BX_CONSTEXPR_FUNC float round(float _a) + { +#if BX_SIMD_SUPPORTED + if (isConstantEvaluated() ) + { + return roundRef(_a); + } + + return roundSimd(_a); +#else + return roundRef(_a); +#endif // BX_SIMD_SUPPORTED + } + + inline BX_CONSTEXPR_FUNC double round(double _a) + { +#if BX_SIMD_SUPPORTED + if (isConstantEvaluated() ) + { + return roundRef(_a); + } + + return roundSimd(_a); +#else + return roundRef(_a); +#endif // BX_SIMD_SUPPORTED + } + + inline BX_CONSTEXPR_FUNC float fract(float _a) + { + return fractT(_a); + } + + inline BX_CONSTEXPR_FUNC double fract(double _a) + { + return fractT(_a); + } + + inline BX_CONSTEXPR_FUNC bool isEqual(double _a, double _b, double _epsilon) + { + const double diff = _a - _b; + const double lhs = abs(diff); + const double absA = abs(_a); + const double absB = abs(_b); + const double scale = max(1.0, absA, absB); + const double rhs = _epsilon * scale; + + return lhs <= rhs; + } + + inline BX_CONSTEXPR_FUNC float ldexp(float _a, int32_t _b) + { + return ldexpT(_a, _b); + } + + inline BX_CONSTEXPR_FUNC double ldexp(double _a, int32_t _b) + { + return ldexpT(_a, _b); + } + + inline constexpr float frexp(float _a, int32_t* _outExp) + { + const uint32_t bits = floatToBits(_a); + const uint32_t magnitude = bits & ~kFloatSignMask; + + if (0 == magnitude + || kFloatExponentMask <= magnitude) + { + *_outExp = 0; + + return _a; + } + + constexpr int32_t kSubnormalShift = int32_t(kFloatMantissaNumBits) + 2; + constexpr uint32_t kSubnormalExp = uint32_t(int32_t(kFloatExponentBias) + kSubnormalShift); + constexpr uint32_t kSubnormalBits = kSubnormalExp << kFloatExponentBitShift; + constexpr uint32_t kNormalSmallest = uint32_t(1) << kFloatExponentBitShift; + constexpr uint32_t kHalfExp = kFloatExponentBias - 1; + constexpr uint32_t kHalfBits = kHalfExp << kFloatExponentBitShift; + + const bool subnormal = magnitude < kNormalSmallest; + const float scale = bitsToFloat(kSubnormalBits); + const float scaled = _a * scale; + const float value = subnormal ? scaled : _a; + const int32_t bias = subnormal ? -kSubnormalShift : 0; + + const uint32_t valueBits = floatToBits(value); + const uint32_t expBits = valueBits & kFloatExponentMask; + const int32_t raw = int32_t(expBits >> kFloatExponentBitShift); + const int32_t exp = raw - int32_t(kFloatExponentBias) + 1; + const uint32_t kept = valueBits & (kFloatSignMask | kFloatMantissaMask); + const uint32_t resBits = kept | kHalfBits; + const float result = bitsToFloat(resBits); + + *_outExp = exp + bias; + + return result; + } + + inline constexpr double frexp(double _a, int32_t* _outExp) + { + const uint64_t bits = doubleToBits(_a); + const uint64_t magnitude = bits & ~kDoubleSignMask; + + if (0 == magnitude + || kDoubleExponentMask <= magnitude) + { + *_outExp = 0; + + return _a; + } + + constexpr int32_t kSubnormalShift = int32_t(kDoubleMantissaNumBits) + 2; + constexpr uint64_t kSubnormalExp = uint64_t(int32_t(kDoubleExponentBias) + kSubnormalShift); + constexpr uint64_t kSubnormalBits = kSubnormalExp << kDoubleExponentShift; + constexpr uint64_t kNormalSmallest = uint64_t(1) << kDoubleExponentShift; + constexpr uint64_t kHalfExp = kDoubleExponentBias - 1; + constexpr uint64_t kHalfBits = kHalfExp << kDoubleExponentShift; + + const bool subnormal = magnitude < kNormalSmallest; + const double scale = bitsToDouble(kSubnormalBits); + const double scaled = _a * scale; + const double value = subnormal ? scaled : _a; + const int32_t bias = subnormal ? -kSubnormalShift : 0; + + const uint64_t valueBits = doubleToBits(value); + const uint64_t expBits = valueBits & kDoubleExponentMask; + const int32_t raw = int32_t(expBits >> kDoubleExponentShift); + const int32_t exp = raw - int32_t(kDoubleExponentBias) + 1; + const uint64_t kept = valueBits & (kDoubleSignMask | kDoubleMantissaMask); + const uint64_t resBits = kept | kHalfBits; + const double result = bitsToDouble(resBits); + + *_outExp = exp + bias; + + return result; + } + + inline constexpr float modf(float _a, float* _outIntegral) + { + if (isInfinite(_a) ) + { + const uint32_t sign = floatToBits(_a) & kFloatSignMask; + + *_outIntegral = _a; + + return bitsToFloat(sign); + } + + const float integral = trunc(_a); + const float fraction = _a - integral; + + *_outIntegral = integral; + + return fraction; + } + + inline constexpr double modf(double _a, double* _outIntegral) + { + if (isInfinite(_a) ) + { + const uint64_t sign = doubleToBits(_a) & kDoubleSignMask; + + *_outIntegral = _a; + + return bitsToDouble(sign); + } + + const double integral = trunc(_a); + const double fraction = _a - integral; + + *_outIntegral = integral; + + return fraction; + } + + inline BX_CONSTEXPR_FUNC float fmod(float _a, float _b) + { + return fmodT(_a, _b); + } + + inline BX_CONSTEXPR_FUNC double fmod(double _a, double _b) + { + return fmodT(_a, _b); + } + + inline BX_CONSTEXPR_FUNC float mod(float _a, float _b) + { + return modT(_a, _b); + } + + inline BX_CONSTEXPR_FUNC double mod(double _a, double _b) + { + return modT(_a, _b); + } + + inline BX_CONSTEXPR_FUNC float log(float _a) + { + return logT(_a); + } + + inline BX_CONSTEXPR_FUNC double log(double _a) + { + return logT(_a); } inline BX_CONSTEXPR_FUNC float log2(float _a) { - return log(_a) * kInvLogNat2; + return log2T(_a); + } + + inline BX_CONSTEXPR_FUNC double log2(double _a) + { + return log2T(_a); + } + + inline BX_CONSTEXPR_FUNC float exp(float _a) + { + return expT(_a); + } + + inline BX_CONSTEXPR_FUNC double exp(double _a) + { + return expT(_a); + } + + inline BX_CONSTEXPR_FUNC float pow(float _a, float _b) + { + return powT(_a, _b); + } + + inline BX_CONSTEXPR_FUNC double pow(double _a, double _b) + { + return powT(_a, _b); + } + + inline BX_CONSTEXPR_FUNC float exp2(float _a) + { + return exp2T(_a); + } + + inline BX_CONSTEXPR_FUNC double exp2(double _a) + { + return exp2T(_a); + } + + template + inline BX_CONSTEXPR_FUNC Ty sqrtT(Ty _a) + { + using Traits = FloatT; + using Bits = typename Traits::Bits; + + if (isNan(_a) ) + { + return _a; + } + + if (_a < Ty(0.0) ) + { + constexpr Bits kNanBits = Traits::kExponentMask | Traits::kMantissaMask; + + return Traits::fromBits(kNanBits); + } + + if (Ty(0.0) == _a) + { + return _a; + } + + if (isInfinite(_a) ) + { + return _a; + } + + int32_t exp = 0; + const Ty frac = frexp(_a, &exp); + + const bool odd = 0 != (exp & 1); + const Ty mant = odd ? frac + frac : frac; + const int32_t rest = odd ? (exp - 1) >> 1 : exp >> 1; + + const Bits mantBits = Traits::toBits(mant); + const Bits seedBits = Traits::kRsqrtSeed - (mantBits >> 1); + + Ty rsq = Traits::fromBits(seedBits); + + for (uint32_t ii = 0; ii < 5; ++ii) + { + const Ty rsqSq = square(rsq); + const Ty scaled = mant * rsqSq; + const Ty halved = Ty(0.5) * scaled; + const Ty corr = Ty(1.5) - halved; + + rsq = rsq * corr; + } + + const Ty approx = mant * rsq; + const Ty quot = mant / approx; + const Ty sum = approx + quot; + const Ty root = Ty(0.5) * sum; + const Ty result = ldexpT(root, rest); + + return result; + } + + template + inline BX_CONSTEXPR_FUNC Ty rsqrtT(Ty _a) + { + using Traits = FloatT; + + if (Ty(0.0) == _a) + { + return Traits::kInfinity; + } + + const Ty root = sqrt(_a); + const Ty result = Ty(1.0) / root; + + return result; + } + + inline BX_CONSTEXPR_FUNC uint64_t mulHi64(uint64_t _a, uint64_t _b) + { + constexpr uint64_t kLoMask = UINT64_C(0xffffffff); + + const uint64_t aLo = _a & kLoMask; + const uint64_t aHi = _a >> 32; + const uint64_t bLo = _b & kLoMask; + const uint64_t bHi = _b >> 32; + + const uint64_t ll = aLo * bLo; + const uint64_t lh = aLo * bHi; + const uint64_t hl = aHi * bLo; + const uint64_t hh = aHi * bHi; + + const uint64_t mid = (ll >> 32) + (lh & kLoMask) + (hl & kLoMask); + const uint64_t result = hh + (lh >> 32) + (hl >> 32) + (mid >> 32); + + return result; + } + + inline BX_CONSTEXPR_FUNC uint64_t twoOverPiWord(int32_t _index) + { + constexpr uint64_t kTwoOverPi[] = + { + UINT64_C(0xa2f9836e4e441529), + UINT64_C(0xfc2757d1f534ddc0), + UINT64_C(0xdb6295993c439041), + UINT64_C(0xfe5163abdebbc561), + UINT64_C(0xb7246e3a424dd2e0), + UINT64_C(0x06492eea09d1921c), + UINT64_C(0xfe1deb1cb129a73e), + UINT64_C(0xe88235f52ebb4484), + UINT64_C(0xe99c7026b45f7e41), + UINT64_C(0x3991d639835339f4), + UINT64_C(0x9c845f8bbdf9283b), + UINT64_C(0x1ff897ffde05980f), + UINT64_C(0xef2f118b5a0a6d1f), + UINT64_C(0x6d367ecf27cb09b7), + UINT64_C(0x4f463f669e5fea2d), + UINT64_C(0x7527bac7ebe5f17b), + UINT64_C(0x3d0739f78a5292ea), + UINT64_C(0x6bfb5fb11f8d5d08), + UINT64_C(0x56033046fc7b6bab), + UINT64_C(0xf0cfbc209af4361d), + UINT64_C(0xa9e391615ee61b08), + UINT64_C(0x6599855f14a06840), + UINT64_C(0x8dffd8804d732731), + UINT64_C(0x06061556ca73a8c9), + }; + constexpr int32_t kNumWords = 24; + + if (0 > _index + || kNumWords <= _index) + { + return 0; + } + + return kTwoOverPi[uint32_t(_index)]; + } + + inline BX_CONSTEXPR_FUNC uint64_t twoOverPiBit(int32_t _bit) + { + if (1 > _bit) + { + return 0; + } + + const int32_t index = (_bit - 1) / 64; + const int32_t shift = 63 - ( (_bit - 1) % 64); + const uint64_t word = twoOverPiWord(index); + + return (word >> shift) & 1; + } + + inline BX_CONSTEXPR_FUNC uint64_t twoOverPiWindow(int32_t _bit, int32_t _which) + { + const int32_t start = _bit + 64*_which; + const int32_t index = (0 < start) ? (start - 1) / 64 : -( (1 - start) / 64 + 1); + const int32_t shift = start - 1 - 64*index; + + const uint64_t lo = twoOverPiWord(index); + const uint64_t hi = twoOverPiWord(index + 1); + + if (0 == shift) + { + return lo; + } + + return (lo << shift) | (hi >> (64 - shift) ); + } + + inline constexpr double reducePiHalfLarge(double _absA, int32_t* _outQuadrant) + { + constexpr double kPiHalfHi = 1.57079632679489655800e+00; + constexpr double kPiHalfLo = 6.12323399573676603587e-17; + + int32_t exp = 0; + + const double frac = frexp(_absA, &exp); + const double mantF = ldexp(frac, 53); + const uint64_t mant = uint64_t(mantF); + const int32_t scale = exp - 53; + + const uint64_t intBit1 = twoOverPiBit(scale - 1); + const uint64_t intBit0 = twoOverPiBit(scale); + const uint64_t intMod4 = 2*intBit1 + intBit0; + + const uint64_t v0 = twoOverPiWindow(scale + 1, 0); + const uint64_t v1 = twoOverPiWindow(scale + 1, 1); + const uint64_t v2 = twoOverPiWindow(scale + 1, 2); + + const uint64_t h0 = mulHi64(mant, v0); + const uint64_t l0 = mant * v0; + const uint64_t h1 = mulHi64(mant, v1); + const uint64_t l1 = mant * v1; + const uint64_t h2 = mulHi64(mant, v2); + + const uint64_t s2 = l1 + h2; + const uint64_t c2 = (s2 < l1) ? 1 : 0; + const uint64_t t1 = l0 + h1; + const uint64_t k1 = (t1 < l0) ? 1 : 0; + const uint64_t s1 = t1 + c2; + const uint64_t k2 = (s1 < t1) ? 1 : 0; + const uint64_t carry = h0 + k1 + k2; + + const uint64_t quad = mant*intMod4 + carry; + + constexpr uint64_t kLowMask = UINT64_C(0x7ff); + + const double hiPart = double(s1 >> 11) * 0x1p-53; + const double loPart = double(s1 & kLowMask) * 0x1p-64 + double(s2 >> 11) * 0x1p-117; + + double fracHi = hiPart; + int32_t quadrant = int32_t(quad & 3); + + if (0.5 <= fracHi + loPart) + { + fracHi -= 1.0; + quadrant += 1; + } + + const double term0 = fracHi * kPiHalfHi; + const double term1 = loPart * kPiHalfHi; + const double term2 = fracHi * kPiHalfLo; + const double result = term0 + (term1 + term2); + + *_outQuadrant = quadrant & 3; + + return result; + } + + inline constexpr double reducePiHalf(double _a, int32_t* _outQuadrant) + { + constexpr double kInvPiHalf = 6.36619772367581382433e-01; + constexpr double kPiHalf1 = 1.57079632673412561417e+00; + constexpr double kPiHalf2 = 6.07710050630396597660e-11; + constexpr double kPiHalf3 = 2.02226624871116645580e-21; + constexpr double kPiHalf3t = 8.47842766036889956997e-32; + constexpr double kLimit = 1647099.0; + + const double absA = abs(_a); + + if (kLimit <= absA) + { + int32_t quadrant = 0; + + const double reduced = reducePiHalfLarge(absA, &quadrant); + + if (signBit(_a) ) + { + *_outQuadrant = (-quadrant) & 3; + + return -reduced; + } + + *_outQuadrant = quadrant; + + return reduced; + } + + const double value = _a; + + const double scaled = value * kInvPiHalf; + const double quadrant = round(scaled); + + const double r0 = value - quadrant*kPiHalf1; + + const double w2 = quadrant*kPiHalf2; + const double r2 = r0 - w2; + const double e2 = (r0 - r2) - w2; + + const double w3 = quadrant*kPiHalf3; + const double r3 = r2 - w3; + const double e3 = (r2 - r3) - w3; + + const double tail = quadrant*kPiHalf3t; + const double acc = e2 + e3; + const double cc = tail - acc; + + const double result = r3 - cc; + + *_outQuadrant = int32_t(quadrant) & 3; + + return result; + } + + inline BX_CONSTEXPR_FUNC double sinPiQuarter(double _a) + { + constexpr double kSinC0 = -1.66666666666666324348e-01; + constexpr double kSinC1 = 8.33333333332248946124e-03; + constexpr double kSinC2 = -1.98412698298579493134e-04; + constexpr double kSinC3 = 2.75573137070700676789e-06; + constexpr double kSinC4 = -2.50507602534068634195e-08; + constexpr double kSinC5 = 1.58969099521155010221e-10; + + const double zz = square(_a); + const double vv = zz * _a; + const double zz2 = square(zz); + const double zz4 = square(zz2); + const double p01 = mad(kSinC1, zz, kSinC0); + const double p23 = mad(kSinC3, zz, kSinC2); + const double p45 = mad(kSinC5, zz, kSinC4); + const double plo = mad(p23, zz2, p01); + const double poly = mad(p45, zz4, plo); + + const double result = mad(vv, poly, _a); + + return result; + } + + inline BX_CONSTEXPR_FUNC double cosPiQuarter(double _a) + { + constexpr double kCosC0 = 4.16666666666666019037e-02; + constexpr double kCosC1 = -1.38888888888741095749e-03; + constexpr double kCosC2 = 2.48015872894767294178e-05; + constexpr double kCosC3 = -2.75573143513906633035e-07; + constexpr double kCosC4 = 2.08757232129817482790e-09; + constexpr double kCosC5 = -1.13596475577881948265e-11; + + const double zz = square(_a); + const double zz2 = square(zz); + const double zz4 = square(zz2); + const double q01 = mad(kCosC1, zz, kCosC0); + const double q23 = mad(kCosC3, zz, kCosC2); + const double q45 = mad(kCosC5, zz, kCosC4); + const double qlo = mad(q23, zz2, q01); + const double qq = mad(q45, zz4, qlo); + const double zq = zz * qq; + const double poly = zz * zq; + + const double halfZz = 0.5 * zz; + const double ww = 1.0 - halfZz; + const double err = (1.0 - ww) - halfZz; + const double corr = err + poly; + + const double result = ww + corr; + + return result; + } + + template + inline BX_CONSTEXPR_FUNC Ty sinT(Ty _a) + { + using Traits = FloatT; + using Bits = typename Traits::Bits; + + if (Ty(0.0) == _a) + { + return _a; + } + + if (isNan(_a) + || isInfinite(_a) ) + { + constexpr Bits kNanBits = Traits::kExponentMask | Traits::kMantissaMask; + + return Traits::fromBits(kNanBits); + } + + int32_t quadrant = 0; + + const double rr = reducePiHalf(double(_a), &quadrant); + + switch (quadrant) + { + case 0: return Ty( sinPiQuarter(rr) ); + case 1: return Ty( cosPiQuarter(rr) ); + case 2: return Ty(-sinPiQuarter(rr) ); + default: return Ty(-cosPiQuarter(rr) ); + } + } + + template + inline BX_CONSTEXPR_FUNC Ty cosT(Ty _a) + { + using Traits = FloatT; + using Bits = typename Traits::Bits; + + if (isNan(_a) + || isInfinite(_a) ) + { + constexpr Bits kNanBits = Traits::kExponentMask | Traits::kMantissaMask; + + return Traits::fromBits(kNanBits); + } + + int32_t quadrant = 0; + + const double rr = reducePiHalf(double(_a), &quadrant); + + switch (quadrant) + { + case 0: return Ty( cosPiQuarter(rr) ); + case 1: return Ty(-sinPiQuarter(rr) ); + case 2: return Ty(-cosPiQuarter(rr) ); + default: return Ty( sinPiQuarter(rr) ); + } + } + + inline BX_CONSTEXPR_FUNC float sin(float _a) + { + return sinT(_a); + } + + inline BX_CONSTEXPR_FUNC double sin(double _a) + { + return sinT(_a); + } + + inline BX_CONSTEXPR_FUNC float cos(float _a) + { + return cosT(_a); + } + + inline BX_CONSTEXPR_FUNC double cos(double _a) + { + return cosT(_a); + } + + inline BX_CONSTEXPR_FUNC float tan(float _a) + { + return tanT(_a); + } + + inline BX_CONSTEXPR_FUNC double tan(double _a) + { + return tanT(_a); + } + + inline BX_CONSTEXPR_FUNC double asinRational(double _t) + { + constexpr double kPs0 = 1.66666666666666657415e-01; + constexpr double kPs1 = -3.25565818622400915405e-01; + constexpr double kPs2 = 2.01212532134862925881e-01; + constexpr double kPs3 = -4.00555345006794114027e-02; + constexpr double kPs4 = 7.91534994289814532176e-04; + constexpr double kPs5 = 3.47933107596021167570e-05; + constexpr double kQs1 = -2.40339491173441421878e+00; + constexpr double kQs2 = 2.02094576023350569471e+00; + constexpr double kQs3 = -6.88283971605453293030e-01; + constexpr double kQs4 = 7.70381505559019352791e-02; + + const double p0 = mad(kPs5, _t, kPs4); + const double p1 = mad(p0, _t, kPs3); + const double p2 = mad(p1, _t, kPs2); + const double p3 = mad(p2, _t, kPs1); + const double p4 = mad(p3, _t, kPs0); + const double pp = _t * p4; + + const double q0 = mad(kQs4, _t, kQs3); + const double q1 = mad(q0, _t, kQs2); + const double q2 = mad(q1, _t, kQs1); + const double qq = mad(q2, _t, 1.0); + + const double result = pp / qq; + + return result; + } + + inline BX_CONSTEXPR_FUNC double clearLowWord(double _a) + { + constexpr uint64_t kHighMask = UINT64_C(0xffffffff00000000); + + const uint64_t bits = doubleToBits(_a); + const uint64_t high = bits & kHighMask; + + return bitsToDouble(high); + } + + template + inline BX_CONSTEXPR_FUNC Ty asinT(Ty _a) + { + using Traits = FloatT; + using Bits = typename Traits::Bits; + + constexpr double kPiHalfHi = 1.57079632679489655800e+00; + constexpr double kPiHalfLo = 6.12323399573676603587e-17; + constexpr double kPiQuartHi = 7.85398163397448278999e-01; + constexpr double kTiny = 7.45058059692383e-09; + constexpr double kNearOne = 0.975; + + if (isNan(_a) ) + { + return _a; + } + + const double aa = double(_a); + const double absA = abs(aa); + + if (absA > 1.0) + { + constexpr Bits kNanBits = Traits::kExponentMask | Traits::kMantissaMask; + + return Traits::fromBits(kNanBits); + } + + if (1.0 == absA) + { + const double hi = aa * kPiHalfHi; + const double lo = aa * kPiHalfLo; + + return Ty(hi + lo); + } + + if (absA < 0.5) + { + if (absA < kTiny) + { + return _a; + } + + const double tt = square(aa); + const double ww = asinRational(tt); + + return Ty(mad(aa, ww, aa) ); + } + + const double ww = 1.0 - absA; + const double tt = 0.5 * ww; + const double rr = asinRational(tt); + const double ss = sqrt(tt); + + double result = 0.0; + + if (absA >= kNearOne) + { + const double sw = mad(ss, rr, ss); + const double tmp0 = 2.0*sw - kPiHalfLo; + + result = kPiHalfHi - tmp0; + } + else + { + const double df = clearLowWord(ss); + const double cc = (tt - df*df) / (ss + df); + const double pp = 2.0*ss*rr - (kPiHalfLo - 2.0*cc); + const double qq = kPiQuartHi - 2.0*df; + + result = kPiQuartHi - (pp - qq); + } + + return aa > 0.0 ? Ty(result) : Ty(-result); + } + + template + inline BX_CONSTEXPR_FUNC Ty acosT(Ty _a) + { + using Traits = FloatT; + using Bits = typename Traits::Bits; + + constexpr double kPiHalfHi = 1.57079632679489655800e+00; + constexpr double kPiHalfLo = 6.12323399573676603587e-17; + constexpr double kPiFull = 3.14159265358979311600e+00; + constexpr double kTiny = 6.93889390390722838e-18; + + if (isNan(_a) ) + { + return _a; + } + + const double aa = double(_a); + const double absA = abs(aa); + + if (absA > 1.0) + { + constexpr Bits kNanBits = Traits::kExponentMask | Traits::kMantissaMask; + + return Traits::fromBits(kNanBits); + } + + if (1.0 == absA) + { + return aa > 0.0 ? Ty(0.0) : Ty(kPiFull); + } + + if (absA < 0.5) + { + if (absA <= kTiny) + { + return Ty(kPiHalfHi + kPiHalfLo); + } + + const double zz = square(aa); + const double rr = asinRational(zz); + const double xr = aa * rr; + + return Ty(kPiHalfHi - (aa - (kPiHalfLo - xr) ) ); + } + + if (aa < 0.0) + { + const double zz = 0.5 * (1.0 + aa); + const double rr = asinRational(zz); + const double ss = sqrt(zz); + const double ww = rr*ss - kPiHalfLo; + const double sw = ss + ww; + + return Ty(kPiFull - 2.0*sw); + } + + const double zz = 0.5 * (1.0 - aa); + const double ss = sqrt(zz); + const double df = clearLowWord(ss); + const double cc = (zz - df*df) / (ss + df); + const double rr = asinRational(zz); + const double ww = mad(rr, ss, cc); + const double dw = df + ww; + + return Ty(2.0 * dw); + } + + inline BX_CONSTEXPR_FUNC float asin(float _a) + { + return asinT(_a); + } + + inline BX_CONSTEXPR_FUNC double asin(double _a) + { + return asinT(_a); + } + + inline BX_CONSTEXPR_FUNC float acos(float _a) + { + return acosT(_a); + } + + inline BX_CONSTEXPR_FUNC double acos(double _a) + { + return acosT(_a); + } + + template + inline BX_CONSTEXPR_FUNC Ty atanT(Ty _a) + { + constexpr double kAtanHi[] = + { + 4.63647609000806093515e-01, + 7.85398163397448278999e-01, + 9.82793723247329054082e-01, + 1.57079632679489655800e+00, + }; + constexpr double kAtanLo[] = + { + 2.26987774529616870924e-17, + 3.06161699786838301793e-17, + 1.39033110312309984516e-17, + 6.12323399573676603587e-17, + }; + + constexpr double kAt0 = 3.33333333333329318027e-01; + constexpr double kAt1 = -1.99999999998764832476e-01; + constexpr double kAt2 = 1.42857142725034663711e-01; + constexpr double kAt3 = -1.11111104054623557880e-01; + constexpr double kAt4 = 9.09088713343650656196e-02; + constexpr double kAt5 = -7.69187620504482999495e-02; + constexpr double kAt6 = 6.66107313738753120669e-02; + constexpr double kAt7 = -5.83357013379057348645e-02; + constexpr double kAt8 = 4.97687799461593236017e-02; + constexpr double kAt9 = -3.65315727442169155270e-02; + constexpr double kAt10 = 1.62858201153657823623e-02; + + constexpr double kHuge = 7.37869762948382e+19; + constexpr double kTiny = 1.86264514923096e-09; + + if (isNan(_a) ) + { + return _a; + } + + const double aa = double(_a); + const double absA = abs(aa); + + if (absA > kHuge) + { + const double big = kAtanHi[3] + kAtanLo[3]; + + return aa > 0.0 ? Ty(big) : Ty(-big); + } + + int32_t id = -1; + double xx = aa; + + if (absA >= 0.4375) + { + if (absA < 1.1875) + { + if (absA < 0.6875) + { + id = 0; + xx = (2.0*absA - 1.0) / (2.0 + absA); + } + else + { + id = 1; + xx = (absA - 1.0) / (absA + 1.0); + } + } + else + { + if (absA < 2.4375) + { + id = 2; + xx = (absA - 1.5) / (1.0 + 1.5*absA); + } + else + { + id = 3; + xx = -1.0 / absA; + } + } + } + else if (absA < kTiny) + { + return _a; + } + + const double zz = square(xx); + const double ww = square(zz); + + const double o0 = mad(kAt10, ww, kAt8); + const double o1 = mad(o0, ww, kAt6); + const double o2 = mad(o1, ww, kAt4); + const double o3 = mad(o2, ww, kAt2); + const double o4 = mad(o3, ww, kAt0); + const double s1 = zz * o4; + + const double e0 = mad(kAt9, ww, kAt7); + const double e1 = mad(e0, ww, kAt5); + const double e2 = mad(e1, ww, kAt3); + const double e3 = mad(e2, ww, kAt1); + const double s2 = ww * e3; + + const double sum = s1 + s2; + + if (id < 0) + { + return Ty(xx - xx*sum); + } + + const double hi = kAtanHi[id]; + const double lo = kAtanLo[id]; + + const double corr = (xx*sum - lo) - xx; + const double result = hi - corr; + + return aa < 0.0 ? Ty(-result) : Ty(result); + } + + template + inline BX_CONSTEXPR_FUNC Ty atan2T(Ty _y, Ty _x) + { + using Traits = FloatT; + using Bits = typename Traits::Bits; + + constexpr double kPiFull = 3.14159265358979311600e+00; + constexpr double kPiLo = 1.22464679914735317722e-16; + constexpr double kHalfPi = 1.57079632679489655800e+00; + constexpr double kQuarterPi = 7.85398163397448278999e-01; + + if (isNan(_x) + || isNan(_y) ) + { + constexpr Bits kNanBits = Traits::kExponentMask | Traits::kMantissaMask; + + return Traits::fromBits(kNanBits); + } + + const bool signY = signBit(_y); + const bool signX = signBit(_x); + + if (Ty(0.0) == _y) + { + if (!signX) + { + return signY ? Ty(-0.0) : Ty(0.0); + } + + return signY ? Ty(-kPiFull) : Ty(kPiFull); + } + + if (Ty(0.0) == _x) + { + return signY ? Ty(-kHalfPi) : Ty(kHalfPi); + } + + if (isInfinite(_x) ) + { + if (isInfinite(_y) ) + { + const double quarter = signX ? 3.0*kQuarterPi : kQuarterPi; + + return signY ? Ty(-quarter) : Ty(quarter); + } + + const double straight = signX ? kPiFull : 0.0; + + return signY ? Ty(-straight) : Ty(straight); + } + + if (isInfinite(_y) ) + { + return signY ? Ty(-kHalfPi) : Ty(kHalfPi); + } + + const double ratio = double(_y) / double(_x); + const double zz = atanT(abs(ratio) ); + + if (!signX) + { + return signY ? Ty(-zz) : Ty(zz); + } + + const double folded = kPiFull - (zz - kPiLo); + + return signY ? Ty(-folded) : Ty(folded); + } + + inline BX_CONSTEXPR_FUNC float atan(float _a) + { + return atanT(_a); + } + + inline BX_CONSTEXPR_FUNC double atan(double _a) + { + return atanT(_a); + } + + inline BX_CONSTEXPR_FUNC float atan2(float _y, float _x) + { + return atan2T(_y, _x); + } + + inline BX_CONSTEXPR_FUNC double atan2(double _y, double _x) + { + return atan2T(_y, _x); + } + + template + inline BX_CONSTEXPR_FUNC Ty asinhT(Ty _a) + { + using Traits = FloatT; + + constexpr Ty kLn2 = Ty(6.93147180559945286227e-01); + + if (isNan(_a) + || isInfinite(_a) ) + { + return _a; + } + + const Ty absA = abs(_a); + + if (absA < Traits::kArcHypTiny) + { + return _a; + } + + Ty magnitude = Ty(0.0); + + if (absA > Traits::kArcHypHuge) + { + const Ty ln = log(absA); + + magnitude = ln + kLn2; + } + else if (absA > Ty(2.0) ) + { + const Ty tt = mad(absA, absA, Ty(1.0) ); + const Ty root = sqrt(tt); + const Ty tail = Ty(1.0) / (root + absA); + const Ty arg = Ty(2.0)*absA + tail; + + magnitude = log(arg); + } + else + { + const Ty tt = square(absA); + const Ty root = sqrt(Ty(1.0) + tt); + const Ty tail = tt / (Ty(1.0) + root); + + magnitude = log1p(absA + tail); + } + + return signBit(_a) ? -magnitude : magnitude; + } + + template + inline BX_CONSTEXPR_FUNC Ty acoshT(Ty _a) + { + using Traits = FloatT; + using Bits = typename Traits::Bits; + + constexpr Ty kLn2 = Ty(6.93147180559945286227e-01); + + if (isNan(_a) ) + { + return _a; + } + + if (Ty(1.0) > _a) + { + constexpr Bits kNanBits = Traits::kExponentMask | Traits::kMantissaMask; + + return Traits::fromBits(kNanBits); + } + + if (Ty(1.0) == _a) + { + return Ty(0.0); + } + + if (isInfinite(_a) ) + { + return _a; + } + + if (_a > Traits::kArcHypHuge) + { + const Ty ln = log(_a); + + return ln + kLn2; + } + + if (_a > Ty(2.0) ) + { + const Ty tt = square(_a); + const Ty root = sqrt(tt - Ty(1.0) ); + const Ty tail = Ty(1.0) / (_a + root); + const Ty arg = Ty(2.0)*_a - tail; + + return log(arg); + } + + const Ty tt = _a - Ty(1.0); + const Ty root = sqrt(Ty(2.0)*tt + square(tt) ); + + return log1p(tt + root); + } + + template + inline BX_CONSTEXPR_FUNC Ty atanhT(Ty _a) + { + using Traits = FloatT; + using Bits = typename Traits::Bits; + + if (isNan(_a) ) + { + return _a; + } + + const Ty absA = abs(_a); + + if (absA > Ty(1.0) ) + { + constexpr Bits kNanBits = Traits::kExponentMask | Traits::kMantissaMask; + + return Traits::fromBits(kNanBits); + } + + if (Ty(1.0) == absA) + { + return signBit(_a) ? -Traits::kInfinity : Traits::kInfinity; + } + + if (absA < Traits::kArcHypTiny) + { + return _a; + } + + Ty magnitude = Ty(0.0); + + if (absA < Ty(0.5) ) + { + const Ty twice = absA + absA; + const Ty tail = twice*absA / (Ty(1.0) - absA); + const Ty half = log1p(twice + tail); + + magnitude = Ty(0.5) * half; + } + else + { + const Ty twice = absA + absA; + const Ty ratio = twice / (Ty(1.0) - absA); + + magnitude = Ty(0.5) * log1p(ratio); + } + + return signBit(_a) ? -magnitude : magnitude; + } + + inline BX_CONSTEXPR_FUNC float log10(float _a) + { + return log10T(_a); + } + + inline BX_CONSTEXPR_FUNC double log10(double _a) + { + return log10T(_a); + } + + inline BX_CONSTEXPR_FUNC float log1p(float _a) + { + return log1pT(_a); + } + + inline BX_CONSTEXPR_FUNC double log1p(double _a) + { + return log1pT(_a); + } + + inline BX_CONSTEXPR_FUNC float asinh(float _a) + { + return asinhT(_a); + } + + inline BX_CONSTEXPR_FUNC double asinh(double _a) + { + return asinhT(_a); + } + + inline BX_CONSTEXPR_FUNC float acosh(float _a) + { + return acoshT(_a); + } + + inline BX_CONSTEXPR_FUNC double acosh(double _a) + { + return acoshT(_a); + } + + inline BX_CONSTEXPR_FUNC float atanh(float _a) + { + return atanhT(_a); + } + inline BX_CONSTEXPR_FUNC double atanh(double _a) + { + return atanhT(_a); } BX_FP_PRECISE_END() @@ -999,12 +2670,7 @@ BX_FP_PRECISE_END() inline BX_CONSTEXPR_FUNC float rsqrtRef(float _a) { - if (_a < kFloatSmallest) - { - return kFloatInfinity; - } - - return pow(_a, -0.5f); + return rsqrtT(_a); } inline BX_CONST_FUNC float rsqrtSimd(float _a) @@ -1025,12 +2691,12 @@ BX_FP_PRECISE_END() inline BX_CONSTEXPR_FUNC float sqrtRef(float _a) { - if (_a < 0.0f) - { - return bitsToFloat(kFloatExponentMask | kFloatMantissaMask); - } + return sqrtT(_a); + } - return _a * pow(_a, -0.5f); + inline BX_CONSTEXPR_FUNC double sqrtRef(double _a) + { + return sqrtT(_a); } inline BX_CONST_FUNC float sqrtSimd(float _a) @@ -1053,6 +2719,22 @@ BX_FP_PRECISE_END() return result; } + inline BX_CONST_FUNC double sqrtSimd(double _a) + { + if (_a < 0.0) + { + return bitsToDouble(kDoubleExponentMask | kDoubleMantissaMask); + } + + const simd128_t aa = simd_splat(_a); + const simd128_t sqrt = simd_f64_sqrt(aa); + + alignas(16) double result[2] = { 0.0, 0.0 }; + simd_st(result, sqrt); + + return result[0]; + } + inline BX_CONSTEXPR_FUNC float rsqrt(float _a) { #if BX_SIMD_SUPPORTED @@ -1067,6 +2749,11 @@ BX_FP_PRECISE_END() #endif // BX_SIMD_SUPPORTED } + inline BX_CONSTEXPR_FUNC double rsqrt(double _a) + { + return rsqrtT(_a); + } + inline BX_CONSTEXPR_FUNC float sqrt(float _a) { #if BX_SIMD_SUPPORTED @@ -1081,6 +2768,20 @@ BX_FP_PRECISE_END() #endif // BX_SIMD_SUPPORTED } + inline BX_CONSTEXPR_FUNC double sqrt(double _a) + { +#if BX_SIMD_SUPPORTED + if (isConstantEvaluated() ) + { + return sqrtRef(_a); + } + + return sqrtSimd(_a); +#else + return sqrtRef(_a); +#endif // BX_SIMD_SUPPORTED + } + inline BX_CONSTEXPR_FUNC bool isEqual(float _a, float _b, float _epsilon) { // Reference(s): diff --git a/include/bx/inline/simd128_ref.inl b/include/bx/inline/simd128_ref.inl index 35c816a..1b986ac 100644 --- a/include/bx/inline/simd128_ref.inl +++ b/include/bx/inline/simd128_ref.inl @@ -17,6 +17,7 @@ namespace bx BX_ALIGN_DECL(16, struct) simd128_i64_ref_t { int64_t i64[2]; }; BX_ALIGN_DECL(16, struct) simd128_u8_ref_t { uint8_t u8[16]; }; BX_ALIGN_DECL(16, struct) simd128_u16_ref_t { uint16_t u16[8]; }; + BX_ALIGN_DECL(16, struct) simd128_u32_ref_t { uint32_t u32[4]; }; BX_ALIGN_DECL(16, struct) simd128_u64_ref_t { uint64_t u64[2]; }; #if BX_SIMD_LANGEXT @@ -672,15 +673,15 @@ namespace bx inline BX_CONSTEXPR_FUNC simd128_ref_t simd128_i32_add(simd128_ref_t _a, simd128_ref_t _b) { #if BX_SIMD_LANGEXT - const simd128_i32_langext_t a = bitCast(_a); - const simd128_i32_langext_t b = bitCast(_b); - const simd128_i32_langext_t sum = a + b; + const simd128_u32_langext_t a = bitCast(_a); + const simd128_u32_langext_t b = bitCast(_b); + const simd128_u32_langext_t sum = a + b; const simd128_ref_t result = bitCast(sum); return result; #else - const simd128_i32_ref_t a = bitCast(_a); - const simd128_i32_ref_t b = bitCast(_b); - const simd128_i32_ref_t result = { { a.i32[0] + b.i32[0], a.i32[1] + b.i32[1], a.i32[2] + b.i32[2], a.i32[3] + b.i32[3] } }; + const simd128_u32_ref_t a = bitCast(_a); + const simd128_u32_ref_t b = bitCast(_b); + const simd128_u32_ref_t result = { { a.u32[0] + b.u32[0], a.u32[1] + b.u32[1], a.u32[2] + b.u32[2], a.u32[3] + b.u32[3] } }; return bitCast(result); #endif // BX_SIMD_LANGEXT } @@ -689,15 +690,15 @@ namespace bx inline BX_CONSTEXPR_FUNC simd128_ref_t simd128_i32_sub(simd128_ref_t _a, simd128_ref_t _b) { #if BX_SIMD_LANGEXT - const simd128_i32_langext_t a = bitCast(_a); - const simd128_i32_langext_t b = bitCast(_b); - const simd128_i32_langext_t diff = a - b; + const simd128_u32_langext_t a = bitCast(_a); + const simd128_u32_langext_t b = bitCast(_b); + const simd128_u32_langext_t diff = a - b; const simd128_ref_t result = bitCast(diff); return result; #else - const simd128_i32_ref_t a = bitCast(_a); - const simd128_i32_ref_t b = bitCast(_b); - const simd128_i32_ref_t result = { { a.i32[0] - b.i32[0], a.i32[1] - b.i32[1], a.i32[2] - b.i32[2], a.i32[3] - b.i32[3] } }; + const simd128_u32_ref_t a = bitCast(_a); + const simd128_u32_ref_t b = bitCast(_b); + const simd128_u32_ref_t result = { { a.u32[0] - b.u32[0], a.u32[1] - b.u32[1], a.u32[2] - b.u32[2], a.u32[3] - b.u32[3] } }; return bitCast(result); #endif // BX_SIMD_LANGEXT } @@ -985,16 +986,16 @@ namespace bx inline BX_CONSTEXPR_FUNC simd128_ref_t simd128_i16_add(simd128_ref_t _a, simd128_ref_t _b) { #if BX_SIMD_LANGEXT - const simd128_i16_langext_t a = bitCast(_a); - const simd128_i16_langext_t b = bitCast(_b); - const simd128_i16_langext_t sum = a + b; + const simd128_u16_langext_t a = bitCast(_a); + const simd128_u16_langext_t b = bitCast(_b); + const simd128_u16_langext_t sum = a + b; const simd128_ref_t result = bitCast(sum); return result; #else - const simd128_i16_ref_t a = bitCast(_a); - const simd128_i16_ref_t b = bitCast(_b); - simd128_i16_ref_t result; - for (int ii = 0; ii < 8; ++ii) { result.i16[ii] = a.i16[ii] + b.i16[ii]; } + const simd128_u16_ref_t a = bitCast(_a); + const simd128_u16_ref_t b = bitCast(_b); + simd128_u16_ref_t result; + for (int ii = 0; ii < 8; ++ii) { result.u16[ii] = a.u16[ii] + b.u16[ii]; } return bitCast(result); #endif // BX_SIMD_LANGEXT } @@ -1003,16 +1004,16 @@ namespace bx inline BX_CONSTEXPR_FUNC simd128_ref_t simd128_i16_sub(simd128_ref_t _a, simd128_ref_t _b) { #if BX_SIMD_LANGEXT - const simd128_i16_langext_t a = bitCast(_a); - const simd128_i16_langext_t b = bitCast(_b); - const simd128_i16_langext_t diff = a - b; + const simd128_u16_langext_t a = bitCast(_a); + const simd128_u16_langext_t b = bitCast(_b); + const simd128_u16_langext_t diff = a - b; const simd128_ref_t result = bitCast(diff); return result; #else - const simd128_i16_ref_t a = bitCast(_a); - const simd128_i16_ref_t b = bitCast(_b); - simd128_i16_ref_t result; - for (int ii = 0; ii < 8; ++ii) { result.i16[ii] = a.i16[ii] - b.i16[ii]; } + const simd128_u16_ref_t a = bitCast(_a); + const simd128_u16_ref_t b = bitCast(_b); + simd128_u16_ref_t result; + for (int ii = 0; ii < 8; ++ii) { result.u16[ii] = a.u16[ii] - b.u16[ii]; } return bitCast(result); #endif // BX_SIMD_LANGEXT } @@ -1021,16 +1022,16 @@ namespace bx inline BX_CONSTEXPR_FUNC simd128_ref_t simd128_i16_mullo(simd128_ref_t _a, simd128_ref_t _b) { #if BX_SIMD_LANGEXT - const simd128_i16_langext_t a = bitCast(_a); - const simd128_i16_langext_t b = bitCast(_b); - const simd128_i16_langext_t prod = a * b; + const simd128_u16_langext_t a = bitCast(_a); + const simd128_u16_langext_t b = bitCast(_b); + const simd128_u16_langext_t prod = a * b; const simd128_ref_t result = bitCast(prod); return result; #else - const simd128_i16_ref_t a = bitCast(_a); - const simd128_i16_ref_t b = bitCast(_b); - simd128_i16_ref_t result; - for (int ii = 0; ii < 8; ++ii) { result.i16[ii] = (int16_t)(a.i16[ii] * b.i16[ii]); } + const simd128_u16_ref_t a = bitCast(_a); + const simd128_u16_ref_t b = bitCast(_b); + simd128_u16_ref_t result; + for (int ii = 0; ii < 8; ++ii) { result.u16[ii] = uint16_t(a.u16[ii] * b.u16[ii]); } return bitCast(result); #endif // BX_SIMD_LANGEXT } @@ -1103,16 +1104,16 @@ namespace bx inline BX_CONSTEXPR_FUNC simd128_ref_t simd128_i8_add(simd128_ref_t _a, simd128_ref_t _b) { #if BX_SIMD_LANGEXT - const simd128_i8_langext_t a = bitCast(_a); - const simd128_i8_langext_t b = bitCast(_b); - const simd128_i8_langext_t sum = a + b; + const simd128_u8_langext_t a = bitCast(_a); + const simd128_u8_langext_t b = bitCast(_b); + const simd128_u8_langext_t sum = a + b; const simd128_ref_t result = bitCast(sum); return result; #else - const simd128_i8_ref_t a = bitCast(_a); - const simd128_i8_ref_t b = bitCast(_b); - simd128_i8_ref_t result; - for (int ii = 0; ii < 16; ++ii) { result.i8[ii] = a.i8[ii] + b.i8[ii]; } + const simd128_u8_ref_t a = bitCast(_a); + const simd128_u8_ref_t b = bitCast(_b); + simd128_u8_ref_t result; + for (int ii = 0; ii < 16; ++ii) { result.u8[ii] = a.u8[ii] + b.u8[ii]; } return bitCast(result); #endif // BX_SIMD_LANGEXT } @@ -1121,16 +1122,16 @@ namespace bx inline BX_CONSTEXPR_FUNC simd128_ref_t simd128_i8_sub(simd128_ref_t _a, simd128_ref_t _b) { #if BX_SIMD_LANGEXT - const simd128_i8_langext_t a = bitCast(_a); - const simd128_i8_langext_t b = bitCast(_b); - const simd128_i8_langext_t diff = a - b; + const simd128_u8_langext_t a = bitCast(_a); + const simd128_u8_langext_t b = bitCast(_b); + const simd128_u8_langext_t diff = a - b; const simd128_ref_t result = bitCast(diff); return result; #else - const simd128_i8_ref_t a = bitCast(_a); - const simd128_i8_ref_t b = bitCast(_b); - simd128_i8_ref_t result; - for (int ii = 0; ii < 16; ++ii) { result.i8[ii] = a.i8[ii] - b.i8[ii]; } + const simd128_u8_ref_t a = bitCast(_a); + const simd128_u8_ref_t b = bitCast(_b); + simd128_u8_ref_t result; + for (int ii = 0; ii < 16; ++ii) { result.u8[ii] = a.u8[ii] - b.u8[ii]; } return bitCast(result); #endif // BX_SIMD_LANGEXT } @@ -1870,15 +1871,15 @@ namespace bx inline BX_CONSTEXPR_FUNC simd128_ref_t simd128_i64_add(simd128_ref_t _a, simd128_ref_t _b) { #if BX_SIMD_LANGEXT - const simd128_i64_langext_t a = bitCast(_a); - const simd128_i64_langext_t b = bitCast(_b); - const simd128_i64_langext_t sum = a + b; + const simd128_u64_langext_t a = bitCast(_a); + const simd128_u64_langext_t b = bitCast(_b); + const simd128_u64_langext_t sum = a + b; const simd128_ref_t result = bitCast(sum); return result; #else - const simd128_i64_ref_t a = bitCast(_a); - const simd128_i64_ref_t b = bitCast(_b); - const simd128_i64_ref_t result = { { a.i64[0] + b.i64[0], a.i64[1] + b.i64[1] } }; + const simd128_u64_ref_t a = bitCast(_a); + const simd128_u64_ref_t b = bitCast(_b); + const simd128_u64_ref_t result = { { a.u64[0] + b.u64[0], a.u64[1] + b.u64[1] } }; return bitCast(result); #endif // BX_SIMD_LANGEXT } @@ -1887,15 +1888,15 @@ namespace bx inline BX_CONSTEXPR_FUNC simd128_ref_t simd128_i64_sub(simd128_ref_t _a, simd128_ref_t _b) { #if BX_SIMD_LANGEXT - const simd128_i64_langext_t a = bitCast(_a); - const simd128_i64_langext_t b = bitCast(_b); - const simd128_i64_langext_t diff = a - b; + const simd128_u64_langext_t a = bitCast(_a); + const simd128_u64_langext_t b = bitCast(_b); + const simd128_u64_langext_t diff = a - b; const simd128_ref_t result = bitCast(diff); return result; #else - const simd128_i64_ref_t a = bitCast(_a); - const simd128_i64_ref_t b = bitCast(_b); - const simd128_i64_ref_t result = { { a.i64[0] - b.i64[0], a.i64[1] - b.i64[1] } }; + const simd128_u64_ref_t a = bitCast(_a); + const simd128_u64_ref_t b = bitCast(_b); + const simd128_u64_ref_t result = { { a.u64[0] - b.u64[0], a.u64[1] - b.u64[1] } }; return bitCast(result); #endif // BX_SIMD_LANGEXT } diff --git a/include/bx/inline/simd256_ref.inl b/include/bx/inline/simd256_ref.inl index 8a79603..54e36a0 100644 --- a/include/bx/inline/simd256_ref.inl +++ b/include/bx/inline/simd256_ref.inl @@ -520,9 +520,9 @@ namespace bx inline BX_CONSTEXPR_FUNC simd256_ref_t simd256_i32_add(simd256_ref_t _a, simd256_ref_t _b) { #if BX_SIMD256_LANGEXT - const simd256_i32_langext_t a = bitCast(_a); - const simd256_i32_langext_t b = bitCast(_b); - const simd256_i32_langext_t sum = a + b; + const simd256_u32_langext_t a = bitCast(_a); + const simd256_u32_langext_t b = bitCast(_b); + const simd256_u32_langext_t sum = a + b; const simd256_ref_t result = bitCast(sum); return result; #else @@ -537,9 +537,9 @@ namespace bx inline BX_CONSTEXPR_FUNC simd256_ref_t simd256_i32_sub(simd256_ref_t _a, simd256_ref_t _b) { #if BX_SIMD256_LANGEXT - const simd256_i32_langext_t a = bitCast(_a); - const simd256_i32_langext_t b = bitCast(_b); - const simd256_i32_langext_t diff = a - b; + const simd256_u32_langext_t a = bitCast(_a); + const simd256_u32_langext_t b = bitCast(_b); + const simd256_u32_langext_t diff = a - b; const simd256_ref_t result = bitCast(diff); return result; #else @@ -554,8 +554,8 @@ namespace bx inline BX_CONST_FUNC simd256_ref_t simd256_i32_neg(simd256_ref_t _a) { #if BX_SIMD256_LANGEXT - const simd256_i32_langext_t a = bitCast(_a); - const simd256_i32_langext_t neg = -a; + const simd256_u32_langext_t a = bitCast(_a); + const simd256_u32_langext_t neg = -a; const simd256_ref_t result = bitCast(neg); return result; #else @@ -1315,9 +1315,9 @@ namespace bx inline BX_CONSTEXPR_FUNC simd256_ref_t simd256_i64_add(simd256_ref_t _a, simd256_ref_t _b) { #if BX_SIMD256_LANGEXT - const simd256_i64_langext_t a = bitCast(_a); - const simd256_i64_langext_t b = bitCast(_b); - const simd256_i64_langext_t sum = a + b; + const simd256_u64_langext_t a = bitCast(_a); + const simd256_u64_langext_t b = bitCast(_b); + const simd256_u64_langext_t sum = a + b; const simd256_ref_t result = bitCast(sum); return result; #else @@ -1332,9 +1332,9 @@ namespace bx inline BX_CONSTEXPR_FUNC simd256_ref_t simd256_i64_sub(simd256_ref_t _a, simd256_ref_t _b) { #if BX_SIMD256_LANGEXT - const simd256_i64_langext_t a = bitCast(_a); - const simd256_i64_langext_t b = bitCast(_b); - const simd256_i64_langext_t diff = a - b; + const simd256_u64_langext_t a = bitCast(_a); + const simd256_u64_langext_t b = bitCast(_b); + const simd256_u64_langext_t diff = a - b; const simd256_ref_t result = bitCast(diff); return result; #else @@ -1383,9 +1383,9 @@ namespace bx inline BX_CONSTEXPR_FUNC simd256_ref_t simd256_i16_add(simd256_ref_t _a, simd256_ref_t _b) { #if BX_SIMD256_LANGEXT - const simd256_i16_langext_t a = bitCast(_a); - const simd256_i16_langext_t b = bitCast(_b); - const simd256_i16_langext_t sum = a + b; + const simd256_u16_langext_t a = bitCast(_a); + const simd256_u16_langext_t b = bitCast(_b); + const simd256_u16_langext_t sum = a + b; const simd256_ref_t result = bitCast(sum); return result; #else @@ -1400,9 +1400,9 @@ namespace bx inline BX_CONSTEXPR_FUNC simd256_ref_t simd256_i16_sub(simd256_ref_t _a, simd256_ref_t _b) { #if BX_SIMD256_LANGEXT - const simd256_i16_langext_t a = bitCast(_a); - const simd256_i16_langext_t b = bitCast(_b); - const simd256_i16_langext_t diff = a - b; + const simd256_u16_langext_t a = bitCast(_a); + const simd256_u16_langext_t b = bitCast(_b); + const simd256_u16_langext_t diff = a - b; const simd256_ref_t result = bitCast(diff); return result; #else @@ -1417,9 +1417,9 @@ namespace bx inline BX_CONSTEXPR_FUNC simd256_ref_t simd256_i16_mullo(simd256_ref_t _a, simd256_ref_t _b) { #if BX_SIMD256_LANGEXT - const simd256_i16_langext_t a = bitCast(_a); - const simd256_i16_langext_t b = bitCast(_b); - const simd256_i16_langext_t prod = a * b; + const simd256_u16_langext_t a = bitCast(_a); + const simd256_u16_langext_t b = bitCast(_b); + const simd256_u16_langext_t prod = a * b; const simd256_ref_t result = bitCast(prod); return result; #else @@ -1483,9 +1483,9 @@ namespace bx inline BX_CONSTEXPR_FUNC simd256_ref_t simd256_i8_add(simd256_ref_t _a, simd256_ref_t _b) { #if BX_SIMD256_LANGEXT - const simd256_i8_langext_t a = bitCast(_a); - const simd256_i8_langext_t b = bitCast(_b); - const simd256_i8_langext_t sum = a + b; + const simd256_u8_langext_t a = bitCast(_a); + const simd256_u8_langext_t b = bitCast(_b); + const simd256_u8_langext_t sum = a + b; const simd256_ref_t result = bitCast(sum); return result; #else @@ -1500,9 +1500,9 @@ namespace bx inline BX_CONSTEXPR_FUNC simd256_ref_t simd256_i8_sub(simd256_ref_t _a, simd256_ref_t _b) { #if BX_SIMD256_LANGEXT - const simd256_i8_langext_t a = bitCast(_a); - const simd256_i8_langext_t b = bitCast(_b); - const simd256_i8_langext_t diff = a - b; + const simd256_u8_langext_t a = bitCast(_a); + const simd256_u8_langext_t b = bitCast(_b); + const simd256_u8_langext_t diff = a - b; const simd256_ref_t result = bitCast(diff); return result; #else diff --git a/include/bx/inline/simd_ni.inl b/include/bx/inline/simd_ni.inl index 82aca59..68b8458 100644 --- a/include/bx/inline/simd_ni.inl +++ b/include/bx/inline/simd_ni.inl @@ -516,22 +516,50 @@ namespace bx template BX_SIMD_INLINE Ty simd_f32_ldexp_ni(Ty _a, Ty _b) { - const Ty signexpmask = simd_splat(uint32_t(kFloatSignMask | kFloatExponentMask) ); - const Ty mantmask = simd_splat(kFloatMantissaMask); - const Ty masked = simd_and(_a, signexpmask); - const Ty expsign0 = simd_x32_sra(masked, kFloatExponentBitShift); - const Ty tmp = simd_i32_add(expsign0, _b); - const Ty expsign1 = simd_x32_sll(tmp, kFloatExponentBitShift); - const Ty mantissa = simd_and(_a, mantmask); - const Ty result = simd_or(mantissa, expsign1); + constexpr int32_t kMinStep = 1 - int32_t(kFloatExponentBias); + constexpr int32_t kMaxStep = int32_t(kFloatExponentBias); + constexpr int32_t kLimit = 2*int32_t(kFloatExponentBias) + int32_t(kFloatMantissaNumBits) + 2; + + const Ty minStep = simd_splat(uint32_t(kMinStep) ); + const Ty maxStep = simd_splat(uint32_t(kMaxStep) ); + const Ty limitHi = simd_splat(uint32_t(kLimit) ); + const Ty limitLo = simd_splat(uint32_t(-kLimit) ); + const Ty bias = simd_splat(uint32_t(kFloatExponentBias) ); + + const Ty totalLo = simd_i32_max(_b, limitLo); + const Ty total = simd_i32_min(totalLo, limitHi); + + const Ty step0Lo = simd_i32_max(total, minStep); + const Ty step0 = simd_i32_min(step0Lo, maxStep); + const Ty rest0 = simd_i32_sub(total, step0); + + const Ty step1Lo = simd_i32_max(rest0, minStep); + const Ty step1 = simd_i32_min(step1Lo, maxStep); + const Ty rest1 = simd_i32_sub(rest0, step1); + + const Ty step2Lo = simd_i32_max(rest1, minStep); + const Ty step2 = simd_i32_min(step2Lo, maxStep); + + const Ty exp0 = simd_i32_add(step0, bias); + const Ty exp1 = simd_i32_add(step1, bias); + const Ty exp2 = simd_i32_add(step2, bias); + + const Ty scale0 = simd_x32_sll(exp0, kFloatExponentBitShift); + const Ty scale1 = simd_x32_sll(exp1, kFloatExponentBitShift); + const Ty scale2 = simd_x32_sll(exp2, kFloatExponentBitShift); + + const Ty tmp0 = simd_f32_mul(_a, scale0); + const Ty tmp1 = simd_f32_mul(tmp0, scale1); + const Ty result = simd_f32_mul(tmp1, scale2); return result; } BX_FP_PRECISE_BEGIN() + /// cos(_a) for _bias 0, sin(_a) for _bias 3, matching bx::sinCosPiHalf. template - BX_SIMD_INLINE Ty simd_f32_cos_ni(Ty _a) + BX_SIMD_INLINE Ty simd_f32_sincos_pi_half_ni(Ty _a, uint32_t _bias) { const Ty two_over_pi = simd_splat(2.0f * kInvPi); const Ty pi_half = simd_splat(kPiHalf); @@ -540,8 +568,10 @@ BX_FP_PRECISE_BEGIN() const Ty tmp_rp = simd_f32_mul(real, pi_half); const Ty xx = simd_f32_sub(_a, tmp_rp); const Ty ireal = simd_f32_ftoi_trunc(real); + const Ty ibias = simd_splat(_bias); + const Ty ibiased = simd_i32_add(ireal, ibias); const Ty three = simd_splat(uint32_t(3) ); - const Ty ibits = simd_and(ireal, three); + const Ty ibits = simd_and(ibiased, three); const Ty ione = simd_splat(uint32_t(1) ); const Ty izero = simd_zero(); const Ty bit0 = simd_and(ibits, ione); @@ -579,13 +609,16 @@ BX_FP_PRECISE_BEGIN() return result; } + template + BX_SIMD_INLINE Ty simd_f32_cos_ni(Ty _a) + { + return simd_f32_sincos_pi_half_ni(_a, 0); + } + template BX_SIMD_INLINE Ty simd_f32_sin_ni(Ty _a) { - const Ty pi_half = simd_splat(kPiHalf); - const Ty shifted = simd_f32_sub(_a, pi_half); - const Ty result = simd_f32_cos(shifted); - return result; + return simd_f32_sincos_pi_half_ni(_a, 3); } template @@ -745,12 +778,52 @@ BX_FP_PRECISE_BEGIN() template BX_SIMD_INLINE Ty simd_f32_sinh_ni(Ty _a) { - const Ty half = simd_splat(0.5f); - const Ty na = simd_f32_neg(_a); - const Ty ea = simd_f32_exp(_a); - const Ty ena = simd_f32_exp(na); - const Ty diff = simd_f32_sub(ea, ena); - const Ty result = simd_f32_mul(half, diff); + const Ty half = simd_splat(0.5f); + const Ty one = simd_splat(1.0f); + const Ty small = simd_splat(0.0625f); + const Ty big = simd_splat(9.0f); + const Ty noFit = simd_splat(88.0f); + const Ty kC3 = simd_splat(1.0f/6.0f); + const Ty kC5 = simd_splat(1.0f/120.0f); + const Ty kC7 = simd_splat(1.0f/5040.0f); + const Ty signmask = simd_splat(uint32_t(kFloatSignMask) ); + + const Ty absA = simd_f32_abs(_a); + const Ty nabsA = simd_f32_neg(absA); + + const Ty zz = simd_f32_mul(absA, absA); + const Ty inn0 = simd_f32_madd(zz, kC7, kC5); + const Ty inn1 = simd_f32_madd(zz, inn0, kC3); + const Ty poly = simd_f32_madd(zz, inn1, one); + const Ty magSmall = simd_f32_mul(absA, poly); + + const Ty ep = simd_f32_exp(absA); + const Ty en = simd_f32_exp(nabsA); + const Ty diff = simd_f32_sub(ep, en); + const Ty magMid = simd_f32_mul(half, diff); + + const Ty magBig = simd_f32_mul(half, ep); + + constexpr int32_t kShift = 235; + constexpr uint32_t kBits = uint32_t(int32_t(kFloatExponentBias) + kShift/2) << kFloatExponentBitShift; + + const Ty kKLn2 = simd_splat(1.62889587431602213e+02f); + const Ty kScale = simd_splat(kBits); + const Ty reduced = simd_f32_sub(absA, kKLn2); + const Ty ee = simd_f32_exp(reduced); + const Ty once = simd_f32_mul(ee, kScale); + const Ty magNoFit = simd_f32_mul(once, kScale); + + const Ty smallMask = simd_f32_cmplt(absA, small); + const Ty bigMask = simd_f32_cmplt(absA, big); + const Ty noFitMask = simd_f32_cmplt(absA, noFit); + const Ty mag0 = simd_selb(noFitMask, magBig, magNoFit); + const Ty mag1 = simd_selb(bigMask, magMid, mag0); + const Ty magnitude = simd_selb(smallMask, magSmall, mag1); + + const Ty sign = simd_and(_a, signmask); + const Ty result = simd_or(magnitude, sign); + return result; } @@ -769,38 +842,82 @@ BX_FP_PRECISE_BEGIN() template BX_SIMD_INLINE Ty simd_f32_tanh_ni(Ty _a) { - const Ty one = simd_splat(1.0f); - const Ty two = simd_splat(2.0f); - const Ty twoA = simd_f32_mul(two, _a); - const Ty tmp0 = simd_f32_exp(twoA); - const Ty tmp1 = simd_f32_sub(tmp0, one); - const Ty tmp2 = simd_f32_add(tmp0, one); - const Ty result = simd_f32_div(tmp1, tmp2); + const Ty one = simd_splat(1.0f); + const Ty two = simd_splat(2.0f); + const Ty small = simd_splat(0.0625f); + const Ty large = simd_splat(10.0f); + const Ty kC3 = simd_splat( -1.0f/3.0f); + const Ty kC5 = simd_splat( 2.0f/15.0f); + const Ty kC7 = simd_splat(-17.0f/315.0f); + const Ty signmask = simd_splat(uint32_t(kFloatSignMask) ); + + const Ty absA = simd_f32_abs(_a); + const Ty zz = simd_f32_mul(absA, absA); + const Ty inn0 = simd_f32_madd(zz, kC7, kC5); + const Ty inn1 = simd_f32_madd(zz, inn0, kC3); + const Ty poly = simd_f32_madd(zz, inn1, one); + const Ty magSmall = simd_f32_mul(absA, poly); + + const Ty twoA = simd_f32_mul(two, absA); + const Ty ee = simd_f32_exp(twoA); + const Ty num = simd_f32_sub(ee, one); + const Ty den = simd_f32_add(ee, one); + const Ty magMid = simd_f32_div(num, den); + + const Ty smallMask = simd_f32_cmplt(absA, small); + const Ty largeMask = simd_f32_cmplt(absA, large); + const Ty mag0 = simd_selb(largeMask, magMid, one); + const Ty magnitude = simd_selb(smallMask, magSmall, mag0); + + const Ty sign = simd_and(_a, signmask); + const Ty signed0 = simd_or(magnitude, sign); + + const Ty self_eq = simd_f32_cmpeq(_a, _a); + const Ty nan_mask = simd_not(self_eq); + const Ty result = simd_selb(nan_mask, _a, signed0); + return result; } - template BX_SIMD_INLINE Ty simd_f32_log_ni(Ty _a) { const Ty zero = simd_zero(); const Ty one = simd_splat(1.0f); + const Ty nanmask = simd_f32_cmpneq(_a, _a); const Ty negmask = simd_f32_cmplt(_a, zero); const Ty zeromask = simd_f32_cmpeq(_a, zero); + const Ty pos_inf = simd_splat(uint32_t(kFloatExponentMask) ); + const Ty infmask = simd_f32_cmpeq(_a, pos_inf); const Ty nan_val = simd_splat(uint32_t(kFloatSignMask | kFloatExponentMask | kFloatMantissaMask) ); const Ty neg_inf = simd_splat(uint32_t(kFloatSignMask | kFloatExponentMask) ); const Ty expmask_c = simd_splat(kFloatExponentMask); const Ty mantmask_c = simd_splat(kFloatMantissaMask); const Ty exp_half = simd_splat(uint32_t(0x3f000000) ); const Ty bias_126 = simd_splat(uint32_t(kFloatExponentBias - 1) ); - const Ty aexp = simd_and(_a, expmask_c); + + constexpr int32_t kSubnormalShift = int32_t(kFloatMantissaNumBits) + 2; + constexpr uint32_t kSubnormalExp = uint32_t(int32_t(kFloatExponentBias) + kSubnormalShift); + constexpr uint32_t kSubnormalBits = kSubnormalExp << kFloatExponentBitShift; + + const Ty smallest = simd_splat(uint32_t(uint32_t(1) << kFloatExponentBitShift) ); + const Ty magnitude = simd_and(_a, simd_splat(uint32_t(~kFloatSignMask) ) ); + const Ty submask = simd_i32_cmplt(magnitude, smallest); + const Ty subscale = simd_splat(kSubnormalBits); + const Ty subscaled = simd_f32_mul(_a, subscale); + const Ty value = simd_selb(submask, subscaled, _a); + const Ty subshift = simd_splat(uint32_t(kSubnormalShift) ); + const Ty subadj = simd_and(submask, subshift); + + const Ty aexp = simd_and(value, expmask_c); const Ty rawexp = simd_x32_srl(aexp, kFloatExponentBitShift); - const Ty exp_i0 = simd_i32_sub(rawexp, bias_126); - const Ty mant = simd_and(_a, mantmask_c); + const Ty exp_n = simd_i32_sub(rawexp, bias_126); + const Ty exp_i0 = simd_i32_sub(exp_n, subadj); + const Ty mant = simd_and(value, mantmask_c); const Ty ff0 = simd_or(mant, exp_half); - const Ty sqrt2_half = simd_splat(kSqrt2 * 0.5f); + const Ty sqrt2_half = simd_splat(7.07106781186547524401e-01f); const Ty adj_mask = simd_f32_cmplt(ff0, sqrt2_half); const Ty two = simd_splat(2.0f); - const Ty ff_dbl = simd_f32_mul(ff0, two); + const Ty ff_dbl = simd_f32_add(ff0, ff0); const Ty ff1 = simd_selb(adj_mask, ff_dbl, ff0); const Ty ione = simd_splat(uint32_t(1) ); const Ty adj_one = simd_and(adj_mask, ione); @@ -812,8 +929,8 @@ BX_FP_PRECISE_BEGIN() const Ty kC4 = simd_splat(1.818357216161805012e-01f); const Ty kC5 = simd_splat(1.531383769920937332e-01f); const Ty kC6 = simd_splat(1.479819860511658591e-01f); - const Ty log2lo = simd_splat(1.90821492927058770002e-10f); - const Ty lognat2 = simd_splat(kLogNat2); + const Ty log2lo = simd_splat(kFloatLn2Lo); + const Ty lognat2 = simd_splat(kFloatLn2Hi); const Ty half = simd_splat(0.5f); const Ty ff = simd_f32_sub(ff1, one); const Ty kk = simd_i32_itof(exp_i); @@ -840,7 +957,9 @@ BX_FP_PRECISE_BEGIN() const Ty inff = simd_f32_sub(inner, ff); const Ty raw = simd_f32_sub(hi, inff); const Ty rz = simd_selb(zeromask, neg_inf, raw); - const Ty result = simd_selb(negmask, nan_val, rz); + const Ty ri = simd_selb(infmask, _a, rz); + const Ty rn = simd_selb(negmask, nan_val, ri); + const Ty result = simd_selb(nanmask, _a, rn); return result; } @@ -850,29 +969,31 @@ BX_FP_PRECISE_BEGIN() { const Ty zero = simd_zero(); const Ty one = simd_splat(1.0f); + const Ty two = simd_splat(2.0f); + const Ty inf = simd_splat(uint32_t(kFloatExponentMask) ); const Ty near_zero = simd_splat(kNearZero); const Ty absa = simd_f32_abs(_a); const Ty near_mask = simd_f32_cmple(absa, near_zero); const Ty near_result = simd_f32_add(_a, one); - const Ty exp_min = simd_splat(-87.33654475f); - const Ty clamp_mask = simd_f32_cmple(_a, exp_min); - const Ty exp_max = simd_splat(88.72283905f); - const Ty inf = simd_splat(kFloatInfinity); - const Ty inf_mask = simd_f32_cmpge(_a, exp_max); + const Ty overflow = simd_splat( 88.7228394f); + const Ty underflow = simd_splat(-103.972084f); + const Ty over_mask = simd_f32_cmpgt(_a, overflow); + const Ty under_mask = simd_f32_cmplt(_a, underflow); + const Ty self_eq = simd_f32_cmpeq(_a, _a); + const Ty nan_mask = simd_not(self_eq); const Ty kC0 = simd_splat( 1.66666666666666019037e-01f); const Ty kC1 = simd_splat(-2.77777777770155933842e-03f); const Ty kC2 = simd_splat( 6.61375632143793436117e-05f); const Ty kC3 = simd_splat(-1.65339022054652515390e-06f); const Ty kC4 = simd_splat( 4.13813679705723846039e-08f); - const Ty log2lo = simd_splat(1.90821492927058770002e-10f); - const Ty lognat2 = simd_splat(kLogNat2); + const Ty ln2Hi = simd_splat(6.9314575195e-01f); + const Ty ln2Lo = simd_splat(1.4286067653e-06f); const Ty invlognat2 = simd_splat(kInvLogNat2); - const Ty two = simd_splat(2.0f); const Ty amul = simd_f32_mul(_a, invlognat2); const Ty kk = simd_f32_round(amul); - const Ty kkln = simd_f32_mul(kk, lognat2); + const Ty kkln = simd_f32_mul(kk, ln2Hi); const Ty hi = simd_f32_sub(_a, kkln); - const Ty lo = simd_f32_mul(kk, log2lo); + const Ty lo = simd_f32_mul(kk, ln2Lo); const Ty hml = simd_f32_sub(hi, lo); const Ty hmlsq = simd_f32_mul(hml, hml); const Ty tmp0 = simd_f32_madd(kC4, hmlsq, kC3); @@ -889,9 +1010,10 @@ BX_FP_PRECISE_BEGIN() const Ty tmp6 = simd_f32_sub(one, tmp5b); const Ty ikk = simd_f32_ftoi_trunc(kk); const Ty raw = simd_f32_ldexp_ni(tmp6, ikk); - const Ty rc = simd_selb(clamp_mask, zero, raw); - const Ty ri = simd_selb(inf_mask, inf, rc); - const Ty result = simd_selb(near_mask, near_result, ri); + const Ty rNear = simd_selb(near_mask, near_result, raw); + const Ty rUnder = simd_selb(under_mask, zero, rNear); + const Ty rOver = simd_selb(over_mask, inf, rUnder); + const Ty result = simd_selb(nan_mask, _a, rOver); return result; } @@ -919,21 +1041,43 @@ BX_FP_PRECISE_BEGIN() { const Ty zero = simd_zero(); const Ty one = simd_splat(1.0f); - const Ty smallest = simd_splat(kFloatSmallest); - const Ty signmask = simd_splat(kFloatSignMask); + const Ty half = simd_splat(0.5f); + const Ty nan = simd_splat(uint32_t(kFloatExponentMask | kFloatMantissaMask) ); + const Ty inf = simd_splat(uint32_t(kFloatExponentMask) ); const Ty absa = simd_f32_abs(_a); - const Ty absb = simd_f32_abs(_b); const Ty loga = simd_f32_log(absa); const Ty bloga = simd_f32_mul(_b, loga); - const Ty pw = simd_f32_exp(bloga); - const Ty asign = simd_and(_a, signmask); - const Ty pwabs = simd_f32_abs(pw); - const Ty result0 = simd_or(pwabs, asign); - const Ty bmask = simd_f32_cmplt(absb, smallest); - const Ty amask = simd_f32_cmplt(absa, smallest); - const Ty result1 = simd_selb(amask, zero, result0); - const Ty result = simd_selb(bmask, one, result1); + const Ty mag = simd_f32_exp(bloga); + + const Ty halfB = simd_f32_mul(half, _b); + const Ty truncB = simd_f32_trunc(_b); + const Ty truncH = simd_f32_trunc(halfB); + const Ty integral = simd_f32_cmpeq(_b, truncB); + const Ty halfEq = simd_f32_cmpeq(halfB, truncH); + const Ty odd = simd_not(halfEq); + + const Ty negA = simd_f32_cmplt(_a, zero); + const Ty negOdd = simd_and(negA, odd); + const Ty negMag = simd_f32_neg(mag); + const Ty signed0 = simd_selb(negOdd, negMag, mag); + + const Ty negNoInt = simd_andc(negA, integral); + const Ty result0 = simd_selb(negNoInt, nan, signed0); + + const Ty zeroA = simd_f32_cmpeq(absa, zero); + const Ty bNeg = simd_f32_cmplt(_b, zero); + const Ty zeroMag = simd_selb(bNeg, inf, zero); + const Ty result1 = simd_selb(zeroA, zeroMag, result0); + + const Ty zeroB = simd_f32_cmpeq(_b, zero); + const Ty result2 = simd_selb(zeroB, one, result1); + + const Ty aEq = simd_f32_cmpeq(_a, _a); + const Ty bEq = simd_f32_cmpeq(_b, _b); + const Ty bothEq = simd_and(aEq, bEq); + const Ty nanIn = simd_not(bothEq); + const Ty result = simd_selb(nanIn, nan, result2); return result; } diff --git a/include/bx/math.h b/include/bx/math.h index acd4393..182b235 100644 --- a/include/bx/math.h +++ b/include/bx/math.h @@ -168,6 +168,86 @@ namespace bx /// BX_CONSTEXPR_FUNC double bitsToDouble(uint64_t _a); + /// Format constants and bit access for a floating-point type. + template + struct FloatT; + + /// `FloatT` for `float`. + template<> + struct FloatT + { + using Type = float; + using Bits = uint32_t; + using Int = int32_t; + + static constexpr int32_t kExponentBias = int32_t(kFloatExponentBias); + static constexpr int32_t kExponentBitShift = int32_t(kFloatExponentBitShift); + static constexpr int32_t kMantissaNumBits = int32_t(kFloatMantissaNumBits); + static constexpr Bits kSignMask = kFloatSignMask; + static constexpr Bits kExponentMask = kFloatExponentMask; + static constexpr Bits kMantissaMask = kFloatMantissaMask; + static constexpr Type kSmallest = kFloatSmallest; + static constexpr Type kLargest = kFloatLargest; + static constexpr Type kInfinity = kFloatInfinity; + static constexpr Type kInvLn2 = kFloatInvLn2; + static constexpr Type kInvLn10 = kFloatInvLn10; + static constexpr Type kLn2Hi = kFloatLn2Hi; + static constexpr Type kLn2Lo = kFloatLn2Lo; + static constexpr Type kExpOverflow = kFloatExpOverflow; + static constexpr Type kExpUnderflow = kFloatExpUnderflow; + static constexpr int32_t kExpHalfShift = kFloatExpHalfShift; + static constexpr Type kExpHalfLn2 = kFloatExpHalfLn2; + static constexpr Type kHypSmall = kFloatHypSmall; + static constexpr Type kHypBig = kFloatHypBig; + static constexpr Type kHypNoFit = kFloatHypNoFit; + static constexpr Type kTanhLarge = kFloatTanhLarge; + static constexpr Type kArcHypTiny = kFloatArcHypTiny; + static constexpr Type kArcHypHuge = kFloatArcHypHuge; + static constexpr Bits kRsqrtSeed = kFloatRsqrtSeed; + static constexpr Type kTanhPoly[] = { -17.0f/315.0f, 2.0f/15.0f, -1.0f/3.0f }; + + static BX_CONSTEXPR_FUNC Bits toBits(Type _a) { return floatToBits(_a); } + static BX_CONSTEXPR_FUNC Type fromBits(Bits _a) { return bitsToFloat(_a); } + }; + + /// `FloatT` for `double`. + template<> + struct FloatT + { + using Type = double; + using Bits = uint64_t; + using Int = int64_t; + + static constexpr int32_t kExponentBias = int32_t(kDoubleExponentBias); + static constexpr int32_t kExponentBitShift = int32_t(kDoubleExponentShift); + static constexpr int32_t kMantissaNumBits = int32_t(kDoubleMantissaNumBits); + static constexpr Bits kSignMask = kDoubleSignMask; + static constexpr Bits kExponentMask = kDoubleExponentMask; + static constexpr Bits kMantissaMask = kDoubleMantissaMask; + static constexpr Type kSmallest = kDoubleSmallest; + static constexpr Type kLargest = kDoubleLargest; + static constexpr Type kInfinity = kDoubleInfinity; + static constexpr Type kInvLn2 = kDoubleInvLn2; + static constexpr Type kInvLn10 = kDoubleInvLn10; + static constexpr Type kLn2Hi = kDoubleLn2Hi; + static constexpr Type kLn2Lo = kDoubleLn2Lo; + static constexpr Type kExpOverflow = kDoubleExpOverflow; + static constexpr Type kExpUnderflow = kDoubleExpUnderflow; + static constexpr int32_t kExpHalfShift = kDoubleExpHalfShift; + static constexpr Type kExpHalfLn2 = kDoubleExpHalfLn2; + static constexpr Type kHypSmall = kDoubleHypSmall; + static constexpr Type kHypBig = kDoubleHypBig; + static constexpr Type kHypNoFit = kDoubleHypNoFit; + static constexpr Type kTanhLarge = kDoubleTanhLarge; + static constexpr Type kArcHypTiny = kDoubleArcHypTiny; + static constexpr Type kArcHypHuge = kDoubleArcHypHuge; + static constexpr Bits kRsqrtSeed = kDoubleRsqrtSeed; + static constexpr Type kTanhPoly[] = { 62.0/2835.0, -17.0/315.0, 2.0/15.0, -1.0/3.0 }; + + static BX_CONSTEXPR_FUNC Bits toBits(Type _a) { return doubleToBits(_a); } + static BX_CONSTEXPR_FUNC Type fromBits(Bits _a) { return bitsToDouble(_a); } + }; + /// Returns sortable floating point value. /// /// @param[in] _value Float bits as uint32_t. @@ -232,6 +312,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float floor(float _f); + /// Returns the largest integer value not greater than _a. + /// + /// @param[in] _a Input value. + /// + /// @returns Floor of _a. + /// + BX_CONSTEXPR_FUNC double floor(double _a); + /// Returns the smallest integer value not less than _f. /// /// @param[in] _f Input value. @@ -240,7 +328,15 @@ namespace bx /// BX_CONSTEXPR_FUNC float ceil(float _f); - /// Returns the nearest integer value to _f, rounding halfway cases away from zero. + /// Returns the smallest integer value not less than _a. + /// + /// @param[in] _a Input value. + /// + /// @returns Ceiling of _a. + /// + BX_CONSTEXPR_FUNC double ceil(double _a); + + /// Returns the nearest integer value to _f, rounding halfway cases to even. /// /// @param[in] _f Input value. /// @@ -248,6 +344,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float round(float _f); + /// Returns the nearest integer value to _a, rounding halfway cases to even. + /// + /// @param[in] _a Input value. + /// + /// @returns Rounded value. + /// + BX_CONSTEXPR_FUNC double round(double _a); + /// Returns linear interpolation between two values _a and _b. /// /// @param[in] _a Start value. @@ -276,6 +380,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float sign(float _a); + /// Returns the sign of _a as -1.0, 0.0 or 1.0. + /// + /// @param[in] _a Input value. + /// + /// @returns Sign of _a. + /// + BX_CONSTEXPR_FUNC double sign(double _a); + /// Returns `true` if the velue `_a` is negative. /// /// @param[in] _a Value. @@ -293,6 +405,15 @@ namespace bx /// BX_CONSTEXPR_FUNC float copySign(float _value, float _sign); + /// Returns the magnitude of _value with the sign of _sign. + /// + /// @param[in] _value Value to take the magnitude from. + /// @param[in] _sign Value to take the sign from. + /// + /// @returns |_value| with the sign of _sign. + /// + BX_CONSTEXPR_FUNC double copySign(double _value, double _sign); + /// Returns the absolute of _a. /// /// @param[in] _a Input value. @@ -301,6 +422,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float abs(float _a); + /// Returns the absolute value of _a. + /// + /// @param[in] _a Input value. + /// + /// @returns |_a|. + /// + BX_CONSTEXPR_FUNC double abs(double _a); + /// Returns the square of _a. /// /// @param[in] _a Input value. @@ -309,6 +438,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float square(float _a); + /// Returns _a squared. + /// + /// @param[in] _a Input value. + /// + /// @returns _a * _a. + /// + BX_CONSTEXPR_FUNC double square(double _a); + /// Returns the both sine and cosine of the argument _a. /// /// @param[out] _outSinApprox Output approximate sine value. @@ -328,6 +465,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float sin(float _a); + /// Returns the sine of _a. + /// + /// @param[in] _a Angle in radians. + /// + /// @returns sin(_a). + /// + BX_CONSTEXPR_FUNC double sin(double _a); + /// Returns hyperbolic sine of the argument _a. /// /// @param[in] _a Input value. @@ -336,6 +481,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float sinh(float _a); + /// Returns the hyperbolic sine of _a. + /// + /// @param[in] _a Input value. + /// + /// @returns sinh(_a). + /// + BX_CONSTEXPR_FUNC double sinh(double _a); + /// Returns radian angle between -pi/2 and +pi/2 whose sine is _a. /// /// @param[in] _a Sine value [-1, 1]. @@ -344,6 +497,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float asin(float _a); + /// Returns the arc sine of _a. + /// + /// @param[in] _a Input value in [-1, 1]. + /// + /// @returns asin(_a) in radians, or NaN if _a is outside [-1, 1]. + /// + BX_CONSTEXPR_FUNC double asin(double _a); + /// Returns the cosine of the argument _a. /// /// @param[in] _a Angle in radians. @@ -352,6 +513,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float cos(float _a); + /// Returns the cosine of _a. + /// + /// @param[in] _a Angle in radians. + /// + /// @returns cos(_a). + /// + BX_CONSTEXPR_FUNC double cos(double _a); + /// Returns hyperbolic cosine of the argument _a. /// /// @param[in] _a Input value. @@ -360,6 +529,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float cosh(float _a); + /// Returns the hyperbolic cosine of _a. + /// + /// @param[in] _a Input value. + /// + /// @returns cosh(_a). + /// + BX_CONSTEXPR_FUNC double cosh(double _a); + /// Returns radian angle between 0 and pi whose cosine is _a. /// /// @param[in] _a Cosine value [-1, 1]. @@ -368,6 +545,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float acos(float _a); + /// Returns the arc cosine of _a. + /// + /// @param[in] _a Input value in [-1, 1]. + /// + /// @returns acos(_a) in radians, or NaN if _a is outside [-1, 1]. + /// + BX_CONSTEXPR_FUNC double acos(double _a); + /// Returns the circular tangent of the radian argument _a. /// /// @param[in] _a Angle in radians. @@ -376,6 +561,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float tan(float _a); + /// Returns the tangent of _a. + /// + /// @param[in] _a Angle in radians. + /// + /// @returns tan(_a). + /// + BX_CONSTEXPR_FUNC double tan(double _a); + /// Returns hyperbolic tangent of the argument _a. /// /// @param[in] _a Input value. @@ -384,6 +577,80 @@ namespace bx /// BX_CONSTEXPR_FUNC float tanh(float _a); + /// Returns the hyperbolic tangent of _a. + /// + /// @param[in] _a Input value. + /// + /// @returns tanh(_a). + /// + BX_CONSTEXPR_FUNC double tanh(double _a); + + /// Returns the natural logarithm of 1 + _a, without the precision loss that + /// log(1.0f + _a) has for a small _a. + /// + /// @param[in] _a Input value (must be >= -1). + /// + /// @returns ln(1 + _a). + /// + BX_CONSTEXPR_FUNC float log1p(float _a); + + /// Returns the natural logarithm of 1 + _a, without the precision loss that + /// log(1.0 + _a) has for a small _a. + /// + /// @param[in] _a Input value (must be >= -1). + /// + /// @returns ln(1 + _a). + /// + BX_CONSTEXPR_FUNC double log1p(double _a); + + /// Returns the inverse hyperbolic sine of _a. + /// + /// @param[in] _a Input value. + /// + /// @returns asinh(_a). + /// + BX_CONSTEXPR_FUNC float asinh(float _a); + + /// Returns the inverse hyperbolic sine of _a. + /// + /// @param[in] _a Input value. + /// + /// @returns asinh(_a). + /// + BX_CONSTEXPR_FUNC double asinh(double _a); + + /// Returns the inverse hyperbolic cosine of _a. + /// + /// @param[in] _a Input value (must be >= 1). + /// + /// @returns acosh(_a), or NaN if _a < 1. + /// + BX_CONSTEXPR_FUNC float acosh(float _a); + + /// Returns the inverse hyperbolic cosine of _a. + /// + /// @param[in] _a Input value (must be >= 1). + /// + /// @returns acosh(_a), or NaN if _a < 1. + /// + BX_CONSTEXPR_FUNC double acosh(double _a); + + /// Returns the inverse hyperbolic tangent of _a. + /// + /// @param[in] _a Input value in [-1, 1]. + /// + /// @returns atanh(_a), or NaN if |_a| > 1. + /// + BX_CONSTEXPR_FUNC float atanh(float _a); + + /// Returns the inverse hyperbolic tangent of _a. + /// + /// @param[in] _a Input value in [-1, 1]. + /// + /// @returns atanh(_a), or NaN if |_a| > 1. + /// + BX_CONSTEXPR_FUNC double atanh(double _a); + /// Returns radian angle between -pi/2 and +pi/2 whose tangent is _a. /// /// @param[in] _a Tangent value. @@ -392,6 +659,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float atan(float _a); + /// Returns the arc tangent of _a. + /// + /// @param[in] _a Input value. + /// + /// @returns atan(_a) in radians. + /// + BX_CONSTEXPR_FUNC double atan(double _a); + /// Returns the inverse tangent of _y/_x. /// /// @param[in] _y Y component. @@ -401,6 +676,15 @@ namespace bx /// BX_CONSTEXPR_FUNC float atan2(float _y, float _x); + /// Returns the arc tangent of _y / _x, using their signs to pick the quadrant. + /// + /// @param[in] _y Numerator. + /// @param[in] _x Denominator. + /// + /// @returns atan2(_y, _x) in radians. + /// + BX_CONSTEXPR_FUNC double atan2(double _y, double _x); + /// Computes _a raised to the _b power. /// /// @param[in] _a Base. @@ -410,6 +694,16 @@ namespace bx /// BX_CONSTEXPR_FUNC float pow(float _a, float _b); + /// Returns _a raised to the _b power. If _a is negative _b must be integral, + /// otherwise the result is NaN. + /// + /// @param[in] _a Base. + /// @param[in] _b Exponent. + /// + /// @returns _a^_b. + /// + BX_CONSTEXPR_FUNC double pow(double _a, double _b); + /// Returns the result of multiplying _a by 2 raised to the power of the exponent `_a * (2^_b)`. /// /// @param[in] _a Significand. @@ -419,6 +713,15 @@ namespace bx /// BX_CONSTEXPR_FUNC float ldexp(float _a, int32_t _b); + /// Returns _a multiplied by 2 raised to the power _b. + /// + /// @param[in] _a Significand. + /// @param[in] _b Exponent. + /// + /// @returns _a * 2^_b. + /// + BX_CONSTEXPR_FUNC double ldexp(double _a, int32_t _b); + /// Returns decomposed given floating point value _a into a normalized fraction and /// an integral power of two. /// @@ -427,7 +730,17 @@ namespace bx /// /// @returns Normalized fraction in [0.5, 1.0). /// - float frexp(float _a, int32_t* _outExp); + constexpr float frexp(float _a, int32_t* _outExp); + + /// Returns decomposed given floating point value _a into a normalized fraction and + /// an integral power of two. + /// + /// @param[in] _a Input value. + /// @param[out] _outExp Output exponent. + /// + /// @returns Normalized fraction in [0.5, 1.0). + /// + constexpr double frexp(double _a, int32_t* _outExp); /// Returns e (2.71828...) raised to the _a power. /// @@ -437,6 +750,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float exp(float _a); + /// Returns e (2.71828...) raised to the _a power. + /// + /// @param[in] _a Exponent. + /// + /// @returns e^_a. + /// + BX_CONSTEXPR_FUNC double exp(double _a); + /// Returns 2 raised to the _a power. /// /// @param[in] _a Exponent. @@ -445,6 +766,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float exp2(float _a); + /// Returns 2 raised to the _a power. + /// + /// @param[in] _a Exponent. + /// + /// @returns 2^_a. + /// + BX_CONSTEXPR_FUNC double exp2(double _a); + /// Returns the base e (2.71828...) logarithm of _a. /// /// @param[in] _a Input value (must be > 0). @@ -453,6 +782,14 @@ namespace bx /// BX_CONSTEXPR_FUNC float log(float _a); + /// Returns the base e (2.71828...) logarithm of _a. + /// + /// @param[in] _a Input value (must be > 0). + /// + /// @returns ln(_a). + /// + BX_CONSTEXPR_FUNC double log(double _a); + /// Returns the base 2 logarithm of _a. /// /// @param[in] _a Input value (must be > 0). @@ -461,6 +798,220 @@ namespace bx /// BX_CONSTEXPR_FUNC float log2(float _a); + /// Returns the base 2 logarithm of _a. + /// + /// @param[in] _a Input value (must be > 0). + /// + /// @returns log2(_a). + /// + BX_CONSTEXPR_FUNC double log2(double _a); + + /// Returns the base 10 logarithm of _a. + /// + /// @param[in] _a Input value (must be > 0). + /// + /// @returns log10(_a). + /// + BX_CONSTEXPR_FUNC float log10(float _a); + + /// Returns the base 10 logarithm of _a. + /// + /// @param[in] _a Input value (must be > 0). + /// + /// @returns log10(_a). + /// + BX_CONSTEXPR_FUNC double log10(double _a); + + /// Returns the nearest integer not greater in magnitude than _a. + /// + /// @param[in] _a Input value. + /// + /// @returns Truncated value. + /// + BX_CONSTEXPR_FUNC float trunc(float _a); + + /// Returns _a truncated towards zero. + /// + /// @param[in] _a Input value. + /// + /// @returns Integer part of _a. + /// + BX_CONSTEXPR_FUNC double trunc(double _a); + + /// Returns the fractional (or decimal) part of _a, which is greater than or equal to 0 + /// and less than 1. + /// + /// @param[in] _a Input value. + /// + /// @returns Fractional part of _a. + /// + BX_CONSTEXPR_FUNC float fract(float _a); + + /// Returns the fractional part of _a, with the sign of _a. + /// + /// @param[in] _a Input value. + /// + /// @returns _a - trunc(_a). + /// + BX_CONSTEXPR_FUNC double fract(double _a); + + /// Returns the absolute of integer _a. + /// + /// @param[in] _a Input value. + /// + /// @returns Absolute value of _a. + /// + template + requires (isInteger() ) + constexpr Ty abs(Ty _a); + + /// Returns true if the sign bit of _a is set. + /// + /// @param[in] _a Input value. + /// + /// @returns True if _a is negative, including negative zero. + /// + BX_CONSTEXPR_FUNC bool signBit(double _a); + + /// Returns true if _a and _b are equal within _epsilon tolerance. + /// + /// @param[in] _a First value. + /// @param[in] _b Second value. + /// @param[in] _epsilon Tolerance. + /// + /// @returns True if _a and _b are within _epsilon of each other. + /// + BX_CONSTEXPR_FUNC bool isEqual(double _a, double _b, double _epsilon); + + /// Returns the fractional part of _a, and stores the integral part in _outIntegral. + /// + /// @param[in] _a Input value. + /// @param[out] _outIntegral Output integral part, truncated towards zero. + /// + /// @returns Fractional part of _a, with the sign of _a. + /// + constexpr float modf(float _a, float* _outIntegral); + + /// Returns the fractional part of _a, and stores the integral part in _outIntegral. + /// + /// @param[in] _a Input value. + /// @param[out] _outIntegral Output integral part, truncated towards zero. + /// + /// @returns Fractional part of _a, with the sign of _a. + /// + constexpr double modf(double _a, double* _outIntegral); + + /// Returns the remainder of _a / _b, truncated towards zero. The result has the + /// sign of _a. Unlike mod, which truncates the quotient towards negative infinity. + /// + /// @param[in] _a Numerator. + /// @param[in] _b Denominator. + /// + /// @returns _a - _b * trunc(_a / _b). + /// + BX_CONSTEXPR_FUNC float fmod(float _a, float _b); + + /// Returns the remainder of _a / _b, truncated towards zero. The result has the + /// sign of _a. Unlike mod, which truncates the quotient towards negative infinity. + /// + /// @param[in] _a Numerator. + /// @param[in] _b Denominator. + /// + /// @returns _a - _b * trunc(_a / _b). + /// + BX_CONSTEXPR_FUNC double fmod(double _a, double _b); + + /// Returns the remainder of _a / _b, with the sign of _b. + /// + /// @param[in] _a Numerator. + /// @param[in] _b Denominator. + /// + /// @returns _a - _b * floor(_a / _b). + /// + BX_CONSTEXPR_FUNC float mod(float _a, float _b); + + /// Returns the remainder of _a / _b, with the sign of _b. + /// + /// @param[in] _a Numerator. + /// @param[in] _b Denominator. + /// + /// @returns _a - _b * floor(_a / _b). + /// + BX_CONSTEXPR_FUNC double mod(double _a, double _b); + + /// Returns the square root of _a. + /// + /// @param[in] _a Input value (must be >= 0). + /// + /// @returns Square root of _a. + /// + BX_CONSTEXPR_FUNC float sqrt(float _a); + + /// Returns the square root of _a. + /// + /// @param[in] _a Input value (must be >= 0). + /// + /// @returns sqrt(_a). + /// + BX_CONSTEXPR_FUNC double sqrt(double _a); + + /// Returns reciprocal square root of _a. + /// + /// @param[in] _a Input value (must be > 0). + /// + /// @returns 1 / sqrt(_a). + /// + BX_CONSTEXPR_FUNC float rsqrt(float _a); + + /// Returns the reciprocal square root of _a. + /// + /// @param[in] _a Input value (must be > 0). + /// + /// @returns 1/sqrt(_a). + /// + BX_CONSTEXPR_FUNC double rsqrt(double _a); + + /// Returns result of multiply and add (_a * _b + _c). + /// + /// @param[in] _a First factor. + /// @param[in] _b Second factor. + /// @param[in] _c Addend. + /// + /// @returns _a * _b + _c. + /// + BX_CONSTEXPR_FUNC float mad(float _a, float _b, float _c); + + /// Returns the result of multiply and add. (_a * _b + _c) + /// + /// @param[in] _a First factor. + /// @param[in] _b Second factor. + /// @param[in] _c Addend. + /// + /// @returns _a * _b + _c. + /// + BX_CONSTEXPR_FUNC double mad(double _a, double _b, double _c); + + /// Returns result of negated multiply-sub operation -(_a * _b - _c) -> _c - _a * _b, + /// rounded like mad. + /// + /// @param[in] _a First factor. + /// @param[in] _b Second factor. + /// @param[in] _c Value to subtract from. + /// + /// @returns _c - _a * _b. + /// + BX_CONSTEXPR_FUNC float nms(float _a, float _b, float _c); + + /// Returns the result of negated multiply and add. (_c - _a * _b) + /// + /// @param[in] _a First factor. + /// @param[in] _b Second factor. + /// @param[in] _c Addend. + /// + /// @returns _c - _a * _b. + /// + BX_CONSTEXPR_FUNC double nms(double _a, double _b, double _c); + /// Count number of bits set. /// /// @param[in] _val Input value. @@ -533,50 +1084,6 @@ namespace bx template BX_CONSTEXPR_FUNC Ty nextPow2(Ty _a); - /// Returns the square root of _a. - /// - /// @param[in] _a Input value (must be >= 0). - /// - /// @returns Square root of _a. - /// - BX_CONSTEXPR_FUNC float sqrt(float _a); - - /// Returns reciprocal square root of _a. - /// - /// @param[in] _a Input value (must be > 0). - /// - /// @returns 1 / sqrt(_a). - /// - BX_CONSTEXPR_FUNC float rsqrt(float _a); - - /// Returns the nearest integer not greater in magnitude than _a. - /// - /// @param[in] _a Input value. - /// - /// @returns Truncated value. - /// - BX_CONSTEXPR_FUNC float trunc(float _a); - - /// Returns the fractional (or decimal) part of _a, which is greater than or equal to 0 - /// and less than 1. - /// - /// @param[in] _a Input value. - /// - /// @returns Fractional part of _a. - /// - BX_CONSTEXPR_FUNC float fract(float _a); - - /// Returns result of negated multiply-sub operation -(_a * _b - _c) -> _c - _a * _b, - /// rounded like mad. - /// - /// @param[in] _a First factor. - /// @param[in] _b Second factor. - /// @param[in] _c Value to subtract from. - /// - /// @returns _c - _a * _b. - /// - BX_CONSTEXPR_FUNC float nms(float _a, float _b, float _c); - /// Returns result of addition (_a + _b). /// /// @param[in] _a First value. @@ -628,16 +1135,6 @@ namespace bx /// BX_CONSTEXPR_FUNC float mul(float _a, float _b); - /// Returns result of multiply and add (_a * _b + _c). - /// - /// @param[in] _a First factor. - /// @param[in] _b Second factor. - /// @param[in] _c Addend. - /// - /// @returns _a * _b + _c. - /// - BX_CONSTEXPR_FUNC float mad(float _a, float _b, float _c); - /// Returns reciprocal of _a. /// /// @param[in] _a Input value. @@ -672,15 +1169,6 @@ namespace bx /// BX_CONSTEXPR_FUNC float divSafe(float _a, float _b); - /// Returns the floating-point remainder of the division operation _a/_b. - /// - /// @param[in] _a Numerator. - /// @param[in] _b Denominator. - /// - /// @returns Remainder of _a / _b. - /// - BX_CONSTEXPR_FUNC float mod(float _a, float _b); - /// Returns true if _a and _b are equal within _epsilon tolerance. /// /// @param[in] _a First value. diff --git a/src/crtnone.cpp b/src/crtnone.cpp index cd67c9e..b783634 100644 --- a/src/crtnone.cpp +++ b/src/crtnone.cpp @@ -151,7 +151,7 @@ extern "C" double fabs(double _x) extern "C" double ldexp(double _x, int _exp) { - return bx::ldexp(float(_x), _exp); + return bx::ldexp(_x, _exp); } extern "C" float expf(float _x) @@ -159,15 +159,29 @@ extern "C" float expf(float _x) return bx::exp(_x); } +extern "C" double exp(double _x) +{ + return bx::exp(_x); +} + extern "C" float logf(float _x) { return bx::log(_x); } +extern "C" double log(double _x) +{ + return bx::log(_x); +} + extern "C" float log10f(float _x) { - BX_UNUSED(_x); - return 0.0f; + return bx::log10(_x); +} + +extern "C" double log10(double _x) +{ + return bx::log10(_x); } extern "C" float powf(float _x, float _y) @@ -175,7 +189,7 @@ extern "C" float powf(float _x, float _y) return bx::pow(_x, _y); } -extern "C" double pow(double _x, float _y) +extern "C" double pow(double _x, double _y) { return bx::pow(_x, _y); } @@ -185,21 +199,41 @@ extern "C" float sinf(float _x) return bx::sin(_x); } +extern "C" double sin(double _x) +{ + return bx::sin(_x); +} + extern "C" float cosf(float _x) { return bx::cos(_x); } +extern "C" double cos(double _x) +{ + return bx::cos(_x); +} + extern "C" float tanf(float _x) { return bx::tan(_x); } +extern "C" double tan(double _x) +{ + return bx::tan(_x); +} + extern "C" float atan2f(float _y, float _x) { return bx::atan2(_y, _x); } +extern "C" double atan2(double _y, double _x) +{ + return bx::atan2(_y, _x); +} + extern "C" float sqrtf(float _x) { return bx::sqrt(_x); @@ -235,9 +269,19 @@ extern "C" float acosf(float _x) return bx::acos(_x); } +extern "C" double acos(double _x) +{ + return bx::acos(_x); +} + extern "C" float fmodf(float _numer, float _denom) { - return bx::mod(_numer, _denom); + return bx::fmod(_numer, _denom); +} + +extern "C" double fmod(double _numer, double _denom) +{ + return bx::fmod(_numer, _denom); } extern "C" int atoi(const char* _str) diff --git a/src/dtoa.cpp b/src/dtoa.cpp index 885c488..aedb33c 100644 --- a/src/dtoa.cpp +++ b/src/dtoa.cpp @@ -562,479 +562,244 @@ namespace bx return toStringUnsigned(_dst, _max, _value, _base, _separator); } - /* - * https://github.com/grzegorz-kraszewski/stringtofloat/ - * - * MIT License - * - * Copyright (c) 2016 Grzegorz Kraszewski - * - * Permission is hereby granted, free of charge, to any person obtaining a copy - * of this software and associated documentation files (the "Software"), to deal - * in the Software without restriction, including without limitation the rights - * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - * copies of the Software, and to permit persons to whom the Software is - * furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included in all - * copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - * SOFTWARE. - */ + constexpr uint32_t kFromStringMaxDigits = 768; - /* - * IMPORTANT - * - * The code works in "round towards zero" mode. This is different from - * GCC standard library strtod(), which uses "round half to even" rule. - * Therefore it cannot be used as a direct drop-in replacement, as in - * some cases results will be different on the least significant bit of - * mantissa. Read more in the README.md file. - */ + constexpr uint32_t kFromStringNumLimbs = 96; -#define DIGITS 18 - -#define DOUBLE_PLUS_ZERO UINT64_C(0x0000000000000000) -#define DOUBLE_MINUS_ZERO UINT64_C(0x8000000000000000) -#define DOUBLE_PLUS_INFINITY UINT64_C(0x7ff0000000000000) -#define DOUBLE_MINUS_INFINITY UINT64_C(0xfff0000000000000) - -#define lsr96(s2, s1, s0, d2, d1, d0) \ - d0 = ( (s0) >> 1) | ( ( (s1) & 1) << 31); \ - d1 = ( (s1) >> 1) | ( ( (s2) & 1) << 31); \ - d2 = (s2) >> 1; - -#define lsl96(s2, s1, s0, d2, d1, d0) \ - d2 = ( (s2) << 1) | ( ( (s1) & (1 << 31) ) >> 31); \ - d1 = ( (s1) << 1) | ( ( (s0) & (1 << 31) ) >> 31); \ - d0 = (s0) << 1; - - /* - * Undefine the below constant if your processor or compiler is slow - * at 64-bit arithmetic. This is a rare case however. 64-bit macros are - * better for deeply pipelined CPUs (no conditional execution), are - * very efficient for 64-bit processors and also fast on 32-bit processors - * featuring extended precision arithmetic (x86, PowerPC_32, M68k and probably - * more). - */ - -#define USE_64BIT_FOR_ADDSUB_MACROS 0 - -#if USE_64BIT_FOR_ADDSUB_MACROS - -#define add96(s2, s1, s0, d2, d1, d0) { \ - uint64_t w; \ - w = (uint64_t)(s0) + (uint64_t)(d0); \ - (s0) = w; \ - w >>= 32; \ - w += (uint64_t)(s1) + (uint64_t)(d1); \ - (s1) = w; \ - w >>= 32; \ - w += (uint64_t)(s2) + (uint64_t)(d2); \ - (s2) = w; } - -#define sub96(s2, s1, s0, d2, d1, d0) { \ - uint64_t w; \ - w = (uint64_t)(s0) - (uint64_t)(d0); \ - (s0) = w; \ - w >>= 32; \ - w += (uint64_t)(s1) - (uint64_t)(d1); \ - (s1) = w; \ - w >>= 32; \ - w += (uint64_t)(s2) - (uint64_t)(d2); \ - (s2) = w; } - -#else - -#define add96(s2, s1, s0, d2, d1, d0) { \ - uint32_t _x, _c; \ - _x = (s0); (s0) += (d0); \ - if ( (s0) < _x) _c = 1; else _c = 0; \ - _x = (s1); (s1) += (d1) + _c; \ - if ( ( (s1) < _x) || ( ( (s1) == _x) && _c) ) _c = 1; else _c = 0; \ - (s2) += (d2) + _c; } - -#define sub96(s2, s1, s0, d2, d1, d0) { \ - uint32_t _x, _c; \ - _x = (s0); (s0) -= (d0); \ - if ( (s0) > _x) _c = 1; else _c = 0; \ - _x = (s1); (s1) -= (d1) + _c; \ - if ( ( (s1) > _x) || ( ( (s1) == _x) && _c) ) _c = 1; else _c = 0; \ - (s2) -= (d2) + _c; } - -#endif /* USE_64BIT_FOR_ADDSUB_MACROS */ - - /* parser state machine states */ - -#define FSM_A 0 -#define FSM_B 1 -#define FSM_C 2 -#define FSM_D 3 -#define FSM_E 4 -#define FSM_F 5 -#define FSM_G 6 -#define FSM_H 7 -#define FSM_I 8 -#define FSM_STOP 9 - - /* The structure is filled by parser, then given to converter. */ - struct PrepNumber + struct FloatBig { - int negative; /* 0 if positive number, 1 if negative */ - int32_t exponent; /* power of 10 exponent */ - uint64_t mantissa; /* integer mantissa */ + uint32_t num; + uint64_t limb[kFromStringNumLimbs]; }; - /* Possible parser return values. */ - -#define PARSER_OK 0 // parser finished OK -#define PARSER_PZERO 1 // no digits or number is smaller than +-2^-1022 -#define PARSER_MZERO 2 // number is negative, module smaller -#define PARSER_PINF 3 // number is higher than +HUGE_VAL -#define PARSER_MINF 4 // number is lower than -HUGE_VAL - - inline char next(const char*& _s, const char* _term) + static bool floatBigMulAdd(FloatBig& _big, uint32_t _mul, uint32_t _add) { - return _s != _term - ? *_s++ - : '\0' - ; + uint64_t carry = _add; + + for (uint32_t ii = 0; ii < _big.num; ++ii) + { + const uint64_t lo = (_big.limb[ii] & UINT32_MAX) * _mul + (carry & UINT32_MAX); + const uint64_t hi = (_big.limb[ii] >> 32) * _mul + (carry >> 32) + (lo >> 32); + + _big.limb[ii] = (hi << 32) | (lo & UINT32_MAX); + carry = hi >> 32; + } + + if (0 != carry) + { + if (kFromStringNumLimbs <= _big.num) + { + return false; + } + + _big.limb[_big.num] = carry; + ++_big.num; + } + + return true; } - static int parser(const char* _s, const char* _term, PrepNumber* _pn) + static bool floatBigShiftLeft(FloatBig& _big, uint32_t _bits) { - int state = FSM_A; - int digx = 0; - char c = ' '; /* initial value for kicking off the state machine */ - int result = PARSER_OK; - int expneg = 0; - int32_t expexp = 0; - - while (state != FSM_STOP) // && _s != _term) + if (0 == _big.num) { - switch (state) + return true; + } + + const uint32_t words = _bits / 64; + const uint32_t bits = _bits % 64; + const uint32_t extra = 0 != bits ? 1 : 0; + + if (kFromStringNumLimbs < _big.num + words + extra) + { + return false; + } + + if (0 != bits) + { + uint64_t carry = 0; + + for (uint32_t ii = 0; ii < _big.num; ++ii) { - case FSM_A: - if (isSpace(c) ) - { - c = next(_s, _term); - } - else - { - state = FSM_B; - } - break; + const uint64_t value = _big.limb[ii]; - case FSM_B: - state = FSM_C; + _big.limb[ii] = (value << bits) | carry; + carry = value >> (64 - bits); + } - if (c == '+') - { - c = next(_s, _term); - } - else if (c == '-') - { - _pn->negative = 1; - c = next(_s, _term); - } - else if (isNumeric(c) ) - { - } - else if (c == '.') - { - } - else - { - state = FSM_STOP; - } - break; - - case FSM_C: - if (c == '0') - { - c = next(_s, _term); - } - else if (c == '.') - { - c = next(_s, _term); - state = FSM_D; - } - else - { - state = FSM_E; - } - break; - - case FSM_D: - if (c == '0') - { - c = next(_s, _term); - if (_pn->exponent > -2147483647) _pn->exponent--; - } - else - { - state = FSM_F; - } - break; - - case FSM_E: - if (isNumeric(c) ) - { - if (digx < DIGITS) - { - _pn->mantissa *= 10; - _pn->mantissa += c - '0'; - digx++; - } - else if (_pn->exponent < 2147483647) - { - _pn->exponent++; - } - - c = next(_s, _term); - } - else if (c == '.') - { - c = next(_s, _term); - state = FSM_F; - } - else - { - state = FSM_F; - } - break; - - case FSM_F: - if (isNumeric(c) ) - { - if (digx < DIGITS) - { - _pn->mantissa *= 10; - _pn->mantissa += c - '0'; - _pn->exponent--; - digx++; - } - - c = next(_s, _term); - } - else if ('e' == toLower(c) ) - { - c = next(_s, _term); - state = FSM_G; - } - else - { - state = FSM_G; - } - break; - - case FSM_G: - if (c == '+') - { - c = next(_s, _term); - } - else if (c == '-') - { - expneg = 1; - c = next(_s, _term); - } - - state = FSM_H; - break; - - case FSM_H: - if (c == '0') - { - c = next(_s, _term); - } - else - { - state = FSM_I; - } - break; - - case FSM_I: - if (isNumeric(c) ) - { - if (expexp < 214748364) - { - expexp *= 10; - expexp += c - '0'; - } - - c = next(_s, _term); - } - else - { - state = FSM_STOP; - } - break; + if (0 != carry) + { + _big.limb[_big.num] = carry; + ++_big.num; } } - if (expneg) + if (0 != words) { - expexp = -expexp; + for (uint32_t ii = _big.num; ii > 0; --ii) + { + _big.limb[ii - 1 + words] = _big.limb[ii - 1]; + } + + for (uint32_t ii = 0; ii < words; ++ii) + { + _big.limb[ii] = 0; + } + + _big.num += words; } - _pn->exponent += expexp; + return true; + } - if (_pn->mantissa == 0) + static uint32_t floatBigDivSmall(FloatBig& _big, uint32_t _div) + { + uint64_t rem = 0; + + for (uint32_t ii = _big.num; ii > 0; --ii) { - if (_pn->negative) + const uint64_t value = _big.limb[ii - 1]; + + const uint64_t hiNum = (rem << 32) | (value >> 32); + const uint64_t hiQ = hiNum / _div; + const uint64_t hiR = hiNum % _div; + + const uint64_t loNum = (hiR << 32) | (value & UINT32_MAX); + const uint64_t loQ = loNum / _div; + + rem = loNum % _div; + + _big.limb[ii - 1] = (hiQ << 32) | loQ; + } + + while (0 != _big.num + && 0 == _big.limb[_big.num - 1]) + { + --_big.num; + } + + return uint32_t(rem); + } + + static uint32_t floatBigBitLength(const FloatBig& _big) + { + if (0 == _big.num) + { + return 0; + } + + return (_big.num - 1)*64 + 64 - countLeadingZeros(_big.limb[_big.num - 1]); + } + + static bool floatBigBit(const FloatBig& _big, uint32_t _bit) + { + const uint32_t word = _bit / 64; + + if (_big.num <= word) + { + return false; + } + + return 0 != ( (_big.limb[word] >> (_bit % 64) ) & 1); + } + + static bool floatBigAnyBitBelow(const FloatBig& _big, uint32_t _bit) + { + const uint32_t word = _bit / 64; + const uint32_t bits = _bit % 64; + + for (uint32_t ii = 0; ii < word && ii < _big.num; ++ii) + { + if (0 != _big.limb[ii]) { - result = PARSER_MZERO; - } - else - { - result = PARSER_PZERO; + return true; } } - else if (_pn->exponent > 309) + + if (word < _big.num + && 0 != bits) { - if (_pn->negative) + const uint64_t mask = (UINT64_C(1) << bits) - 1; + + if (0 != (_big.limb[word] & mask) ) { - result = PARSER_MINF; - } - else - { - result = PARSER_PINF; + return true; } } - else if (_pn->exponent < -328) + + return false; + } + + static uint64_t floatBigBits(const FloatBig& _big, uint32_t _from, uint32_t _count) + { + uint64_t result = 0; + + for (uint32_t ii = 0; ii < _count; ++ii) { - if (_pn->negative) + if (floatBigBit(_big, _from + ii) ) { - result = PARSER_MZERO; - } - else - { - result = PARSER_PZERO; + result |= UINT64_C(1) << ii; } } return result; } - static double converter(PrepNumber* _pn) + static double floatBigRound(const FloatBig& _big, int32_t _binExp, bool _sticky, bool _negative, int32_t _mantBits, int32_t _minExp) { - int binexp = 92; - uint32_t s2, s1, s0; /* 96-bit precision integer */ - uint32_t q2, q1, q0; /* 96-bit precision integer */ - uint32_t r2, r1, r0; /* 96-bit precision integer */ - uint32_t mask28 = UINT32_C(0xf) << 28; + const uint32_t length = floatBigBitLength(_big); - uint64_t hdu = 0; - - s0 = (uint32_t)(_pn->mantissa & UINT32_MAX); - s1 = (uint32_t)(_pn->mantissa >> 32); - s2 = 0; - - while (_pn->exponent > 0) + if (0 == length) { - lsl96(s2, s1, s0, q2, q1, q0); // q = p << 1 - lsl96(q2, q1, q0, r2, r1, r0); // r = p << 2 - lsl96(r2, r1, r0, s2, s1, s0); // p = p << 3 - add96(s2, s1, s0, q2, q1, q0); // p = (p << 3) + (p << 1) - - _pn->exponent--; - - while (s2 & mask28) - { - lsr96(s2, s1, s0, q2, q1, q0); - binexp++; - s2 = q2; - s1 = q1; - s0 = q0; - } + return _negative ? -0.0 : 0.0; } - while (_pn->exponent < 0) + const int32_t unbiased = int32_t(length) - 1 + _binExp; + const int32_t smallest = _minExp - _mantBits + 1; + + int32_t keep = _mantBits; + + if (_minExp > unbiased) { - while (!(s2 & (1 << 31) ) ) - { - lsl96(s2, s1, s0, q2, q1, q0); - binexp--; - s2 = q2; - s1 = q1; - s0 = q0; - } - - q2 = s2 / 10; - r1 = s2 % 10; - r2 = (s1 >> 8) | (r1 << 24); - q1 = r2 / 10; - r1 = r2 % 10; - r2 = ( (s1 & 0xFF) << 16) | (s0 >> 16) | (r1 << 24); - r0 = r2 / 10; - r1 = r2 % 10; - q1 = (q1 << 8) | ( (r0 & 0x00FF0000) >> 16); - q0 = r0 << 16; - r2 = (s0 & UINT16_MAX) | (r1 << 16); - q0 |= r2 / 10; - s2 = q2; - s1 = q1; - s0 = q0; - - _pn->exponent++; + keep = _mantBits + (unbiased - _minExp); } - if (s2 || s1 || s0) + if (0 >= keep) { - while (!(s2 & mask28) ) - { - lsl96(s2, s1, s0, q2, q1, q0); - binexp--; - s2 = q2; - s1 = q1; - s0 = q0; - } + const bool half = 0 == keep && floatBigBit(_big, length - 1); + const bool more = _sticky || floatBigAnyBitBelow(_big, length - 1); + const double tiny = (half && more) ? ldexp(1.0, smallest) : 0.0; + + return _negative ? -tiny : tiny; } - binexp += 1023; + const int32_t shift = int32_t(length) - keep; - if (binexp > 2046) + uint64_t mantissa = 0; + bool roundBit = false; + bool tail = _sticky; + + if (0 < shift) { - if (_pn->negative) - { - hdu = DOUBLE_MINUS_INFINITY; - } - else - { - hdu = DOUBLE_PLUS_INFINITY; - } + mantissa = floatBigBits(_big, uint32_t(shift), uint32_t(keep) ); + roundBit = floatBigBit(_big, uint32_t(shift) - 1); + tail = tail || (1 < shift && floatBigAnyBitBelow(_big, uint32_t(shift) - 1) ); } - else if (binexp < 1) + else { - if (_pn->negative) - { - hdu = DOUBLE_MINUS_ZERO; - } - } - else if (s2) - { - uint64_t q; - uint64_t binexs2 = (uint64_t)binexp; - - binexs2 <<= 52; - q = ( (uint64_t)(s2 & ~mask28) << 24) - | ( ( (uint64_t)s1 + 128) >> 8) | binexs2; - - if (_pn->negative) - { - q |= (1ULL << 63); - } - - hdu = q; + mantissa = floatBigBits(_big, 0, length); + mantissa <<= uint32_t(-shift); } - return bitCast(hdu); + if (roundBit + && (tail || 0 != (mantissa & 1) ) ) + { + ++mantissa; + } + + const double scaled = ldexp(double(mantissa), _binExp + shift); + + return _negative ? -scaled : scaled; } int32_t toString(char* _out, int32_t _max, bool _value) @@ -1051,14 +816,6 @@ namespace bx return 0 != _str.getLength(); } - bool fromString(float* _out, const StringView& _str) - { - double dbl; - bool result = fromString(&dbl, _str); - *_out = float(dbl); - return result; - } - static bool toHexDigit(uint32_t* _out, char _ch) { if (_ch >= '0' && _ch <= '9') @@ -1082,130 +839,256 @@ namespace bx return false; } - static bool fromStringHex(double* _out, const StringView& _str) + static bool fromStringFloating(double* _out, const StringView& _str, int32_t _mantBits, int32_t _minExp) { - const char* ptr = _str.getPtr(); - const char* end = _str.getTerm(); + const char* ptr = _str.getPtr(); + const char* term = _str.getTerm(); + + while (ptr != term + && isSpace(*ptr) ) + { + ++ptr; + } bool negative = false; - if (ptr < end - && ('-' == *ptr || '+' == *ptr) ) + if (ptr != term + && ('+' == *ptr || '-' == *ptr) ) { - negative = ('-' == *ptr); + negative = '-' == *ptr; ++ptr; } - if ( (end - ptr) < 2 - || '0' != ptr[0] - || ('x' != ptr[1] && 'X' != ptr[1]) ) + const bool hex = (term - ptr) >= 2 + && '0' == ptr[0] + && ('x' == ptr[1] || 'X' == ptr[1]) + ; + + if (hex) { - return false; + ptr += 2; } - ptr += 2; + const uint32_t radix = hex ? 16 : 10; + const int32_t bitsEach = hex ? 4 : 1; - double mantissa = 0.0; - bool anyDigit = false; - uint32_t digit; + FloatBig big; + big.num = 0; + big.limb[0] = 0; - for (; ptr < end && toHexDigit(&digit, *ptr); ++ptr) + uint32_t numDigits = 0; + int32_t pointExp = 0; + int32_t dropped = 0; + bool sticky = false; + bool anyDigit = false; + bool fraction = false; + + for (; ptr != term; ++ptr) { - mantissa = mantissa * 16.0 + double(digit); - anyDigit = true; - } + const char ch = *ptr; - if (ptr < end && '.' == *ptr) - { - ++ptr; - double scale = 1.0 / 16.0; - for (; ptr < end && toHexDigit(&digit, *ptr); ++ptr) + if ('.' == ch) { - mantissa += double(digit) * scale; - scale *= 1.0 / 16.0; - anyDigit = true; + if (fraction) + { + break; + } + + fraction = true; + continue; + } + + uint32_t digit = 0; + + if (hex) + { + if (!toHexDigit(&digit, ch) ) + { + break; + } + } + else + { + if (!isNumeric(ch) ) + { + break; + } + + digit = uint32_t(ch - '0'); + } + + anyDigit = true; + + if (0 == numDigits + && 0 == digit) + { + if (fraction) + { + pointExp -= bitsEach; + } + + continue; + } + + if (kFromStringMaxDigits > numDigits) + { + if (!floatBigMulAdd(big, radix, digit) ) + { + return false; + } + + ++numDigits; + + if (fraction) + { + pointExp -= bitsEach; + } + } + else + { + sticky = sticky || (0 != digit); + + if (!fraction) + { + dropped += bitsEach; + } } } if (!anyDigit) { - return false; // "0x" with no hex digits + return false; } - int32_t exponent = 0; - if (ptr < end - && ('p' == *ptr || 'P' == *ptr) ) + int32_t expPart = 0; + + const bool expMark = ptr != term + && ( hex + ? ('p' == *ptr || 'P' == *ptr) + : ('e' == *ptr || 'E' == *ptr) + ) + ; + + if (expMark) { + const char* save = ptr; + ++ptr; - int32_t expSign = 1; - if (ptr < end && ('+' == *ptr || '-' == *ptr) ) + bool expNeg = false; + + if (ptr != term + && ('+' == *ptr || '-' == *ptr) ) { - expSign = ('-' == *ptr) ? -1 : 1; ++ptr; + expNeg = '-' == *ptr; + ++ptr; } - int32_t magnitude = 0; - for (; ptr < end && *ptr >= '0' && *ptr <= '9'; ++ptr) + if (ptr == term + || !isNumeric(*ptr) ) { - magnitude = magnitude * 10 + int32_t(*ptr - '0'); + ptr = save; + } + else + { + for (; ptr != term && isNumeric(*ptr); ++ptr) + { + if (100000 > expPart) + { + expPart = expPart*10 + (*ptr - '0'); + } + } + + if (expNeg) + { + expPart = -expPart; + } + } + } + + if (0 == big.num) + { + *_out = negative ? -0.0 : 0.0; + return true; + } + + if (hex) + { + *_out = floatBigRound(big, expPart + pointExp + dropped, sticky, negative, _mantBits, _minExp); + return true; + } + + const int32_t decExp = expPart + pointExp + dropped; + + const int32_t magnitude = decExp + int32_t(numDigits); + const int32_t maxBinExp = 1 - _minExp; + const int32_t minBinExp = _minExp - _mantBits + 1; + + if (magnitude > (maxBinExp + 1)*10/33 + 2) + { + *_out = negative ? -kDoubleInfinity : kDoubleInfinity; + return true; + } + + if (magnitude < minBinExp*10/33 - 2) + { + *_out = negative ? -0.0 : 0.0; + return true; + } + + int32_t binExp = decExp; + + if (0 <= decExp) + { + for (int32_t ii = 0; ii < decExp; ++ii) + { + if (!floatBigMulAdd(big, 5, 0) ) + { + *_out = negative ? -kDoubleInfinity : kDoubleInfinity; + return true; + } + } + } + else + { + const int32_t count = -decExp; + + const uint32_t prescale = uint32_t(count)*7/3 + 70; + + if (!floatBigShiftLeft(big, prescale) ) + { + return false; } - exponent = expSign * magnitude; + binExp -= int32_t(prescale); + + for (int32_t ii = 0; ii < count; ++ii) + { + const uint32_t rem = floatBigDivSmall(big, 5); + + sticky = sticky || (0 != rem); + } } - double value = mantissa; - - for (int32_t ii = 0; ii < exponent; ++ii) - { - value *= 2.0; - } - - for (int32_t ii = 0; ii > exponent; --ii) - { - value *= 0.5; - } - - *_out = negative ? -value : value; + *_out = floatBigRound(big, binExp, sticky, negative, _mantBits, _minExp); return true; } + bool fromString(float* _out, const StringView& _str) + { + double value; + + const bool result = fromStringFloating(&value, _str, 24, -126); + + *_out = float(value); + + return result; + } + bool fromString(double* _out, const StringView& _str) { - if (fromStringHex(_out, _str) ) - { - return true; - } - - PrepNumber pn; - pn.mantissa = 0; - pn.negative = 0; - pn.exponent = 0; - - switch (parser(_str.getPtr(), _str.getTerm(), &pn) ) - { - case PARSER_OK: - *_out = converter(&pn); - break; - - case PARSER_PZERO: - *_out = bitCast(DOUBLE_PLUS_ZERO); - break; - - case PARSER_MZERO: - *_out = bitCast(DOUBLE_MINUS_ZERO); - break; - - case PARSER_PINF: - *_out = bitCast(DOUBLE_PLUS_INFINITY); - break; - - case PARSER_MINF: - *_out = bitCast(DOUBLE_MINUS_INFINITY); - break; - } - - return true; + return fromStringFloating(_out, _str, 53, -1022); } bool fromString(long long* _out, const StringView& _str) diff --git a/src/hash.cpp b/src/hash.cpp index bd43e47..1ac1118 100644 --- a/src/hash.cpp +++ b/src/hash.cpp @@ -162,7 +162,7 @@ void mixTail32(Ty& _self, const uint8_t*& _data, int32_t& _len) if (4 == _self.m_count) { - uint32_t kk = *( (uint32_t*)&_self.m_tail[0]); + uint32_t kk = loadUnaligned(&_self.m_tail[0]); _self.mix(kk); _self.m_count = 0; } @@ -217,8 +217,8 @@ void mixTail128(Ty& _self, const uint8_t*& _data, int32_t& _len) if (16 == _self.m_count) { - const uint64_t kk1 = *( (uint64_t*)&_self.m_tail[0]); - const uint64_t kk2 = *( (uint64_t*)&_self.m_tail[8]); + const uint64_t kk1 = loadUnaligned(&_self.m_tail[0]); + const uint64_t kk2 = loadUnaligned(&_self.m_tail[8]); _self.mix(kk1, kk2); _self.m_count = 0; } diff --git a/src/math.cpp b/src/math.cpp index 4c7e742..9825fe9 100644 --- a/src/math.cpp +++ b/src/math.cpp @@ -8,24 +8,6 @@ namespace bx { - float frexp(float _a, int32_t* _outExp) - { - const simd32_t expMask = simd32_splat(kFloatExponentMask); - const simd32_t ftob = simd32_splat(_a); - const simd32_t masked0 = simd_and(ftob, expMask); - const simd32_t exp0 = simd_x32_srl(masked0, kFloatExponentBitShift); - - const simd32_t sMantMask = simd32_splat(kFloatSignMask | kFloatMantissaMask); - const simd32_t masked1 = simd_and(ftob, sMantMask); - const simd32_t half = simd32_splat(0x3f000000u); - const simd32_t bits = simd_or(masked1, half); - const float result = bitsToFloat(bits.u32); - - *_outExp = int32_t(exp0.u32 - 0x7e); - - return result; - } - void mtxLookAt(float* _result, const Vec3& _eye, const Vec3& _at, const Vec3& _up, Handedness::Enum _handedness) { const Vec3 view = normalize( diff --git a/src/string.cpp b/src/string.cpp index 5a174ef..a21ebac 100644 --- a/src/string.cpp +++ b/src/string.cpp @@ -949,12 +949,18 @@ namespace bx char* fracBegin = &str[dot - str + min(prec + _param.spec, 1)]; const int32_t curPrec = int32_t(fracEnd - fracBegin); - const int32_t outPrec = 'g' == _param.fmt ? min(prec, curPrec) : prec; + const int32_t exponentLen = int32_t(strEnd - fracEnd); + const int32_t maxPrec = int32_t(sizeof(str) ) - 1 + - int32_t(fracBegin - str) + - exponentLen + ; + + const int32_t reqPrec = 'g' == _param.fmt ? min(prec, curPrec) : prec; + const int32_t outPrec = min(reqPrec, maxPrec); // Move exponent to its final location after trimming or adding extra 0s. if (fracEnd != strEnd) { - const int32_t exponentLen = int32_t(strEnd - fracEnd); char* finalExponentPtr = &fracBegin[outPrec]; memMove(finalExponentPtr, fracEnd, exponentLen); diff --git a/tests/hash_test.cpp b/tests/hash_test.cpp index 8ed3d37..d37ed4d 100644 --- a/tests/hash_test.cpp +++ b/tests/hash_test.cpp @@ -115,7 +115,7 @@ uint32_t MurmurHash2A(const void * key, int len, uint32_t seed = 0) while(len >= 4) { - uint32_t k = *(uint32_t*)data; + uint32_t k = bx::loadUnaligned(data); mmix(h,k); @@ -206,11 +206,11 @@ uint32_t MurmurHash3_x86_32(const void * key, int len, uint32_t seed) //---------- // body - const uint32_t * blocks = (const uint32_t *)(data + nblocks*4); + const uint8_t * blocks = data + nblocks*4; for(int i = -nblocks; i; i++) { - uint32_t k1 = blocks[i]; + uint32_t k1 = bx::loadUnaligned(blocks + i*4); k1 *= c1; k1 = rotl32(k1,15); diff --git a/tests/math_test.cpp b/tests/math_test.cpp index d132020..847fa74 100644 --- a/tests/math_test.cpp +++ b/tests/math_test.cpp @@ -735,9 +735,9 @@ TEST_CASE("pow", "[math][libm]") STATIC_REQUIRE(0.0f == bx::pow(0.0f, 2.0f) ); STATIC_REQUIRE( 4.0f == bx::pow( 2.0f, 2.0f) ); - STATIC_REQUIRE( -4.0f == bx::pow(-2.0f, 2.0f) ); + STATIC_REQUIRE( 4.0f == bx::pow(-2.0f, 2.0f) ); STATIC_REQUIRE( 0.25f == bx::pow( 2.0f, -2.0f) ); - STATIC_REQUIRE( -0.25f == bx::pow(-2.0f, -2.0f) ); + STATIC_REQUIRE( 0.25f == bx::pow(-2.0f, -2.0f) ); STATIC_REQUIRE( 8.0f == bx::pow( 2.0f, 3.0f) ); STATIC_REQUIRE( -8.0f == bx::pow(-2.0f, 3.0f) ); STATIC_REQUIRE( 0.125f == bx::pow( 2.0f, -3.0f) ); @@ -968,6 +968,622 @@ TEST_CASE("copySign", "[math][libm]") STATIC_REQUIRE(-bx::kFloatInfinity == bx::copySign(bx::kFloatInfinity, -1389) ); } +typedef double (*MathDoubleFn)(double); + +template +static void testMathFunc1Double(double _value) +{ + REQUIRE(CrtT(_value) == BxT(_value) ); +} + +TEST_CASE("floor, ceil, trunc, round double", "[math][libm]") +{ + STATIC_REQUIRE( 13.0 == bx::floor( 13.89) ); + STATIC_REQUIRE(-14.0 == bx::floor(-13.89) ); + STATIC_REQUIRE( 14.0 == bx::ceil( 13.89) ); + STATIC_REQUIRE(-13.0 == bx::ceil( -13.89) ); + STATIC_REQUIRE( 13.0 == bx::trunc( 13.89) ); + STATIC_REQUIRE(-13.0 == bx::trunc(-13.89) ); + + STATIC_REQUIRE( 2.0 == bx::round( 2.5) ); + STATIC_REQUIRE( 4.0 == bx::round( 3.5) ); + STATIC_REQUIRE( -2.0 == bx::round( -2.5) ); + STATIC_REQUIRE( 0.0 == bx::round( 0.49999999999999994) ); + + STATIC_REQUIRE( 1.0e300 == bx::trunc( 1.0e300) ); + STATIC_REQUIRE(-1.0e300 == bx::floor(-1.0e300) ); + STATIC_REQUIRE( 0x1p52 == bx::round(0x1p52) ); + + STATIC_REQUIRE(bx::kDoubleSignMask == bx::bitCast(bx::trunc(-0.5) ) ); + STATIC_REQUIRE(bx::kDoubleSignMask == bx::bitCast(bx::ceil( -0.5) ) ); + STATIC_REQUIRE(bx::kDoubleSignMask == bx::bitCast(bx::round(-0.4) ) ); + STATIC_REQUIRE(bx::kDoubleExponentMask == bx::bitCast(bx::floor(bx::kDoubleInfinity) ) ); + + testMathFunc1Double( 13.89); + testMathFunc1Double(-13.89); + testMathFunc1Double( 13.89); + testMathFunc1Double(-13.89); + testMathFunc1Double( 13.89); + testMathFunc1Double(-13.89); + + bx::RngMwc rng; + + for (uint32_t ii = 0; ii < 100000; ++ii) + { + const uint64_t hi = uint64_t(rng.gen() ); + const uint64_t lo = uint64_t(rng.gen() ); + const uint64_t bits = (hi<<32) | lo; + const double vv = bx::bitsToDouble(bits); + + if (bx::isNan(vv) ) + { + continue; + } + + REQUIRE(::floor(vv) == bx::floor(vv) ); + REQUIRE(::ceil(vv) == bx::ceil(vv) ); + REQUIRE(::trunc(vv) == bx::trunc(vv) ); + REQUIRE(::nearbyint(vv) == bx::round(vv) ); + } +} + +TEST_CASE("abs, sign, copySign double", "[math][libm]") +{ + STATIC_REQUIRE(13.89 == bx::abs(-13.89) ); + STATIC_REQUIRE(13.89 == bx::abs( 13.89) ); + STATIC_REQUIRE(0.0 == bx::abs(-0.0) ); + STATIC_REQUIRE(!bx::signBit(bx::abs(-0.0) ) ); + + STATIC_REQUIRE(-1 == bx::sign(-0.1389) ); + STATIC_REQUIRE( 0 == bx::sign( 0.0000) ); + STATIC_REQUIRE( 1 == bx::sign( 0.1389) ); + + STATIC_REQUIRE( bx::signBit(-0.0) ); + STATIC_REQUIRE(!bx::signBit( 0.0) ); + STATIC_REQUIRE( bx::signBit(-bx::kDoubleInfinity) ); + + STATIC_REQUIRE(-13.89 == bx::copySign(13.89, -1.0) ); + STATIC_REQUIRE( 13.89 == bx::copySign(13.89, 1.0) ); + STATIC_REQUIRE(bx::signBit(bx::copySign(0.0, -1.0) ) ); +} + +TEST_CASE("ldexp double", "[math][libm]") +{ + STATIC_REQUIRE(48.0 == bx::ldexp(3.0, 4) ); + STATIC_REQUIRE(0.4375 == bx::ldexp(7.0, -4) ); + STATIC_REQUIRE(0x1p1023 == bx::ldexp(1.0, 1023) ); + STATIC_REQUIRE(0x1p-1022 == bx::ldexp(1.0, -1022) ); + + STATIC_REQUIRE(0x1p-1074 == bx::ldexp(1.0, -1074) ); + STATIC_REQUIRE(0.0 == bx::ldexp(1.0, -1200) ); + STATIC_REQUIRE(1.0 == bx::ldexp(0x1p-1074, 1074) ); + REQUIRE(bx::kDoubleInfinity == bx::ldexp(1.0, 1200) ); + STATIC_REQUIRE(bx::signBit(bx::ldexp(-0.0, 10) ) ); + + for (int32_t yy = -1100; yy < 1100; ++yy) + { + for (double xx = -100.0; xx < 100.0; xx += 7.3) + { + REQUIRE(::ldexp(xx, yy) == bx::ldexp(xx, yy) ); + } + } +} + +TEST_CASE("frexp", "[math][libm]") +{ + constexpr auto frexpConst = []() + { + int32_t exp = 0; + const double frac = bx::frexp(12.0, &exp); + + return 0.75 == frac && 4 == exp; + }(); + STATIC_REQUIRE(frexpConst); + + bx::RngMwc rng; + + for (uint32_t ii = 0; ii < 100000; ++ii) + { + const uint64_t hi = uint64_t(rng.gen() ); + const uint64_t lo = uint64_t(rng.gen() ); + const uint64_t bits = (hi<<32) | lo; + const double vv = bx::bitsToDouble(bits); + + if (bx::isNan(vv) ) + { + continue; + } + + int32_t mine = 0; + int theirs = 0; + const double fa = bx::frexp(vv, &mine); + const double fb = ::frexp(vv, &theirs); + + REQUIRE(fa == fb); + REQUIRE(mine == theirs); + } + + for (uint32_t ii = 0; ii < 100000; ++ii) + { + const float vv = bx::bitsToFloat(rng.gen() ); + + if (bx::isNan(vv) ) + { + continue; + } + + int32_t mine = 0; + int theirs = 0; + const float fa = bx::frexp(vv, &mine); + const float fb = ::frexpf(vv, &theirs); + + REQUIRE(fa == fb); + REQUIRE(mine == theirs); + } +} + +TEST_CASE("modf, fmod, mod double", "[math][libm]") +{ + constexpr auto modfConst = []() + { + double integral = 0.0; + const double frac = bx::modf(-2.5, &integral); + + return -2.0 == integral && -0.5 == frac; + }(); + STATIC_REQUIRE(modfConst); + + STATIC_REQUIRE(-1.0 == bx::fmod(-7.0, 3.0) ); + STATIC_REQUIRE( 2.0 == bx::mod( -7.0, 3.0) ); + STATIC_REQUIRE( 1.0 == bx::fmod( 7.0, -3.0) ); + STATIC_REQUIRE( 1.0 == bx::fmod( 7.0, 3.0) ); + STATIC_REQUIRE( 1.0 == bx::mod( 7.0, 3.0) ); + + STATIC_REQUIRE(-1.0f == bx::fmod(-7.0f, 3.0f) ); + STATIC_REQUIRE( 2.0f == bx::mod( -7.0f, 3.0f) ); +} + +TEST_CASE("sqrt, log, exp double", "[math][libm]") +{ + STATIC_REQUIRE(4.0 == bx::sqrt(16.0) ); + STATIC_REQUIRE(0.0 == bx::sqrt(0.0) ); + STATIC_REQUIRE(bx::signBit(bx::sqrt(-0.0) ) ); + STATIC_REQUIRE(bx::isNan(bx::sqrt(-1.0) ) ); + STATIC_REQUIRE(bx::isEqual(1.4142135623730951, bx::sqrt(2.0), 1.0e-15) ); + + STATIC_REQUIRE(bx::isEqual(0.0, bx::log(1.0), 1.0e-15) ); + STATIC_REQUIRE(bx::isEqual(1.0, bx::log(2.71828182845904524), 1.0e-15) ); + STATIC_REQUIRE(bx::isEqual(0.69314718055994531, bx::log(2.0), 1.0e-15) ); + STATIC_REQUIRE(-bx::kDoubleInfinity == bx::log(0.0) ); + STATIC_REQUIRE(bx::isNan(bx::log(-1.0) ) ); + + STATIC_REQUIRE( 8.0 == bx::log2(256.0) ); + STATIC_REQUIRE(10.0 == bx::log2(1024.0) ); + + STATIC_REQUIRE(bx::isEqual(1.0, bx::exp(0.0), 1.0e-15) ); + STATIC_REQUIRE(bx::isEqual(2.71828182845904524, bx::exp(1.0), 1.0e-15) ); + STATIC_REQUIRE(0.0 == bx::exp(-800.0) ); + STATIC_REQUIRE(bx::kDoubleInfinity == bx::exp(800.0) ); + + STATIC_REQUIRE(1024.0 == bx::exp2(10.0) ); + + STATIC_REQUIRE( 81.0 == bx::pow( 3.0, 4.0) ); + STATIC_REQUIRE( -8.0 == bx::pow(-2.0, 3.0) ); + STATIC_REQUIRE( 4.0 == bx::pow(-2.0, 2.0) ); + STATIC_REQUIRE(0.0625 == bx::pow( 2.0, -4.0) ); + STATIC_REQUIRE( 1.0 == bx::pow( 7.0, 0.0) ); + STATIC_REQUIRE(bx::isNan(bx::pow(-2.0, 0.5) ) ); + STATIC_REQUIRE(bx::isEqual(1.4142135623730951, bx::pow(2.0, 0.5), 1.0e-15) ); + + for (double xx = 0.125; xx < 1.0e12; xx *= 1.7) + { + REQUIRE(bx::isEqual(::sqrt(xx), bx::sqrt(xx), 1.0e-15*xx) ); + REQUIRE(bx::isEqual(::log(xx), bx::log(xx), 1.0e-14) ); + } + + for (double xx = -700.0; xx < 700.0; xx += 3.1) + { + const double ref = ::exp(xx); + + REQUIRE(bx::isEqual(ref, bx::exp(xx), 1.0e-14*ref) ); + } +} + +TEST_CASE("sin, cos, tan double", "[math][libm]") +{ + STATIC_REQUIRE(0.0 == bx::sin(0.0) ); + STATIC_REQUIRE(1.0 == bx::cos(0.0) ); + STATIC_REQUIRE(0.0 == bx::tan(0.0) ); + + for (double xx = -100.0; xx < 100.0; xx += 0.37) + { + REQUIRE(bx::isEqual(::sin(xx), bx::sin(xx), 1.0e-14) ); + REQUIRE(bx::isEqual(::cos(xx), bx::cos(xx), 1.0e-14) ); + } + + REQUIRE(bx::isEqual(::sin(1.0e6), bx::sin(1.0e6), 1.0e-12) ); + REQUIRE(bx::isEqual(::cos(1.0e6), bx::cos(1.0e6), 1.0e-12) ); +} + +TEST_CASE("asin, acos, atan double", "[math][libm]") +{ + STATIC_REQUIRE(0.0 == bx::asin(0.0) ); + STATIC_REQUIRE(0.0 == bx::atan(0.0) ); + STATIC_REQUIRE(0.0 == bx::acos(1.0) ); + STATIC_REQUIRE(bx::isNan(bx::asin(1.5) ) ); + STATIC_REQUIRE(bx::isNan(bx::acos(1.5) ) ); + + STATIC_REQUIRE(bx::isEqual(1.57079632679489662, bx::asin(1.0), 1.0e-15) ); + STATIC_REQUIRE(bx::isEqual(3.14159265358979324, bx::acos(-1.0), 1.0e-15) ); + STATIC_REQUIRE(bx::isEqual(0.78539816339744831, bx::atan(1.0), 1.0e-15) ); + STATIC_REQUIRE(bx::isEqual(0.78539816339744831, bx::atan2(1.0,1.0), 1.0e-15) ); + + for (double xx = -1.0; xx < 1.0; xx += 0.013) + { + REQUIRE(bx::isEqual(::asin(xx), bx::asin(xx), 1.0e-14) ); + REQUIRE(bx::isEqual(::acos(xx), bx::acos(xx), 1.0e-14) ); + } + + for (double xx = -1000.0; xx < 1000.0; xx += 7.3) + { + REQUIRE(bx::isEqual(::atan(xx), bx::atan(xx), 1.0e-14) ); + } + + for (double yy = -10.0; yy < 10.0; yy += 1.3) + { + for (double xx = -10.0; xx < 10.0; xx += 1.3) + { + REQUIRE(bx::isEqual(::atan2(yy, xx), bx::atan2(yy, xx), 1.0e-14) ); + } + } + + REQUIRE(::atan2( 0.0, 1.0) == bx::atan2( 0.0, 1.0) ); + REQUIRE(::atan2(-0.0, 1.0) == bx::atan2(-0.0, 1.0) ); + REQUIRE(::atan2( 0.0, -1.0) == bx::atan2( 0.0, -1.0) ); + REQUIRE(::atan2(-0.0, -1.0) == bx::atan2(-0.0, -1.0) ); + REQUIRE(::atan2( 1.0, 0.0) == bx::atan2( 1.0, 0.0) ); + REQUIRE(::atan2(-1.0, 0.0) == bx::atan2(-1.0, 0.0) ); +} + +TEST_CASE("sinh, cosh, tanh double", "[math][libm]") +{ + STATIC_REQUIRE(0.0 == bx::sinh(0.0) ); + STATIC_REQUIRE(1.0 == bx::cosh(0.0) ); + STATIC_REQUIRE(0.0 == bx::tanh(0.0) ); + + for (double xx = -20.0; xx < 20.0; xx += 0.17) + { + const double sh = ::sinh(xx); + const double ch = ::cosh(xx); + + REQUIRE(bx::isEqual(sh, bx::sinh(xx), 1.0e-13*bx::abs(sh) ) ); + REQUIRE(bx::isEqual(ch, bx::cosh(xx), 1.0e-13*ch) ); + REQUIRE(bx::isEqual(::tanh(xx), bx::tanh(xx), 1.0e-14) ); + } + + for (double xx = 1.0e-9; xx < 0.1; xx *= 1.3) + { + REQUIRE(bx::isEqual(::sinh(xx), bx::sinh(xx), 1.0e-13*xx) ); + REQUIRE(bx::isEqual(::tanh(xx), bx::tanh(xx), 1.0e-13*xx) ); + } +} + +TEST_CASE("log1p, asinh, acosh, atanh", "[math][libm]") +{ + STATIC_REQUIRE(0.0 == bx::log1p(0.0) ); + STATIC_REQUIRE(0.0 == bx::asinh(0.0) ); + STATIC_REQUIRE(0.0 == bx::acosh(1.0) ); + STATIC_REQUIRE(0.0 == bx::atanh(0.0) ); + STATIC_REQUIRE(-bx::kDoubleInfinity == bx::log1p(-1.0) ); + STATIC_REQUIRE(bx::kDoubleInfinity == bx::atanh(1.0) ); + STATIC_REQUIRE(bx::isNan(bx::acosh(0.5) ) ); + STATIC_REQUIRE(bx::isNan(bx::atanh(1.5) ) ); + + STATIC_REQUIRE(bx::isEqual(0.88137358701954303, bx::asinh(1.0), 1.0e-15) ); + STATIC_REQUIRE(bx::isEqual(1.31695789692481670, bx::acosh(2.0), 1.0e-15) ); + STATIC_REQUIRE(bx::isEqual(0.54930614433405489, bx::atanh(0.5), 1.0e-15) ); + + for (double xx = 1.0e-9; xx < 1.0e9; xx *= 1.7) + { + REQUIRE(bx::isEqual(::log1p(xx), bx::log1p(xx), 1.0e-14*bx::abs(::log1p(xx) ) ) ); + REQUIRE(bx::isEqual(::asinh(xx), bx::asinh(xx), 1.0e-14*::asinh(xx) ) ); + REQUIRE(bx::isEqual(::asinh(-xx), bx::asinh(-xx), 1.0e-14*::asinh(xx) ) ); + } + + for (double xx = 0.0; xx < 1.0; xx += 0.0131) + { + REQUIRE(bx::isEqual(::atanh(xx), bx::atanh(xx), 1.0e-14) ); + } + + for (double tt = 0x1p-40; tt < 0x1p-18; tt *= 2.0) + { + const double xx = 1.0 + tt; + const double poly = 1.0 - tt/12.0 + 3.0*tt*tt/160.0; + const double ref = bx::sqrt(2.0*tt) * poly; + + REQUIRE(bx::isEqual(ref, bx::acosh(xx), 1.0e-14*ref) ); + } +} + +TEST_CASE("ldexp float range", "[math][libm]") +{ + STATIC_REQUIRE(0x1p127f == bx::ldexp(1.0f, 127) ); + STATIC_REQUIRE(0x1p-126f == bx::ldexp(1.0f, -126) ); + STATIC_REQUIRE(0x1p-149f == bx::ldexp(1.0f, -149) ); + STATIC_REQUIRE(0.0f == bx::ldexp(1.0f, -150) ); + STATIC_REQUIRE(0.0f == bx::ldexp(1.0f, -400) ); + STATIC_REQUIRE(1.0f == bx::ldexp(0x1p-149f, 149) ); + STATIC_REQUIRE(bx::signBit(bx::ldexp(-0.0f, 10) ) ); + + REQUIRE(bx::isInfinite(bx::ldexp(1.0f, 128) ) ); + REQUIRE(bx::isInfinite(bx::ldexp(1.0f, 400) ) ); + + for (int32_t yy = -300; yy < 300; ++yy) + { + for (float xx = -100.0f; xx < 100.0f; xx += 3.7f) + { + REQUIRE(::ldexpf(xx, yy) == bx::ldexp(xx, yy) ); + } + } +} + +TEST_CASE("exp float overflow", "[math][libm]") +{ + REQUIRE(bx::isInfinite(bx::exp( 89.0f) ) ); + REQUIRE(bx::isInfinite(bx::exp(100.0f) ) ); + REQUIRE(bx::isInfinite(bx::exp(500.0f) ) ); + REQUIRE(0.0f == bx::exp(-500.0f) ); + + REQUIRE(bx::isInfinite(bx::exp( 800.0) ) ); + REQUIRE(0.0 == bx::exp(-800.0) ); +} + +TEST_CASE("sinh, cosh, tanh saturation", "[math][libm]") +{ + REQUIRE( 1.0f == bx::tanh( 90.0f) ); + REQUIRE(-1.0f == bx::tanh(-90.0f) ); + REQUIRE( 1.0 == bx::tanh( 9876.0) ); + REQUIRE(-1.0 == bx::tanh(-9876.0) ); + + REQUIRE(bx::isFinite(bx::cosh(88.5f) ) ); + REQUIRE(bx::isFinite(bx::sinh(88.5f) ) ); + REQUIRE(bx::isFinite(bx::cosh(710.0) ) ); + REQUIRE(bx::isFinite(bx::sinh(710.0) ) ); + REQUIRE(bx::isInfinite(bx::cosh(720.0) ) ); + + for (float xx = 1.0e-7f; xx < 1.0f; xx *= 1.3f) + { + const float ref = float(bx::tanh(double(xx) ) ); + + REQUIRE(bx::isEqual(ref, bx::tanh(xx), 1.0e-6f*ref) ); + } +} + +TEST_CASE("log1p, asinh, acosh, atanh float", "[math][libm]") +{ + STATIC_REQUIRE(0.0f == bx::log1p(0.0f) ); + STATIC_REQUIRE(0.0f == bx::asinh(0.0f) ); + STATIC_REQUIRE(0.0f == bx::acosh(1.0f) ); + STATIC_REQUIRE(0.0f == bx::atanh(0.0f) ); + STATIC_REQUIRE(-bx::kFloatInfinity == bx::log1p(-1.0f) ); + STATIC_REQUIRE( bx::kFloatInfinity == bx::atanh(1.0f) ); + STATIC_REQUIRE(bx::isNan(bx::acosh(0.5f) ) ); + STATIC_REQUIRE(bx::isNan(bx::atanh(1.5f) ) ); + + for (float xx = 1.0e-6f; xx < 1.0e6f; xx *= 1.7f) + { + const float rl = float(bx::log1p(double(xx) ) ); + const float ra = float(bx::asinh(double(xx) ) ); + + REQUIRE(bx::isEqual(rl, bx::log1p(xx), 1.0e-6f*bx::abs(rl) ) ); + REQUIRE(bx::isEqual(ra, bx::asinh(xx), 1.0e-6f*ra) ); + REQUIRE(bx::isEqual(-ra, bx::asinh(-xx), 1.0e-6f*ra) ); + } + + for (float tt = 1.0e-6f; tt < 100.0f; tt *= 1.7f) + { + const float xx = 1.0f + tt; + const float ref = float(bx::acosh(double(xx) ) ); + + REQUIRE(bx::isEqual(ref, bx::acosh(xx), 1.0e-6f*ref) ); + } + + for (float xx = 0.0f; xx < 1.0f; xx += 0.013f) + { + const float ref = float(bx::atanh(double(xx) ) ); + + REQUIRE(bx::isEqual(ref, bx::atanh(xx), 1.0e-6f) ); + } +} + +TEST_CASE("odd and even symmetry", "[math][libm]") +{ + bx::RngMwc rng; + + for (uint32_t ii = 0; ii < 20000; ++ii) + { + const double uu = double(rng.gen() ) / 4294967296.0; + const double xx = bx::pow(10.0, uu*14.0 - 7.0); + + REQUIRE(bx::sin(-xx) == -bx::sin(xx) ); + REQUIRE(bx::cos(-xx) == bx::cos(xx) ); + REQUIRE(bx::tan(-xx) == -bx::tan(xx) ); + REQUIRE(bx::sinh(-xx) == -bx::sinh(xx) ); + REQUIRE(bx::cosh(-xx) == bx::cosh(xx) ); + REQUIRE(bx::tanh(-xx) == -bx::tanh(xx) ); + REQUIRE(bx::atan(-xx) == -bx::atan(xx) ); + REQUIRE(bx::asinh(-xx) == -bx::asinh(xx) ); + REQUIRE(bx::atan2(-xx, 3.0) == -bx::atan2(xx, 3.0) ); + REQUIRE(bx::atan2(-xx,-3.0) == -bx::atan2(xx,-3.0) ); + + const double ss = uu*2.0 - 1.0; + + REQUIRE(bx::asin(-ss) == -bx::asin(ss) ); + REQUIRE(bx::atanh(-ss*0.999) == -bx::atanh(ss*0.999) ); + } + + STATIC_REQUIRE(bx::signBit(bx::sin(-0.0) ) ); + STATIC_REQUIRE(bx::signBit(bx::tan(-0.0) ) ); + STATIC_REQUIRE(bx::signBit(bx::asin(-0.0) ) ); + STATIC_REQUIRE(bx::signBit(bx::atan(-0.0) ) ); + STATIC_REQUIRE(bx::signBit(bx::sinh(-0.0) ) ); + STATIC_REQUIRE(bx::signBit(bx::tanh(-0.0) ) ); + STATIC_REQUIRE(bx::signBit(bx::asinh(-0.0) ) ); + STATIC_REQUIRE(bx::signBit(bx::atanh(-0.0) ) ); + STATIC_REQUIRE(bx::signBit(bx::log1p(-0.0) ) ); + STATIC_REQUIRE(bx::signBit(bx::log1p(-0.0f) ) ); + + for (double xx = 1.0e-9; xx < 0.4375; xx *= 1.3) + { + REQUIRE(bx::atan(-xx) == -bx::atan(xx) ); + REQUIRE(bx::isEqual(::atan(-xx), bx::atan(-xx), 1.0e-14) ); + } +} + +TEST_CASE("pow negative base", "[math][libm]") +{ + STATIC_REQUIRE( 4.0f == bx::pow(-2.0f, 2.0f) ); + STATIC_REQUIRE( -8.0f == bx::pow(-2.0f, 3.0f) ); + STATIC_REQUIRE( 16.0f == bx::pow(-2.0f, 4.0f) ); + STATIC_REQUIRE(-0.125f == bx::pow(-2.0f, -3.0f) ); + STATIC_REQUIRE(bx::isNan(bx::pow(-2.0f, 0.5f) ) ); + STATIC_REQUIRE(bx::isNan(bx::pow(-2.0f, 1.5f) ) ); + + STATIC_REQUIRE( 4.0 == bx::pow(-2.0, 2.0) ); + STATIC_REQUIRE( -8.0 == bx::pow(-2.0, 3.0) ); + STATIC_REQUIRE(bx::isNan(bx::pow(-2.0, 0.5) ) ); + + STATIC_REQUIRE( 81.0f == bx::pow(3.0f, 4.0f) ); + STATIC_REQUIRE( 1024.0f == bx::pow(2.0f, 10.0f) ); + STATIC_REQUIRE( 0.0625f == bx::pow(2.0f, -4.0f) ); + STATIC_REQUIRE( 1024.0f == bx::exp2(10.0f) ); + STATIC_REQUIRE( 1.0f == bx::pow(7.0f, 0.0f) ); + STATIC_REQUIRE( 4.0f == bx::sqrt(16.0f) ); + + STATIC_REQUIRE(bx::kFloatInfinity == bx::pow(0.0f, -2.0f) ); + STATIC_REQUIRE(0.0f == bx::pow(0.0f, 2.0f) ); +} + +TEST_CASE("exp float accuracy", "[math][libm]") +{ + for (float xx = -103.0f; xx < 88.0f; xx += 0.37f) + { + const float ref = float(bx::exp(double(xx) ) ); + + if (ref >= bx::kFloatSmallest) + { + REQUIRE(bx::isEqual(ref, bx::exp(xx), 1.0e-6f*ref) ); + } + } + + REQUIRE(::expf(-100.0f) == bx::exp(-100.0f) ); + REQUIRE(::expf(-103.9f) == bx::exp(-103.9f) ); + REQUIRE(0.0f == bx::exp(-104.1f) ); + REQUIRE(bx::isFinite(bx::exp(88.7f) ) ); + REQUIRE(bx::isInfinite(bx::exp(88.8f) ) ); +} + +TEST_CASE("sinh, cosh near overflow", "[math][libm]") +{ + for (float xx = 88.0f; xx < 89.4f; xx += 0.013f) + { + const float rs = float(bx::sinh(double(xx) ) ); + const float rc = float(bx::cosh(double(xx) ) ); + + if (bx::isFinite(rs) ) + { + REQUIRE(bx::isEqual(rs, bx::sinh(xx), 1.0e-6f*rs) ); + REQUIRE(bx::isEqual(rc, bx::cosh(xx), 1.0e-6f*rc) ); + } + } + + REQUIRE(bx::isFinite(bx::sinh(709.5) ) ); + REQUIRE(bx::isFinite(bx::cosh(709.5) ) ); + REQUIRE(bx::isInfinite(bx::cosh(711.0) ) ); + REQUIRE(bx::signBit(bx::sinh(-709.5) ) ); +} + +TEST_CASE("trig large argument reduction", "[math][libm]") +{ + REQUIRE(bx::isFinite(bx::sin(1.0e308) ) ); + REQUIRE(bx::isFinite(bx::cos(1.0e308) ) ); + REQUIRE(bx::isFinite(bx::tan(1.0e308) ) ); + REQUIRE(bx::isFinite(bx::sin(7.8145278949113378e+202) ) ); + REQUIRE(bx::isFinite(bx::cos(7.8145278949113378e+202) ) ); + + { + int32_t quadrant = 0; + const double reduced = bx::reducePiHalfLarge(1.0e16, &quadrant); + + REQUIRE(1 == quadrant); + REQUIRE(bx::isEqual(0.4307553505349006*1.5707963267948966, reduced, 1.0e-16) ); + } + + { + int32_t quadrant = 0; + const double reduced = bx::reducePiHalfLarge(1.0e22, &quadrant); + + REQUIRE(3 == quadrant); + REQUIRE(bx::isEqual(0.35053490057448137839*1.5707963267948966, reduced, 1.0e-16) ); + } + + for (double xx = 1.0e3; xx < 1.6e6; xx *= 1.3) + { + int32_t quadrant = 0; + + const double large = bx::reducePiHalfLarge(xx, &quadrant); + const double small = bx::sin(xx); + + switch (quadrant) + { + case 0: REQUIRE(bx::isEqual(small, bx::sinPiQuarter(large), 1.0e-15) ); break; + case 1: REQUIRE(bx::isEqual(small, bx::cosPiQuarter(large), 1.0e-15) ); break; + case 2: REQUIRE(bx::isEqual(small, -bx::sinPiQuarter(large), 1.0e-15) ); break; + default: REQUIRE(bx::isEqual(small, -bx::cosPiQuarter(large), 1.0e-15) ); break; + } + } + + for (int32_t pp = 21; pp <= 52; ++pp) + { + const double xx = bx::ldexp(1.0, pp) + 1.0; + const double dd = bx::ldexp(1.0, -8); + + if (xx + dd - xx != dd) + { + continue; + } + + const double lhs = bx::sin(xx + dd); + const double rhs = bx::sin(xx)*bx::cos(dd) + bx::cos(xx)*bx::sin(dd); + + REQUIRE(bx::isEqual(lhs, rhs, 1.0e-15) ); + } + + bx::RngMwc rng; + + for (uint32_t ii = 0; ii < 50000; ++ii) + { + const uint64_t hi = uint64_t(rng.gen() ); + const uint64_t lo = uint64_t(rng.gen() ); + const uint64_t bits = (hi<<32) | lo; + const double xx = bx::bitsToDouble(bits); + + if (bx::isNan(xx) || bx::isInfinite(xx) ) + { + continue; + } + + const double ss = bx::sin(xx); + const double cc = bx::cos(xx); + + REQUIRE(bx::isFinite(ss) ); + REQUIRE(bx::isFinite(cc) ); + REQUIRE(bx::abs(ss) <= 1.0); + REQUIRE(bx::abs(cc) <= 1.0); + REQUIRE(bx::isEqual(1.0, ss*ss + cc*cc, 1.0e-15) ); + } +} + TEST_CASE("bitsToFloat, floatToBits, bitsToDouble, doubleToBits", "[math]") { STATIC_REQUIRE(0x12345678u == bx::floatToBits( bx::bitsToFloat (0x12345678u) ) ); diff --git a/tests/simd_test.cpp b/tests/simd_test.cpp index f92e604..03ac8b7 100644 --- a/tests/simd_test.cpp +++ b/tests/simd_test.cpp @@ -710,6 +710,116 @@ static void parity2(const char* _name, SimdFn _simd, ScalarFn _scalar, float _lo CHECK(0 == bad); } +static float accuracyInput(bx::RngMwc& _rng, float _lo, float _hi) +{ + const uint32_t r = _rng.gen(); + + switch (r & 31) + { + case 0: return 0.0f; + case 1: return clamp( 1.0f, _lo, _hi); + case 2: return clamp(-1.0f, _lo, _hi); + case 3: return clamp( 0.5f, _lo, _hi); + case 4: return _lo; + case 5: return _hi; + default: break; + } + + const float t = float(_rng.gen() & 0xffffff) * (1.0f/16777216.0f); + + return _lo + (_hi - _lo) * t; +} + +template +static void accuracy1(const char* _name, SimdFn _simd, RefFn _ref, ScaleFn _scale, double _tol, float _lo, float _hi) +{ + constexpr uint32_t kNumLanes = sizeof(Ty)/sizeof(float); + bx::RngMwc rng; + double maxErr = 0.0; + float badIn = 0.0f; + + for (uint32_t ii = 0; ii < 2048; ++ii) + { + alignas(sizeof(Ty) ) float in[kNumLanes]; + alignas(sizeof(Ty) ) float out[kNumLanes]; + + for (uint32_t jj = 0; jj < kNumLanes; ++jj) + { + in[jj] = accuracyInput(rng, _lo, _hi); + } + + simd_st(out, _simd(simd_ld(in) ) ); + + for (uint32_t jj = 0; jj < kNumLanes; ++jj) + { + const double ref = _ref(double(in[jj]) ); + + if (bx::abs(ref) > double(bx::kFloatLargest) ) + { + continue; + } + + const double err = bx::abs(double(out[jj]) - ref) / _scale(ref); + + if (err > maxErr) + { + maxErr = err; + badIn = in[jj]; + } + } + } + + INFO(_name << " lanes=" << kNumLanes << ": worst " << maxErr << " at " << badIn << ", budget " << _tol); + CHECK(maxErr <= _tol); +} + +template +static void accuracy2(const char* _name, SimdFn _simd, RefFn _ref, ScaleFn _scale, double _tol, float _lo0, float _hi0, float _lo1, float _hi1) +{ + constexpr uint32_t kNumLanes = sizeof(Ty)/sizeof(float); + bx::RngMwc rng; + double maxErr = 0.0; + float badIn0 = 0.0f; + float badIn1 = 0.0f; + + for (uint32_t ii = 0; ii < 2048; ++ii) + { + alignas(sizeof(Ty) ) float in0[kNumLanes]; + alignas(sizeof(Ty) ) float in1[kNumLanes]; + alignas(sizeof(Ty) ) float out[kNumLanes]; + + for (uint32_t jj = 0; jj < kNumLanes; ++jj) + { + in0[jj] = accuracyInput(rng, _lo0, _hi0); + in1[jj] = accuracyInput(rng, _lo1, _hi1); + } + + simd_st(out, _simd(simd_ld(in0), simd_ld(in1) ) ); + + for (uint32_t jj = 0; jj < kNumLanes; ++jj) + { + const double ref = _ref(double(in0[jj]), double(in1[jj]) ); + + if (bx::abs(ref) > double(bx::kFloatLargest) ) + { + continue; + } + + const double err = bx::abs(double(out[jj]) - ref) / _scale(ref); + + if (err > maxErr) + { + maxErr = err; + badIn0 = in0[jj]; + badIn1 = in1[jj]; + } + } + } + + INFO(_name << " lanes=" << kNumLanes << ": worst " << maxErr << " at " << badIn0 << ", " << badIn1 << ", budget " << _tol); + CHECK(maxErr <= _tol); +} + template static void parityAll() { @@ -721,21 +831,21 @@ static void parityAll() parity1("fract", [](Ty a ) { return simd_f32_fract(a); }, [](float a ) { return bx::fract(a); }, -1.0e6f, 1.0e6f); parity1("sign", [](Ty a ) { return simd_f32_sign(a); }, [](float a ) { return bx::sign(a); }, -10.0f, 10.0f, false); parity1("smoothstep", [](Ty a ) { return simd_f32_smoothstep(a); }, [](float a ) { return bx::smoothStep(a); }, -2.0f, 2.0f); - parity1("cos", [](Ty a ) { return simd_f32_cos(a); }, [](float a ) { return bx::cos(a); }, -100.0f, 100.0f); - parity1("sin", [](Ty a ) { return simd_f32_sin(a); }, [](float a ) { return bx::sin(a); }, -100.0f, 100.0f); - parity1("tan", [](Ty a ) { return simd_f32_tan(a); }, [](float a ) { return bx::tan(a); }, -100.0f, 100.0f); + accuracy1("cos", [](Ty a) { return simd_f32_cos(a); }, [](double a) { return bx::cos(a); }, [](double ) { return 1.0; }, 1.0e-5, -100.0f, 100.0f); + accuracy1("sin", [](Ty a) { return simd_f32_sin(a); }, [](double a) { return bx::sin(a); }, [](double ) { return 1.0; }, 1.0e-5, -100.0f, 100.0f); + accuracy1("tan", [](Ty a) { return simd_f32_tan(a); }, [](double a) { return bx::tan(a); }, [](double _r ) { return 1.0 + _r*_r; }, 1.0e-5, -100.0f, 100.0f); parity1("exp", [](Ty a ) { return simd_f32_exp(a); }, [](float a ) { return bx::exp(a); }, -90.0f, 90.0f); parity1("log", [](Ty a ) { return simd_f32_log(a); }, [](float a ) { return bx::log(a); }, 1.0e-30f, 1.0e30f, false); parity1("exp2", [](Ty a ) { return simd_f32_exp2(a); }, [](float a ) { return bx::exp2(a); }, -120.0f, 120.0f, false); parity1("log2", [](Ty a ) { return simd_f32_log2(a); }, [](float a ) { return bx::log2(a); }, 1.0e-30f, 1.0e30f, false); - parity1("acos", [](Ty a ) { return simd_f32_acos(a); }, [](float a ) { return bx::acos(a); }, -1.5f, 1.5f); - parity1("asin", [](Ty a ) { return simd_f32_asin(a); }, [](float a ) { return bx::asin(a); }, -1.5f, 1.5f); - parity1("atan", [](Ty a ) { return simd_f32_atan(a); }, [](float a ) { return bx::atan(a); }, -1000.0f, 1000.0f, false); + accuracy1("acos", [](Ty a) { return simd_f32_acos(a); }, [](double a) { return bx::acos(a); }, [](double ) { return 1.0; }, 1.0e-4, -1.5f, 1.5f); + accuracy1("asin", [](Ty a) { return simd_f32_asin(a); }, [](double a) { return bx::asin(a); }, [](double ) { return 1.0; }, 1.0e-4, -1.5f, 1.5f); + accuracy1("atan", [](Ty a) { return simd_f32_atan(a); }, [](double a) { return bx::atan(a); }, [](double ) { return 1.0; }, 5.0e-6, -1000.0f, 1000.0f); parity1("sinh", [](Ty a ) { return simd_f32_sinh(a); }, [](float a ) { return bx::sinh(a); }, -30.0f, 30.0f); parity1("cosh", [](Ty a ) { return simd_f32_cosh(a); }, [](float a ) { return bx::cosh(a); }, -30.0f, 30.0f); parity1("tanh", [](Ty a ) { return simd_f32_tanh(a); }, [](float a ) { return bx::tanh(a); }, -30.0f, 30.0f); - parity2("pow", [](Ty a, Ty b) { return simd_f32_pow(a, b); }, [](float a, float b) { return bx::pow(a, b); }, 0.01f, 100.0f, -10.0f, 10.0f, false); - parity2("atan2", [](Ty y, Ty x) { return simd_f32_atan2(y, x); }, [](float y, float x) { return bx::atan2(y, x); }, -100.0f, 100.0f, -100.0f, 100.0f, false); + accuracy2("pow", [](Ty a, Ty b) { return simd_f32_pow(a, b); }, [](double a, double b) { return bx::pow(a, b); }, [](double _r) { return bx::abs(_r); }, 1.0e-5, 0.01f, 100.0f, -10.0f, 10.0f); + accuracy2("atan2", [](Ty y, Ty x) { return simd_f32_atan2(y, x); }, [](double y, double x) { return bx::atan2(y, x); }, [](double ) { return 1.0; }, 5.0e-6, -100.0f, 100.0f, -100.0f, 100.0f); parity2("step", [](Ty e, Ty a) { return simd_f32_step(e, a); }, [](float e, float a) { return bx::step(e, a); }, -10.0f, 10.0f, -10.0f, 10.0f); parity2("mod", [](Ty a, Ty b) { return simd_f32_mod(a, b); }, [](float a, float b) { return bx::mod(a, b); }, -100.0f, 100.0f, -10.0f, 10.0f); } diff --git a/tests/string_test.cpp b/tests/string_test.cpp index f0c456a..31f19dd 100644 --- a/tests/string_test.cpp +++ b/tests/string_test.cpp @@ -514,6 +514,156 @@ TEST_CASE("fromString hex float", "[string]") REQUIRE(f == bx::kFloatSmallest); } +TEST_CASE("fromString round half to even", "[string]") +{ + double d = 0.0; + + REQUIRE(bx::fromString(&d, "1e23") ); + REQUIRE(bx::doubleToBits(d) == UINT64_C(0x44b52d02c7e14af6) ); + + REQUIRE(bx::fromString(&d, "9007199254740993") ); + REQUIRE(d == 9007199254740992.0); + + REQUIRE(bx::fromString(&d, "7.2057594037927933e16") ); + REQUIRE(d == 72057594037927936.0); + + REQUIRE(bx::fromString(&d, "9007199254740992.5") ); + REQUIRE(d == 9007199254740992.0); + + REQUIRE(bx::fromString(&d, "9007199254740993.5") ); + REQUIRE(d == 9007199254740994.0); + + REQUIRE(bx::fromString(&d, "1.00000000000000011102230246251565404236316680908203125") ); + REQUIRE(d == 1.0); + + REQUIRE(bx::fromString(&d, "1.000000000000000111022302462515654042363166809082031250000001") ); + REQUIRE(d == 1.0000000000000002); + + REQUIRE(bx::fromString(&d, "1.00000000000000011102230246251565404236316680908203124999999") ); + REQUIRE(d == 1.0); +} + +TEST_CASE("fromString subnormal", "[string]") +{ + double d = 1.0; + + REQUIRE(bx::fromString(&d, "5e-324") ); + REQUIRE(d == 0x1p-1074); + + REQUIRE(bx::fromString(&d, "2.5e-324") ); + REQUIRE(d == 0x1p-1074); + + REQUIRE(bx::fromString(&d, "2e-324") ); + REQUIRE(d == 0.0); + + REQUIRE(bx::fromString(&d, "1e-320") ); + REQUIRE(bx::doubleToBits(d) == UINT64_C(0x00000000000007e8) ); + + REQUIRE(bx::fromString(&d, "1e-308") ); + REQUIRE(bx::doubleToBits(d) == UINT64_C(0x000730d67819e8d2) ); + + for (int32_t exp = -1022; exp >= -1074; --exp) + { + char tmp[1200]; + const double value = bx::ldexp(1.0, exp); + + bx::toString(tmp, BX_COUNTOF(tmp), value); + + double read = 0.0; + REQUIRE(bx::fromString(&read, tmp) ); + REQUIRE(read == value); + } + + float f = 1.0f; + + REQUIRE(bx::fromString(&f, "1e-45") ); + REQUIRE(f == 0x1p-149f); + + REQUIRE(bx::fromString(&f, "1e-46") ); + REQUIRE(f == 0.0f); +} + +TEST_CASE("fromString rejects malformed", "[string]") +{ + const char* bad[] = { "", "abc", "+", "-", ".", "e10", " ", "0x", "-0x", "+.e5" }; + + for (const char* str : bad) + { + double d = 1234.5; + float f = 1234.5f; + + REQUIRE(!bx::fromString(&d, str) ); + REQUIRE(!bx::fromString(&f, str) ); + } + + double d = 0.0; + REQUIRE(bx::fromString(&d, "12abc") ); + REQUIRE(d == 12.0); +} + +TEST_CASE("fromString float rounds once", "[string]") +{ + float f = 0.0f; + + REQUIRE(bx::fromString(&f, "1e23") ); + REQUIRE(bx::floatToBits(f) == UINT32_C(0x65a96816) ); + + REQUIRE(bx::fromString(&f, "0.1") ); + REQUIRE(f == 0.1f); + + REQUIRE(bx::fromString(&f, "16777217") ); + REQUIRE(f == 16777216.0f); + + REQUIRE(bx::fromString(&f, "3.4028234663852886e+38") ); + REQUIRE(f == bx::kFloatLargest); + + REQUIRE(bx::fromString(&f, "1e39") ); + REQUIRE(bx::isInfinite(f) ); +} + +TEST_CASE("fromString long digit strings", "[string]") +{ + const double values[] = + { + 1.0, 0.1, 1.0e-300, 1.0e300, 0x1p-1074, 2.2250738585072014e-308, + 1.7976931348623157e308, 3.141592653589793, -79.39773355813419, + }; + + for (double value : values) + { + char tmp[1200]; + bx::snprintf(tmp, BX_COUNTOF(tmp), "%.760e", value); + + double read = 0.0; + REQUIRE(bx::fromString(&read, tmp) ); + REQUIRE(read == value); + } + + double d = 0.0; + REQUIRE(bx::fromString(&d, "1.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001") ); + REQUIRE(d == 1.0); +} + +TEST_CASE("fromString hex wide mantissa", "[string]") +{ + double d = 0.0; + + REQUIRE(bx::fromString(&d, "0x1fffffffffffffp0") ); + REQUIRE(d == 9007199254740991.0); + + REQUIRE(bx::fromString(&d, "0x1ffffffffffffff80p0") ); + REQUIRE(d == 36893488147419103232.0); + + REQUIRE(bx::fromString(&d, "0x1p-1074") ); + REQUIRE(d == 0x1p-1074); + + REQUIRE(bx::fromString(&d, "0x1p-1075") ); + REQUIRE(d == 0.0); + + REQUIRE(bx::fromString(&d, "0x3p-1075") ); + REQUIRE(d == 0x1p-1073); +} + static bool testFromString(int32_t _value, const char* _input) { char tmp[1024];