36725334879

Committed 30 Sep 2026 01:56PM UTC coverage: 91.621% (-4.5%) from 96.121%

Build # 36725334879

Build Type

Pull #385

github

Committed by

web-flow

Commit Message

Merge 07c2a2318 into d53a78873

Pull Request Pull Request #385: Updating Tuning Information

Coverage Stats

167 of 184 new or added lines in 5 files covered. (90.76%)

318 existing lines in 4 files now uncovered.

5992 of 6540 relevant lines covered (91.62%)

6484122.34 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

90.48

/src/mul_dispatch.cpp

// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
// SPDX-License-Identifier: BSL-1.0

#include <beman/big_int/detail/mul_impl.hpp>

#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <optional>
#include <span>

#include <beman/big_int/detail/config.hpp>
#include <beman/big_int/detail/multiply_long_runtime.hpp>
#include <beman/big_int/detail/scratch_allocator.hpp>
#include <beman/big_int/detail/span_ops.hpp>
#include <beman/big_int/detail/square_long_runtime.hpp>

// The runtime multiplication tier ladders, compiled once. The header
// dispatchers (multiply_dispatch / square_dispatch) keep the constexpr
// small-operand and constant-evaluation paths and forward every runtime
// multi-limb product here; kernel workspaces come from the type-erased heap
// hooks, so a single compiled definition serves every allocator.

namespace beman::big_int::detail {

std::size_t square_runtime(const std::span<uint_multiprecision_t>       result,
                           const std::span<const uint_multiprecision_t> a,
                           const scratch_heap_source&                   heap) {
    BEMAN_BIG_INT_DEBUG_ASSERT(a.size() >= 2);
    BEMAN_BIG_INT_DEBUG_ASSERT(a.back() != 0);
    BEMAN_BIG_INT_DEBUG_ASSERT(result.size() >= 2 * a.size());
    BEMAN_BIG_INT_DEBUG_ASSERT(result.data() != a.data());

    const std::size_t n            = a.size();
    const std::size_t result_total = 2 * n;

    // Tiny squares: plain schoolbook beats the squaring basecase.
    if (n < square_long_cutoff) {
        if BEMAN_BIG_INT_IS_NOT_CONSTEVAL {
            ::beman_big_int_multiply_long_runtime(
                result.first(result_total).data(), a.data(), a.size(), a.data(), a.size());
        } else {
            multiply_long(result.first(result_total), a, a);
        }
        return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});
    }

    // (2^k)^2 = 2^(2k): a shifted copy beats any squaring kernel.
    if (is_power_of_two_span(a)) {
        return multiply_power_of_two(result, a, a);
    }

    if (n < square_karatsuba_cutoff) {
        ::beman_big_int_square_long_runtime(result.first(result_total).data(), a.data(), n);
        return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});
    }

    // The FFT kernel packs into 64-bit words, so it is gated to 64-bit limbs (the branch is discarded otherwise).
    if constexpr (width_v<uint_multiprecision_t> == 64) {
        if (square_fft_worthwhile(n)) {
#if defined(BEMAN_BIG_INT_SIMD_MUL)
            scratch_heap_array<double>        fp_ws(heap, square_fft_fp_storage_size(n));
            scratch_heap_array<std::uint64_t> int_ws(heap, square_fft_int_storage_size(n));
            square_fft(result.first(result_total), a, fp_ws.span(), int_ws.span());
#else
            scratch_heap_array<std::uint64_t> ws(heap, square_fft_storage_size(n));
            square_fft(result.first(result_total), a, ws.span());
#endif
            return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});
        }
    }

    const auto in_heap_scratch = [&](const std::size_t limbs, auto&& kernel) {
        scratch_heap_array<uint_multiprecision_t> buf(heap, limbs);
        scratch_allocator_base                    scratch(buf.data(), limbs);
        kernel(scratch);
    };

    if (n < square_toom_cook_3_cutoff) {
        in_heap_scratch(karatsuba_storage_size(n), [&](scratch_allocator_base& scratch) {
            square_karatsuba(result.first(result_total), a, scratch);
        });
    } else if (n < square_toom_cook_4_cutoff) {
        in_heap_scratch(toom_cook_3_storage_size(n), [&](scratch_allocator_base& scratch) {
            square_toom_cook_3(result.first(result_total), a, scratch);
        });
    } else if (n < square_toom_cook_6_5_cutoff) {
        in_heap_scratch(toom_cook_4_storage_size(n), [&](scratch_allocator_base& scratch) {
            square_toom_cook_4(result.first(result_total), a, scratch);
        });
    } else if (n < square_toom_cook_8_5_cutoff) {
        in_heap_scratch(toom_cook_6_5_storage_size(n), [&](scratch_allocator_base& scratch) {
            square_toom_cook_6_5(result.first(result_total), a, scratch);
        });
    } else {
        in_heap_scratch(toom_cook_8_5_storage_size(n), [&](scratch_allocator_base& scratch) {
            square_toom_cook_8_5(result.first(result_total), a, scratch);
        });
    }

    return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});
}

