Add lagmma

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~~

hwy/contrib/math/math-inl.h

@@ -384,6 +384,21 @@ HWY_NOINLINE V CallTgamma(const D d, VecArg<V> 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 <class D, class V>
+HWY_INLINE V LogGamma(D d, V x);
+template <class D, class V>
+HWY_NOINLINE V CallLogGamma(const D d, VecArg<V> 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 <class FloatOrDouble>
 struct GammaImpl {};
 template <class FloatOrDouble>
+struct LgammaImpl {};
+template <class FloatOrDouble>
 struct LogImpl {};
 template <class FloatOrDouble>
 struct ExtPrecLog2ForPowImpl;
@@ -1007,6 +1024,47 @@ struct GammaImpl<float> {
   }
 };
 
+template <>
+struct LgammaImpl<float> {
+  // Use the Stirling asymptotic path for w >= StirlingLimit().
+  template <class D, class V = VFromD<D>, 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 <class D, class V = VFromD<D>, 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 <class D, class V = VFromD<D>, 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<double> {
   }
 };
 
+template <>
+struct LgammaImpl<double> {
+  // Use the Stirling asymptotic path for w >= StirlingLimit().
+  template <class D, class V = VFromD<D>, 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 <class D, class V = VFromD<D>, 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 <class D, class V = VFromD<D>, 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<decltype(d)> 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 D, class V = VFromD<D>, class M = MFromD<D>>
+HWY_INLINE V Lgamma(D d, V a) {
+  using T = TFromD<D>;
+  static_assert(IsFloat<T>(), "Only makes sense for floating-point");
+  LgammaImpl<T> impl;
+
+  const V kHalf = Set(d, static_cast<T>(0.5));
+  const V kOne = Set(d, static_cast<T>(1.0));
+  const V kTwo = Set(d, static_cast<T>(2.0));
+  const V kZero = Zero(d);
+  const V kPi = Set(d, static_cast<T>(+3.14159265358979323846264));
+  const V kLogPi = Set(d, static_cast<T>(+1.1447298858494001741434));
+  const V kLowCenter = Set(d, static_cast<T>(0.75));
+  const V kMidCenter = Set(d, static_cast<T>(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<T>::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 <class D, class V>
+HWY_INLINE V LogGamma(const D d, V x) {
+  const V kZero = Zero(d);
+  V result = impl::Lgamma(d, x);
+
+  const MFromD<D> 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

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