Document the new DPF surfaces in one command set, and test the field, half-tree, and multipoint edges.

Co-authored-by: Cursor <cursoragent@cursor.com>

include/grotto/range_lut.hpp

@@ -10,6 +10,11 @@
 ///          `sinh` / `cosh` / `tanh` / `sech` share the hyperbolic addition;
 ///          `coth` and `csch` use their principal small-argument tables;
 ///          `sqrt` / `inv` / `rsqrt` / `invsq` are dyadic lifts of `[1/2, 1]`.
+///          `expm1` and `log1p` use those same reductions. On `|x| < ln 2`
+///          and `|x| <= 1/2` they sum the Taylor series in extra bits so the
+///          cancellation in `exp(x)-1` and `ln(1+x)` is not rounded away.
+///          Outside that, `expm1` rebuilds `2^n exp(r) - 1` and `log1p`
+///          calls the positive logarithm on the exact fixed-point `1+x`.
 /// @copyright Copyright (c) 2019-2026 Ryan Henry and [others](@ref authors)
 /// @license Released under a GNU General Public v2.0 (GPLv2) license.
 
@@ -50,6 +55,8 @@ enum class reduced : unsigned
     inv,
     rsqrt,
     invsq,
+    expm1,
+    log1p,
 };
 
 HEDLEY_WARN_UNUSED_RESULT
@@ -60,12 +67,14 @@ namespace range_detail
 
 using u128 = unsigned __int128;
 
-constexpr u128 words(std::uint64_t hi, std::uint64_t lo)
+HEDLEY_CONST
+HEDLEY_NO_THROW
+constexpr u128 words(std::uint64_t hi, std::uint64_t lo) noexcept
 {
     return (u128{hi} << 64) | lo;
 }
 
-/// `value * 2^64`, rounded half away from zero. Values above `2^64` keep the
+/// @brief `value * 2^64`, rounded half away from zero. Values above `2^64` keep the
 /// high limb so the constant is not truncated.
 inline constexpr u128 ln2_64 = words(0, 12786308645202655660ULL);
 inline constexpr u128 inv_ln2_64 = words(1, 8166282121979093367ULL);
@@ -78,7 +87,7 @@ inline constexpr u128 two_over_pi_64 = words(0, 11743562013128004906ULL);
 inline constexpr u128 four_over_pi_64 = words(1, 5040379952546458196ULL);
 inline constexpr u128 pi_over_4_64 = words(0, 14488038916154245685ULL);
 
