diff --git a/g3doc/op_wishlist.md b/g3doc/op_wishlist.md index a2cdb30ee0..9f0e460863 100644 --- a/g3doc/op_wishlist.md +++ b/g3doc/op_wishlist.md @@ -37,7 +37,7 @@ _mm512_getmant (f32/f64) High-precision! Consider copying from SLEEF. See #1650. -fmod, ilogb, lgamma, logb, modf, nextafter, nexttoward, scalbn +fmod, ilogb, logb, modf, nextafter, nexttoward, scalbn ### Remaining STL functions for hwy/contrib/algo @@ -195,3 +195,4 @@ For SVE (svld1sb_u32)+WASM? Compiler can probably already fuse. * ~~pow~~ * ~~Lookup32~~ * ~~tgamma~~ +* ~~lgamma~~ diff --git a/hwy/contrib/math/math-inl.h b/hwy/contrib/math/math-inl.h index 493fecff6f..b3d9b86de3 100644 --- a/hwy/contrib/math/math-inl.h +++ b/hwy/contrib/math/math-inl.h @@ -384,6 +384,21 @@ HWY_NOINLINE V CallTgamma(const D d, VecArg x) { return Tgamma(d, x); } +/** + * Highway SIMD version of std::lgamma(x). + * + * Valid Lane Types: float32, float64 + * Max Error: ULP = 6 (float32), 10 (float64) + * Valid Range: float32(0, +FLT_MAX], float64(0, +DBL_MAX] + * @return natural log of the absolute value of the gamma function of 'x' + */ +template +HWY_INLINE V LogGamma(D d, V x); +template +HWY_NOINLINE V CallLogGamma(const D d, VecArg x) { + return LogGamma(d, x); +} + /** * Highway SIMD version of SinCos. * Compute the sine and cosine at the same time @@ -734,6 +749,8 @@ struct ExpImpl {}; template struct GammaImpl {}; template +struct LgammaImpl {}; +template struct LogImpl {}; template struct ExtPrecLog2ForPowImpl; @@ -1007,6 +1024,47 @@ struct GammaImpl { } }; +template <> +struct LgammaImpl { + // Use the Stirling asymptotic path for w >= StirlingLimit(). + template , HWY_IF_F32_D(D)> + HWY_INLINE V StirlingLimit(D d) { + return Set(d, +4.0f); + } + // Recurrence steps to reduce w into [1, 2): ceil(StirlingLimit - 2). + static constexpr int kReduceSteps = 2; + + // logGamma(x) = (x-1)*(x-2)*MidPoly(t) on [1, 2], argument t = x - 1.5. + template , HWY_IF_F32_D(D)> + HWY_INLINE V MidPoly(D d, V t) { + const V c0 = Set(d, +0.00238831134f); + const V c1 = Set(d, -0.00397213376f); + const V c2 = Set(d, +0.00535511656f); + const V c3 = Set(d, -0.00927653644f); + const V c4 = Set(d, +0.016803567f); + const V c5 = Set(d, -0.0313190425f); + const V c6 = Set(d, +0.0629112789f); + const V c7 = Set(d, -0.145959697f); + const V c8 = Set(d, +0.483128951f); + return Estrin(t, c8, c7, c6, c5, c4, c3, c2, c1, c0); + } + + // logGamma(x) = (x-1)*LowPoly(t) on [0.5, 1), argument t = x - 0.75. + template , HWY_IF_F32_D(D)> + HWY_INLINE V LowPoly(D d, V t) { + const V c0 = Set(d, -1.40863339f); + const V c1 = Set(d, +1.18704979f); + const V c2 = Set(d, -0.81872104f); + const V c3 = Set(d, +0.723862927f); + const V c4 = Set(d, -0.673221395f); + const V c5 = Set(d, +0.655416733f); + const V c6 = Set(d, -0.719965008f); + const V c7 = Set(d, +1.09094739f); + const V c8 = Set(d, -0.813123806f); + return Estrin(t, c8, c7, c6, c5, c4, c3, c2, c1, c0); + } +}; + #if HWY_HAVE_FLOAT64 && HWY_HAVE_INTEGER64 template <> @@ -1056,6 +1114,70 @@ struct GammaImpl { } }; +template <> +struct LgammaImpl { + // Use the Stirling asymptotic path for w >= StirlingLimit(). + template , HWY_IF_F64_D(D)> + HWY_INLINE V StirlingLimit(D d) { + return Set(d, +12.0); + } + + // Recurrence steps to reduce w into [1, 2): ceil(StirlingLimit - 2). + static constexpr int kReduceSteps = 10; + + // logGamma(x) = (x-1)*(x-2)*MidPoly(t) on [1, 2], argument t = x - 1.5. + template , HWY_IF_F64_D(D)> + HWY_INLINE V MidPoly(D d, V t) { + const V c0 = Set(d, +2.7601736760716535e-05); + const V c1 = Set(d, -4.3458333543198137e-05); + const V c2 = Set(d, +3.581888030972756e-05); + const V c3 = Set(d, -5.6969476167985029e-05); + const V c4 = Set(d, +0.0001073270909194498); + const V c5 = Set(d, -0.0001715654244759322); + const V c6 = Set(d, +0.00027099300921864671); + const V c7 = Set(d, -0.00043754969398502409); + const V c8 = Set(d, +0.00071155062268775133); + const V c9 = Set(d, -0.0011644589445235255); + const V c10 = Set(d, +0.0019229144861301855); + const V c11 = Set(d, -0.0032120714838724603); + const V c12 = Set(d, +0.0054464578966048312); + const V c13 = Set(d, -0.0094256225300762171); + const V c14 = Set(d, +0.016797098630444426); + const V c15 = Set(d, -0.031308487499847056); + const V c16 = Set(d, +0.062911401074568038); + const V c17 = Set(d, -0.14595989591431094); + const V c18 = Set(d, +0.48312895054098087); + return Estrin(t, c18, c17, c16, c15, c14, c13, c12, c11, c10, c9, c8, c7, + c6, c5, c4, c3, c2, c1, c0); + } + + // logGamma(x) = (x-1)*LowPoly(t) on [0.5, 1), argument t = x - 0.75. + template , HWY_IF_F64_D(D)> + HWY_INLINE V LowPoly(D d, V t) { + const V c0 = Set(d, -15.548294260608042); + const V c1 = Set(d, +12.284843340929337); + const V c2 = Set(d, -5.118578439889296); + const V c3 = Set(d, +4.0915134400148023); + const V c4 = Set(d, -3.8614319065200999); + const V c5 = Set(d, +3.1052563349956928); + const V c6 = Set(d, -2.4708509229390843); + const V c7 = Set(d, +2.0105874995065292); + const V c8 = Set(d, -1.6495725282916478); + const V c9 = Set(d, +1.3638192653661416); + const V c10 = Set(d, -1.1396020445425534); + const V c11 = Set(d, +0.96515178838859494); + const V c12 = Set(d, -0.83191158526198028); + const V c13 = Set(d, +0.73473508727089576); + const V c14 = Set(d, -0.6729879446836875); + const V c15 = Set(d, +0.65522436942156825); + const V c16 = Set(d, -0.71996611036258851); + const V c17 = Set(d, +1.0909482962451835); + const V c18 = Set(d, -0.8131238057251815); + return Estrin(t, c18, c17, c16, c15, c14, c13, c12, c11, c10, c9, c8, c7, + c6, c5, c4, c3, c2, c1, c0); + } +}; + #endif template <> @@ -1798,12 +1920,79 @@ HWY_INLINE V Gamma(D d, V a) { const RebindToSigned di; const M odd = RebindMask(d, Ne(And(ConvertTo(di, ra), Set(di, 1)), Zero(di))); - const V s_signed = IfThenElse(odd, Neg(s_mag), s_mag); + const V s_signed = MaskedXorOr(s_mag, odd, s_mag, SignBit(d)); const V refl = Div(kPi, Mul(s_signed, gamma_w)); return IfThenElse(neg, refl, gamma_w); } +template , class M = MFromD> +HWY_INLINE V Lgamma(D d, V a) { + using T = TFromD; + static_assert(IsFloat(), "Only makes sense for floating-point"); + LgammaImpl impl; + + const V kHalf = Set(d, static_cast(0.5)); + const V kOne = Set(d, static_cast(1.0)); + const V kTwo = Set(d, static_cast(2.0)); + const V kZero = Zero(d); + const V kPi = Set(d, static_cast(+3.14159265358979323846264)); + const V kLogPi = Set(d, static_cast(+1.1447298858494001741434)); + const V kLowCenter = Set(d, static_cast(0.75)); + const V kMidCenter = Set(d, static_cast(1.5)); + const V kStirlingLimit = impl.StirlingLimit(d); + + // Reduce to w >= 0.5: w = 1 - a when a < 0.5. + const M neg = Lt(a, kHalf); + const V w = MaskedSubOr(a, neg, kOne, a); + + // [0.5, 1): logGamma = (w-1)*LowPoly(w - 0.75). + const V low_poly = impl.LowPoly(d, Sub(w, kLowCenter)); + V low; + if constexpr (HWY_NATIVE_FMA) { + low = MulSub(w, low_poly, low_poly); + } else { + low = Mul(Sub(w, kOne), low_poly); + } + + // Shift w into [1, 2) via logGamma(y) = logGamma(y-1) + log(y-1). + V y = w; + V acc = kZero; + for (int i = 0; i < LgammaImpl::kReduceSteps; ++i) { + const M down = Ge(y, kTwo); + y = MaskedSubOr(y, down, y, kOne); + acc = MaskedAddOr(acc, down, acc, impl::Log(d, y)); + } + + // [1, 2): logGamma = (y-1)*(y-2)*MidPoly(y - 1.5) + acc. + const V t = Sub(y, kMidCenter); + V zero_factors; + if constexpr (HWY_NATIVE_FMA) { + const V q = Sub(kTwo, y); + zero_factors = NegMulAdd(y, q, q); + } else { + zero_factors = Mul(Sub(y, kOne), Sub(y, kTwo)); + } + const V mid = MulAdd(zero_factors, impl.MidPoly(d, t), acc); + + // w >= StirlingLimit: Stirling series in double-double. + V stir_lo; + const V stir_hi = StirlingLogGamma(d, w, stir_lo); + const V stir = Add(stir_hi, stir_lo); + + const M is_low = Lt(w, kOne); + const M is_stir = Ge(w, kStirlingLimit); + V loggamma_w = IfThenElse(is_low, low, mid); + loggamma_w = IfThenElse(is_stir, stir, loggamma_w); + + // For a < 0.5: logGamma(a) = log(pi) - log|sin(pi*a)| - logGamma(1 - a). + const V frac = Sub(a, Round(a)); + const V s_mag = Abs(Sin(d, Mul(kPi, frac))); + const V refl = Sub(Sub(kLogPi, impl::Log(d, s_mag)), loggamma_w); + + return IfThenElse(neg, refl, loggamma_w); +} + // SinCos // Based on "sse_mathfun.h", by Julien Pommier // http://gruntthepeon.free.fr/ssemath/ @@ -2853,6 +3042,18 @@ HWY_INLINE V Tgamma(const D d, V x) { return result; } +template +HWY_INLINE V LogGamma(const D d, V x) { + const V kZero = Zero(d); + V result = impl::Lgamma(d, x); + + const MFromD is_pole = And(Eq(x, Round(x)), Le(x, kZero)); + result = IfThenElse(is_pole, Inf(d), result); + result = IfThenElse(IsInf(x), Inf(d), result); + result = IfThenElse(IsNaN(x), x, result); + return result; +} + // NOLINTNEXTLINE(google-readability-namespace-comments) } // namespace HWY_NAMESPACE } // namespace hwy diff --git a/hwy/contrib/math/math_test.cc b/hwy/contrib/math/math_test.cc index 8f96cd7e1d..45d3999428 100644 --- a/hwy/contrib/math/math_test.cc +++ b/hwy/contrib/math/math_test.cc @@ -72,6 +72,9 @@ DEFINE_MATH_TEST(Cbrt, DEFINE_MATH_TEST(Tgamma, std::tgamma, CallTgamma, +0.5f, +35.0f, 6, std::tgamma, CallTgamma, +0.5, +171.6, 8) +DEFINE_MATH_TEST(LogGamma, + std::lgamma, CallLogGamma, +0.5f, +1000.0f, 6, + std::lgamma, CallLogGamma, +0.5, +1000.0, 10) // clang-format on struct TestPow { @@ -208,6 +211,7 @@ HWY_EXPORT_AND_TEST_P(HwyMathTest, TestAllLog1p); HWY_EXPORT_AND_TEST_P(HwyMathTest, TestAllLog2); HWY_EXPORT_AND_TEST_P(HwyMathTest, TestAllCbrt); HWY_EXPORT_AND_TEST_P(HwyMathTest, TestAllTgamma); +HWY_EXPORT_AND_TEST_P(HwyMathTest, TestAllLogGamma); HWY_EXPORT_AND_TEST_P(HwyMathTest, TestAllPow); HWY_AFTER_TEST(); } // namespace