namespace {

enum class slice_mode { automatic, forced, disabled };

// Karatsuba over a stack buffer. Kept out of line so the buffer never sits in
// the frames of the dispatch and slicing recursion.
BEMAN_BIG_INT_NOINLINE void karatsuba_on_stack(const std::span<uint_multiprecision_t>       result,
                                               const std::span<const uint_multiprecision_t> a,
                                               const std::span<const uint_multiprecision_t> b) noexcept {
    uint_multiprecision_t  stack_buf[karatsuba_stack_threshold];
    scratch_allocator_base scratch(stack_buf, karatsuba_stack_threshold);
    multiply_karatsuba(result, a, b, scratch);
}

// Karatsuba / Toom-3 / Toom-4 / Toom-6.5 / Toom-8.5 chosen by the shorter
// operand. `result` is pre-zeroed with exactly a.size() + b.size() limbs and
// `scratch` holds toom_ladder_storage_size(min, max) limbs.
void toom_ladder(const std::span<uint_multiprecision_t>       result,
                 const std::span<const uint_multiprecision_t> a,
                 const std::span<const uint_multiprecision_t> b,
                 scratch_allocator_base&                      scratch) noexcept {
    const std::size_t min_size = std::min(a.size(), b.size());
    if (min_size < toom_cook_3_cutoff) {
        multiply_karatsuba(result, a, b, scratch);
    } else if (min_size < toom_cook_4_cutoff) {
        multiply_toom_cook_3(result, a, b, scratch);
    } else if (min_size < toom_cook_6_5_cutoff) {
        multiply_toom_cook_4(result, a, b, scratch);
    } else if (min_size < toom_cook_8_5_cutoff) {
        multiply_toom_cook_6_5(result, a, b, scratch);
    } else {
        multiply_toom_cook_8_5(result, a, b, scratch);
    }
}

void run_ladder(const std::span<uint_multiprecision_t>       result,
                const std::span<const uint_multiprecision_t> a,
                const std::span<const uint_multiprecision_t> b,
                const scratch_heap_source&                   heap) {
    const std::size_t min_size = std::min(a.size(), b.size());
    const std::size_t s        = std::max(a.size(), b.size());
    if (min_size < toom_cook_3_cutoff && karatsuba_storage_size(s) <= karatsuba_stack_threshold) {
        karatsuba_on_stack(result, a, b);
        return;
    }
    const std::size_t                         limbs = toom_ladder_storage_size(min_size, s);
    scratch_heap_array<uint_multiprecision_t> buf(heap, limbs);
    scratch_allocator_base                    scratch(buf.data(), limbs);
    toom_ladder(result, a, b, scratch);
}

std::size_t multiply_runtime_impl(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,
                                  slice_mode                             mode);

// Euclidean slicing: cuts the longer operand into pieces of the shorter one's
// length m and accumulates piece * short at each offset. In automatic mode
// pieces are cut while the remainder is at least 2m or still enters slicing
// (mul_should_slice), so the last piece is below 2m for every ratio; the
// forced mode cuts until at most m remains. Full pieces (at least m limbs) run
// the FFT when the cost-model gate takes the m x piece product (automatic mode
// only; 64-bit limbs) and the Karatsuba/Toom ladder otherwise, over scratch
// allocated up front (the FFT workspace on first use, sized for the largest
// piece); a real
// tail or a piece that trims below m re-dispatches through multiply_runtime_any
// (it may slice the other way or take the FFT) and allocates its own scratch,
// possibly after piece 0 has been written. On x86 (model off, gate = floor on
// min) slicing is only entered when min = m is below the floor, so no piece
// takes the FFT and this path is a no-op there. Callers pre-zero the result and
// discard it if an allocation throws. Requires trimmed operands,
// min >= karatsuba_cutoff and a pre-zeroed result.
std::size_t multiply_sliced(const std::span<uint_multiprecision_t>       result,
                            const std::span<const uint_multiprecision_t> a,
                            const std::span<const uint_multiprecision_t> b,
                            const scratch_heap_source&                   heap,
                            const bool                                   force) {
    const auto        lng = a.size() >= b.size() ? a : b;
    const auto        sht = a.size() >= b.size() ? b : a;
    const std::size_t m   = sht.size();
    const std::size_t n   = lng.size();

    std::size_t pieces_before_last = 0;
    std::size_t rem                = n;
    while (force ? rem > m : (rem >= 2 * m || mul_should_slice(m, rem))) {
        ++pieces_before_last;
        rem -= m;
    }

    const std::size_t s_max    = std::max(m, rem);
    const std::size_t tmp_size = m + s_max;
    const std::size_t ws_size  = toom_ladder_storage_size(m, s_max);

    scratch_heap_array<uint_multiprecision_t> buf(heap, tmp_size + ws_size);
    const std::span<uint_multiprecision_t>    tmp(buf.data(), tmp_size);
#if defined(BEMAN_BIG_INT_SIMD_MUL)
    [[maybe_unused]] std::optional<scratch_heap_array<double>>        fft_fp_ws;
    [[maybe_unused]] std::optional<scratch_heap_array<std::uint64_t>> fft_int_ws;
#else
    [[maybe_unused]] std::optional<scratch_heap_array<std::uint64_t>> fft_ws;
#endif

    for (std::size_t i = 0; i <= pieces_before_last; ++i) {
        const std::size_t off   = i * m;
        const std::size_t len   = (i < pieces_before_last) ? m : rem;
        const auto        piece = lng.subspan(off, len);
        const std::size_t pn    = trimmed_size_span(piece);
        if (pn == 1 && piece[0] == 0) {
            continue;
        }
        const auto p   = piece.first(pn);
        const auto dst = (i == 0 ? result : tmp).first(pn + m);
        if (i != 0) {
            std::ranges::fill(dst, uint_multiprecision_t{0});
        }

        if (pn < m) {
            multiply_runtime_any(dst, p, sht, heap);
        } else {
            bool piece_done = false;
            // Full pieces (m x pn, pn >= m) consult the FFT gate first, except in the forced mode.
            if constexpr (width_v<uint_multiprecision_t> == 64) {
                if (!force && fft_mul_worthwhile(m, pn)) {
#if defined(BEMAN_BIG_INT_SIMD_MUL)
                    if (!fft_fp_ws) {
                        fft_fp_ws.emplace(heap, fft_mul_fp_storage_size(m, s_max));
                        fft_int_ws.emplace(heap, fft_mul_int_storage_size(m, s_max));
                    }
                    multiply_fft(dst, p, sht, fft_fp_ws->span(), fft_int_ws->span());
#else
                    if (!fft_ws) {
                        fft_ws.emplace(heap, fft_mul_storage_size(m, s_max));
                    }
                    multiply_fft(dst, p, sht, fft_ws->span());
#endif
                    piece_done = true;
                }
            }
            if (!piece_done) {
                scratch_allocator_base ws(buf.data() + tmp_size, ws_size);
                toom_ladder(dst, p, sht, ws);
            }
        }

        if (i != 0) {
            add_shifted(result.first(off + len + m), off, dst);
        }
    }
    return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), n + m});
}

