Implement Item 1 from gmp_gap_analysis.md

include/beman/big_int/charconv.hpp

@@ -423,6 +423,8 @@ from_chars(const char* const begin, const char* const end, basic_big_int<b, L, A
         // In any other case, we have the guarantee that at least one digit can be parsed.
         out.set_zero();
         out.grow(limbs_needed);
+        // The digit blocks are OR-ed in, so a reused destination must not keep stale limbs.
+        std::fill_n(out.limb_ptr(), limbs_needed, uint_multiprecision_t{0});
         while (true) {
             const auto        digit_block_length = std::min(current_end - current_begin, max_digits_per_iteration);
             const char* const digit_block_begin  = current_end - digit_block_length;

include/beman/big_int/detail/config.hpp

@@ -1080,7 +1080,13 @@ BEMAN_BIG_INT_END_NAMESPACE
 // GCOVR_EXCL_STOP
 // LCOV_EXCL_STOP
 
+// BEMAN_BIG_INT_DEBUG is defined when NDEBUG is not. Debug-only code, such as
+// BEMAN_BIG_INT_DEBUG_ASSERT, tests it instead of NDEBUG.
 #ifndef NDEBUG
+    #define BEMAN_BIG_INT_DEBUG 1
+#endif
+
+#ifdef BEMAN_BIG_INT_DEBUG
     #define BEMAN_BIG_INT_DEBUG_ASSERT(...) BEMAN_BIG_INT_ASSERT(__VA_ARGS__)
 #else
     // The requires expression makes sure that we still check for expression validity,

include/beman/big_int/detail/div_impl.hpp

@@ -52,7 +52,10 @@ constexpr std::size_t divide_unsigned_storage_size(const std::size_t dividend_li
 // Preconditions: trimmed operands, divisor.size() >= 2, dividend.size() >=
 // divisor.size(), quotient.size() >= dividend.size() - divisor.size() + 1,
 // remainder.size() >= dividend.size() + 1, no aliasing, and `scratch`
-// provides at least divide_unsigned_storage_size(...) limbs.
+// provides at least divide_unsigned_storage_size(...) limbs. The one
+// exception to no aliasing: `quotient` may start at the dividend's first limb
+// (same data()), because the dividend is copied whole into scratch before the
+// first quotient limb is written and never read again.
 // ---------------------------------------------------------------------------
 constexpr void divide_unsigned(const std::span<uint_multiprecision_t>       quotient,
                                const std::span<uint_multiprecision_t>       remainder,
@@ -66,7 +69,6 @@ constexpr void divide_unsigned(const std::span<uint_multiprecision_t>       quot
     BEMAN_BIG_INT_DEBUG_ASSERT(dividend.size() >= divisor.size());
     BEMAN_BIG_INT_DEBUG_ASSERT(quotient.size() >= dividend.size() - divisor.size() + 1);
     BEMAN_BIG_INT_DEBUG_ASSERT(remainder.size() >= dividend.size() + 1);
-    BEMAN_BIG_INT_DEBUG_ASSERT(quotient.data() != dividend.data());
     BEMAN_BIG_INT_DEBUG_ASSERT(quotient.data() != divisor.data());
     BEMAN_BIG_INT_DEBUG_ASSERT(remainder.data() != dividend.data());
     BEMAN_BIG_INT_DEBUG_ASSERT(remainder.data() != divisor.data());
@@ -908,6 +910,37 @@ void divide_barrett_preinv(std::span<uint_multiprecision_t>       quotient,
                            std::span<const uint_multiprecision_t> inv,
                            scratch_allocator_base&                scratch);
 
+// True when the runtime dispatchers send an m-limb by s-limb division to the
+// Barrett tier (see the gate rules above).
+[[nodiscard]] constexpr bool divide_takes_barrett(const std::size_t m, const std::size_t s) noexcept {
+    const bool barrett_march    = s >= barrett_march_cutoff && m / 16 >= s;
+    const bool barrett_march8   = s >= barrett_march8_cutoff && m / 8 >= s;
+    const bool barrett_quarter  = m >= barrett_quarter_cutoff && m / 4 >= s;
+    const bool barrett_balanced = m >= barrett_balanced_cutoff && m - s >= s;
+    return barrett_march || barrett_march8 || barrett_quarter || barrett_balanced;
+}
+
+// True when the runtime dispatchers take the schoolbook kernel (divide_unsigned)
+// for an m-limb by s-limb division: neither the Barrett nor the
+// divide-and-conquer gate fires. Constant evaluation always takes it.
+[[nodiscard]] constexpr bool divide_takes_schoolbook(const std::size_t m, const std::size_t s) noexcept {
+    if (divide_takes_barrett(m, s)) {
+        return false;
+    }
+    return !(s >= burnikel_ziegler_cutoff && m - s >= burnikel_ziegler_offset);
+}
+
+// Scratch limbs the schoolbook path of divide_dispatch (remainder_in_scratch =
+// false) or divide_dispatch_q (true: the discarded remainder also lives in
+// scratch) needs for an m-limb by s-limb division.
+[[nodiscard]] constexpr std::size_t
+divide_schoolbook_storage_size(const std::size_t m, const std::size_t s, const bool remainder_in_scratch) noexcept {
+    return divide_unsigned_storage_size(m, s) + (remainder_in_scratch ? m + 1 : 0);
+}
+
+// Scratch limbs up to which a caller may keep division workspace on the stack.
+inline constexpr std::size_t small_division_stack_limbs = 128;
+
 // ---------------------------------------------------------------------------
 // Top-level division dispatcher (counterpart of multiply_dispatch): the Barrett
 // gates above are checked first, then the divide-and-conquer path, which needs
@@ -930,16 +963,12 @@ constexpr void divide_dispatch(const std::span<uint_multiprecision_t>       quot
         const std::size_t m = dividend.size();
         const std::size_t s = divisor.size();
 
-        const bool barrett_march    = s >= barrett_march_cutoff && m / 16 >= s;
-        const bool barrett_march8   = s >= barrett_march8_cutoff && m / 8 >= s;
-        const bool barrett_quarter  = m >= barrett_quarter_cutoff && m / 4 >= s;
-        const bool barrett_balanced = m >= barrett_balanced_cutoff && m - s >= s;
-        if (barrett_march || barrett_march8 || barrett_quarter || barrett_balanced) {
-            divide_barrett(quotient, remainder, dividend, divisor, alloc);
-            return;
-        }
-        if (s >= burnikel_ziegler_cutoff && m - s >= burnikel_ziegler_offset) {
-            divide_burnikel_ziegler(quotient, remainder, dividend, divisor, alloc);
+        if (!divide_takes_schoolbook(m, s)) {
+            if (divide_takes_barrett(m, s)) {
+                divide_barrett(quotient, remainder, dividend, divisor, alloc);
+            } else {
+                divide_burnikel_ziegler(quotient, remainder, dividend, divisor, alloc);
+            }
             return;
         }
     }
@@ -967,17 +996,13 @@ constexpr void divide_dispatch_q(const std::span<uint_multiprecision_t>       qu
     const std::size_t s = divisor.size();
 
     if BEMAN_BIG_INT_IS_NOT_CONSTEVAL {
-        const bool barrett_march    = s >= barrett_march_cutoff && m / 16 >= s;
-        const bool barrett_march8   = s >= barrett_march8_cutoff && m / 8 >= s;
-        const bool barrett_quarter  = m >= barrett_quarter_cutoff && m / 4 >= s;
-        const bool barrett_balanced = m >= barrett_balanced_cutoff && m - s >= s;
-        if (barrett_march || barrett_march8 || barrett_quarter || barrett_balanced) {
+        if (divide_takes_barrett(m, s)) {
             const std::span<uint_multiprecision_t> r = scratch.allocate(m + 1);
             divide_barrett(quotient, r, dividend, divisor, alloc);
             scratch.deallocate(m + 1);
             return;
         }
-        if (s >= burnikel_ziegler_cutoff && m - s >= burnikel_ziegler_offset) {
+        if (!divide_takes_schoolbook(m, s)) {
             divide_quotient(quotient, dividend, divisor, alloc);
             return;
         }

include/beman/big_int/detail/mul_impl.hpp

@@ -1307,18 +1307,37 @@ static_assert(mul_slice_zones_valid(), "invalid slicing zone table (see mul_slic
 // the shape-aware dispatch names above.
 #define BEMAN_BIG_INT_HAS_SHAPE_DISPATCH 1
 
+// ---------------------------------------------------------------------------
+// multiply_single_limb that also writes result[a.size()] when the product has
+// no carry limb, so all a.size() + 1 limbs are defined. Same contract
+// otherwise: `result` has space for a.size() + 1 limbs, no aliasing. Returns
+// the number of significant result limbs.
+// ---------------------------------------------------------------------------
+constexpr std::size_t multiply_by_limb_full(const std::span<uint_multiprecision_t>       result,
+                                            const std::span<const uint_multiprecision_t> a,
+                                            const uint_multiprecision_t                  val) noexcept {
+    const std::size_t n = multiply_single_limb(result, a, val);
+    if (n == a.size()) {
+        result[n] = 0;
+    }
+    return n;
+}
+
 // ---------------------------------------------------------------------------
 // Runtime multiplication tier ladders (src/mul_dispatch.cpp): operands must
-// be trimmed with at least two limbs; `result` pre-zeroed, sized for the
-// full product, non-aliasing. Kernel workspaces come from the type-erased
-// heap hooks. Returns the trimmed significant size.
+// be trimmed with at least two limbs; `result` sized for the full product,
+// non-aliasing, and need NOT be pre-zeroed: limbs [0, a.size() + b.size()) are
+// all written (tiers that accumulate zero-fill them first), limbs beyond are
+// untouched. Kernel workspaces come from the type-erased heap hooks. Returns
+// the trimmed significant size.
 // ---------------------------------------------------------------------------
 std::size_t multiply_runtime(std::span<uint_multiprecision_t>       result,
                              std::span<const uint_multiprecision_t> a,
                              std::span<const uint_multiprecision_t> b,
                              const scratch_heap_source&             heap);
 
-// Benchmark/test-only variants of multiply_runtime, same contract. The
+// Benchmark/test-only variants of multiply_runtime, same contract (a result
+// that is not pre-zeroed is fine). The
 // schoolbook (min < karatsuba_cutoff) and power-of-two shortcuts still apply in
 // both. _sliced forces slicing of the top-level product whenever max > min,
 // ignoring the ratio test, and skips the FFT gate even where it would take the
@@ -1341,18 +1360,37 @@ std::size_t square_runtime(std::span<uint_multiprecision_t>       result,
 
 // Untrimmed/short-operand-tolerant runtime product for the src-side callers
 // (the division tiers and mulmod): trims, takes the single-limb shortcuts,
-// then runs the tier ladder. The runtime mirror of multiply_dispatch.
+// then runs the tier ladder. The runtime mirror of multiply_dispatch, with the
+// same non-zeroed-result contract.
 std::size_t multiply_runtime_any(std::span<uint_multiprecision_t>       result,
                                  std::span<const uint_multiprecision_t> a_untrimmed,
                                  std::span<const uint_multiprecision_t> b_untrimmed,
                                  const scratch_heap_source&             heap);
 
+// Master switch for the header-side basecase shortcut in multiply_dispatch.
+inline constexpr bool mul_header_basecase_enabled = true;
+
+// Operand length below which multiply_dispatch calls multiply_basecase_runtime
+// instead of building the scratch hooks and entering the ladder. One constant
+// for every configuration; src/mul_dispatch.cpp static_asserts that it is at
+// most every Karatsuba cutoff of the compiled configuration.
+inline constexpr std::size_t mul_header_basecase_limbs = 32;
+
+// Product (or square, when a and b are the same span) of trimmed operands with
+// at least two limbs each, min(a.size(), b.size()) < mul_header_basecase_limbs,
+// via the basecase kernels only: no scratch. Same result contract as
+// multiply_runtime (no pre-zeroing needed). Returns the trimmed size.
+std::size_t multiply_basecase_runtime(std::span<uint_multiprecision_t>       result,
+                                      std::span<const uint_multiprecision_t> a,
+                                      std::span<const uint_multiprecision_t> b) noexcept;
+
 // ---------------------------------------------------------------------------
 // result <- a * p2 where p2 is a trimmed power of two: a shifted copy of `a`
 // placed at limb offset p2.size() - 1, bit-shifted by countr_zero(p2.back()).
 // O(n) with no scratch, versus a full kernel dispatch (see GMP mpz_mul_2exp).
-// `result` must be pre-zeroed, have space for a.size() + p2.size() limbs, and
-// must NOT alias `a`. Returns the number of significant result limbs.
+// Writes all a.size() + p2.size() limbs of `result` (it need NOT be pre-zeroed),
+// which must have space for them and must NOT alias `a`. Returns the number of
+// significant result limbs.
 // ---------------------------------------------------------------------------
 constexpr std::size_t multiply_power_of_two(const std::span<uint_multiprecision_t>       result,
                                             const std::span<const uint_multiprecision_t> a,
@@ -1367,14 +1405,20 @@ constexpr std::size_t multiply_power_of_two(const std::span<uint_multiprecision_
     const auto        bit_shift   = static_cast<unsigned>(std::countr_zero(p2.back()));
     const auto        dst         = result.subspan(limb_offset);
 
+    std::ranges::fill(result.first(limb_offset), uint_multiprecision_t{0});
     std::ranges::copy(a, dst.begin());
-    return limb_offset + shift_left_n(dst, a.size(), bit_shift);
+    const std::size_t shifted = shift_left_n(dst, a.size(), bit_shift);
+    if (shifted == a.size()) {
+        dst[shifted] = 0;
+    }
+    return limb_offset + shifted;
 }
 
 // ---------------------------------------------------------------------------
 // Squaring dispatcher: result <- a * a using kernels that exploit the
 // symmetric cross products. Same contract as multiply_dispatch: `result` must
-// be pre-zeroed with space for 2 * a.size() limbs and must NOT alias `a`.
+// have space for 2 * a.size() limbs (it need NOT be pre-zeroed) and must NOT
+// alias `a`.
 // `a` must be trimmed with at least two limbs. Returns the number of
 // significant result limbs.
 // ---------------------------------------------------------------------------
@@ -1393,7 +1437,10 @@ std::size_t square_dispatch(const std::span<uint_multiprecision_t>       result,
 
 // ---------------------------------------------------------------------------
 // Top-level multiplication dispatcher.
-// `result` must be pre-zeroed and have space for a.size() + b.size() limbs.
+// `result` must have space for a.size() + b.size() limbs; it need NOT be
+// pre-zeroed. On return limbs [0, trimmed a.size() + trimmed b.size()) hold the
+// product (the high ones zero when the product is shorter); limbs beyond are
+// untouched, so a caller that wants zeros there fills them itself.
 // `result` must NOT alias `a` or `b`.
 // Returns the number of significant result limbs (trimmed).
 // ---------------------------------------------------------------------------
@@ -1419,10 +1466,10 @@ constexpr std::size_t multiply_dispatch(const std::span<uint_multiprecision_t>
         return hi != 0 ? 2 : 1;
     }
     if (a.size() == 1) {
-        return multiply_single_limb(result, b, a[0]);
+        return multiply_by_limb_full(result, b, a[0]);
     }
     if (b.size() == 1) {
-        return multiply_single_limb(result, a, b[0]);
+        return multiply_by_limb_full(result, a, b[0]);
     }
 
     // This check has 0 runtime cost, but could speed up/reduce depth of constant evaluation
@@ -1435,6 +1482,7 @@ constexpr std::size_t multiply_dispatch(const std::span<uint_multiprecision_t>
         }
         // Squaring halves the widening muls, and so the consteval step count.
         if (a.data() == b.data() && a.size() == b.size()) {
+            std::ranges::fill(result.first(2 * a.size()), uint_multiprecision_t{0});
             square_long(result, a);
             return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), 2 * a.size()});
         }
@@ -1444,6 +1492,11 @@ constexpr std::size_t multiply_dispatch(const std::span<uint_multiprecision_t>
     // keeps long multiplication (the recursive tiers could blow up consteval
     // step limits, and the kernels are runtime-compiled).
     if BEMAN_BIG_INT_IS_NOT_CONSTEVAL {
+        if constexpr (mul_header_basecase_enabled) {
+            if (std::min(a.size(), b.size()) < mul_header_basecase_limbs) {
+                return multiply_basecase_runtime(result, a, b);
+            }
+        }
         const scratch_allocator<Allocator> hooks(alloc);
         return multiply_runtime(result, a, b, hooks.heap());
     } else {