-/// `exp(2^{i-13}) * 2^64`.
+/// @brief `exp(2^{i-13}) * 2^64`.
 inline constexpr u128 exp_chunk_64[13] = {
     words(1, 2251937258231296ULL),
     words(1, 4504149427926357ULL),
@@ -266,7 +275,9 @@ inline std::int64_t eval_exp_at_scale(unsigned fractional_bits, std::int64_t raw
     return shift_pow2(exp_s, static_cast<int>(n_bin));
 }
 
-inline std::int64_t fractional_raw(std::int64_t raw, unsigned fractional_bits, std::int64_t & whole)
+HEDLEY_NO_THROW
+constexpr std::int64_t fractional_raw(std::int64_t raw, unsigned fractional_bits,
+    std::int64_t & whole) noexcept
 {
     const std::int64_t one = one_raw(fractional_bits);
     std::int64_t q = raw / one;
@@ -303,8 +314,13 @@ struct angle
     std::int64_t frac_raw;
 };
 
-/// `{ |x| * multiplier }` at this precision, with the integer part reduced
+/// @brief `{ |x| * multiplier }` at this precision, with the integer part reduced
 /// only as far as the low bits the quadrant logic reads.
+/// @param fractional_bits the number of fractional bits
+/// @param raw the underlying integer
+/// @param multiplier_64 the `multiplier_64`
+/// @return `{ |x| * multiplier }` at this precision, with the integer part reduced only as far as
+/// the low bits the quadrant logic reads
 inline angle reduce_positive(unsigned fractional_bits, std::int64_t raw, u128 multiplier_64)
 {
     const u128 scaled = magnitude_of(raw) * multiplier_64;
@@ -480,7 +496,10 @@ inline hyp sinh_cosh(unsigned fractional_bits, std::int64_t raw)
     return hyp{sh, ch};
 }
 
-/// `ln(2^{k+1} ± 1) / 2`, the saturation threshold used by `tanh` and `coth`.
+/// @brief `ln(2^{k+1} ± 1) / 2`, the saturation threshold used by `tanh` and `coth`.
+/// @param fractional_bits the number of fractional bits
+/// @param plus the `plus`
+/// @return `ln(2^{k+1} ± 1) / 2`, the saturation threshold used by `tanh` and `coth`
 inline std::int64_t beta_raw(unsigned fractional_bits, bool plus)
 {
     u128 ln = u128{fractional_bits + 1} * ln2_64;
@@ -492,11 +511,141 @@ inline std::int64_t beta_raw(unsigned fractional_bits, bool plus)
     return round_mag(ln, 65u - fractional_bits, false);
 }
 
-inline int half_pow_of(int power)
+HEDLEY_CONST
+HEDLEY_NO_THROW
+constexpr int half_pow_of(int power) noexcept
 {
     return (power & 1) != 0 ? (power - 1) / 2 : power / 2;
 }
 
+inline __int128 div_round_i128(__int128 num, int den)
+{
+    const bool neg = num < 0;
+    const auto mag = static_cast<u128>(neg ? -num : num);
+    const auto d = static_cast<u128>(den);
+    const u128 quot = (mag + d / 2) / d;
+    return neg ? -static_cast<__int128>(quot) : static_cast<__int128>(quot);
+}
+
+inline __int128 shr_round_i128(__int128 num, unsigned shift)
+{
+    if (shift == 0)
+        return num;
+    const bool neg = num < 0;
+    auto mag = static_cast<u128>(neg ? -num : num);
+    mag = (mag + (u128{1} << (shift - 1))) >> shift;
+    return neg ? -static_cast<__int128>(mag) : static_cast<__int128>(mag);
+}
+
+/// @brief `expm1` on `|x| < ln 2`, summed at `k+48` fractional bits.
+/// @param fractional_bits the number of fractional bits
+/// @param raw the underlying integer
+/// @return `expm1` on `|x| < ln 2`, summed at `k+48` fractional bits
+inline std::int64_t expm1_series(unsigned fractional_bits, std::int64_t raw)
+{
+    constexpr unsigned extra = 48;
+    __int128 power = static_cast<__int128>(raw) << extra;
+    __int128 acc = 0;
+    for (int n = 1; n <= 24; ++n)
+    {
+        const __int128 term = div_round_i128(power, n);
+        acc += term;
+        power = shr_round_i128(term * static_cast<__int128>(raw), fractional_bits);
+        if (power == 0)
+            break;
+    }
+    return round_i128(acc, extra);
+}
+
+inline std::int64_t eval_expm1(unsigned fractional_bits, std::int64_t raw)
+{
+    if (raw == 0)
+        return 0;
+    // Match `exp`: precisions below the 2^{-13} reduction evaluate one
+    // scale up and round once, so the power-of-two lift is not rounded early.
+    if (fractional_bits < 13)
+    {
+        const int lift = static_cast<int>(16u - fractional_bits);
+        const __int128 lifted_arg = static_cast<__int128>(raw) << lift;
+        if (lifted_arg > INT64_MAX || lifted_arg < INT64_MIN)
+            throw std::overflow_error("range lut: exponent overflow");
+        const std::int64_t lifted = eval_expm1(16, static_cast<std::int64_t>(lifted_arg));
+        return round_i128(lifted, static_cast<unsigned>(lift));
+    }
+    const std::int64_t ln2 = ln2_raw(fractional_bits);
+    if (raw > -ln2 && raw < ln2)
+        return expm1_series(fractional_bits, raw);
+
+    std::int64_t n_bin = raw / ln2;
+    std::int64_t remainder = raw - n_bin * ln2;
+    if (remainder < 0)
+    {
+        remainder += ln2;
+        --n_bin;
+    }
+    while (remainder >= ln2)
+    {
+        remainder -= ln2;
+        ++n_bin;
+    }
+    const std::int64_t exp_r = eval_exp_at_scale(fractional_bits, remainder);
+    const std::int64_t one = one_raw(fractional_bits);
+    if (n_bin >= 0)
+    {
+        const __int128 wide = static_cast<__int128>(shift_pow2(exp_r, static_cast<int>(n_bin))) - one;
+        if (wide > INT64_MAX || wide < INT64_MIN)
+            throw std::overflow_error("range lut: exponent overflow");
+        return static_cast<std::int64_t>(wide);
+    }
+    const int places = static_cast<int>(-n_bin);
+    if (places > static_cast<int>(fractional_bits) + 1)
+        return -one;
+    return shift_pow2(exp_r, -places) - one;
+}
+
+/// @brief `log1p` on `|x| <= 1/2`. Every term of a negative argument is negative.
+/// @param fractional_bits the number of fractional bits
+/// @param raw the underlying integer
+/// @return `log1p` on `|x| <= 1/2`
+inline std::int64_t log1p_series(unsigned fractional_bits, std::int64_t raw)
+{
+    constexpr unsigned extra = 48;
+    const bool xneg = raw < 0;
+    const std::int64_t mag_raw = xneg ? -raw : raw;
+    __int128 power = static_cast<__int128>(mag_raw) << extra;
+    __int128 acc = 0;
+    for (int n = 1; n <= 80; ++n)
+    {
+        const __int128 term = div_round_i128(power, n);
+        const bool neg = xneg || (n % 2 == 0);
+        acc += neg ? -term : term;
+        power = shr_round_i128(power * static_cast<__int128>(mag_raw), fractional_bits);
+        if (power == 0)
+            break;
+    }
+    return round_i128(acc, extra);
+}
+
+inline std::int64_t eval_log1p(unsigned fractional_bits, std::int64_t raw)
+{
+    const std::int64_t one = one_raw(fractional_bits);
+    if (raw == 0)
+        return 0;
+    if (raw <= -one)
+        throw std::domain_error("range lut: log1p argument is <= -1");
+    const std::int64_t half = one >> 1;
+    if (raw >= -half && raw <= half)
+        return log1p_series(fractional_bits, raw);
+    if (raw > INT64_MAX - one)
+    {
+        const std::int64_t ln_x = eval_ln_positive(fractional_bits, raw);
+        const u128 num = u128{1} << (2u * fractional_bits);
+        const auto corr = static_cast<std::int64_t>((num + static_cast<u128>(raw) / 2) / static_cast<u128>(raw));
+        return ln_x + corr;
+    }
+    return eval_ln_positive(fractional_bits, one + raw);
+}
+
 }  // namespace range_detail
 
 HEDLEY_WARN_UNUSED_RESULT
@@ -667,6 +816,10 @@ inline std::int64_t eval_reduced(reduced which, unsigned fractional_bits, std::i
             principal::invsq, fractional_bits, part.mantissa_raw);
         return shift_pow2(square, -2 * part.power);
     }
+    case reduced::expm1:
+        return eval_expm1(fractional_bits, raw);
+    case reduced::log1p:
+        return eval_log1p(fractional_bits, raw);
     }
     throw std::invalid_argument("range lut: unknown map");
 }