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

83.87

/src/mulmod_bnm1.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>

// Wraparound multiplication a * b mod (B^w - 1), compiled once; the header
// keeps the constexpr fold helpers, the size chooser, and the storage model.
// `scratch` must carry the type-erased heap hooks (every
// scratch_allocator<Allocator> installs them) for the internal products and
// the cyclic NTT workspaces.

namespace beman::big_int::detail {

namespace {

// r = a * b mod (B^h + 1) with h + 1 = r.size(); a and b canonical in
// [0, B^h] as produced by fold_mod_bnp1. One full h x h product plus a
// signed fold; operands equal to B^h itself (== -1) shortcut to a negation.
void multiply_mod_bnp1(const std::span<uint_multiprecision_t>       r,
                       const std::span<const uint_multiprecision_t> a,
                       const std::span<const uint_multiprecision_t> b,
                       scratch_allocator_base&                      scratch) {
    const std::size_t h = r.size() - 1;
    BEMAN_BIG_INT_DEBUG_ASSERT(a.size() == h + 1);
    BEMAN_BIG_INT_DEBUG_ASSERT(b.size() == h + 1);

    const auto negate_into = [&](const std::span<const uint_multiprecision_t> x) {
        // r = (B^h + 1) - x for x in (0, B^h], r = 0 for x == 0.
        if (is_span_zero(x)) {
            std::ranges::fill(r, uint_multiprecision_t{0});
            return;
        }
        std::ranges::fill(r, uint_multiprecision_t{0});
        r[0] = 1;
        r[h] = 1;
        subtract_unsigned_spans(r, r, x);
    };

    if (a[h] != 0) {
        BEMAN_BIG_INT_DEBUG_ASSERT(is_span_zero(a.first(h)));
        negate_into(b);
        return;
    }
    if (b[h] != 0) {
        BEMAN_BIG_INT_DEBUG_ASSERT(is_span_zero(b.first(h)));
        negate_into(a);
        return;
    }

    const std::span<uint_multiprecision_t> prod = scratch.allocate(2 * h);
    std::ranges::fill(prod, uint_multiprecision_t{0});
    multiply_runtime_any(prod, a.first(h), b.first(h), scratch.heap());
    fold_mod_bnp1(r, prod);
    scratch.deallocate(2 * h);
}

} // namespace

// ---------------------------------------------------------------------------
// r = a * b mod (B^w - 1) with w = r.size(), semi-canonical (all-ones means
// zero). a.size() and b.size() must be at most w (fold larger operands
// first); r must not alias the inputs. `scratch` provides
// multiply_mod_bnm1_storage_size(w) limbs.
// Odd wrap sizes fall back to the plain product (size via
// multiply_mod_bnm1_next_size to keep the recursion even).
// `cutoff_override` is a test-only escape hatch forcing deep recursion.
// ---------------------------------------------------------------------------
void multiply_mod_bnm1(const std::span<uint_multiprecision_t>       r,
                       const std::span<const uint_multiprecision_t> a,
                       const std::span<const uint_multiprecision_t> b,
                       scratch_allocator_base&                      scratch,
                       const std::size_t                            cutoff_override) {
    const std::size_t w = r.size();
    BEMAN_BIG_INT_DEBUG_ASSERT(w >= 1);
    BEMAN_BIG_INT_DEBUG_ASSERT(!a.empty());
    BEMAN_BIG_INT_DEBUG_ASSERT(!b.empty());
    BEMAN_BIG_INT_DEBUG_ASSERT(a.size() <= w);
    BEMAN_BIG_INT_DEBUG_ASSERT(b.size() <= w);
    BEMAN_BIG_INT_DEBUG_ASSERT(r.data() != a.data());
    BEMAN_BIG_INT_DEBUG_ASSERT(r.data() != b.data());

    const std::size_t cutoff = cutoff_override != 0 ? cutoff_override : multiply_mod_bnm1_cutoff;

    // Plain product when the wrap cannot engage (no wraparound, odd size, or
    // too small to be worth the CRT split).
    if (w <= cutoff || (w % 2) != 0 || a.size() + b.size() <= w) {
        const std::size_t                      p_len = a.size() + b.size();
        const std::span<uint_multiprecision_t> prod  = scratch.allocate(p_len);
        std::ranges::fill(prod, uint_multiprecision_t{0});
        multiply_runtime_any(prod, a, b, scratch.heap());
        fold_mod_bnm1(r, std::span<const uint_multiprecision_t>{prod.data(), p_len});
        scratch.deallocate(p_len);
        return;
    }

    // Cyclic NTT tier: one length-L transform set computes the wrapped
    // product directly when w is a chooser size (next_size produces exactly
    // these above the cutoff). Transform workspaces live on the heap like
    // multiply_dispatch's FFT branch, so the scratch model is untouched.
    // The test-only override keeps forcing the CRT recursion.
    if constexpr (width_v<uint_multiprecision_t> == 64) {
        if (cutoff_override == 0 && w >= fft_cyclic_cutoff) {
            const fft_cyclic_params params = multiply_fft_cyclic_next_size(w);
            if (params.wrap_limbs == w) {
                if (is_span_zero(a) || is_span_zero(b)) {
                    std::ranges::fill(r, uint_multiprecision_t{0});
                    return;
                }
#if defined(BEMAN_BIG_INT_SIMD_MUL)
                scratch_heap_array<double>        fp_ws(scratch.heap(), fft_cyclic_fp_storage_size(params));
                scratch_heap_array<std::uint64_t> int_ws(scratch.heap(), fft_cyclic_int_storage_size(params));
                multiply_fft_cyclic(r, a, b, params, fp_ws.span(), int_ws.span());
#else
                scratch_heap_array<std::uint64_t> ws(scratch.heap(), fft_cyclic_storage_size(params));
                multiply_fft_cyclic(r, a, b, params, ws.span());
#endif
                return;
            }
        }
    }

    const std::size_t h = w / 2;

    // Half 1 (recursive): rm1 = a*b mod (B^h - 1), built into r's low half.
    {
        const std::span<uint_multiprecision_t> am1 = scratch.allocate(h);
        const std::span<uint_multiprecision_t> bm1 = scratch.allocate(h);
        fold_mod_bnm1(am1, a);
        fold_mod_bnm1(bm1, b);
        multiply_mod_bnm1(r.first(h),
                          std::span<const uint_multiprecision_t>{am1.data(), h},
                          std::span<const uint_multiprecision_t>{bm1.data(), h},
                          scratch,
                          cutoff_override);
        scratch.deallocate(h);
        scratch.deallocate(h);
    }

    // Half 2: rp1 = a*b mod (B^h + 1).
    const std::span<uint_multiprecision_t> ap1 = scratch.allocate(h + 1);
    const std::span<uint_multiprecision_t> bp1 = scratch.allocate(h + 1);
    const std::span<uint_multiprecision_t> rp1 = scratch.allocate(h + 1);
    fold_mod_bnp1(ap1, a);
    fold_mod_bnp1(bp1, b);
    multiply_mod_bnp1(rp1,
                      std::span<const uint_multiprecision_t>{ap1.data(), h + 1},
                      std::span<const uint_multiprecision_t>{bp1.data(), h + 1},
                      scratch);

    // CRT: r = rm1 + t * (B^h - 1) with t = (rm1 - rp1) / 2 mod (B^h + 1).
    // Reuse ap1's buffer for t.
    const std::span<uint_multiprecision_t> t = ap1;
    {
        std::ranges::copy(r.first(h), t.begin());
        t[h] = 0;
        if (compare_unsigned_spans(t.first(h), rp1) == std::strong_ordering::less) {
            // t = rm1 + (B^h + 1) before the subtraction.
            t[h] = increment_span(t.first(h)) ? 2 : 1;
        }
        subtract_unsigned_spans(t, t, rp1);
        if ((t[0] & 1u) != 0) {
            // Make the value even by adding B^h + 1 once more before halving.
            const bool wrapped = increment_span(t.first(h));
            t[h]               = t[h] + uint_multiprecision_t{1} + uint_multiprecision_t{wrapped};
        }
        const uint_multiprecision_t dropped = shift_right_n(t, 1u);
        BEMAN_BIG_INT_DEBUG_ASSERT(dropped == 0);
        BEMAN_BIG_INT_DEBUG_ASSERT(t[h] <= 1);
    }

    // Assemble in place: r = [rm1 | t_low] (+ 1 if t's top limb carries the
    // B^w == 1 wrap), then a modular subtraction of t.
    std::ranges::copy(t.first(h), r.begin() + static_cast<std::ptrdiff_t>(h));
    if (t[h] != 0) {
        if (increment_span(r)) {
            r[0] = 1;
        }
    }
    if (subtract_unsigned_spans_borrow_out(r, r, std::span<const uint_multiprecision_t>{t.data(), h + 1})) {
        // Wrapped past zero: -B^w == -1 (mod B^w - 1).
        [[maybe_unused]] const bool all_zero = decrement_span(r);
    }

    scratch.deallocate(h + 1);
    scratch.deallocate(h + 1);
    scratch.deallocate(h + 1);
}

} // 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	// Wraparound multiplication a * b mod (B^w - 1), compiled once; the header
16	// keeps the constexpr fold helpers, the size chooser, and the storage model.
17	// `scratch` must carry the type-erased heap hooks (every
18	// scratch_allocator<Allocator> installs them) for the internal products and
19	// the cyclic NTT workspaces.
20
21	namespace beman::big_int::detail {
22
23	namespace {
24
25	// r = a * b mod (B^h + 1) with h + 1 = r.size(); a and b canonical in
26	// [0, B^h] as produced by fold_mod_bnp1. One full h x h product plus a
27	// signed fold; operands equal to B^h itself (== -1) shortcut to a negation.
28	void multiply_mod_bnp1(const std::span<uint_multiprecision_t> r,	4,221✔
29	const std::span<const uint_multiprecision_t> a,
30	const std::span<const uint_multiprecision_t> b,
31	scratch_allocator_base& scratch) {
32	const std::size_t h = r.size() - 1;	4,221✔
33	BEMAN_BIG_INT_DEBUG_ASSERT(a.size() == h + 1);	4,221✔
34	BEMAN_BIG_INT_DEBUG_ASSERT(b.size() == h + 1);	4,221✔
35
NEW 36	const auto negate_into = [&](const std::span<const uint_multiprecision_t> x) {	×
37	// r = (B^h + 1) - x for x in (0, B^h], r = 0 for x == 0.
NEW 38	if (is_span_zero(x)) {	×
NEW 39	std::ranges::fill(r, uint_multiprecision_t{0});	×
NEW 40	return;	×
41	}
NEW 42	std::ranges::fill(r, uint_multiprecision_t{0});	×
NEW 43	r[0] = 1;	×
NEW 44	r[h] = 1;	×
NEW 45	subtract_unsigned_spans(r, r, x);	×
46	};	4,221✔
47
48	if (a[h] != 0) {	4,221✔
NEW 49	BEMAN_BIG_INT_DEBUG_ASSERT(is_span_zero(a.first(h)));	×
NEW 50	negate_into(b);	×
NEW 51	return;	×
52	}
53	if (b[h] != 0) {	4,221✔
NEW 54	BEMAN_BIG_INT_DEBUG_ASSERT(is_span_zero(b.first(h)));	×
NEW 55	negate_into(a);	×
NEW 56	return;	×
57	}
58
59	const std::span<uint_multiprecision_t> prod = scratch.allocate(2 * h);	4,221✔
60	std::ranges::fill(prod, uint_multiprecision_t{0});	4,221✔
61	multiply_runtime_any(prod, a.first(h), b.first(h), scratch.heap());	4,221✔
62	fold_mod_bnp1(r, prod);	4,221✔
63	scratch.deallocate(2 * h);	4,221✔
64	}
65
66	} // namespace
67
68	// ---------------------------------------------------------------------------
69	// r = a * b mod (B^w - 1) with w = r.size(), semi-canonical (all-ones means
70	// zero). a.size() and b.size() must be at most w (fold larger operands
71	// first); r must not alias the inputs. `scratch` provides
72	// multiply_mod_bnm1_storage_size(w) limbs.
73	// Odd wrap sizes fall back to the plain product (size via
74	// multiply_mod_bnm1_next_size to keep the recursion even).
75	// `cutoff_override` is a test-only escape hatch forcing deep recursion.
76	// ---------------------------------------------------------------------------
77	void multiply_mod_bnm1(const std::span<uint_multiprecision_t> r,	10,173✔
78	const std::span<const uint_multiprecision_t> a,
79	const std::span<const uint_multiprecision_t> b,
80	scratch_allocator_base& scratch,
81	const std::size_t cutoff_override) {
82	const std::size_t w = r.size();	10,173✔
83	BEMAN_BIG_INT_DEBUG_ASSERT(w >= 1);	10,173✔
84	BEMAN_BIG_INT_DEBUG_ASSERT(!a.empty());	10,173✔
85	BEMAN_BIG_INT_DEBUG_ASSERT(!b.empty());	10,173✔
86	BEMAN_BIG_INT_DEBUG_ASSERT(a.size() <= w);	10,173✔
87	BEMAN_BIG_INT_DEBUG_ASSERT(b.size() <= w);	10,173✔
88	BEMAN_BIG_INT_DEBUG_ASSERT(r.data() != a.data());	10,173✔
89	BEMAN_BIG_INT_DEBUG_ASSERT(r.data() != b.data());	10,173✔
90
91	const std::size_t cutoff = cutoff_override != 0 ? cutoff_override : multiply_mod_bnm1_cutoff;	10,173✔
92
93	// Plain product when the wrap cannot engage (no wraparound, odd size, or
94	// too small to be worth the CRT split).
95	if (w <= cutoff \|\| (w % 2) != 0 \|\| a.size() + b.size() <= w) {	10,173✔
96	const std::size_t p_len = a.size() + b.size();	5,942✔
97	const std::span<uint_multiprecision_t> prod = scratch.allocate(p_len);	5,942✔
98	std::ranges::fill(prod, uint_multiprecision_t{0});	5,942✔
99	multiply_runtime_any(prod, a, b, scratch.heap());	5,942✔
100	fold_mod_bnm1(r, std::span<const uint_multiprecision_t>{prod.data(), p_len});	5,942✔
101	scratch.deallocate(p_len);	5,942✔
102	return;	5,942✔
103	}
104
105	// Cyclic NTT tier: one length-L transform set computes the wrapped
106	// product directly when w is a chooser size (next_size produces exactly
107	// these above the cutoff). Transform workspaces live on the heap like
108	// multiply_dispatch's FFT branch, so the scratch model is untouched.
109	// The test-only override keeps forcing the CRT recursion.
110	if constexpr (width_v<uint_multiprecision_t> == 64) {
111	if (cutoff_override == 0 && w >= fft_cyclic_cutoff) {	4,231✔
112	const fft_cyclic_params params = multiply_fft_cyclic_next_size(w);	11✔
113	if (params.wrap_limbs == w) {	11✔
114	if (is_span_zero(a) \|\| is_span_zero(b)) {	10✔
115	std::ranges::fill(r, uint_multiprecision_t{0});	1✔
116	return;	1✔
117	}
118	#if defined(BEMAN_BIG_INT_SIMD_MUL)
119	scratch_heap_array<double> fp_ws(scratch.heap(), fft_cyclic_fp_storage_size(params));
120	scratch_heap_array<std::uint64_t> int_ws(scratch.heap(), fft_cyclic_int_storage_size(params));
121	multiply_fft_cyclic(r, a, b, params, fp_ws.span(), int_ws.span());
122	#else
123	scratch_heap_array<std::uint64_t> ws(scratch.heap(), fft_cyclic_storage_size(params));	9✔
124	multiply_fft_cyclic(r, a, b, params, ws.span());	9✔
125	#endif
126	return;	9✔
127	}	9✔
128	}
129	}
130
131	const std::size_t h = w / 2;	4,221✔
132
133	// Half 1 (recursive): rm1 = a*b mod (B^h - 1), built into r's low half.
134	{
135	const std::span<uint_multiprecision_t> am1 = scratch.allocate(h);	4,221✔
136	const std::span<uint_multiprecision_t> bm1 = scratch.allocate(h);	4,221✔
137	fold_mod_bnm1(am1, a);	4,221✔
138	fold_mod_bnm1(bm1, b);	4,221✔
139	multiply_mod_bnm1(r.first(h),	4,221✔
140	std::span<const uint_multiprecision_t>{am1.data(), h},	4,221✔
141	std::span<const uint_multiprecision_t>{bm1.data(), h},	4,221✔
142	scratch,
143	cutoff_override);
144	scratch.deallocate(h);	4,221✔
145	scratch.deallocate(h);	4,221✔
146	}
147
148	// Half 2: rp1 = a*b mod (B^h + 1).
149	const std::span<uint_multiprecision_t> ap1 = scratch.allocate(h + 1);	4,221✔
150	const std::span<uint_multiprecision_t> bp1 = scratch.allocate(h + 1);	4,221✔
151	const std::span<uint_multiprecision_t> rp1 = scratch.allocate(h + 1);	4,221✔
152	fold_mod_bnp1(ap1, a);	4,221✔
153	fold_mod_bnp1(bp1, b);	4,221✔
154	multiply_mod_bnp1(rp1,	4,221✔
155	std::span<const uint_multiprecision_t>{ap1.data(), h + 1},	4,221✔
156	std::span<const uint_multiprecision_t>{bp1.data(), h + 1},	4,221✔
157	scratch);
158
159	// CRT: r = rm1 + t * (B^h - 1) with t = (rm1 - rp1) / 2 mod (B^h + 1).
160	// Reuse ap1's buffer for t.
161	const std::span<uint_multiprecision_t> t = ap1;	4,221✔
162	{
163	std::ranges::copy(r.first(h), t.begin());	4,221✔
164	t[h] = 0;	4,221✔
165	if (compare_unsigned_spans(t.first(h), rp1) == std::strong_ordering::less) {	4,221✔
166	// t = rm1 + (B^h + 1) before the subtraction.
167	t[h] = increment_span(t.first(h)) ? 2 : 1;	1,999✔
168	}
169	subtract_unsigned_spans(t, t, rp1);	4,221✔
170	if ((t[0] & 1u) != 0) {	4,221✔
171	// Make the value even by adding B^h + 1 once more before halving.
172	const bool wrapped = increment_span(t.first(h));	2,141✔
173	t[h] = t[h] + uint_multiprecision_t{1} + uint_multiprecision_t{wrapped};	2,141✔
174	}
175	const uint_multiprecision_t dropped = shift_right_n(t, 1u);	4,221✔
176	BEMAN_BIG_INT_DEBUG_ASSERT(dropped == 0);	4,221✔
177	BEMAN_BIG_INT_DEBUG_ASSERT(t[h] <= 1);	4,221✔
178	}
179
180	// Assemble in place: r = [rm1 \| t_low] (+ 1 if t's top limb carries the
181	// B^w == 1 wrap), then a modular subtraction of t.
182	std::ranges::copy(t.first(h), r.begin() + static_cast<std::ptrdiff_t>(h));	8,442✔
183	if (t[h] != 0) {	4,221✔
184	if (increment_span(r)) {	64✔
NEW 185	r[0] = 1;	×
186	}
187	}
188	if (subtract_unsigned_spans_borrow_out(r, r, std::span<const uint_multiprecision_t>{t.data(), h + 1})) {	4,221✔
189	// Wrapped past zero: -B^w == -1 (mod B^w - 1).
190	[[maybe_unused]] const bool all_zero = decrement_span(r);	6✔
191	}
192
193	scratch.deallocate(h + 1);	4,221✔
194	scratch.deallocate(h + 1);	4,221✔
195	scratch.deallocate(h + 1);	4,221✔
196	}
197
198	} // 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