std::size_t multiply_runtime_impl(const std::span<uint_multiprecision_t>       result,
                                  const std::span<const uint_multiprecision_t> a,
                                  const std::span<const uint_multiprecision_t> b,
                                  const scratch_heap_source&                   heap,
                                  const slice_mode                             mode) {
    BEMAN_BIG_INT_DEBUG_ASSERT(a.size() >= 2);
    BEMAN_BIG_INT_DEBUG_ASSERT(b.size() >= 2);
    BEMAN_BIG_INT_DEBUG_ASSERT(a.back() != 0);
    BEMAN_BIG_INT_DEBUG_ASSERT(b.back() != 0);
    BEMAN_BIG_INT_DEBUG_ASSERT(result.size() >= a.size() + b.size());

    // x * x and x *= x pass the same span twice, so squaring detection is
    // a pointer compare that almost always fails fast for ordinary mul.
    if (a.data() == b.data() && a.size() == b.size()) {
        return square_runtime(result, a, heap);
    }

    const std::size_t min_size     = std::min(a.size(), b.size());
    const std::size_t max_size     = std::max(a.size(), b.size());
    const std::size_t result_total = a.size() + b.size();

    if (min_size < karatsuba_cutoff) {
        // Schoolbook long multiplication runtime fallback (known to be on the runtime path).
        ::beman_big_int_multiply_long_runtime(result.data(), a.data(), a.size(), b.data(), b.size());
        return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});
    }

    // Power-of-two operands reduce to a shifted copy of the other operand.
    // This is only worth checking if we're about to do a big number mul anyway.
    if (is_power_of_two_span(b)) {
        return multiply_power_of_two(result, a, b);
    }
    if (is_power_of_two_span(a)) {
        return multiply_power_of_two(result, b, a);
    }

    // The FFT kernel packs into 64-bit words, so it is gated to 64-bit limbs.
    if constexpr (width_v<uint_multiprecision_t> == 64) {
        const bool use_fft = mode != slice_mode::forced && fft_mul_worthwhile(min_size, max_size);
        if (use_fft) {
#if defined(BEMAN_BIG_INT_SIMD_MUL)
            scratch_heap_array<double>        fp_ws(heap, fft_mul_fp_storage_size(a.size(), b.size()));
            scratch_heap_array<std::uint64_t> int_ws(heap, fft_mul_int_storage_size(a.size(), b.size()));
            multiply_fft(result.first(result_total), a, b, fp_ws.span(), int_ws.span());
#else
            scratch_heap_array<std::uint64_t> ws(heap, fft_mul_storage_size(a.size(), b.size()));
            multiply_fft(result.first(result_total), a, b, ws.span());
#endif
            return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});
        }
    }

    const bool slice = mode == slice_mode::forced      ? max_size > min_size
                       : mode == slice_mode::automatic ? mul_should_slice(min_size, max_size)
                                                       : false;
    if (slice) {
        return multiply_sliced(result, a, b, heap, mode == slice_mode::forced);
    }

    run_ladder(result.first(result_total), a, b, heap);
    return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});
}

} // namespace

