27377983558

Committed 11 Jun 2026 09:16PM UTC coverage: 96.536% (-0.2%) from 96.715%

Build # 27377983558

Build Type

push

github

Committed by

web-flow

Commit Message

Division From Burnikel-Ziegler (#247)

Coverage Stats

1122 of 1173 new or added lines in 10 files covered. (95.65%)

42 existing lines in 2 files now uncovered.

4961 of 5139 relevant lines covered (96.54%)

6840874.82 hits per line

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

94.87

/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 <span>

#include <beman/big_int/detail/config.hpp>
#include <beman/big_int/detail/scratch_allocator.hpp>
#include <beman/big_int/detail/span_ops.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 three-pass squaring basecase.
    if (n < square_long_cutoff) {
        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) {
        square_long(result, a);
        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 {
        // n >= square_toom_cook_6_5_cutoff: Toom-6.5 / Toom-8.5, or FFT once
        // it overtakes at square_fft_cutoff. The FFT kernel packs into 64-bit
        // words, so it is gated to 64-bit limbs; on a 32-bit build the branch
        // is discarded and execution falls through to the Toom chain.
        bool used_fft = false;
        if constexpr (width_v<uint_multiprecision_t> == 64) {
            if (n >= square_fft_cutoff) {
#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
                used_fft = true;
            }
        }
        if (!used_fft) {
            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});
}

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) {
    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());
    if (min_size >= karatsuba_cutoff) {
        // 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);
        }

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

        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 (min_size < toom_cook_3_cutoff) {
            const std::size_t storage_size = karatsuba_storage_size(s);
            if (storage_size <= karatsuba_stack_threshold) {
                uint_multiprecision_t  stack_buf[karatsuba_stack_threshold];
                scratch_allocator_base scratch(stack_buf, karatsuba_stack_threshold);
                multiply_karatsuba(result.first(result_total), a, b, scratch);
            } else {
                in_heap_scratch(storage_size, [&](scratch_allocator_base& scratch) {
                    multiply_karatsuba(result.first(result_total), a, b, scratch);
                });
            }
        } else if (min_size < toom_cook_4_cutoff) {
            in_heap_scratch(toom_cook_3_storage_size(s), [&](scratch_allocator_base& scratch) {
                multiply_toom_cook_3(result.first(result_total), a, b, scratch);
            });
        } else if (min_size < toom_cook_6_5_cutoff) {
            in_heap_scratch(toom_cook_4_storage_size(s), [&](scratch_allocator_base& scratch) {
                multiply_toom_cook_4(result.first(result_total), a, b, scratch);
            });
        } else {
            // min_size >= toom_cook_6_5_cutoff: Toom-6.5 / Toom-8.5, or FFT
            // once it overtakes at fft_mul_cutoff. Gated to 64-bit limbs (the
            // FFT kernel packs into 64-bit words); on a 32-bit build the
            // branch is discarded and execution falls to the Toom chain.
            bool used_fft = false;
            if constexpr (width_v<uint_multiprecision_t> == 64) {
                if (min_size >= fft_mul_cutoff) {
#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
                    used_fft = true;
                }
            }
            if (!used_fft) {
                if (min_size < toom_cook_8_5_cutoff) {
                    in_heap_scratch(toom_cook_6_5_storage_size(s), [&](scratch_allocator_base& scratch) {
                        multiply_toom_cook_6_5(result.first(result_total), a, b, scratch);
                    });
                } else {
                    in_heap_scratch(toom_cook_8_5_storage_size(s), [&](scratch_allocator_base& scratch) {
                        multiply_toom_cook_8_5(result.first(result_total), a, b, scratch);
                    });
                }
            }
        }
        return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});
    }

    // Long multiplication fallback.
    multiply_long(result, a, b);
    return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), a.size() + b.size()});
}

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 <span>
10
11	#include <beman/big_int/detail/config.hpp>
12	#include <beman/big_int/detail/scratch_allocator.hpp>
13	#include <beman/big_int/detail/span_ops.hpp>
14
15	// The runtime multiplication tier ladders, compiled once. The header
16	// dispatchers (multiply_dispatch / square_dispatch) keep the constexpr
17	// small-operand and constant-evaluation paths and forward every runtime
18	// multi-limb product here; kernel workspaces come from the type-erased heap
19	// hooks, so a single compiled definition serves every allocator.
20
21	namespace beman::big_int::detail {
22
23	std::size_t square_runtime(const std::span<uint_multiprecision_t> result,	949,851✔
24	const std::span<const uint_multiprecision_t> a,
25	const scratch_heap_source& heap) {
26	BEMAN_BIG_INT_DEBUG_ASSERT(a.size() >= 2);	949,851✔
27	BEMAN_BIG_INT_DEBUG_ASSERT(a.back() != 0);	949,851✔
28	BEMAN_BIG_INT_DEBUG_ASSERT(result.size() >= 2 * a.size());	949,851✔
29	BEMAN_BIG_INT_DEBUG_ASSERT(result.data() != a.data());	949,851✔
30
31	const std::size_t n = a.size();	949,851✔
32	const std::size_t result_total = 2 * n;	949,851✔
33
34	// Tiny squares: plain schoolbook beats the three-pass squaring basecase.
35	if (n < square_long_cutoff) {	949,851✔
36	multiply_long(result.first(result_total), a, a);	3,176✔
37	return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});	3,176✔
38	}
39
40	// (2^k)^2 = 2^(2k): a shifted copy beats any squaring kernel.
41	if (is_power_of_two_span(a)) {	946,675✔
42	return multiply_power_of_two(result, a, a);	139✔
43	}
44
45	if (n < square_karatsuba_cutoff) {	946,536✔
46	square_long(result, a);	946,524✔
47	return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});	946,524✔
48	}
49
50	const auto in_heap_scratch = [&](const std::size_t limbs, auto&& kernel) {	9✔
51	scratch_heap_array<uint_multiprecision_t> buf(heap, limbs);	9✔
52	scratch_allocator_base scratch(buf.data(), limbs);	9✔
53	kernel(scratch);	9✔
54	};	21✔
55
56	if (n < square_toom_cook_3_cutoff) {	12✔
57	in_heap_scratch(karatsuba_storage_size(n), [&](scratch_allocator_base& scratch) {	4✔
58	square_karatsuba(result.first(result_total), a, scratch);	4✔
59	});	4✔
60	} else if (n < square_toom_cook_4_cutoff) {	8✔
61	in_heap_scratch(toom_cook_3_storage_size(n), [&](scratch_allocator_base& scratch) {	2✔
62	square_toom_cook_3(result.first(result_total), a, scratch);	2✔
63	});	2✔
64	} else if (n < square_toom_cook_6_5_cutoff) {	6✔
NEW 65	in_heap_scratch(toom_cook_4_storage_size(n), [&](scratch_allocator_base& scratch) {	×
NEW 66	square_toom_cook_4(result.first(result_total), a, scratch);	×
NEW 67	});	×
68	} else {
69	// n >= square_toom_cook_6_5_cutoff: Toom-6.5 / Toom-8.5, or FFT once
70	// it overtakes at square_fft_cutoff. The FFT kernel packs into 64-bit
71	// words, so it is gated to 64-bit limbs; on a 32-bit build the branch
72	// is discarded and execution falls through to the Toom chain.
73	bool used_fft = false;	6✔
74	if constexpr (width_v<uint_multiprecision_t> == 64) {
75	if (n >= square_fft_cutoff) {	6✔
76	#if defined(BEMAN_BIG_INT_SIMD_MUL)
77	scratch_heap_array<double> fp_ws(heap, square_fft_fp_storage_size(n));
78	scratch_heap_array<std::uint64_t> int_ws(heap, square_fft_int_storage_size(n));
79	square_fft(result.first(result_total), a, fp_ws.span(), int_ws.span());
80	#else
81	scratch_heap_array<std::uint64_t> ws(heap, square_fft_storage_size(n));	3✔
82	square_fft(result.first(result_total), a, ws.span());	3✔
83	#endif
84	used_fft = true;	3✔
85	}	3✔
86	}
87	if (!used_fft) {	6✔
88	if (n < square_toom_cook_8_5_cutoff) {	3✔
89	in_heap_scratch(toom_cook_6_5_storage_size(n), [&](scratch_allocator_base& scratch) {	3✔
90	square_toom_cook_6_5(result.first(result_total), a, scratch);	3✔
91	});	3✔
92	} else {
NEW 93	in_heap_scratch(toom_cook_8_5_storage_size(n), [&](scratch_allocator_base& scratch) {	×
NEW 94	square_toom_cook_8_5(result.first(result_total), a, scratch);	×
NEW 95	});	×
96	}
97	}
98	}
99
100	return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});	12✔
101	}
102
103	std::size_t multiply_runtime(const std::span<uint_multiprecision_t> result,	1,560,445✔
104	const std::span<const uint_multiprecision_t> a,
105	const std::span<const uint_multiprecision_t> b,
106	const scratch_heap_source& heap) {
107	BEMAN_BIG_INT_DEBUG_ASSERT(a.size() >= 2);	1,560,445✔
108	BEMAN_BIG_INT_DEBUG_ASSERT(b.size() >= 2);	1,560,445✔
109	BEMAN_BIG_INT_DEBUG_ASSERT(a.back() != 0);	1,560,445✔
110	BEMAN_BIG_INT_DEBUG_ASSERT(b.back() != 0);	1,560,445✔
111	BEMAN_BIG_INT_DEBUG_ASSERT(result.size() >= a.size() + b.size());	1,560,445✔
112
113	// x * x and x *= x pass the same span twice, so squaring detection is
114	// a pointer compare that almost always fails fast for ordinary mul.
115	if (a.data() == b.data() && a.size() == b.size()) {	1,560,445✔
116	return square_runtime(result, a, heap);	949,851✔
117	}
118
119	const std::size_t min_size = std::min(a.size(), b.size());	610,594✔
120	if (min_size >= karatsuba_cutoff) {	610,594✔
121	// Power-of-two operands reduce to a shifted copy of the other
122	// operand. This is only worth checking if we're about to do a big
123	// number mul anyway.
124	if (is_power_of_two_span(b)) {	6,204✔
125	return multiply_power_of_two(result, a, b);	49✔
126	}
127	if (is_power_of_two_span(a)) {	6,155✔
128	return multiply_power_of_two(result, b, a);	20✔
129	}
130
131	const std::size_t s = std::max(a.size(), b.size());	6,135✔
132	const std::size_t result_total = a.size() + b.size();	6,135✔
133
134	const auto in_heap_scratch = [&](const std::size_t limbs, auto&& kernel) {	5,591✔
135	scratch_heap_array<uint_multiprecision_t> buf(heap, limbs);	5,591✔
136	scratch_allocator_base scratch(buf.data(), limbs);	5,591✔
137	kernel(scratch);	5,591✔
138	};	11,726✔
139
140	if (min_size < toom_cook_3_cutoff) {	6,135✔
141	const std::size_t storage_size = karatsuba_storage_size(s);	5,156✔
142	if (storage_size <= karatsuba_stack_threshold) {	5,156✔
143	uint_multiprecision_t stack_buf[karatsuba_stack_threshold];
144	scratch_allocator_base scratch(stack_buf, karatsuba_stack_threshold);	533✔
145	multiply_karatsuba(result.first(result_total), a, b, scratch);	533✔
146	} else {
147	in_heap_scratch(storage_size, [&](scratch_allocator_base& scratch) {	4,623✔
148	multiply_karatsuba(result.first(result_total), a, b, scratch);	4,623✔
149	});	4,623✔
150	}
151	} else if (min_size < toom_cook_4_cutoff) {	979✔
152	in_heap_scratch(toom_cook_3_storage_size(s), [&](scratch_allocator_base& scratch) {	763✔
153	multiply_toom_cook_3(result.first(result_total), a, b, scratch);	763✔
154	});	763✔
155	} else if (min_size < toom_cook_6_5_cutoff) {	216✔
156	in_heap_scratch(toom_cook_4_storage_size(s), [&](scratch_allocator_base& scratch) {	99✔
157	multiply_toom_cook_4(result.first(result_total), a, b, scratch);	99✔
158	});	99✔
159	} else {
160	// min_size >= toom_cook_6_5_cutoff: Toom-6.5 / Toom-8.5, or FFT
161	// once it overtakes at fft_mul_cutoff. Gated to 64-bit limbs (the
162	// FFT kernel packs into 64-bit words); on a 32-bit build the
163	// branch is discarded and execution falls to the Toom chain.
164	bool used_fft = false;	117✔
165	if constexpr (width_v<uint_multiprecision_t> == 64) {
166	if (min_size >= fft_mul_cutoff) {	117✔
167	#if defined(BEMAN_BIG_INT_SIMD_MUL)
168	scratch_heap_array<double> fp_ws(heap, fft_mul_fp_storage_size(a.size(), b.size()));
169	scratch_heap_array<std::uint64_t> int_ws(heap, fft_mul_int_storage_size(a.size(), b.size()));
170	multiply_fft(result.first(result_total), a, b, fp_ws.span(), int_ws.span());
171	#else
172	scratch_heap_array<std::uint64_t> ws(heap, fft_mul_storage_size(a.size(), b.size()));	11✔
173	multiply_fft(result.first(result_total), a, b, ws.span());	11✔
174	#endif
175	used_fft = true;	11✔
176	}	11✔
177	}
178	if (!used_fft) {	117✔
179	if (min_size < toom_cook_8_5_cutoff) {	106✔
180	in_heap_scratch(toom_cook_6_5_storage_size(s), [&](scratch_allocator_base& scratch) {	95✔
181	multiply_toom_cook_6_5(result.first(result_total), a, b, scratch);	95✔
182	});	95✔
183	} else {
184	in_heap_scratch(toom_cook_8_5_storage_size(s), [&](scratch_allocator_base& scratch) {	11✔
185	multiply_toom_cook_8_5(result.first(result_total), a, b, scratch);	11✔
186	});	11✔
187	}
188	}
189	}
190	return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), result_total});	6,135✔
191	}
192
193	// Long multiplication fallback.
194	multiply_long(result, a, b);	604,390✔
195	return trimmed_size_span(std::span<const uint_multiprecision_t>{result.data(), a.size() + b.size()});	604,390✔
196	}
197
198	std::size_t multiply_runtime_any(const std::span<uint_multiprecision_t> result,	246,962✔
199	const std::span<const uint_multiprecision_t> a_untrimmed,
200	const std::span<const uint_multiprecision_t> b_untrimmed,
201	const scratch_heap_source& heap) {
202	BEMAN_BIG_INT_DEBUG_ASSERT(!a_untrimmed.empty());	246,962✔
203	BEMAN_BIG_INT_DEBUG_ASSERT(!b_untrimmed.empty());	246,962✔
204	BEMAN_BIG_INT_DEBUG_ASSERT(result.size() >= a_untrimmed.size() + b_untrimmed.size());	246,962✔
205
206	const auto a = a_untrimmed.first(trimmed_size_span(a_untrimmed));	246,962✔
207	const auto b = b_untrimmed.first(trimmed_size_span(b_untrimmed));	246,962✔
208
209	if (a.size() == 1 && b.size() == 1) {	246,962✔
210	const auto [lo, hi] = widening_mul(a[0], b[0]);	76,298✔
211	result[0] = lo;	76,298✔
212	result[1] = hi;	76,298✔
213	return hi != 0 ? 2 : 1;	76,298✔
214	}
215	if (a.size() == 1) {	170,664✔
216	return multiply_single_limb(result, b, a[0]);	49,132✔
217	}
218	if (b.size() == 1) {	121,532✔
219	return multiply_single_limb(result, a, b[0]);	919✔
220	}
221
222	return multiply_runtime(result, a, b, heap);	120,613✔
223	}
224
225	} // namespace beman::big_int::detail

eisenwave / std-big-int / 27377983558

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