std::size_t multiply_runtime(const std::span<uint_multiprecision_t>       result,
                             const std::span<const uint_multiprecision_t> a,
                             const std::span<const uint_multiprecision_t> b,
                             const scratch_heap_source&                   heap) {
    return multiply_runtime_impl(result, a, b, heap, slice_mode::automatic);
}

std::size_t multiply_runtime_sliced(const std::span<uint_multiprecision_t>       result,
                                    const std::span<const uint_multiprecision_t> a,
                                    const std::span<const uint_multiprecision_t> b,
                                    const scratch_heap_source&                   heap) {
    return multiply_runtime_impl(result, a, b, heap, slice_mode::forced);
}

std::size_t multiply_runtime_unsliced(const std::span<uint_multiprecision_t>       result,
                                      const std::span<const uint_multiprecision_t> a,
                                      const std::span<const uint_multiprecision_t> b,
                                      const scratch_heap_source&                   heap) {
    return multiply_runtime_impl(result, a, b, heap, slice_mode::disabled);
}

std::size_t multiply_runtime_any(const std::span<uint_multiprecision_t>       result,
                                 const std::span<const uint_multiprecision_t> a_untrimmed,
                                 const std::span<const uint_multiprecision_t> b_untrimmed,
                                 const scratch_heap_source&                   heap) {
    BEMAN_BIG_INT_DEBUG_ASSERT(!a_untrimmed.empty());
    BEMAN_BIG_INT_DEBUG_ASSERT(!b_untrimmed.empty());
    BEMAN_BIG_INT_DEBUG_ASSERT(result.size() >= a_untrimmed.size() + b_untrimmed.size());

    const auto a = a_untrimmed.first(trimmed_size_span(a_untrimmed));
    const auto b = b_untrimmed.first(trimmed_size_span(b_untrimmed));

    if (a.size() == 1 && b.size() == 1) {
        const auto [lo, hi] = widening_mul(a[0], b[0]);
        result[0]           = lo;
        result[1]           = hi;
        return hi != 0 ? 2 : 1;
    }
    if (a.size() == 1) {
        return multiply_single_limb(result, b, a[0]);
    }
    if (b.size() == 1) {
        return multiply_single_limb(result, a, b[0]);
    }

    return multiply_runtime(result, a, b, heap);
}

} // namespace beman::big_int::detail

1	// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
2	// SPDX-License-Identifier: BSL-1.0
3
4	#include <beman/big_int/detail/mul_impl.hpp>
5
6	#include <algorithm>
7	#include <cstddef>
8	#include <cstdint>
9	#include <optional>
10	#include <span>
11
12	#include <beman/big_int/detail/config.hpp>
13	#include <beman/big_int/detail/multiply_long_runtime.hpp>
14	#include <beman/big_int/detail/scratch_allocator.hpp>
15	#include <beman/big_int/detail/span_ops.hpp>
16	#include <beman/big_int/detail/square_long_runtime.hpp>
17
18	// The runtime multiplication tier ladders, compiled once. The header
19	// dispatchers (multiply_dispatch / square_dispatch) keep the constexpr
20	// small-operand and constant-evaluation paths and forward every runtime
21	// multi-limb product here; kernel workspaces come from the type-erased heap
22	// hooks, so a single compiled definition serves every allocator.
23
24	namespace beman::big_int::detail {
25
26	std::size_t square_runtime(const std::span<uint_multiprecision_t> result,	958,222 ✔
27	const std::span<const uint_multiprecision_t> a,
28	const scratch_heap_source& heap) {
29	BEMAN_BIG_INT_DEBUG_ASSERT(a.size() >= 2);	958,222 ✔
30	BEMAN_BIG_INT_DEBUG_ASSERT(a.back() != 0);	958,222 ✔
31	BEMAN_BIG_INT_DEBUG_ASSERT(result.size() >= 2 * a.size());	958,222 ✔
32	BEMAN_BIG_INT_DEBUG_ASSERT(result.data() != a.data());	958,222 ✔
33
34	const std::size_t n = a.size();	958,222 ✔
35	const std::size_t result_total = 2 * n;	958,222 ✔
36
37	// Tiny squares: plain schoolbook beats the squaring basecase.
38	if (n < square_long_cutoff) {	958,222 ✔
39	if BEMAN_BIG_INT_IS_NOT_CONSTEVAL {
40	::beman_big_int_multiply_long_runtime(	14,652 ✔
41	result.first(result_total).data(), a.data(), a.size(), a.data(), a.size());	7,326 ✔
42	} else {
43	multiply_long(result.first(result_total), a, a);
44	}
45	return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});	7,326 ✔
46	}
47
48	// (2^k)^2 = 2^(2k): a shifted copy beats any squaring kernel.
49	if (is_power_of_two_span(a)) {	950,896 ✔
50	return multiply_power_of_two(result, a, a);	149 ✔
51	}
52
53	if (n < square_karatsuba_cutoff) {	950,747 ✔
54	::beman_big_int_square_long_runtime(result.first(result_total).data(), a.data(), n);	950,082 ✔
55	return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});	950,082 ✔
56	}
57
58	// The FFT kernel packs into 64-bit words, so it is gated to 64-bit limbs (the branch is discarded otherwise).
59	if constexpr (width_v<uint_multiprecision_t> == 64) {
60	if (square_fft_worthwhile(n)) {	665 ✔
61	#if defined(BEMAN_BIG_INT_SIMD_MUL)
62	scratch_heap_array<double> fp_ws(heap, square_fft_fp_storage_size(n));
63	scratch_heap_array<std::uint64_t> int_ws(heap, square_fft_int_storage_size(n));
64	square_fft(result.first(result_total), a, fp_ws.span(), int_ws.span());
65	#else
NEW 66	scratch_heap_array<std::uint64_t> ws(heap, square_fft_storage_size(n));	×
NEW 67	square_fft(result.first(result_total), a, ws.span());	×
68	#endif
NEW 69	return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});	×
NEW 70	}	×
71	}
72
73	const auto in_heap_scratch = [&](const std::size_t limbs, auto&& kernel) {	665 ✔
74	scratch_heap_array<uint_multiprecision_t> buf(heap, limbs);	665 ✔
75	scratch_allocator_base scratch(buf.data(), limbs);	665 ✔
76	kernel(scratch);	665 ✔
77	};	1,330 ✔
78
79	if (n < square_toom_cook_3_cutoff) {	665 ✔
80	in_heap_scratch(karatsuba_storage_size(n), [&](scratch_allocator_base& scratch) {	430 ✔
81	square_karatsuba(result.first(result_total), a, scratch);	430 ✔
82	});	430 ✔
83	} else if (n < square_toom_cook_4_cutoff) {	235 ✔
84	in_heap_scratch(toom_cook_3_storage_size(n), [&](scratch_allocator_base& scratch) {	217 ✔
85	square_toom_cook_3(result.first(result_total), a, scratch);	217 ✔
86	});	217 ✔
87	} else if (n < square_toom_cook_6_5_cutoff) {	18 ✔
88	in_heap_scratch(toom_cook_4_storage_size(n), [&](scratch_allocator_base& scratch) {	×
89	square_toom_cook_4(result.first(result_total), a, scratch);	×
90	});	×
91	} else if (n < square_toom_cook_8_5_cutoff) {	18 ✔
NEW 92	in_heap_scratch(toom_cook_6_5_storage_size(n), [&](scratch_allocator_base& scratch) {	×
NEW 93	square_toom_cook_6_5(result.first(result_total), a, scratch);	×
NEW 94	});	×
95	} else {
96	in_heap_scratch(toom_cook_8_5_storage_size(n), [&](scratch_allocator_base& scratch) {	18 ✔
97	square_toom_cook_8_5(result.first(result_total), a, scratch);	18 ✔
98	});	18 ✔
99	}
100
101	return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});	665 ✔
102	}
103
104	namespace {
105
106	enum class slice_mode { automatic, forced, disabled };
107
108	// Karatsuba over a stack buffer. Kept out of line so the buffer never sits in
109	// the frames of the dispatch and slicing recursion.
110	BEMAN_BIG_INT_NOINLINE void karatsuba_on_stack(const std::span<uint_multiprecision_t> result,	595 ✔
111	const std::span<const uint_multiprecision_t> a,
112	const std::span<const uint_multiprecision_t> b) noexcept {
113	uint_multiprecision_t stack_buf[karatsuba_stack_threshold];
114	scratch_allocator_base scratch(stack_buf, karatsuba_stack_threshold);	595 ✔
115	multiply_karatsuba(result, a, b, scratch);	595 ✔
116	}	595 ✔
117
118	// Karatsuba / Toom-3 / Toom-4 / Toom-6.5 / Toom-8.5 chosen by the shorter
119	// operand. `result` is pre-zeroed with exactly a.size() + b.size() limbs and
120	// `scratch` holds toom_ladder_storage_size(min, max) limbs.
121	void toom_ladder(const std::span<uint_multiprecision_t> result,	25,639 ✔
122	const std::span<const uint_multiprecision_t> a,
123	const std::span<const uint_multiprecision_t> b,
124	scratch_allocator_base& scratch) noexcept {
125	const std::size_t min_size = std::min(a.size(), b.size());	25,639 ✔
126	if (min_size < toom_cook_3_cutoff) {	25,639 ✔
127	multiply_karatsuba(result, a, b, scratch);	19,300 ✔
128	} else if (min_size < toom_cook_4_cutoff) {	6,339 ✔
129	multiply_toom_cook_3(result, a, b, scratch);	5,188 ✔
130	} else if (min_size < toom_cook_6_5_cutoff) {	1,151 ✔
NEW 131	multiply_toom_cook_4(result, a, b, scratch);	×
132	} else if (min_size < toom_cook_8_5_cutoff) {	1,151 ✔
NEW 133	multiply_toom_cook_6_5(result, a, b, scratch);	×
134	} else {
135	multiply_toom_cook_8_5(result, a, b, scratch);	1,151 ✔
136	}
137	}	25,639 ✔
138
139	void run_ladder(const std::span<uint_multiprecision_t> result,	17,908 ✔
140	const std::span<const uint_multiprecision_t> a,
141	const std::span<const uint_multiprecision_t> b,
142	const scratch_heap_source& heap) {
143	const std::size_t min_size = std::min(a.size(), b.size());	17,908 ✔
144	const std::size_t s = std::max(a.size(), b.size());	17,908 ✔
145	if (min_size < toom_cook_3_cutoff && karatsuba_storage_size(s) <= karatsuba_stack_threshold) {	17,908 ✔
146	karatsuba_on_stack(result, a, b);	595 ✔
147	return;	595 ✔
148	}
149	const std::size_t limbs = toom_ladder_storage_size(min_size, s);	17,313 ✔
150	scratch_heap_array<uint_multiprecision_t> buf(heap, limbs);	17,313 ✔
151	scratch_allocator_base scratch(buf.data(), limbs);	17,313 ✔
152	toom_ladder(result, a, b, scratch);	17,313 ✔
153	}	17,313 ✔
154
155	std::size_t multiply_runtime_impl(std::span<uint_multiprecision_t> result,
156	std::span<const uint_multiprecision_t> a,
157	std::span<const uint_multiprecision_t> b,
158	const scratch_heap_source& heap,
159	slice_mode mode);
160
161	// Euclidean slicing: cuts the longer operand into pieces of the shorter one's
162	// length m and accumulates piece * short at each offset. In automatic mode
163	// pieces are cut while the remainder is at least 2m or still enters slicing
164	// (mul_should_slice), so the last piece is below 2m for every ratio; the
165	// forced mode cuts until at most m remains. Full pieces (at least m limbs) run
166	// the FFT when the cost-model gate takes the m x piece product (automatic mode
167	// only; 64-bit limbs) and the Karatsuba/Toom ladder otherwise, over scratch
168	// allocated up front (the FFT workspace on first use, sized for the largest
169	// piece); a real
170	// tail or a piece that trims below m re-dispatches through multiply_runtime_any
171	// (it may slice the other way or take the FFT) and allocates its own scratch,
172	// possibly after piece 0 has been written. On x86 (model off, gate = floor on
173	// min) slicing is only entered when min = m is below the floor, so no piece
174	// takes the FFT and this path is a no-op there. Callers pre-zero the result and
175	// discard it if an allocation throws. Requires trimmed operands,
176	// min >= karatsuba_cutoff and a pre-zeroed result.
177	std::size_t multiply_sliced(const std::span<uint_multiprecision_t> result,	2,214 ✔
178	const std::span<const uint_multiprecision_t> a,
179	const std::span<const uint_multiprecision_t> b,
180	const scratch_heap_source& heap,
181	const bool force) {
182	const auto lng = a.size() >= b.size() ? a : b;	2,214 ✔
183	const auto sht = a.size() >= b.size() ? b : a;	2,214 ✔
184	const std::size_t m = sht.size();	2,214 ✔
185	const std::size_t n = lng.size();	2,214 ✔
186
187	std::size_t pieces_before_last = 0;	2,214 ✔
188	std::size_t rem = n;	2,214 ✔
189	while (force ? rem > m : (rem >= 2 * m \|\| mul_should_slice(m, rem))) {	9,622 ✔
190	++pieces_before_last;	7,408 ✔
191	rem -= m;	7,408 ✔
192	}
193
194	const std::size_t s_max = std::max(m, rem);	2,214 ✔
195	const std::size_t tmp_size = m + s_max;	2,214 ✔
196	const std::size_t ws_size = toom_ladder_storage_size(m, s_max);	2,214 ✔
197
198	scratch_heap_array<uint_multiprecision_t> buf(heap, tmp_size + ws_size);	2,214 ✔
199	const std::span<uint_multiprecision_t> tmp(buf.data(), tmp_size);	2,214 ✔
200	#if defined(BEMAN_BIG_INT_SIMD_MUL)
201	[[maybe_unused]] std::optional<scratch_heap_array<double>> fft_fp_ws;
202	[[maybe_unused]] std::optional<scratch_heap_array<std::uint64_t>> fft_int_ws;
203	#else
204	[[maybe_unused]] std::optional<scratch_heap_array<std::uint64_t>> fft_ws;	2,214 ✔
205	#endif
206
207	for (std::size_t i = 0; i <= pieces_before_last; ++i) {	11,836 ✔
208	const std::size_t off = i * m;	9,622 ✔
209	const std::size_t len = (i < pieces_before_last) ? m : rem;	9,622 ✔
210	const auto piece = lng.subspan(off, len);	9,622 ✔
211	const std::size_t pn = trimmed_size_span(piece);	9,622 ✔
212	if (pn == 1 && piece[0] == 0) {	9,622 ✔
213	continue;	59 ✔
214	}
215	const auto p = piece.first(pn);	9,563 ✔
216	const auto dst = (i == 0 ? result : tmp).first(pn + m);	9,563 ✔
217	if (i != 0) {	9,563 ✔
218	std::ranges::fill(dst, uint_multiprecision_t{0});	7,355 ✔
219	}
220
221	if (pn < m) {	9,563 ✔
222	multiply_runtime_any(dst, p, sht, heap);	1,237 ✔
223	} else {
224	bool piece_done = false;	8,326 ✔
225	// Full pieces (m x pn, pn >= m) consult the FFT gate first, except in the forced mode.
226	if constexpr (width_v<uint_multiprecision_t> == 64) {
227	if (!force && fft_mul_worthwhile(m, pn)) {	8,326 ✔
228	#if defined(BEMAN_BIG_INT_SIMD_MUL)
229	if (!fft_fp_ws) {
230	fft_fp_ws.emplace(heap, fft_mul_fp_storage_size(m, s_max));
231	fft_int_ws.emplace(heap, fft_mul_int_storage_size(m, s_max));
232	}
233	multiply_fft(dst, p, sht, fft_fp_ws->span(), fft_int_ws->span());
234	#else
NEW 235	if (!fft_ws) {	×
NEW 236	fft_ws.emplace(heap, fft_mul_storage_size(m, s_max));	×
237	}
NEW 238	multiply_fft(dst, p, sht, fft_ws->span());	×
239	#endif
NEW 240	piece_done = true;	×
241	}
242	}
243	if (!piece_done) {	8,326 ✔
244	scratch_allocator_base ws(buf.data() + tmp_size, ws_size);	8,326 ✔
245	toom_ladder(dst, p, sht, ws);	8,326 ✔
246	}
247	}
248
249	if (i != 0) {	9,563 ✔
250	add_shifted(result.first(off + len + m), off, dst);	7,355 ✔
251	}
252	}
253	return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), n + m});	4,428 ✔
254	}	2,214 ✔
255
256	std::size_t multiply_runtime_impl(const std::span<uint_multiprecision_t> result,	1,594,531 ✔
257	const std::span<const uint_multiprecision_t> a,
258	const std::span<const uint_multiprecision_t> b,
259	const scratch_heap_source& heap,
260	const slice_mode mode) {
261	BEMAN_BIG_INT_DEBUG_ASSERT(a.size() >= 2);	1,594,531 ✔
262	BEMAN_BIG_INT_DEBUG_ASSERT(b.size() >= 2);	1,594,531 ✔
263	BEMAN_BIG_INT_DEBUG_ASSERT(a.back() != 0);	1,594,531 ✔
264	BEMAN_BIG_INT_DEBUG_ASSERT(b.back() != 0);	1,594,531 ✔
265	BEMAN_BIG_INT_DEBUG_ASSERT(result.size() >= a.size() + b.size());	1,594,531 ✔
266
267	// x * x and x *= x pass the same span twice, so squaring detection is
268	// a pointer compare that almost always fails fast for ordinary mul.
269	if (a.data() == b.data() && a.size() == b.size()) {	1,594,531 ✔
270	return square_runtime(result, a, heap);	949,852 ✔
271	}
272
273	const std::size_t min_size = std::min(a.size(), b.size());	644,679 ✔
274	const std::size_t max_size = std::max(a.size(), b.size());	644,679 ✔
275	const std::size_t result_total = a.size() + b.size();	644,679 ✔
276
277	if (min_size < karatsuba_cutoff) {	644,679 ✔
278	// Schoolbook long multiplication runtime fallback (known to be on the runtime path).
279	::beman_big_int_multiply_long_runtime(result.data(), a.data(), a.size(), b.data(), b.size());	624,408 ✔
280	return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});	624,408 ✔
281	}
282
283	// Power-of-two operands reduce to a shifted copy of the other operand.
284	// This is only worth checking if we're about to do a big number mul anyway.
285	if (is_power_of_two_span(b)) {	20,271 ✔
286	return multiply_power_of_two(result, a, b);	69 ✔
287	}
288	if (is_power_of_two_span(a)) {	20,202 ✔
289	return multiply_power_of_two(result, b, a);	72 ✔
290	}
291
292	// The FFT kernel packs into 64-bit words, so it is gated to 64-bit limbs.
293	if constexpr (width_v<uint_multiprecision_t> == 64) {
294	const bool use_fft = mode != slice_mode::forced && fft_mul_worthwhile(min_size, max_size);	20,130 ✔
295	if (use_fft) {	20,130 ✔
296	#if defined(BEMAN_BIG_INT_SIMD_MUL)
297	scratch_heap_array<double> fp_ws(heap, fft_mul_fp_storage_size(a.size(), b.size()));
298	scratch_heap_array<std::uint64_t> int_ws(heap, fft_mul_int_storage_size(a.size(), b.size()));
299	multiply_fft(result.first(result_total), a, b, fp_ws.span(), int_ws.span());
300	#else
301	scratch_heap_array<std::uint64_t> ws(heap, fft_mul_storage_size(a.size(), b.size()));	8 ✔
302	multiply_fft(result.first(result_total), a, b, ws.span());	8 ✔
303	#endif
304	return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});	8 ✔
305	}	8 ✔
306	}
307
308	const bool slice = mode == slice_mode::forced ? max_size > min_size	20,122 ✔
309	: mode == slice_mode::automatic ? mul_should_slice(min_size, max_size)	19,592 ✔
310	: false;	20,122 ✔
311	if (slice) {	20,122 ✔
312	return multiply_sliced(result, a, b, heap, mode == slice_mode::forced);	2,214 ✔
313	}
314
315	run_ladder(result.first(result_total), a, b, heap);	17,908 ✔
316	return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});	17,908 ✔
317	}
318
319	} // namespace
320
321	std::size_t multiply_runtime(const std::span<uint_multiprecision_t> result,	1,593,883 ✔
322	const std::span<const uint_multiprecision_t> a,
323	const std::span<const uint_multiprecision_t> b,
324	const scratch_heap_source& heap) {
325	return multiply_runtime_impl(result, a, b, heap, slice_mode::automatic);	1,593,883 ✔
326	}
327
328	std::size_t multiply_runtime_sliced(const std::span<uint_multiprecision_t> result,	558 ✔
329	const std::span<const uint_multiprecision_t> a,
330	const std::span<const uint_multiprecision_t> b,
331	const scratch_heap_source& heap) {
332	return multiply_runtime_impl(result, a, b, heap, slice_mode::forced);	558 ✔
333	}
334
335	std::size_t multiply_runtime_unsliced(const std::span<uint_multiprecision_t> result,	90 ✔
336	const std::span<const uint_multiprecision_t> a,
337	const std::span<const uint_multiprecision_t> b,
338	const scratch_heap_source& heap) {
339	return multiply_runtime_impl(result, a, b, heap, slice_mode::disabled);	90 ✔
340	}
341
342	std::size_t multiply_runtime_any(const std::span<uint_multiprecision_t> result,	268,086 ✔
343	const std::span<const uint_multiprecision_t> a_untrimmed,
344	const std::span<const uint_multiprecision_t> b_untrimmed,
345	const scratch_heap_source& heap) {
346	BEMAN_BIG_INT_DEBUG_ASSERT(!a_untrimmed.empty());	268,086 ✔
347	BEMAN_BIG_INT_DEBUG_ASSERT(!b_untrimmed.empty());	268,086 ✔
348	BEMAN_BIG_INT_DEBUG_ASSERT(result.size() >= a_untrimmed.size() + b_untrimmed.size());	268,086 ✔
349
350	const auto a = a_untrimmed.first(trimmed_size_span(a_untrimmed));	268,086 ✔
351	const auto b = b_untrimmed.first(trimmed_size_span(b_untrimmed));	268,086 ✔
352
353	if (a.size() == 1 && b.size() == 1) {	268,086 ✔
354	const auto [lo, hi] = widening_mul(a[0], b[0]);	76,314 ✔
355	result[0] = lo;	76,314 ✔
356	result[1] = hi;	76,314 ✔
357	return hi != 0 ? 2 : 1;	76,314 ✔
358	}
359	if (a.size() == 1) {	191,772 ✔
360	return multiply_single_limb(result, b, a[0]);	51,840 ✔
361	}
362	if (b.size() == 1) {	139,932 ✔
363	return multiply_single_limb(result, a, b[0]);	1,380 ✔
364	}
365
366	return multiply_runtime(result, a, b, heap);	138,552 ✔
367	}
368
369	} // namespace beman::big_int::detail

bemanproject / big_int / 36725334879

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous