From 69005cbf81246139359a67e27201595b8db7abf0 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Fri, 7 Aug 2026 18:47:01 +0200 Subject: [PATCH 1/3] fft_small-negacyclic pointwise products in SS fft convolution --- doc/source/fft.rst | 30 ++ doc/source/fft_small.rst | 64 +++ src/fft.h | 26 + src/fft/convolution.c | 329 +++++++++++++ src/fft/convolution_precache.c | 81 +++ src/fft/fft_mfa_truncate_sqrt2_inner.c | 86 +++- src/fft_small.h | 50 ++ src/fft_small/mpn_mul.c | 56 ++- src/fft_small/negacyclic.c | 464 ++++++++++++++++++ src/fft_small/test/main.c | 4 + src/fft_small/test/t-mpn_ctx_to_from_ffts.c | 129 +++++ .../test/t-sd_fft_mpn_mulmod_2expp1.c | 213 ++++++++ src/fmpz_poly.h | 3 + src/fmpz_poly/mul.c | 3 +- src/fmpz_poly/mullow.c | 3 +- src/fmpz_poly/mullow_SS_precache.c | 35 +- src/fmpz_poly/mulmid.c | 3 +- src/fmpz_poly/mulmid_SS.c | 54 +- 18 files changed, 1588 insertions(+), 45 deletions(-) create mode 100644 src/fft_small/negacyclic.c create mode 100644 src/fft_small/test/t-mpn_ctx_to_from_ffts.c create mode 100644 src/fft_small/test/t-sd_fft_mpn_mulmod_2expp1.c diff --git a/doc/source/fft.rst b/doc/source/fft.rst index f726621f96..26d67086c7 100644 --- a/doc/source/fft.rst +++ b/doc/source/fft.rst @@ -578,3 +578,33 @@ FFT Precaching As per ``fft_convolution`` except that it is assumed ``fft_precache`` has been called on ``jj`` with the same parameters. This will then run faster than if ``fft_convolution`` had been run with the original ``jj``. +.. function:: void fft_convolution_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, mp_limb_t ** ii, mp_limb_t ** jj, slong depth, slong limbs, slong trunc, mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** s1, mp_limb_t ** tt) + + As ``fft_convolution_basic``, except that when ``sdctx`` is not + ``NULL`` the pointwise multiplications modulo `2^N + 1` are + performed by ``sd_fft_mpn_mulmod_2expp1`` with the given context instead of + ``fft_mulmod_2expp1``. Passing ``NULL`` reproduces the classic + behavior exactly. The context must have been initialized for + `N` equal to ``limbs * FLINT_BITS``. + +.. function:: void fft_convolution_precache_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, mp_limb_t ** ii, mp_limb_t ** jj, slong depth, slong limbs, slong trunc, mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** s1, mp_limb_t ** tt) + + As ``fft_convolution_precache`` with the same optional negacyclic + pointwise context as ``fft_convolution_sd_fft``. + +.. function:: slong fft_adjust_limbs_sd_fft(slong tight, slong n_outer, slong * m_out) + + Negacyclic-aware ring selection for Schoenhage-Strassen + convolutions. The instantiated digit sizes of the negacyclic + pointwise engine form a lattice of rings `N' = mb` with `m` a + power of two; among candidates whose growth over + ``fft_adjust_limbs(tight)`` stays within a per-prime-count cap, + the one minimizing an estimated cost is selected. On success + ``*m_out`` receives the digit count and the returned limb count + defines `N'`; otherwise ``*m_out`` is zero and the standard + Nussbaumer rounding is returned, so the caller keeps + ``fft_mulmod_2expp1``. The caps and cost weights are tuning + parameters, calibrated conservatively; on targets with tuned + assembly the measured advantage of the engine is larger and + admits looser caps. + diff --git a/doc/source/fft_small.rst b/doc/source/fft_small.rst index 3d0aa60d9f..ceb488658c 100644 --- a/doc/source/fft_small.rst +++ b/doc/source/fft_small.rst @@ -152,3 +152,67 @@ Transform plans and operands than exhaust memory; the matrix multiplication drivers respond by multiplying in blocks instead. Mutable for tuning; the default is 4 GiB on 64-bit machines. + + +Negacyclic pointwise multiplication +------------------------------------------------------------------------------- + +These functions multiply residues modulo `2^N + 1` by splitting the +operands into `m` digits of `b = N/m` bits, so that the product is a +length-`m` negacyclic convolution of the digit sequences, computed via +the small-prime FFT machinery of this module and Chinese remaindering. +They serve as the pointwise stage of Schoenhage-Strassen convolutions: +see :func:`fft_convolution_sd_fft` and the ``fmpz_poly`` multiplication +functions, which engage this engine for large coefficient sizes. + +Operands and results occupy `n + 1 = N/\mathtt{FLINT\_BITS} + 1` limbs. +An *operand* may have top limb `0` or `1` with arbitrary remaining +limbs; its value is interpreted modulo `2^N + 1` (in particular the +representation need not be reduced, and both representatives of a +residue with top limb set are accepted). The *result* is returned +fully reduced: the unique representative in `[0, 2^N]`, with top limb +`1` exactly when the value is `2^N` and the remaining limbs then +zero. + +.. function:: void sd_fft_mpn_mulmod_2expp1_ctx_init(sd_fft_mpn_mulmod_2expp1_ctx_struct * C, mpn_ctx_struct * R, slong N, slong m) + + Initializes a context for multiplication modulo `2^N + 1` with + digit count *m*, requiring `N = mb` for an integer digit size + *b*, *m* a power of two at least 16, and *N* a multiple of + ``FLINT_BITS``. Passing `m = 0` chooses the digit count + internally by an estimated-cost comparison over the instantiated + digit sizes dividing *N*. The number of FFT primes is chosen as + the smallest count, between 3 and 8, whose exact prime product + exceeds the coefficient bound `2^{2b + \log_2 m + 1}` and for + which a digit-splitting input pass is instantiated. The primes + and transform contexts are taken from *R*, normally + ``get_default_mpn_ctx()``. + + An exception is raised when no supported configuration exists. + With `m = 0` and the default context, initialization succeeds + precisely when `N = 2^k b` for some `k \geq 4` and an + instantiated digit size + `b \in \{64, 84, 88, 92, 112, 116, 120, 126, 128, 136, 140, + 144, 160, 164, 168, 184, 188, 192\}`, subject to prime + capacity: unlimited in practice (`k \geq 40` supported) for + `b \leq 168`, and `k \leq 28`, `20`, `12` for + `b = 184`, `188`, `192` respectively, since those digit sizes + admit no larger prime count to grow into. In particular every + power of two `N \geq 1024` is supported through `b = 64` and + `b = 128`. + +.. function:: void sd_fft_mpn_mulmod_2expp1(sd_fft_mpn_mulmod_2expp1_ctx_struct * C, nn_ptr z, nn_srcptr x, nn_srcptr y, sd_fft_mpn_mulmod_2expp1_scratch_struct * S) + void sd_fft_mpn_mulmod_2expp1_scratch_init(sd_fft_mpn_mulmod_2expp1_scratch_struct * S, const sd_fft_mpn_mulmod_2expp1_ctx_struct * C) + void sd_fft_mpn_mulmod_2expp1_scratch_clear(sd_fft_mpn_mulmod_2expp1_scratch_struct * S) + + The context is read-only during multiplication; all mutable state + lives in a scratch object, so a single context may serve several + threads concurrently, each with its own scratch. + The pointwise stage of the + ``fmpz_poly`` Schoenhage-Strassen convolutions multiplies its + coefficients in parallel over the available threads this way. + +.. function:: void sd_fft_mpn_mulmod_2expp1_ctx_clear(sd_fft_mpn_mulmod_2expp1_ctx_struct * C) + + Releases the memory held by the context. + diff --git a/src/fft.h b/src/fft.h index 09959b9f14..fb11c90180 100644 --- a/src/fft.h +++ b/src/fft.h @@ -183,6 +183,12 @@ void fft_mfa_truncate_sqrt2_inner(mp_limb_t ** ii, mp_limb_t ** jj, mp_size_t n, flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** FLINT_UNUSED(temp), mp_size_t n1, mp_size_t trunc, mp_limb_t ** tt); +struct sd_fft_mpn_mulmod_2expp1_ctx_struct; +void fft_mfa_truncate_sqrt2_inner_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, mp_size_t n, + flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** s1, + mp_size_t sqrt, mp_size_t trunc, mp_limb_t ** tt); + void ifft_mfa_truncate_sqrt2_outer(mp_limb_t ** ii, mp_size_t n, flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** temp, mp_size_t n1, mp_size_t trunc); @@ -211,6 +217,22 @@ void fft_convolution_basic(mp_limb_t ** ii, mp_limb_t ** jj, slong depth, slong limbs, slong trunc, mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** s1, mp_limb_t ** tt); +struct sd_fft_mpn_mulmod_2expp1_ctx_struct; +#define FLINT_FFT_SD_FFT_MIN_RING 16384 /* bits */ +slong fft_adjust_limbs_sd_fft(slong tight, slong n_outer, slong * m_out); +void _fft_pointwise_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, slong j0, slong j1, + slong limbs); + +void _fft_pointwise_rows_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, slong base, + slong srt, slong nrows, slong rvbits, slong limbs); + +void fft_convolution_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, slong depth, + slong limbs, slong trunc, mp_limb_t ** t1, + mp_limb_t ** t2, mp_limb_t ** s1, mp_limb_t ** tt); + void fft_convolution(mp_limb_t ** ii, mp_limb_t ** jj, slong depth, slong limbs, slong trunc, mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** s1, mp_limb_t ** tt); @@ -223,6 +245,10 @@ void fft_precache(mp_limb_t ** jj, slong depth, slong limbs, void fft_convolution_precache(mp_limb_t ** ii, mp_limb_t ** jj, slong depth, slong limbs, slong trunc, mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** s1, mp_limb_t ** tt); +void fft_convolution_precache_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, slong depth, + slong limbs, slong trunc, mp_limb_t ** t1, + mp_limb_t ** t2, mp_limb_t ** s1, mp_limb_t ** tt); #ifdef __cplusplus } diff --git a/src/fft/convolution.c b/src/fft/convolution.c index 4eddd13ba5..e45beb6503 100644 --- a/src/fft/convolution.c +++ b/src/fft/convolution.c @@ -1,5 +1,8 @@ /* Copyright (C) 2008-2011, 2020 William Hart + Copyright (C) 2026 Fredrik Johansson + + Developed using Claude Fable 5 This file is part of FLINT. @@ -11,6 +14,9 @@ #include "ulong_extras.h" #include "fft.h" +#include "fft_small.h" +#include "thread_pool.h" +#include "thread_support.h" void fft_convolution(mp_limb_t ** ii, mp_limb_t ** jj, slong depth, slong limbs, slong trunc, mp_limb_t ** t1, @@ -134,3 +140,326 @@ void fft_convolution_basic(mp_limb_t ** ii, mp_limb_t ** jj, slong depth, } } } + +#if FLINT_HAVE_FFT_SMALL +/* Threaded pointwise stage of the negacyclic convolutions. The + coefficient work is embarrassingly parallel: the context is + read-only during multiplication and each worker gets its own + scratch. Two range shapes are needed: a contiguous [j0, j1), and + the truncated-MFA row set j = base + n_revbin(s, rvbits)*srt + u + for s in [s0, s1), u in [0, srt). */ +typedef struct +{ + sd_fft_mpn_mulmod_2expp1_ctx_struct * C; + mp_limb_t ** ii; + mp_limb_t ** jj; + slong j0; + slong j1; + slong base; + slong srt; + slong rvbits; + slong limbs; + sd_fft_mpn_mulmod_2expp1_scratch_struct * scr; + int rows; +} +_sd_fft_pw_arg_struct; + +static void +_sd_fft_pw_one(_sd_fft_pw_arg_struct * X, slong j) +{ + mpn_normmod_2expp1(X->ii[j], X->limbs); + if (X->ii != X->jj) + mpn_normmod_2expp1(X->jj[j], X->limbs); + sd_fft_mpn_mulmod_2expp1(X->C, X->ii[j], X->ii[j], X->jj[j], X->scr); +} + +static void +_sd_fft_pw_worker(void * varg) +{ + _sd_fft_pw_arg_struct * X = (_sd_fft_pw_arg_struct *) varg; + slong s, u, j; + + if (!X->rows) + { + for (j = X->j0; j < X->j1; j++) + _sd_fft_pw_one(X, j); + } + else + { + for (s = X->j0; s < X->j1; s++) + { + slong t = (slong) n_revbin((ulong) s, (ulong) X->rvbits); + for (u = 0; u < X->srt; u++) + _sd_fft_pw_one(X, X->base + t*X->srt + u); + } + } +} + +static void +_sd_fft_pw_run(_sd_fft_pw_arg_struct * proto, slong count) +{ + thread_pool_handle * handles; + slong nhandles, nthreads, i; + _sd_fft_pw_arg_struct * args; + sd_fft_mpn_mulmod_2expp1_scratch_struct * scr; + slong j0 = proto->j0; + + if (count <= 0) + return; + + nhandles = flint_request_threads(&handles, FLINT_MIN(count, + flint_get_num_threads())); + nthreads = nhandles + 1; + + if (nhandles == 0) + { + sd_fft_mpn_mulmod_2expp1_scratch_t S; + sd_fft_mpn_mulmod_2expp1_scratch_init(S, proto->C); + proto->scr = S; + _sd_fft_pw_worker(proto); + sd_fft_mpn_mulmod_2expp1_scratch_clear(S); + flint_give_back_threads(handles, nhandles); + return; + } + + args = flint_malloc(nthreads * sizeof(_sd_fft_pw_arg_struct)); + scr = flint_malloc(nthreads + * sizeof(sd_fft_mpn_mulmod_2expp1_scratch_struct)); + for (i = 0; i < nthreads; i++) + sd_fft_mpn_mulmod_2expp1_scratch_init(scr + i, proto->C); + + for (i = 0; i < nthreads; i++) + { + args[i] = *proto; + args[i].j0 = j0 + (i * count) / nthreads; + args[i].j1 = j0 + ((i + 1) * count) / nthreads; + args[i].scr = scr + i; + } + + for (i = nhandles; i > 0; i--) + thread_pool_wake(global_thread_pool, handles[i - 1], 0, + _sd_fft_pw_worker, args + i); + _sd_fft_pw_worker(args + 0); + for (i = nhandles; i > 0; i--) + thread_pool_wait(global_thread_pool, handles[i - 1]); + + for (i = 0; i < nthreads; i++) + sd_fft_mpn_mulmod_2expp1_scratch_clear(scr + i); + flint_free(scr); + flint_free(args); + flint_give_back_threads(handles, nhandles); +} + +void +_fft_pointwise_sd_fft(sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, mp_limb_t ** ii, + mp_limb_t ** jj, slong j0, slong j1, slong limbs) +{ + _sd_fft_pw_arg_struct X = { sdctx, ii, jj, j0, j1, + 0, 0, 0, limbs, NULL, 0 }; + _sd_fft_pw_run(&X, j1 - j0); +} + +void +_fft_pointwise_rows_sd_fft(sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, slong base, + slong srt, slong nrows, slong rvbits, slong limbs) +{ + _sd_fft_pw_arg_struct X = { sdctx, ii, jj, 0, nrows, + base, srt, rvbits, limbs, NULL, 1 }; + _sd_fft_pw_run(&X, nrows); +} +#else +/* without the fft_small engine no context can exist; these are + never reached and exist only to satisfy references */ +void +_fft_pointwise_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, slong j0, + slong j1, slong limbs) +{ + FLINT_ASSERT(sdctx == NULL); + (void) sdctx; (void) ii; (void) jj; + (void) j0; (void) j1; (void) limbs; +} + +void +_fft_pointwise_rows_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, slong base, + slong srt, slong nrows, slong rvbits, slong limbs) +{ + FLINT_ASSERT(sdctx == NULL); + (void) sdctx; (void) ii; (void) jj; (void) base; + (void) srt; (void) nrows; (void) rvbits; (void) limbs; +} +#endif + +/* fft_convolution_basic with an optional negacyclic pointwise + engine; sdctx == NULL reproduces the classic behavior */ +void fft_convolution_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, slong depth, + slong limbs, slong trunc, mp_limb_t ** t1, + mp_limb_t ** t2, mp_limb_t ** s1, mp_limb_t ** tt) +{ + slong n = (WORD(1)< cap[np - 3]) + continue; + cost = (double) np * (double) m + * ((double) depth + 4.0 + (double) b / 16.0); + if (bestm == 0 || cost < bestc) + { + bestc = cost; + bestl = (m * b) / FLINT_BITS; + bestm = m; + } + } + *m_out = bestm; + return bestl; +#endif +} diff --git a/src/fft/convolution_precache.c b/src/fft/convolution_precache.c index 8269c2ac37..e18ea77179 100644 --- a/src/fft/convolution_precache.c +++ b/src/fft/convolution_precache.c @@ -1,5 +1,8 @@ /* Copyright (C) 2008-2011, 2020 William Hart + Copyright (C) 2026 Fredrik Johansson + + Developed using Claude Fable 5 This file is part of FLINT. @@ -11,6 +14,8 @@ #include "ulong_extras.h" #include "fft.h" +#include "fft_small.h" + void fft_precache(mp_limb_t ** jj, slong depth, slong limbs, slong trunc, mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** s1) @@ -117,3 +122,79 @@ void fft_convolution_precache(mp_limb_t ** ii, mp_limb_t ** jj, slong depth, } } } + +/* fft_convolution_precache with an optional negacyclic pointwise + engine; sdctx == NULL reproduces the classic behavior */ +void fft_convolution_precache_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, slong depth, + slong limbs, slong trunc, mp_limb_t ** t1, + mp_limb_t ** t2, mp_limb_t ** s1, mp_limb_t ** tt) +{ + slong n = (WORD(1)<. */ +#include "fft_small.h" #include "thread_pool.h" + +struct sd_fft_mpn_mulmod_2expp1_ctx_struct; +struct sd_fft_mpn_mulmod_2expp1_scratch_struct; #include "thread_support.h" #include "ulong_extras.h" #include "fft.h" @@ -29,6 +36,8 @@ typedef struct mp_limb_t ** t1; mp_limb_t ** t2; mp_limb_t * tt; + struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx; + struct sd_fft_mpn_mulmod_2expp1_scratch_struct * nscr; #if FLINT_USES_PTHREAD pthread_mutex_t * mutex; #endif @@ -78,7 +87,13 @@ _fft_inner1_worker(void * arg_ptr) mp_size_t t = i*n1 + j; mpn_normmod_2expp1(ii[t], limbs); if (ii != jj) mpn_normmod_2expp1(jj[t], limbs); - fft_mulmod_2expp1(ii[t], ii[t], jj[t], n, w, tt); + #if FLINT_HAVE_FFT_SMALL + if (arg.sdctx != NULL) + sd_fft_mpn_mulmod_2expp1(arg.sdctx, ii[t], ii[t], jj[t], + arg.nscr); + else +#endif + fft_mulmod_2expp1(ii[t], ii[t], jj[t], n, w, tt); } ifft_radix2(ii + i*n1, n1/2, w*n2, t1, t2); @@ -126,7 +141,13 @@ _fft_inner2_worker(void * arg_ptr) mp_size_t t = i*n1 + j; mpn_normmod_2expp1(ii[t], limbs); if (ii != jj) mpn_normmod_2expp1(jj[t], limbs); - fft_mulmod_2expp1(ii[t], ii[t], jj[t], n, w, tt); + #if FLINT_HAVE_FFT_SMALL + if (arg.sdctx != NULL) + sd_fft_mpn_mulmod_2expp1(arg.sdctx, ii[t], ii[t], jj[t], + arg.nscr); + else +#endif + fft_mulmod_2expp1(ii[t], ii[t], jj[t], n, w, tt); } ifft_radix2(ii + i*n1, n1/2, w*n2, t1, t2); @@ -134,7 +155,9 @@ _fft_inner2_worker(void * arg_ptr) } } -void fft_mfa_truncate_sqrt2_inner(mp_limb_t ** ii, mp_limb_t ** jj, mp_size_t n, +static void +_fft_mfa_truncate_sqrt2_inner_impl(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, mp_size_t n, flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** FLINT_UNUSED(temp), mp_size_t n1, mp_size_t trunc, mp_limb_t ** tt) { @@ -161,8 +184,17 @@ void fft_mfa_truncate_sqrt2_inner(mp_limb_t ** ii, mp_limb_t ** jj, mp_size_t n, /* convolutions on relevant rows */ + /* Phase 1 has 2n/n1 columns and phase 2 has trunc2 rows of + work-stolen units; each unit is heavy at large limb counts + (n1 column transforms plus n1 pointwise multiplications), so + for the negacyclic engine one unit per thread is worthwhile. + The classic cap is kept unchanged. */ num_threads = flint_request_threads(&threads, - FLINT_MIN(flint_get_num_threads(), (trunc2 + 15)/16)); + sdctx != NULL + ? FLINT_MIN(flint_get_num_threads(), + FLINT_MAX((2*n)/n1, trunc2)) + : FLINT_MIN(flint_get_num_threads(), + (trunc2 + 15)/16)); args = (fft_inner_arg_t *) flint_malloc(sizeof(fft_inner_arg_t)*(num_threads + 1)); @@ -182,11 +214,25 @@ void fft_mfa_truncate_sqrt2_inner(mp_limb_t ** ii, mp_limb_t ** jj, mp_size_t n, args[i].t1 = t1 + i; args[i].t2 = t2 + i; args[i].tt = tt[i]; + args[i].sdctx = sdctx; + args[i].nscr = NULL; #if FLINT_USES_PTHREAD args[i].mutex = &mutex; #endif } +#if FLINT_HAVE_FFT_SMALL + if (sdctx != NULL) + { + for (i = 0; i < num_threads + 1; i++) + { + args[i].nscr = flint_malloc( + sizeof(sd_fft_mpn_mulmod_2expp1_scratch_struct)); + sd_fft_mpn_mulmod_2expp1_scratch_init(args[i].nscr, sdctx); + } + } +#endif + for (i = 0; i < num_threads; i++) thread_pool_wake(global_thread_pool, threads[i], 0, _fft_inner1_worker, &args[i]); @@ -220,9 +266,41 @@ void fft_mfa_truncate_sqrt2_inner(mp_limb_t ** ii, mp_limb_t ** jj, mp_size_t n, flint_give_back_threads(threads, num_threads); +#if FLINT_HAVE_FFT_SMALL + if (sdctx != NULL) + { + for (i = 0; i < num_threads + 1; i++) + { + sd_fft_mpn_mulmod_2expp1_scratch_clear(args[i].nscr); + flint_free(args[i].nscr); + } + } +#endif flint_free(args); #if FLINT_USES_PTHREAD pthread_mutex_destroy(&mutex); #endif } + +void fft_mfa_truncate_sqrt2_inner(mp_limb_t ** ii, mp_limb_t ** jj, mp_size_t n, + flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** s1, + mp_size_t n1, mp_size_t trunc, mp_limb_t ** tt) +{ + _fft_mfa_truncate_sqrt2_inner_impl(NULL, ii, jj, n, w, t1, t2, s1, + n1, trunc, tt); +} + +#if FLINT_HAVE_FFT_SMALL +/* as fft_mfa_truncate_sqrt2_inner with the pointwise multiplications + performed by the negacyclic engine; each worker gets its own + scratch, the context being read-only during multiplication */ +void fft_mfa_truncate_sqrt2_inner_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, mp_size_t n, + flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** s1, + mp_size_t n1, mp_size_t trunc, mp_limb_t ** tt) +{ + _fft_mfa_truncate_sqrt2_inner_impl(sdctx, ii, jj, n, w, t1, t2, s1, + n1, trunc, tt); +} +#endif diff --git a/src/fft_small.h b/src/fft_small.h index a48ef9da07..e43ec92537 100644 --- a/src/fft_small.h +++ b/src/fft_small.h @@ -788,6 +788,56 @@ int _fmpz_poly_mul_mid_default_mpn_ctx( const fmpz * a, slong an, const fmpz * b, slong bn); + +/* negacyclic pointwise engine (negacyclic.c) */ +typedef void (*from_ffts_func)( + ulong* z, ulong lo, ulong hi, ulong c_lo, ulong clen, + sd_fft_ctx_struct* Rffts, double* d, ulong dstride, + crt_data_struct* Rcrts, + ulong bits, + ulong start_easy, ulong stop_easy, + ulong* overhang, ulong* boundbuf); + +to_ffts_func _mpn_ctx_to_ffts_func(ulong np, ulong bits); +from_ffts_func _mpn_ctx_from_ffts_func(ulong np); + +typedef struct sd_fft_mpn_mulmod_2expp1_scratch_struct +{ + double * buf; /* transform buffers, 2 np data blocks */ + nn_ptr w; /* recomposition window */ + double * digs; /* digit doubles for b <= 50, else NULL */ +} sd_fft_mpn_mulmod_2expp1_scratch_struct; + +typedef sd_fft_mpn_mulmod_2expp1_scratch_struct sd_fft_mpn_mulmod_2expp1_scratch_t[1]; + +typedef struct sd_fft_mpn_mulmod_2expp1_ctx_struct +{ + slong N, b, m, depth, np; + nn_ptr biaspoly; /* precomputed bias polynomial (np == 4) */ + slong whi; /* recomposition limb count (np == 4) */ + ulong biasres[8]; /* bias mod p_k */ + ulong p[8]; + nmod_t mod[8]; + double * tw[8]; /* w^i, twisted input scale */ + double * itw[8]; /* w^-i * m^-1, output scale */ + slong * crt2; /* scratch: two-prime corrections */ + fmpz_t Pbig, Phalf; + mpn_ctx_struct * R; + slong dsz; +} sd_fft_mpn_mulmod_2expp1_ctx_struct; + +typedef sd_fft_mpn_mulmod_2expp1_ctx_struct sd_fft_mpn_mulmod_2expp1_ctx_t[1]; + +slong sd_fft_mpn_mulmod_2expp1_choose_np(mpn_ctx_struct * R, slong b, slong depth); +void sd_fft_mpn_mulmod_2expp1_ctx_init(sd_fft_mpn_mulmod_2expp1_ctx_struct * C, + mpn_ctx_struct * R, slong N, slong m); +void sd_fft_mpn_mulmod_2expp1_ctx_clear(sd_fft_mpn_mulmod_2expp1_ctx_struct * C); +void sd_fft_mpn_mulmod_2expp1(sd_fft_mpn_mulmod_2expp1_ctx_struct * C, nn_ptr z, nn_srcptr x, + nn_srcptr y, sd_fft_mpn_mulmod_2expp1_scratch_struct * S); +void sd_fft_mpn_mulmod_2expp1_scratch_init(sd_fft_mpn_mulmod_2expp1_scratch_struct * S, + const sd_fft_mpn_mulmod_2expp1_ctx_struct * C); +void sd_fft_mpn_mulmod_2expp1_scratch_clear(sd_fft_mpn_mulmod_2expp1_scratch_struct * S); + #ifdef __cplusplus } #endif diff --git a/src/fft_small/mpn_mul.c b/src/fft_small/mpn_mul.c index e89c5a6e15..14622f6e95 100644 --- a/src/fft_small/mpn_mul.c +++ b/src/fft_small/mpn_mul.c @@ -1,6 +1,9 @@ /* Copyright (C) 2022 Daniel Schultz Copyright (C) 2024 Fredrik Johansson + Copyright (C) 2026 Fredrik Johansson + + Developed using Claude Fable 5 This file is part of FLINT. @@ -543,6 +546,7 @@ FLINT_STATIC_NOINLINE void CAT(mpn_to_ffts_hard, NP)( \ (d + l*dstride)[i] = 0.0; \ } +DEFINE_IT(3) DEFINE_IT(4) DEFINE_IT(5) DEFINE_IT(6) @@ -656,18 +660,36 @@ static void CAT3(mpn_to_ffts, NP, BITS)( \ start_hard, stop_hard, bits); \ } +DEFINE_IT(3, 64) +DEFINE_IT(4, 64) DEFINE_IT(4, 84) DEFINE_IT(4, 88) DEFINE_IT(4, 92) +DEFINE_IT(5, 84) +DEFINE_IT(5, 88) +DEFINE_IT(5, 92) DEFINE_IT(5,112) DEFINE_IT(5,116) DEFINE_IT(5,120) +DEFINE_IT(6,112) +DEFINE_IT(6,116) +DEFINE_IT(6,120) +DEFINE_IT(6,126) +DEFINE_IT(6,128) DEFINE_IT(6,136) DEFINE_IT(6,140) DEFINE_IT(6,144) +DEFINE_IT(7,126) +DEFINE_IT(7,128) +DEFINE_IT(7,136) +DEFINE_IT(7,140) +DEFINE_IT(7,144) DEFINE_IT(7,160) DEFINE_IT(7,164) DEFINE_IT(7,168) +DEFINE_IT(8,160) +DEFINE_IT(8,164) +DEFINE_IT(8,168) DEFINE_IT(8,184) DEFINE_IT(8,188) DEFINE_IT(8,192) @@ -729,13 +751,7 @@ DEFINE_IT(6, 7) DEFINE_IT(7, 8) #undef DEFINE_IT -typedef void (*from_ffts_func)( - ulong* z, ulong lo, ulong hi, ulong c_lo, ulong clen, - sd_fft_ctx_struct* Rffts, double* d, ulong dstride, - crt_data_struct* Rcrts, - ulong bits, - ulong start_easy, ulong stop_easy, - ulong* overhang, ulong* boundbuf); +/* from_ffts_func typedef now in fft_small.h */ /* The "n" here is the limb count Rcrts[np-1].coeff_len, which is big enough @@ -2194,3 +2210,29 @@ void fft_small_export_mpn(ulong* z, ulong zn, const fft_small_op_t X, flint_give_back_threads(handles, nworkers); } + +/* exported dispatch for the negacyclic engine (negacyclic.c) */ +to_ffts_func _mpn_ctx_to_ffts_func(ulong np, ulong bits) +{ +#define P(NP, B) if (np == NP && bits == B) return CAT3(mpn_to_ffts, NP, B); + P(3, 64) P(4, 64) P(4, 84) P(4, 88) P(4, 92) + P(5, 84) P(5, 88) P(5, 92) + P(5, 112) P(5, 116) P(5, 120) + P(6, 112) P(6, 116) P(6, 120) + P(6, 126) P(6, 128) P(6, 136) P(6, 140) P(6, 144) + P(7, 126) P(7, 128) P(7, 136) P(7, 140) P(7, 144) + P(7, 160) P(7, 164) P(7, 168) + P(8, 160) P(8, 164) P(8, 168) + P(8, 184) P(8, 188) P(8, 192) +#undef P + return NULL; +} + +from_ffts_func _mpn_ctx_from_ffts_func(ulong np) +{ + static from_ffts_func tab[8-3+1] = {_mpn_from_ffts_3, + _mpn_from_ffts_4, _mpn_from_ffts_5, _mpn_from_ffts_6, + _mpn_from_ffts_7, _mpn_from_ffts_8}; + FLINT_ASSERT(3 <= np && np <= 8); + return tab[np - 3]; +} diff --git a/src/fft_small/negacyclic.c b/src/fft_small/negacyclic.c new file mode 100644 index 0000000000..e8b846a363 --- /dev/null +++ b/src/fft_small/negacyclic.c @@ -0,0 +1,464 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + Developed using Claude Fable 5 + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +/* + Negacyclic pointwise multiplication: z = x*y mod 2^N+1 on + top-limb-0-or-1 (n+1)-limb operands, via CRT over sd_fft primes + with b-bit digits, N = m*b, m a power of two. Pointwise stage + for Schoenhage-Strassen convolutions (fmpz_poly mul_SS). +*/ + +#include +#include "flint.h" +#include "ulong_extras.h" +#include "nmod.h" +#include "fmpz.h" +#include "mpn_extras.h" +#include "fft_small.h" +#include "crt_helpers.h" + +#define LG_BLK 7 /* LG_BLK_SZ */ + + + +/* Smallest np whose exact prime product covers the biased + coefficient bound 2^(2b + depth + 1), among those with an + instantiated input pass. Least np is fastest: every stage + (input, transforms, pointwise, crt) scales with np. */ +static slong _mulmod_2expp1_choose_m(mpn_ctx_struct * R, + slong N); + +slong +sd_fft_mpn_mulmod_2expp1_choose_np(mpn_ctx_struct * R, slong b, slong depth) +{ + fmpz_t P; + slong np, need = 2 * b + depth + 1; + + fmpz_init_set_ui(P, R->ffts[0].p); + fmpz_mul_ui(P, P, R->ffts[1].p); + for (np = 3; np <= 8; np++) + { + fmpz_mul_ui(P, P, R->ffts[np - 1].p); + if ((slong) fmpz_bits(P) - 1 >= need + && (b <= 50 + || _mpn_ctx_to_ffts_func((ulong) np, (ulong) b) != NULL)) + { + fmpz_clear(P); + return np; + } + } + fmpz_clear(P); + return -1; +} + +void +sd_fft_mpn_mulmod_2expp1_ctx_init(sd_fft_mpn_mulmod_2expp1_ctx_struct * C, mpn_ctx_struct * R, + slong N, slong m) +{ + slong np; + + if (m == 0) + { + m = _mulmod_2expp1_choose_m(R, N); + if (m == 0) + flint_throw(FLINT_ERROR, + "sd_fft_mpn_mulmod_2expp1_ctx_init: no supported " + "digit splitting for N = %wd\n", N); + } + if (N % m != 0 || m < 16 || (m & (m - 1)) != 0) + flint_throw(FLINT_ERROR, "sd_fft_mpn_mulmod_2expp1_ctx_init: need N = m*b " + "with m a power of two >= 16\n"); + np = sd_fft_mpn_mulmod_2expp1_choose_np(R, N / m, FLINT_BIT_COUNT(m) - 1); + if (np < 0) + flint_throw(FLINT_ERROR, "sd_fft_mpn_mulmod_2expp1_ctx_init: no supported " + "prime configuration for b = %wd, m = %wd\n", + N / m, m); + C->R = R; + C->N = N; + C->m = m; + C->b = N / m; + C->np = np; + C->depth = FLINT_BIT_COUNT(m) - 1; + C->dsz = (slong) sd_fft_ctx_data_size(C->depth); + C->biaspoly = NULL; + for (slong k = 0; k < np; k++) + { + C->p[k] = C->R->ffts[k].p; + nmod_init(&C->mod[k], C->p[k]); + sd_fft_ctx_fit_depth(C->R->ffts + k, C->depth); + } + + for (slong k = 0; k < np; k++) + { + /* 2m-th root: g^((p-1)/2m) */ + ulong g = n_primitive_root_prime(C->p[k]); + ulong w = nmod_pow_ui(g, (C->p[k] - 1) / (2 * (ulong) m), + C->mod[k]); + ulong minv = nmod_inv(m, C->mod[k]); + /* the from_ffts contract (determined empirically with a + planted-coefficient probe): stored residues carry the + inverse of the reduced crt cofactor */ + minv = nmod_mul(minv, nmod_inv( + *crt_data_co_prime_red(C->R->crts + C->np - 1, k) % C->p[k], + C->mod[k]), C->mod[k]); + ulong winv = nmod_inv(w, C->mod[k]); + C->tw[k] = flint_aligned_alloc(FLINT_FFT_SMALL_ALIGNMENT, + n_round_up(m * sizeof(double), FLINT_FFT_SMALL_ALIGNMENT)); + C->itw[k] = flint_aligned_alloc(FLINT_FFT_SMALL_ALIGNMENT, + n_round_up(m * sizeof(double), FLINT_FFT_SMALL_ALIGNMENT)); + ulong wp = 1, wip = minv; + for (slong i = 0; i < m; i++) + { + C->tw[k][i] = (double)(slong)(wp <= C->p[k]/2 ? wp + : wp - C->p[k]); + C->itw[k][i] = (double)(slong)(wip <= C->p[k]/2 ? wip + : wip - C->p[k]); + wp = nmod_mul(wp, w, C->mod[k]); + wip = nmod_mul(wip, winv, C->mod[k]); + } + } + fmpz_init(C->Pbig); + fmpz_init(C->Phalf); + fmpz_one(C->Pbig); + for (slong k = 0; k < np; k++) + fmpz_mul_ui(C->Pbig, C->Pbig, C->p[k]); + fmpz_fdiv_q_2exp(C->Phalf, C->Pbig, 1); + { + /* bias B = 2^(2b + depth + 2) exceeds any coefficient + magnitude; biased coefficients stay below the four-prime + product. biaspoly = B * (2^(b m) - 1)/(2^b - 1). */ + slong Bexp = 2 * C->b + C->depth; + slong nc2 = (slong) C->R->crts[np - 1].coeff_len; + C->whi = (C->b * (m - 1)) / FLINT_BITS + nc2 + 2; + (void) Bexp; + for (slong k = 0; k < np; k++) + { + ulong r = nmod_pow_ui(2, (ulong) Bexp, C->mod[k]); + r = nmod_mul(r, nmod_inv( + *crt_data_co_prime_red(C->R->crts + np - 1, k) % C->p[k], + C->mod[k]), C->mod[k]); + C->biasres[k] = r; + } + C->biaspoly = flint_calloc(C->whi + 1, sizeof(ulong)); + for (slong i = 0; i < m; i++) + { + slong bit = C->b * i + Bexp, l = bit / FLINT_BITS, o = bit % FLINT_BITS; + ulong cy = mpn_add_1(C->biaspoly + l, C->biaspoly + l, + C->whi + 1 - l, UWORD(1) << o); + (void) cy; + } + } + C->crt2 = flint_malloc(m * sizeof(slong)); +} + + +/* For a caller-fixed ring 2^N + 1, choose the digit count m over + the instantiated digit sizes dividing N (m = N/b a power of + two), minimizing np * m * (depth + 4 + b/16). For power-of-two N + this switches from (3, 64) to (6, 128) at half the transform + length once three-prime capacity runs out, flattening the large-N + cost staircase. */ +static slong +_mulmod_2expp1_choose_m(mpn_ctx_struct * R, slong N) +{ + static const slong bs[] = { 64, 84, 88, 92, 112, 116, 120, 126, + 128, 136, 140, 144, 160, 164, 168, + 184, 188, 192 }; + slong ii, bestm = 0; + double bestc = 0.0; + + for (ii = 0; ii < (slong) (sizeof(bs) / sizeof(bs[0])); ii++) + { + slong b = bs[ii], m, depth, np; + double cost; + + if (N % b != 0) + continue; + m = N / b; + if (m < 16 || (m & (m - 1)) != 0) + continue; + depth = FLINT_BIT_COUNT(m) - 1; + np = sd_fft_mpn_mulmod_2expp1_choose_np(R, b, depth); + if (np < 0) + continue; + cost = (double) np * (double) m + * ((double) depth + 4.0 + (double) b / 16.0); + if (bestm == 0 || cost < bestc) + { + bestc = cost; + bestm = m; + } + } + return bestm; +} + +void +sd_fft_mpn_mulmod_2expp1_ctx_clear(sd_fft_mpn_mulmod_2expp1_ctx_struct * C) +{ + for (slong k = 0; k < C->np; k++) + { + flint_aligned_free(C->tw[k]); + flint_aligned_free(C->itw[k]); + } + flint_free(C->crt2); + if (C->biaspoly) flint_free(C->biaspoly); + fmpz_clear(C->Pbig); fmpz_clear(C->Phalf); +} + +/* digits of the N-bit residue x (n limbs), b-bit fields, as doubles + ready for the vectorized twist. Vectorized in the style of the + fft_small input pass, across output digits: byte-granular unaligned + loads put each field within one 64-bit word (b <= 56 keeps the field + plus the sub-byte offset inside it), a per-lane variable shift + aligns it, a mask isolates it, and the limited conversion (values + below 2^52) produces the doubles. The scalar tail also covers the + final fields, whose byte loads could run past the n+1 limbs. */ +static void +digits_of_d(double * dig, nn_srcptr x, slong n, slong b, + slong m) +{ + const unsigned char * bytes = (const unsigned char *) x; + ulong mask = (UWORD(1) << b) - 1; + slong i = 0; +#if defined(__AVX2__) + if (b <= 56) + { + /* last i whose 8-byte load ends within the n+1 limbs */ + slong safe = (FLINT_BITS * (slong)(n + 1) - 71) / b - 3; + __m256i vmask = _mm256_set1_epi64x((slong) mask); + for (; i + 4 <= m && i + 4 <= safe; i += 4) + { + ulong b0 = (ulong)(b*(i+0)), b1 = (ulong)(b*(i+1)), + b2 = (ulong)(b*(i+2)), b3 = (ulong)(b*(i+3)); + ulong w0, w1, w2, w3; + memcpy(&w0, bytes + b0/8, 8); + memcpy(&w1, bytes + b1/8, 8); + memcpy(&w2, bytes + b2/8, 8); + memcpy(&w3, bytes + b3/8, 8); + __m256i w = _mm256_setr_epi64x((slong) w0, (slong) w1, + (slong) w2, (slong) w3); + __m256i sh = _mm256_setr_epi64x((slong)(b0 % 8), + (slong)(b1 % 8), (slong)(b2 % 8), (slong)(b3 % 8)); + w = _mm256_and_si256(_mm256_srlv_epi64(w, sh), vmask); + vec4d_store(dig + i, vec4n_convert_limited_vec4d(w)); + } + } +#endif + for (; i < m; i++) + { + slong bit = i * b, l = bit / FLINT_BITS, o = bit % FLINT_BITS; + ulong lo = x[l] >> o; + if (o + b > FLINT_BITS && l + 1 < n) + lo |= x[l+1] << (64 - o); + dig[i] = (double)(slong)(lo & mask); + } +} + + +/* z = x*y mod 2^N+1; x, y in [0, 2^N] (top limb 0 or the value 2^N), + all n = N/64 limbs + 1 top limb; z likewise */ +void +sd_fft_mpn_mulmod_2expp1(sd_fft_mpn_mulmod_2expp1_ctx_struct * C, nn_ptr z, nn_srcptr x, + nn_srcptr y, sd_fft_mpn_mulmod_2expp1_scratch_struct * S) +{ + slong m = C->m, b = C->b, n = C->N / 64, np = C->np; + + slong dsz = (slong) sd_fft_ctx_data_size(C->depth); + double * abuf = S->buf, * bbuf = S->buf + np * dsz; + ulong bits = (ulong) b; + if (b > 50) + { + /* the templated wide-field input pass for the instantiated + (np, bits) pairs; two_pow table by vector-group count */ + ulong stop_easy = n_min((ulong) m, (FLINT_BITS*(ulong)n - 33)/bits); + /* easy-pass alignment contract of the templated input + passes: 4, 8 or 16 coefficients per unrolled iteration + for bits divisible by 8, 4 or 2 respectively */ + stop_easy &= -(ulong) (bits % 8 == 0 ? 4 : bits % 4 == 0 ? 8 + : 16); + ulong stop_hard = n_min((ulong) m, + (FLINT_BITS*(ulong)n + bits - 1)/bits); + const vec4d * tp = C->R->vec_two_pow_tab[(np + 3)/4 - 1]; + to_ffts_func tf = _mpn_ctx_to_ffts_func((ulong) np, bits); + FLINT_ASSERT(tf != NULL); /* guaranteed by ctx_init */ + tf(C->R->ffts, abuf, dsz, x, n, (ulong) m, tp, + 0, stop_easy, stop_easy, stop_hard); + tf(C->R->ffts, bbuf, dsz, y, n, (ulong) m, tp, + 0, stop_easy, stop_easy, stop_hard); + } + else + { + /* fields fit doubles: the vectorized extraction, shared + across primes, twist folded in below */ + digits_of_d(S->digs, x, n, b, m); + S->digs[0] -= (double)(x[n] != 0); + digits_of_d(S->digs + m, y, n, b, m); + S->digs[m] -= (double)(y[n] != 0); + } + for (slong k = 0; k < np; k++) + { + sd_fft_ctx_struct * Q = C->R->ffts + k; + double * X = abuf + k * dsz, * Y = bbuf + k * dsz; + vec1d pd = Q->p, pinv = Q->pinv; + { + vec8d n8 = vec8d_set_d(pd), ninv8 = vec8d_set_d(pinv); + if (b > 50) + { + /* wide fields arrived per prime; top bit folds as + -1 into digit zero */ + if (x[n]) X[0] -= 1.0; + if (y[n]) Y[0] -= 1.0; + for (slong i = 0; i + 8 <= m; i += 8) + { + vec8d tw8 = vec8d_load(C->tw[k] + i); + vec8d_store(X + i, + vec8d_mulmod(vec8d_load(X + i), tw8, + n8, ninv8)); + vec8d_store(Y + i, + vec8d_mulmod(vec8d_load(Y + i), tw8, + n8, ninv8)); + } + } + else + { + const double * dx = S->digs, * dy = S->digs + m; + for (slong i = 0; i + 8 <= m; i += 8) + { + vec8d tw8 = vec8d_load(C->tw[k] + i); + vec8d_store(X + i, + vec8d_mulmod(vec8d_load(dx + i), tw8, + n8, ninv8)); + vec8d_store(Y + i, + vec8d_mulmod(vec8d_load(dy + i), tw8, + n8, ninv8)); + } + } + } + sd_fft_trunc(Q, X, C->depth, m, m); + sd_fft_trunc(Q, Y, C->depth, m, m); + { + vec8d n8 = vec8d_set_d(pd), ninv8 = vec8d_set_d(pinv); + for (slong i = 0; i + 8 <= m; i += 8) + { + vec8d t = vec8d_reduce_to_pm1n(vec8d_load(X + i), + n8, ninv8); + vec8d_store(X + i, vec8d_mulmod(t, + vec8d_load(Y + i), n8, ninv8)); + } + } + sd_ifft_trunc(Q, X, C->depth, m); + { + /* untwist, then bias into nonnegative residues for the + unsigned from_ffts reconstruction */ + vec8d n8 = vec8d_set_d(pd), ninv8 = vec8d_set_d(pinv); + vec1d br = (double)(slong)(C->biasres[k] <= C->p[k]/2 + ? C->biasres[k] + : C->biasres[k] - C->p[k]); + vec8d b8 = vec8d_set_d(br); + for (slong i = 0; i + 8 <= m; i += 8) + { + vec8d t = vec8d_mulmod(vec8d_load(X + i), + vec8d_load(C->itw[k] + i), n8, ninv8); + t = vec8d_add(t, b8); + t = vec8d_reduce_to_pm1n(t, n8, ninv8); + vec8d_store(X + i, t); + } + } + } + { + nn_ptr w = S->w; + slong ncoef = (slong) C->R->crts[np - 1].coeff_len; + ulong E0 = 0; + ulong E1 = (((ulong) C->whi >= (ulong) ncoef + 1) + ? (ulong) C->whi - ((ulong) ncoef + 1) : UWORD(0)) + * FLINT_BITS / (ulong) b; + E1 &= -(ulong) BLK_SZ; + if (E1 > (ulong) m) E1 = (ulong) m & -(ulong) BLK_SZ; + /* from_ffts zeroes its own window */ + { + /* CRT recomposition via mpn_mul.c instances */ + _mpn_ctx_from_ffts_func((ulong) np)(w, 0, (ulong) C->whi, 0, (ulong) m, + C->R->ffts, abuf, dsz, C->R->crts, (ulong) b, + E0, E1, NULL, NULL); + } + mpn_sub_n(w, w, C->biaspoly, C->whi); + /* fold: same shape as the small path; w[n..] is the + signed excess */ + { + /* the signed excess spans b + depth + 2 bits above 2^N: + take the full window, sign-extending past whi */ + slong elen = (C->b + C->depth + 2) / FLINT_BITS + 2; + ulong ex[8]; + ulong sign; + slong t; + + FLINT_ASSERT(elen <= 8); + for (t = 0; t < elen; t++) + ex[t] = (n + t < C->whi) ? w[n + t] + : (ulong)((slong) ex[t - 1] >> (FLINT_BITS - 1)); + sign = (slong) ex[elen - 1] < 0; + if (sign) + { + for (t = 0; t < elen; t++) + ex[t] = ~ex[t]; + mpn_add_1(ex, ex, elen, 1); + } + z[n] = 0; + if (!sign) + { + if (mpn_sub(z, w, n, ex, elen)) + { + if (mpn_add_1(z, z, n, 1)) + { + memset(z, 0, n * sizeof(ulong)); + z[n] = 1; + } + } + } + else + { + if (mpn_add(z, w, n, ex, elen)) + { + if (mpn_sub_1(z, z, n, 1)) + { + memset(z, 0, n * sizeof(ulong)); + z[n] = 1; + } + } + } + } + } +} + +/* per-thread mutable state for sd_fft_mpn_mulmod_2expp1; the context itself is + read-only during multiplication and may be shared */ +void +sd_fft_mpn_mulmod_2expp1_scratch_init(sd_fft_mpn_mulmod_2expp1_scratch_struct * S, + const sd_fft_mpn_mulmod_2expp1_ctx_struct * C) +{ + S->buf = flint_aligned_alloc(FLINT_FFT_SMALL_ALIGNMENT, + n_round_up(2 * C->np * C->dsz * sizeof(double), + FLINT_FFT_SMALL_ALIGNMENT)); + S->w = flint_malloc((C->N/FLINT_BITS + 6) * sizeof(ulong)); + S->digs = (C->b <= 50) + ? flint_aligned_alloc(64, 2 * C->m * sizeof(double)) : NULL; +} + +void +sd_fft_mpn_mulmod_2expp1_scratch_clear(sd_fft_mpn_mulmod_2expp1_scratch_struct * S) +{ + flint_aligned_free(S->buf); + flint_free(S->w); + if (S->digs != NULL) + flint_aligned_free(S->digs); +} diff --git a/src/fft_small/test/main.c b/src/fft_small/test/main.c index 15f94b26a0..71e5428ee1 100644 --- a/src/fft_small/test/main.c +++ b/src/fft_small/test/main.c @@ -15,6 +15,8 @@ #include "t-gr_transformed_nmod_poly.c" #include "t-mpn_add_inplace_c.c" #include "t-mpn_mulmid.c" +#include "t-mpn_ctx_to_from_ffts.c" +#include "t-sd_fft_mpn_mulmod_2expp1.c" #include "t-mul.c" #include "t-nmod_poly_divrem.c" #include "t-nmod_poly_mul.c" @@ -32,6 +34,8 @@ test_struct tests[] = TEST_FUNCTION(flint_mpn_add_inplace_c), TEST_FUNCTION(_mpn_ctx_mpn_mul_range), TEST_FUNCTION(mpn_ctx_mpn_mul), + TEST_FUNCTION(mpn_ctx_to_from_ffts), + TEST_FUNCTION(sd_fft_mpn_mulmod_2expp1), TEST_FUNCTION(_nmod_poly_divrem_mpn_ctx), TEST_FUNCTION(_nmod_poly_mul_mid_mpn_ctx), TEST_FUNCTION(fft_small_op_mpn), diff --git a/src/fft_small/test/t-mpn_ctx_to_from_ffts.c b/src/fft_small/test/t-mpn_ctx_to_from_ffts.c new file mode 100644 index 0000000000..40f5e1c145 --- /dev/null +++ b/src/fft_small/test/t-mpn_ctx_to_from_ffts.c @@ -0,0 +1,129 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + Developed using Claude Fable 5 + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include +#include +#include "test_helpers.h" +#include "ulong_extras.h" +#include "fft_small.h" + +/* Stress the templated digit-splitting input passes directly: for + every instantiated (np, bits) pair, extract digits from random + mpn data through the vectorized pass and compare each residue + against an independently computed digit modulo each prime. */ + +static ulong +ref_digit(nn_srcptr x, slong n, slong bits, slong i, nmod_t mod) +{ + /* value of the bits-wide window at bit offset bits*i, mod p */ + slong lo = bits * i, hi = lo + bits; + ulong r = 0, pow = 1; + slong t; + + for (t = lo; t < hi && t < FLINT_BITS * n; t++) + { + if ((x[t / FLINT_BITS] >> (t % FLINT_BITS)) & 1) + r = nmod_add(r, pow, mod); + pow = nmod_add(pow, pow, mod); + } + return r; +} + +TEST_FUNCTION_START(mpn_ctx_to_from_ffts, state) +{ + static const slong cfg[][2] = { + { 3, 64 }, { 4, 64 }, { 4, 84 }, { 4, 88 }, { 4, 92 }, + { 5, 84 }, { 5, 88 }, { 5, 92 }, + { 5, 112 }, { 5, 116 }, { 5, 120 }, + { 6, 112 }, { 6, 116 }, { 6, 120 }, + { 6, 126 }, { 6, 128 }, { 6, 136 }, { 6, 140 }, { 6, 144 }, + { 7, 126 }, { 7, 128 }, { 7, 136 }, { 7, 140 }, { 7, 144 }, + { 7, 160 }, { 7, 164 }, { 7, 168 }, + { 8, 160 }, { 8, 164 }, { 8, 168 }, + { 8, 184 }, { 8, 188 }, { 8, 192 } + }; + mpn_ctx_struct * R = get_default_mpn_ctx(); + slong iter; + + for (iter = 0; iter < 40 * flint_test_multiplier(); iter++) + { + slong ci = n_randint(state, sizeof(cfg) / sizeof(cfg[0])); + slong np = cfg[ci][0], bits = cfg[ci][1]; + slong m = WORD(32) << n_randint(state, 3); + slong n = (bits * m) / FLINT_BITS + 1, i, k; + to_ffts_func tf = _mpn_ctx_to_ffts_func((ulong) np, + (ulong) bits); + ulong stop_easy, stop_hard, align; + slong depth = FLINT_BIT_COUNT(m) - 1, dsz; + nn_ptr x; + double * buf; + const vec4d * tp; + + if (tf == NULL) + TEST_FUNCTION_FAIL("missing pass np = %wd, bits = %wd\n", + np, bits); + + for (k = 0; k < np; k++) + sd_fft_ctx_fit_depth(R->ffts + k, depth); + dsz = (slong) sd_fft_ctx_data_size(depth); + + x = flint_malloc(n * sizeof(ulong)); + for (i = 0; i < n; i++) + x[i] = n_randlimb(state); + if (n_randint(state, 4) == 0) + x[n - 1] = 0; + + buf = flint_aligned_alloc(FLINT_FFT_SMALL_ALIGNMENT, + n_round_up(np * dsz * sizeof(double), + FLINT_FFT_SMALL_ALIGNMENT)); + + align = bits % 8 == 0 ? 4 : bits % 4 == 0 ? 8 : 16; + stop_easy = ((FLINT_BITS*(ulong)(n - 1) - 33)/bits) + & -(ulong) align; + stop_easy = FLINT_MIN(stop_easy, (ulong) m & -(ulong) align); + stop_hard = FLINT_MIN((ulong) m, + (FLINT_BITS*(ulong)(n - 1) + bits - 1)/bits); + tp = R->vec_two_pow_tab[(np + 3)/4 - 1]; + + tf(R->ffts, buf, dsz, x, n - 1, (ulong) m, tp, + 0, stop_easy, stop_easy, stop_hard); + + for (k = 0; k < np; k++) + { + nmod_t mod; + nmod_init(&mod, R->ffts[k].p); + for (i = 0; i < (slong) stop_hard; i++) + { + double v = (buf + k * dsz)[i]; + slong iv = (slong) v; + ulong got, want; + + if (v != floor(v)) + TEST_FUNCTION_FAIL("nonintegral residue " + "np = %wd bits = %wd k = %wd i = %wd\n", + np, bits, k, i); + got = iv < 0 ? (ulong)(iv + (slong) mod.n) + : (ulong) iv; + got = n_mod2_preinv(got, mod.n, mod.ninv); + want = ref_digit(x, n - 1, bits, i, mod); + if (got != want) + TEST_FUNCTION_FAIL("np = %wd bits = %wd k = %wd " + "i = %wd (m = %wd)\n", np, bits, k, i, m); + } + } + flint_aligned_free(buf); + flint_free(x); + } + + TEST_FUNCTION_END(state); +} diff --git a/src/fft_small/test/t-sd_fft_mpn_mulmod_2expp1.c b/src/fft_small/test/t-sd_fft_mpn_mulmod_2expp1.c new file mode 100644 index 0000000000..0ce4c183cd --- /dev/null +++ b/src/fft_small/test/t-sd_fft_mpn_mulmod_2expp1.c @@ -0,0 +1,213 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + Developed using Claude Fable 5 + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include +#include "test_helpers.h" +#include "ulong_extras.h" +#include "fmpz.h" +#include "mpn_extras.h" +#include "fft_small.h" + +/* reference: z = x*y mod 2^N+1 by a plain product and two folds, + v = f0 - f1 + f2 with f0, f1 the low/next n limbs and f2 the top + two limbs of the (2n+2)-limb full product */ +static void +negmul_ref(nn_ptr z, nn_srcptr x, nn_srcptr y, slong n) +{ + nn_ptr full = flint_malloc((2 * n + 2) * sizeof(ulong)); + ulong cy, br; + + flint_mpn_mul(full, x, n + 1, y, n + 1); + cy = mpn_add(z, full, n, full + 2 * n, 2); + br = mpn_sub_n(z, z, full + n, n); + z[n] = 0; + if (cy > br) + { + if (mpn_sub_1(z, z, n, 1)) + { + memset(z, 0, n * sizeof(ulong)); + z[n] = 1; + } + } + else if (br > cy) + { + if (mpn_add_1(z, z, n, 1)) + { + memset(z, 0, n * sizeof(ulong)); + z[n] = 1; + } + } + flint_free(full); +} + +TEST_FUNCTION_START(sd_fft_mpn_mulmod_2expp1, state) +{ + /* every instantiated wide digit size, plus small-b digit + shapes; np is chosen internally by the context */ + static const slong cfg[] = { + 64, 84, 88, 92, 112, 116, 120, 126, 128, 136, 140, 144, + 160, 164, 168, 184, 188, 192, 34, 40, 46, 50 + }; + slong iter; + + for (iter = 0; iter < 300 * flint_test_multiplier(); iter++) + { + slong ci = n_randint(state, sizeof(cfg) / sizeof(cfg[0])); + slong b = cfg[ci]; + slong m = WORD(16) << n_randint(state, 4); + slong N = m * b, n = N / FLINT_BITS; + sd_fft_mpn_mulmod_2expp1_ctx_t C; + nn_ptr x, y, z, r; + slong i; + + if (N % FLINT_BITS != 0) + continue; + + sd_fft_mpn_mulmod_2expp1_ctx_init(C, get_default_mpn_ctx(), N, + (b > 50 && n_randint(state, 4) == 0) ? 0 : m); + if (C->m != m) /* default splitting may differ; adopt it */ + { + sd_fft_mpn_mulmod_2expp1_ctx_clear(C); + sd_fft_mpn_mulmod_2expp1_ctx_init(C, get_default_mpn_ctx(), + N, m); + } + sd_fft_mpn_mulmod_2expp1_scratch_t S; + sd_fft_mpn_mulmod_2expp1_scratch_init(S, C); + + /* the context must pick the least workable prime count */ + if (C->np > 3) + { + fmpz_t P; + slong k; + fmpz_init_set_ui(P, UWORD(1)); + for (k = 0; k < C->np - 1; k++) + fmpz_mul_ui(P, P, get_default_mpn_ctx()->ffts[k].p); + if ((slong) fmpz_bits(P) - 1 >= 2 * b + + FLINT_BIT_COUNT(m) - 1 + 1 + && (b <= 50 + || _mpn_ctx_to_ffts_func((ulong)(C->np - 1), + (ulong) b) != NULL)) + TEST_FUNCTION_FAIL("np = %wd not minimal for b = %wd, " + "m = %wd\n", C->np, b, m); + fmpz_clear(P); + } + + x = flint_malloc((n + 1) * sizeof(ulong)); + y = flint_malloc((n + 1) * sizeof(ulong)); + z = flint_malloc((n + 1) * sizeof(ulong)); + r = flint_malloc((n + 1) * sizeof(ulong)); + + for (i = 0; i < n; i++) + { + x[i] = n_randlimb(state); + y[i] = n_randlimb(state); + } + x[n] = y[n] = 0; + + switch (n_randint(state, 9)) + { + case 8: /* unreduced: top limb set, low limbs kept */ + x[n] = 1; + break; + case 0: /* x = 2^N */ + memset(x, 0, n * sizeof(ulong)); + x[n] = 1; + break; + case 1: /* y = 2^N */ + memset(y, 0, n * sizeof(ulong)); + y[n] = 1; + break; + case 2: /* squaring */ + memcpy(y, x, (n + 1) * sizeof(ulong)); + break; + case 3: /* sparse operand */ + memset(x, 0, n * sizeof(ulong)); + x[n_randint(state, n)] = n_randlimb(state); + break; + default: + break; + } + + negmul_ref(r, x, y, n); + sd_fft_mpn_mulmod_2expp1(C, z, x, y, S); + + if (z[n] != r[n] || mpn_cmp(z, r, n) != 0) + TEST_FUNCTION_FAIL("b = %wd, m = %wd, N = %wd, np = %wd\n", + b, m, N, C->np); + + /* in-place: z aliases x, as the convolution pointwise does */ + memcpy(z, x, (n + 1) * sizeof(ulong)); + sd_fft_mpn_mulmod_2expp1(C, z, z, y, S); + if (z[n] != r[n] || mpn_cmp(z, r, n) != 0) + TEST_FUNCTION_FAIL("aliased: b = %wd, m = %wd\n", + b, m); + + flint_free(x); + flint_free(y); + flint_free(z); + flint_free(r); + sd_fft_mpn_mulmod_2expp1_scratch_clear(S); + sd_fft_mpn_mulmod_2expp1_ctx_clear(C); + } + + /* second-family shapes: depths where the family prime count + runs out of capacity and np is bumped (e.g. b = 92 at large + depth needs five primes) */ + { + static const slong cfg2[][2] = { + { 92, 4096 }, { 128, 8192 }, { 168, 2048 } + }; + slong ci; + for (ci = 0; ci < 3; ci++) + { + slong b = cfg2[ci][0], m = cfg2[ci][1]; + slong N = m * b, n = N / FLINT_BITS, i; + sd_fft_mpn_mulmod_2expp1_ctx_t C; + nn_ptr x, y, z, r; + + sd_fft_mpn_mulmod_2expp1_ctx_init(C, get_default_mpn_ctx(), N, + (b > 50 && n_randint(state, 4) == 0) ? 0 : m); + if (C->m != m) /* default splitting may differ; adopt it */ + { + sd_fft_mpn_mulmod_2expp1_ctx_clear(C); + sd_fft_mpn_mulmod_2expp1_ctx_init(C, get_default_mpn_ctx(), + N, m); + } + sd_fft_mpn_mulmod_2expp1_scratch_t S; + sd_fft_mpn_mulmod_2expp1_scratch_init(S, C); + x = flint_malloc((n + 1) * sizeof(ulong)); + y = flint_malloc((n + 1) * sizeof(ulong)); + z = flint_malloc((n + 1) * sizeof(ulong)); + r = flint_malloc((n + 1) * sizeof(ulong)); + for (i = 0; i < n; i++) + { + x[i] = n_randlimb(state); + y[i] = n_randlimb(state); + } + x[n] = y[n] = 0; + negmul_ref(r, x, y, n); + sd_fft_mpn_mulmod_2expp1(C, z, x, y, S); + if (z[n] != r[n] || mpn_cmp(z, r, n) != 0) + TEST_FUNCTION_FAIL("deep: b = %wd, m = %wd, np = %wd\n", + b, m, C->np); + flint_free(x); + flint_free(y); + flint_free(z); + flint_free(r); + sd_fft_mpn_mulmod_2expp1_scratch_clear(S); + sd_fft_mpn_mulmod_2expp1_ctx_clear(C); + } + } + + TEST_FUNCTION_END(state); +} diff --git a/src/fmpz_poly.h b/src/fmpz_poly.h index aeeea12d11..fda13afafd 100644 --- a/src/fmpz_poly.h +++ b/src/fmpz_poly.h @@ -51,10 +51,13 @@ typedef struct slong bits2; slong limbs; fmpz_poly_t poly2; + struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx; /* negacyclic pointwise, or NULL */ } fmpz_poly_mul_precache_struct; typedef fmpz_poly_mul_precache_struct fmpz_poly_mul_precache_t[1]; +struct sd_fft_mpn_mulmod_2expp1_ctx_struct; + /* Memory management ********************************************************/ void fmpz_poly_init(fmpz_poly_t poly); diff --git a/src/fmpz_poly/mul.c b/src/fmpz_poly/mul.c index 3610a7d16e..daaa1949c5 100644 --- a/src/fmpz_poly/mul.c +++ b/src/fmpz_poly/mul.c @@ -153,7 +153,8 @@ _fmpz_poly_mul(fmpz * res, const fmpz * poly1, only wins in a specific medium-size region and for huge products when using many threads. */ ((len2 >= 8 && len2 <= 75 && bits1 + bits2 >= 800 && bits1 + bits2 <= 4000) || - (len1 + len2 >= 5000 && bits1 + bits2 >= 5000 + (len1 + len2) / 10 && flint_get_num_threads() >= 4)) + (len1 + len2 >= 5000 && bits1 + bits2 >= 5000 + (len1 + len2) / 10 && flint_get_num_threads() >= 4) || + (bits1 + bits2 >= 80000 && (len1 + len2) * (bits1 + bits2) > 4e9)) _fmpz_poly_mul_SS(res, poly1, len1, poly2, len2); else _fmpz_poly_mul_KS(res, poly1, len1, poly2, len2); diff --git a/src/fmpz_poly/mullow.c b/src/fmpz_poly/mullow.c index 2b5dac82e4..c19a5f2b2e 100644 --- a/src/fmpz_poly/mullow.c +++ b/src/fmpz_poly/mullow.c @@ -165,7 +165,8 @@ _fmpz_poly_mullow(fmpz * res, const fmpz * poly1, slong len1, _fmpz_poly_mullow_karatsuba(res, poly1, len1, poly2, len2, n); else if ((len2 >= 8 && len2 <= 75 && bits1 + bits2 >= 800 && bits1 + bits2 <= 4000) || - (len1 + len2 >= 5000 && bits1 + bits2 >= 5000 + (len1 + len2) / 10 && flint_get_num_threads() >= 4)) + (len1 + len2 >= 5000 && bits1 + bits2 >= 5000 + (len1 + len2) / 10 && flint_get_num_threads() >= 4) || + (bits1 + bits2 >= 80000 && (len1 + len2) * (bits1 + bits2) > 4e9)) _fmpz_poly_mullow_SS(res, poly1, len1, poly2, len2, n); else _fmpz_poly_mullow_KS(res, poly1, len1, poly2, len2, n); diff --git a/src/fmpz_poly/mullow_SS_precache.c b/src/fmpz_poly/mullow_SS_precache.c index d70bedc334..46ca423c53 100644 --- a/src/fmpz_poly/mullow_SS_precache.c +++ b/src/fmpz_poly/mullow_SS_precache.c @@ -1,5 +1,8 @@ /* Copyright (C) 2008-2011, 2020 William Hart + Copyright (C) 2026 Fredrik Johansson + + Developed using Claude Fable 5 This file is part of FLINT. @@ -14,10 +17,12 @@ #include "fmpz_vec.h" #include "fmpz_poly.h" #include "fft.h" +#include "fft_small.h" void fmpz_poly_mul_SS_precache_init(fmpz_poly_mul_precache_t pre, slong len1, slong bits1, const fmpz_poly_t poly2) { + slong negm; slong i, len_out, loglen2; slong output_bits, size; ulong size1, size2; @@ -43,10 +48,20 @@ void fmpz_poly_mul_SS_precache_init(fmpz_poly_mul_precache_t pre, /* round up for sqrt2 trick */ output_bits = (((output_bits - 1) >> (pre->loglen - 2)) + 1) << (pre->loglen - 2); - pre->limbs = (output_bits - 1) / FLINT_BITS + 1; /* initial size of FFT coeffs */ - - if (pre->limbs > FFT_MULMOD_2EXPP1_CUTOFF) /* can't be worse than next power of 2 limbs */ - pre->limbs = (WORD(1) << FLINT_CLOG2(pre->limbs)); + /* final ring selection happens before any allocation */ + pre->limbs = (output_bits - 1) / FLINT_BITS + 1; + pre->limbs = fft_adjust_limbs_sd_fft(pre->limbs, pre->n, &negm); + pre->sdctx = NULL; +#if FLINT_HAVE_FFT_SMALL + if (negm != 0) + { + pre->sdctx = flint_malloc(sizeof(sd_fft_mpn_mulmod_2expp1_ctx_struct)); + sd_fft_mpn_mulmod_2expp1_ctx_init(pre->sdctx, get_default_mpn_ctx(), + FLINT_BITS*pre->limbs, negm); + } +#else + (void) negm; +#endif size = pre->limbs + 1; /* allocate space for ffts */ @@ -83,8 +98,6 @@ void fmpz_poly_mul_SS_precache_init(fmpz_poly_mul_precache_t pre, output_bits = (((output_bits - 1) >> (pre->loglen - 2)) + 1) << (pre->loglen - 2); pre->limbs = (output_bits - 1) / FLINT_BITS + 1; - pre->limbs = fft_adjust_limbs(pre->limbs); /* round up limbs for Nussbaumer */ - fft_precache(pre->jj, pre->loglen - 2, pre->limbs, len_out, t1, t2, s1); fmpz_poly_init(pre->poly2); @@ -93,6 +106,14 @@ void fmpz_poly_mul_SS_precache_init(fmpz_poly_mul_precache_t pre, void fmpz_poly_mul_precache_clear(fmpz_poly_mul_precache_t pre) { +#if FLINT_HAVE_FFT_SMALL + if (pre->sdctx != NULL) + { + sd_fft_mpn_mulmod_2expp1_ctx_clear(pre->sdctx); + flint_free(pre->sdctx); + pre->sdctx = NULL; + } +#endif flint_free(pre->jj); fmpz_poly_clear(pre->poly2); } @@ -140,7 +161,7 @@ void _fmpz_poly_mullow_SS_precache(fmpz * output, const fmpz * input1, for (i = len1; i < 4*pre->n; i++) flint_mpn_zero(ii[i], size); - fft_convolution_precache(ii, pre->jj, pre->loglen - 2, pre->limbs, + fft_convolution_precache_sd_fft(pre->sdctx, ii, pre->jj, pre->loglen - 2, pre->limbs, len_out, t1, t2, s1, tt); _fmpz_vec_set_fft(output, trunc, ii, pre->limbs, 1); /* write output */ diff --git a/src/fmpz_poly/mulmid.c b/src/fmpz_poly/mulmid.c index 6360791f5b..a22bce4caa 100644 --- a/src/fmpz_poly/mulmid.c +++ b/src/fmpz_poly/mulmid.c @@ -231,7 +231,8 @@ _fmpz_poly_mulmid(fmpz * res, const fmpz * poly1, slong len1, _fmpz_poly_mullow_karatsuba(res, poly1, len1, poly2, len2, nhi); else if ((len2 >= 8 && len2 <= 75 && bits1 + bits2 >= 800 && bits1 + bits2 <= 4000) || - (len1 + len2 >= 5000 && bits1 + bits2 >= 5000 + (len1 + len2) / 10 && flint_get_num_threads() >= 4)) + (len1 + len2 >= 5000 && bits1 + bits2 >= 5000 + (len1 + len2) / 10 && flint_get_num_threads() >= 4) || + (bits1 + bits2 >= 80000 && (len1 + len2) * (bits1 + bits2) > 4e9)) _fmpz_poly_mulmid_SS(res, poly1, len1, poly2, len2, nlo, nhi); else _fmpz_poly_mulmid_KS(res, poly1, len1, poly2, len2, nlo, nhi); diff --git a/src/fmpz_poly/mulmid_SS.c b/src/fmpz_poly/mulmid_SS.c index 8ff8b99ac6..9a70755662 100644 --- a/src/fmpz_poly/mulmid_SS.c +++ b/src/fmpz_poly/mulmid_SS.c @@ -15,15 +15,17 @@ #include "fmpz_vec.h" #include "fmpz_poly.h" #include "fft.h" +#include "fft_small.h" + void _fmpz_poly_mulmid_SS(fmpz * res, const fmpz * poly1, slong len1, const fmpz * poly2, slong len2, slong nlo, slong nhi) { + slong negm; slong len_out, loglen, loglen2, n; slong res_bits, limbs, size, i; ulong * ptr, ** t1, ** t2, ** tt, ** s1, ** ii, ** jj; slong bits1, bits2; - ulong size1, size2; int sign = 0; int N; TMP_INIT; @@ -101,18 +103,23 @@ void _fmpz_poly_mulmid_SS(fmpz * res, const fmpz * poly1, slong len1, else bits2 = _fmpz_vec_max_bits(poly2, len2); - size1 = (FLINT_ABS(bits1) + FLINT_BITS - 1) / FLINT_BITS; - size2 = (FLINT_ABS(bits2) + FLINT_BITS - 1) / FLINT_BITS; - /* Start with an upper bound on the number of bits needed */ - res_bits = FLINT_BITS * (size1 + size2) + loglen2 + 1; + if (bits1 < WORD(0) || bits2 < WORD(0)) + { + sign = 1; + bits1 = FLINT_ABS(bits1); + bits2 = FLINT_ABS(bits2); + } - /* round up for sqrt2 trick */ + /* exact bit bound, rounded up for the sqrt2 trick */ + res_bits = bits1 + bits2 + loglen2 + sign; res_bits = (((res_bits - 1) >> (loglen - 2)) + 1) << (loglen - 2); - limbs = (res_bits - 1) / FLINT_BITS + 1; /* initial size of FFT coeffs */ - if (limbs > FFT_MULMOD_2EXPP1_CUTOFF) /* can't be worse than next power of 2 limbs */ - limbs = (WORD(1) << FLINT_CLOG2(limbs)); + /* final ring selection happens before any allocation; without + the fft_small engine the chooser degenerates to the standard + Nussbaumer rounding with negm = 0 */ + limbs = (res_bits - 1) / FLINT_BITS + 1; + limbs = fft_adjust_limbs_sd_fft(limbs, n, &negm); size = limbs + 1; /* allocate space for ffts */ @@ -158,24 +165,23 @@ void _fmpz_poly_mulmid_SS(fmpz * res, const fmpz * poly1, slong len1, flint_mpn_zero(jj[i], limbs + 1); } - if (bits1 < WORD(0) || bits2 < WORD(0)) { - sign = 1; - bits1 = FLINT_ABS(bits1); - bits2 = FLINT_ABS(bits2); +#if FLINT_HAVE_FFT_SMALL + if (negm != 0) + { + sd_fft_mpn_mulmod_2expp1_ctx_t sdctx; + sd_fft_mpn_mulmod_2expp1_ctx_init(sdctx, get_default_mpn_ctx(), + FLINT_BITS*limbs, negm); + fft_convolution_sd_fft(sdctx, ii, jj, loglen - 2, limbs, + len_out, t1, t2, s1, tt); + sd_fft_mpn_mulmod_2expp1_ctx_clear(sdctx); + } + else +#endif + fft_convolution(ii, jj, loglen - 2, limbs, len_out, + t1, t2, s1, tt); } - /* Recompute the number of bits/limbs now that we know how large everything is */ - res_bits = bits1 + bits2 + loglen2 + sign; - - /* round up res bits for sqrt2 */ - res_bits = (((res_bits - 1) >> (loglen - 2)) + 1) << (loglen - 2); - - limbs = (res_bits - 1) / FLINT_BITS + 1; - limbs = fft_adjust_limbs(limbs); /* round up limbs for Nussbaumer */ - - fft_convolution(ii, jj, loglen - 2, limbs, len_out, t1, t2, s1, tt); - _fmpz_vec_set_fft(res, trunc - nlo, ii + nlo, limbs, sign); /* write res */ flint_free(ii); From 8ce753065dc6134e283af8c2555c7f25c43fb9c9 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Sat, 8 Aug 2026 01:23:31 +0200 Subject: [PATCH 2/3] Improve test coverage and fix bug in the new mullow_SS_precache --- src/fft.h | 7 + src/fft/convolution.c | 45 ++++- src/fft/convolution_precache.c | 6 +- src/fft/test/main.c | 4 + src/fft/test/t-convolution_precache_sd_fft.c | 156 ++++++++++++++++ src/fft/test/t-convolution_sd_fft.c | 169 ++++++++++++++++++ .../test/t-sd_fft_mpn_mulmod_2expp1.c | 64 +++++++ src/fmpz_poly/mullow_SS_precache.c | 20 +-- src/fmpz_poly/test/t-mullow_SS_precache.c | 43 +++++ 9 files changed, 494 insertions(+), 20 deletions(-) create mode 100644 src/fft/test/t-convolution_precache_sd_fft.c create mode 100644 src/fft/test/t-convolution_sd_fft.c diff --git a/src/fft.h b/src/fft.h index fb11c90180..22e125c3e7 100644 --- a/src/fft.h +++ b/src/fft.h @@ -228,6 +228,13 @@ void _fft_pointwise_rows_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdc mp_limb_t ** ii, mp_limb_t ** jj, slong base, slong srt, slong nrows, slong rvbits, slong limbs); +void _fft_pointwise_sd_fft_ro(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, slong j0, slong j1, + slong limbs); +void _fft_pointwise_rows_sd_fft_ro(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, slong base, + slong srt, slong nrows, slong rvbits, slong limbs); + void fft_convolution_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, mp_limb_t ** ii, mp_limb_t ** jj, slong depth, slong limbs, slong trunc, mp_limb_t ** t1, diff --git a/src/fft/convolution.c b/src/fft/convolution.c index e45beb6503..d37c70779a 100644 --- a/src/fft/convolution.c +++ b/src/fft/convolution.c @@ -161,6 +161,7 @@ typedef struct slong limbs; sd_fft_mpn_mulmod_2expp1_scratch_struct * scr; int rows; + int normjj; } _sd_fft_pw_arg_struct; @@ -168,7 +169,7 @@ static void _sd_fft_pw_one(_sd_fft_pw_arg_struct * X, slong j) { mpn_normmod_2expp1(X->ii[j], X->limbs); - if (X->ii != X->jj) + if (X->normjj && X->ii != X->jj) mpn_normmod_2expp1(X->jj[j], X->limbs); sd_fft_mpn_mulmod_2expp1(X->C, X->ii[j], X->ii[j], X->jj[j], X->scr); } @@ -255,7 +256,18 @@ _fft_pointwise_sd_fft(sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, mp_limb_t ** mp_limb_t ** jj, slong j0, slong j1, slong limbs) { _sd_fft_pw_arg_struct X = { sdctx, ii, jj, j0, j1, - 0, 0, 0, limbs, NULL, 0 }; + 0, 0, 0, limbs, NULL, 0, 1 }; + _sd_fft_pw_run(&X, j1 - j0); +} + +/* as above with the second operand read-only, for the precached + convolutions whose transformed operand must not be modified */ +void +_fft_pointwise_sd_fft_ro(sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, mp_limb_t ** ii, + mp_limb_t ** jj, slong j0, slong j1, slong limbs) +{ + _sd_fft_pw_arg_struct X = { sdctx, ii, jj, j0, j1, + 0, 0, 0, limbs, NULL, 0, 0 }; _sd_fft_pw_run(&X, j1 - j0); } @@ -265,7 +277,17 @@ _fft_pointwise_rows_sd_fft(sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, slong srt, slong nrows, slong rvbits, slong limbs) { _sd_fft_pw_arg_struct X = { sdctx, ii, jj, 0, nrows, - base, srt, rvbits, limbs, NULL, 1 }; + base, srt, rvbits, limbs, NULL, 1, 1 }; + _sd_fft_pw_run(&X, nrows); +} + +void +_fft_pointwise_rows_sd_fft_ro(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, slong base, + slong srt, slong nrows, slong rvbits, slong limbs) +{ + _sd_fft_pw_arg_struct X = { sdctx, ii, jj, 0, nrows, + base, srt, rvbits, limbs, NULL, 1, 0 }; _sd_fft_pw_run(&X, nrows); } #else @@ -290,6 +312,23 @@ _fft_pointwise_rows_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, (void) sdctx; (void) ii; (void) jj; (void) base; (void) srt; (void) nrows; (void) rvbits; (void) limbs; } + +void +_fft_pointwise_sd_fft_ro(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, slong j0, + slong j1, slong limbs) +{ + _fft_pointwise_sd_fft(sdctx, ii, jj, j0, j1, limbs); +} + +void +_fft_pointwise_rows_sd_fft_ro(struct sd_fft_mpn_mulmod_2expp1_ctx_struct * sdctx, + mp_limb_t ** ii, mp_limb_t ** jj, slong base, + slong srt, slong nrows, slong rvbits, slong limbs) +{ + _fft_pointwise_rows_sd_fft(sdctx, ii, jj, base, srt, nrows, + rvbits, limbs); +} #endif /* fft_convolution_basic with an optional negacyclic pointwise diff --git a/src/fft/convolution_precache.c b/src/fft/convolution_precache.c index e18ea77179..6875dacc2d 100644 --- a/src/fft/convolution_precache.c +++ b/src/fft/convolution_precache.c @@ -141,7 +141,7 @@ void fft_convolution_precache_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct fft_truncate_sqrt2(ii, n, w, t1, t2, s1, trunc); if (sdctx != NULL) - _fft_pointwise_sd_fft(sdctx, ii, jj, 0, trunc, limbs); + _fft_pointwise_sd_fft_ro(sdctx, ii, jj, 0, trunc, limbs); else for (j = 0; j < trunc; j++) { mpn_normmod_2expp1(ii[j], limbs); @@ -163,7 +163,7 @@ void fft_convolution_precache_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct fft_mfa_truncate_sqrt2(ii, n, w, t1, t2, s1, sqrt, trunc); if (sdctx != NULL) - _fft_pointwise_sd_fft(sdctx, ii, jj, 0, 2*n, limbs); + _fft_pointwise_sd_fft_ro(sdctx, ii, jj, 0, 2*n, limbs); else for (j = 0; j < 2*n; j++) { mpn_normmod_2expp1(ii[j], limbs); @@ -174,7 +174,7 @@ void fft_convolution_precache_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct trunc2 = (trunc - 2*n)/sqrt; if (sdctx != NULL) - _fft_pointwise_rows_sd_fft(sdctx, ii, jj, 2*n, sqrt, trunc2, + _fft_pointwise_rows_sd_fft_ro(sdctx, ii, jj, 2*n, sqrt, trunc2, depth - depth/2 + 1, limbs); else for (s = 0; s < trunc2; s++) { diff --git a/src/fft/test/main.c b/src/fft/test/main.c index 30d36f769b..b4403fce2d 100644 --- a/src/fft/test/main.c +++ b/src/fft/test/main.c @@ -20,6 +20,8 @@ #include "t-butterfly_twiddle.c" #include "t-convolution.c" #include "t-convolution_precache.c" +#include "t-convolution_precache_sd_fft.c" +#include "t-convolution_sd_fft.c" #include "t-div_2expmod_2expp1.c" #include "t-fft_ifft_mfa_truncate_sqrt2.c" #include "t-fft_ifft_negacyclic.c" @@ -48,6 +50,8 @@ test_struct tests[] = TEST_FUNCTION(fft_ifft_butterfly_twiddle), TEST_FUNCTION(fft_convolution), TEST_FUNCTION(fft_convolution_precache), + TEST_FUNCTION(fft_convolution_precache_sd_fft), + TEST_FUNCTION(fft_convolution_sd_fft), TEST_FUNCTION(mpn_div_2expmod_2expp1), TEST_FUNCTION(fft_ifft_mfa_truncate_sqrt2), TEST_FUNCTION(fft_ifft_negacyclic), diff --git a/src/fft/test/t-convolution_precache_sd_fft.c b/src/fft/test/t-convolution_precache_sd_fft.c new file mode 100644 index 0000000000..1ab277b6dc --- /dev/null +++ b/src/fft/test/t-convolution_precache_sd_fft.c @@ -0,0 +1,156 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + Developed using Claude Fable 5 + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include +#include "test_helpers.h" +#include "mpn_extras.h" +#include "thread_support.h" +#include "fft.h" +#include "fft_small.h" + +/* fft_convolution_precache_sd_fft against classic + fft_convolution_precache on identical inputs with a shared + precached operand, across depths and thread counts. */ + +TEST_FUNCTION_START(fft_convolution_precache_sd_fft, state) +{ +#if FLINT_HAVE_FFT_SMALL + slong depth, ti, iter; + static const int threads[] = { 1, 2, 4 }; + slong save_threads = flint_get_num_threads(); + + for (depth = 3; depth <= 8; depth++) + for (ti = 0; ti < 3; ti++) + for (iter = 0; iter < 2; iter++) + { + mp_size_t n = (UWORD(1) << depth); + mp_size_t n1 = (UWORD(1) << (depth/2)); + slong b, mm, N = 0; + mp_size_t limbs, size, trunc, i, len1, len2; + mp_limb_t * ptr; + mp_limb_t ** ii, ** jj, ** ii2, ** jj2; + mp_limb_t ** t1, ** t2, ** s1, ** tt; + sd_fft_mpn_mulmod_2expp1_ctx_t C; + slong nt; + int squaring; + + /* an engine-compatible ring that the outer transform + accepts: N a multiple of n with an even w */ + for (b = 64; N == 0 && b <= 128; b += 64) + for (mm = 16; mm <= 512; mm *= 2) + if ((mm * b) % n == 0 && ((mm * b) / n) % 2 == 0 + && mm * b >= 24 * n) + { + N = mm * b; + break; + } + if (N == 0) + continue; + + flint_set_num_threads(threads[ti]); + nt = flint_get_num_threads(); + + limbs = N / FLINT_BITS; + size = limbs + 1; + trunc = 2*n + n_randint(state, 2*n) + 1; + trunc = 2*n1*((trunc + 2*n1 - 1)/(2*n1)); + squaring = 0; + len1 = 1 + n_randint(state, trunc); + len2 = squaring ? len1 : 1 + n_randint(state, trunc); + + sd_fft_mpn_mulmod_2expp1_ctx_init(C, get_default_mpn_ctx(), + N, mm); + + ii = flint_malloc((4*(n + n*size) + 5*size*nt)*sizeof(mp_limb_t) + + (4*n + 5*nt)*sizeof(mp_limb_t *)); + for (i = 0, ptr = (mp_limb_t *)(ii + 4*n + 5*nt); i < 4*n; + i++, ptr += size) + { + ii[i] = ptr; + if (i < len1) + random_fermat(ii[i], state, limbs); + else + flint_mpn_zero(ii[i], size); + mpn_normmod_2expp1(ii[i], limbs); + } + t1 = (mp_limb_t **)(ii + 4*n); + t2 = t1 + nt; + s1 = t2 + nt; + tt = s1 + nt; + for (i = 0; i < nt; i++) + { + t1[i] = ptr; ptr += size; + t2[i] = ptr; ptr += size; + s1[i] = ptr; ptr += size; + tt[i] = ptr; ptr += 2*size; + } + + if (squaring) + jj = ii; + else + { + jj = flint_malloc(4*(n + n*size)*sizeof(mp_limb_t) + + 4*n*sizeof(mp_limb_t *)); + for (i = 0, ptr = (mp_limb_t *)(jj + 4*n); i < 4*n; + i++, ptr += size) + { + jj[i] = ptr; + if (i < len2) + random_fermat(jj[i], state, limbs); + else + flint_mpn_zero(jj[i], size); + mpn_normmod_2expp1(jj[i], limbs); + } + } + + fft_precache(jj, depth, limbs, trunc, t1, t2, s1); + + /* classic copies (the precached operand is shared) */ + ii2 = flint_malloc(4*(n + n*size)*sizeof(mp_limb_t) + + 4*n*sizeof(mp_limb_t *)); + for (i = 0, ptr = (mp_limb_t *)(ii2 + 4*n); i < 4*n; + i++, ptr += size) + { + ii2[i] = ptr; + flint_mpn_copyi(ii2[i], ii[i], size); + } + jj2 = jj; + + fft_convolution_precache_sd_fft(C, ii, jj, depth, limbs, trunc, + t1, t2, s1, tt); + fft_convolution_precache(ii2, jj2, depth, limbs, trunc, + t1, t2, s1, tt); + + for (i = 0; i < trunc; i++) + { + mpn_normmod_2expp1(ii[i], limbs); + mpn_normmod_2expp1(ii2[i], limbs); + if (mpn_cmp(ii[i], ii2[i], size) != 0) + TEST_FUNCTION_FAIL("depth = %wd, N = %wd, m = %wd, " + "threads = %wd, trunc = %wd, i = %wd, sq = %d\n", + depth, N, mm, nt, (slong) trunc, (slong) i, + squaring); + } + + sd_fft_mpn_mulmod_2expp1_ctx_clear(C); + flint_free(ii); + if (!squaring) + flint_free(jj); + flint_free(ii2); + } + + flint_set_num_threads(save_threads); +#endif + + TEST_FUNCTION_END(state); +} diff --git a/src/fft/test/t-convolution_sd_fft.c b/src/fft/test/t-convolution_sd_fft.c new file mode 100644 index 0000000000..57f4a567c7 --- /dev/null +++ b/src/fft/test/t-convolution_sd_fft.c @@ -0,0 +1,169 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + Developed using Claude Fable 5 + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include +#include "test_helpers.h" +#include "mpn_extras.h" +#include "thread_support.h" +#include "fft.h" +#include "fft_small.h" + +/* fft_convolution_sd_fft against classic fft_convolution on + identical inputs, across transform depths, thread counts (to + exercise the threaded pointwise helpers and both sides of the + fused-versus-flat dispatch), truncation, and squaring. */ + +TEST_FUNCTION_START(fft_convolution_sd_fft, state) +{ +#if FLINT_HAVE_FFT_SMALL + slong depth, ti, iter; + static const int threads[] = { 1, 2, 4 }; + slong save_threads = flint_get_num_threads(); + + for (depth = 3; depth <= 8; depth++) + for (ti = 0; ti < 3; ti++) + for (iter = 0; iter < 2; iter++) + { + mp_size_t n = (UWORD(1) << depth); + mp_size_t n1 = (UWORD(1) << (depth/2)); + slong b, mm, N = 0; + mp_size_t limbs, size, trunc, i, len1, len2; + mp_limb_t * ptr; + mp_limb_t ** ii, ** jj, ** ii2, ** jj2; + mp_limb_t ** t1, ** t2, ** s1, ** tt; + sd_fft_mpn_mulmod_2expp1_ctx_t C; + slong nt; + int squaring; + + /* an engine-compatible ring that the outer transform + accepts: N a multiple of n with an even w */ + for (b = 64; N == 0 && b <= 128; b += 64) + for (mm = 16; mm <= 512; mm *= 2) + if ((mm * b) % n == 0 && ((mm * b) / n) % 2 == 0 + && mm * b >= 24 * n) + { + N = mm * b; + break; + } + if (N == 0) + continue; + + flint_set_num_threads(threads[ti]); + nt = flint_get_num_threads(); + + limbs = N / FLINT_BITS; + size = limbs + 1; + trunc = 2*n + n_randint(state, 2*n) + 1; + trunc = 2*n1*((trunc + 2*n1 - 1)/(2*n1)); + squaring = (n_randint(state, 4) == 0); + len1 = 1 + n_randint(state, trunc); + len2 = squaring ? len1 : 1 + n_randint(state, trunc); + + sd_fft_mpn_mulmod_2expp1_ctx_init(C, get_default_mpn_ctx(), + N, mm); + + ii = flint_malloc((4*(n + n*size) + 5*size*nt)*sizeof(mp_limb_t) + + (4*n + 5*nt)*sizeof(mp_limb_t *)); + for (i = 0, ptr = (mp_limb_t *)(ii + 4*n + 5*nt); i < 4*n; + i++, ptr += size) + { + ii[i] = ptr; + if (i < len1) + random_fermat(ii[i], state, limbs); + else + flint_mpn_zero(ii[i], size); + mpn_normmod_2expp1(ii[i], limbs); + } + t1 = (mp_limb_t **)(ii + 4*n); + t2 = t1 + nt; + s1 = t2 + nt; + tt = s1 + nt; + for (i = 0; i < nt; i++) + { + t1[i] = ptr; ptr += size; + t2[i] = ptr; ptr += size; + s1[i] = ptr; ptr += size; + tt[i] = ptr; ptr += 2*size; + } + + if (squaring) + jj = ii; + else + { + jj = flint_malloc(4*(n + n*size)*sizeof(mp_limb_t) + + 4*n*sizeof(mp_limb_t *)); + for (i = 0, ptr = (mp_limb_t *)(jj + 4*n); i < 4*n; + i++, ptr += size) + { + jj[i] = ptr; + if (i < len2) + random_fermat(jj[i], state, limbs); + else + flint_mpn_zero(jj[i], size); + mpn_normmod_2expp1(jj[i], limbs); + } + } + + /* classic copies */ + ii2 = flint_malloc(4*(n + n*size)*sizeof(mp_limb_t) + + 4*n*sizeof(mp_limb_t *)); + for (i = 0, ptr = (mp_limb_t *)(ii2 + 4*n); i < 4*n; + i++, ptr += size) + { + ii2[i] = ptr; + flint_mpn_copyi(ii2[i], ii[i], size); + } + if (squaring) + jj2 = ii2; + else + { + jj2 = flint_malloc(4*(n + n*size)*sizeof(mp_limb_t) + + 4*n*sizeof(mp_limb_t *)); + for (i = 0, ptr = (mp_limb_t *)(jj2 + 4*n); i < 4*n; + i++, ptr += size) + { + jj2[i] = ptr; + flint_mpn_copyi(jj2[i], jj[i], size); + } + } + + fft_convolution_sd_fft(C, ii, jj, depth, limbs, trunc, + t1, t2, s1, tt); + fft_convolution(ii2, jj2, depth, limbs, trunc, + t1, t2, s1, tt); + + for (i = 0; i < trunc; i++) + { + mpn_normmod_2expp1(ii[i], limbs); + mpn_normmod_2expp1(ii2[i], limbs); + if (mpn_cmp(ii[i], ii2[i], size) != 0) + TEST_FUNCTION_FAIL("depth = %wd, N = %wd, m = %wd, " + "threads = %wd, trunc = %wd, i = %wd, sq = %d\n", + depth, N, mm, nt, (slong) trunc, (slong) i, + squaring); + } + + sd_fft_mpn_mulmod_2expp1_ctx_clear(C); + flint_free(ii); + if (!squaring) + flint_free(jj); + flint_free(ii2); + if (!squaring) + flint_free(jj2); + } + + flint_set_num_threads(save_threads); +#endif + + TEST_FUNCTION_END(state); +} diff --git a/src/fft_small/test/t-sd_fft_mpn_mulmod_2expp1.c b/src/fft_small/test/t-sd_fft_mpn_mulmod_2expp1.c index 0ce4c183cd..7564de2a85 100644 --- a/src/fft_small/test/t-sd_fft_mpn_mulmod_2expp1.c +++ b/src/fft_small/test/t-sd_fft_mpn_mulmod_2expp1.c @@ -209,5 +209,69 @@ TEST_FUNCTION_START(sd_fft_mpn_mulmod_2expp1, state) } } + /* fold carry cascades: products at the wrap boundary, where + the final subtraction borrows or the final addition carries + all the way into a 2^N result */ + { + static const slong cfg3[][2] = { { 64, 16 }, { 92, 32 }, + { 64, 64 } }; + slong ci; + for (ci = 0; ci < 3; ci++) + { + slong b = cfg3[ci][0], m = cfg3[ci][1]; + slong N = m * b, n = N / FLINT_BITS, cs; + sd_fft_mpn_mulmod_2expp1_ctx_t C; + sd_fft_mpn_mulmod_2expp1_scratch_t S; + nn_ptr x, y, z, r; + + if (N % (2 * FLINT_BITS) != 0) + continue; + + sd_fft_mpn_mulmod_2expp1_ctx_init(C, get_default_mpn_ctx(), + N, m); + sd_fft_mpn_mulmod_2expp1_scratch_init(S, C); + x = flint_calloc(n + 1, sizeof(ulong)); + y = flint_calloc(n + 1, sizeof(ulong)); + z = flint_calloc(n + 1, sizeof(ulong)); + r = flint_calloc(n + 1, sizeof(ulong)); + + for (cs = 0; cs < 4; cs++) + { + memset(x, 0, (n + 1) * sizeof(ulong)); + memset(y, 0, (n + 1) * sizeof(ulong)); + switch (cs) + { + case 0: /* (2^(N/2))^2 = 2^N */ + x[n/2] = y[n/2] = 1; + break; + case 1: /* (2^N)^2 = 1 */ + x[n] = y[n] = 1; + break; + case 2: /* (2^(N/2)+1)(2^(N/2)-1) = 2^N - 1 */ + x[n/2] = 1; x[0] = 1; + y[n/2] = 1; + mpn_sub_1(y, y, n, 1); + break; + case 3: /* 2^N * 2^(N/2) = -2^(N/2) */ + x[n] = 1; + y[n/2] = 1; + break; + } + negmul_ref(r, x, y, n); + sd_fft_mpn_mulmod_2expp1(C, z, x, y, S); + if (z[n] != r[n] || mpn_cmp(z, r, n) != 0) + TEST_FUNCTION_FAIL("wrap case %wd: b = %wd, " + "m = %wd\n", cs, b, m); + } + + flint_free(x); + flint_free(y); + flint_free(z); + flint_free(r); + sd_fft_mpn_mulmod_2expp1_scratch_clear(S); + sd_fft_mpn_mulmod_2expp1_ctx_clear(C); + } + } + TEST_FUNCTION_END(state); } diff --git a/src/fmpz_poly/mullow_SS_precache.c b/src/fmpz_poly/mullow_SS_precache.c index 46ca423c53..6aa85a2ac8 100644 --- a/src/fmpz_poly/mullow_SS_precache.c +++ b/src/fmpz_poly/mullow_SS_precache.c @@ -25,7 +25,6 @@ void fmpz_poly_mul_SS_precache_init(fmpz_poly_mul_precache_t pre, slong negm; slong i, len_out, loglen2; slong output_bits, size; - ulong size1, size2; ulong * ptr; ulong ** t1, ** t2, ** s1; int N; @@ -38,17 +37,17 @@ void fmpz_poly_mul_SS_precache_init(fmpz_poly_mul_precache_t pre, pre->loglen = FLINT_CLOG2(len_out); loglen2 = FLINT_CLOG2(FLINT_MIN(len1, pre->len2)); pre->n = (WORD(1) << (pre->loglen - 2)); - size1 = FLINT_ABS(bits1); - size1 = (size1 + FLINT_BITS - 1)/FLINT_BITS; - size2 = (pre->bits2 + FLINT_BITS - 1)/FLINT_BITS; - /* Start with an upper bound on the number of bits needed */ - output_bits = FLINT_BITS*(size1 + size2) + loglen2 + 1; + /* The transform ring is fixed once, up front: the output bit + bound uses the actual bit counts (both available here), and + ring selection happens before any allocation so that the + buffers, the packed operand, the transforms and the optional + sd_fft pointwise context all live in the same ring. */ + output_bits = bits1 + pre->bits2 + loglen2 + 1; /* round up for sqrt2 trick */ output_bits = (((output_bits - 1) >> (pre->loglen - 2)) + 1) << (pre->loglen - 2); - /* final ring selection happens before any allocation */ pre->limbs = (output_bits - 1) / FLINT_BITS + 1; pre->limbs = fft_adjust_limbs_sd_fft(pre->limbs, pre->n, &negm); pre->sdctx = NULL; @@ -91,13 +90,6 @@ void fmpz_poly_mul_SS_precache_init(fmpz_poly_mul_precache_t pre, for (i = pre->len2; i < 4*pre->n; i++) flint_mpn_zero(pre->jj[i], size); - /* Recompute the number of bits/limbs now that we know how large everything is */ - output_bits = bits1 + pre->bits2 + loglen2 + 1; - - /* round up output bits for sqrt2 */ - output_bits = (((output_bits - 1) >> (pre->loglen - 2)) + 1) << (pre->loglen - 2); - - pre->limbs = (output_bits - 1) / FLINT_BITS + 1; fft_precache(pre->jj, pre->loglen - 2, pre->limbs, len_out, t1, t2, s1); fmpz_poly_init(pre->poly2); diff --git a/src/fmpz_poly/test/t-mullow_SS_precache.c b/src/fmpz_poly/test/t-mullow_SS_precache.c index d35e2fe33d..a83bdd5745 100644 --- a/src/fmpz_poly/test/t-mullow_SS_precache.c +++ b/src/fmpz_poly/test/t-mullow_SS_precache.c @@ -102,5 +102,48 @@ TEST_FUNCTION_START(fmpz_poly_mullow_SS_precache, state) fmpz_poly_clear(d); } + /* engine-triggering coefficient sizes across thread counts: + exercises the sd_fft context in the precache structure and + the mfa branch of the precached convolution */ + { + slong iter2; + slong save_threads = flint_get_num_threads(); + for (iter2 = 0; iter2 < 4; iter2++) + { + fmpz_poly_t a, b, c, d; + fmpz_poly_mul_precache_t pre; + slong len1 = 100 + n_randint(state, 220); + slong len2 = 100 + n_randint(state, 220); + slong bits = 14000 + n_randint(state, 12000); + slong tr; + + flint_set_num_threads(1 + n_randint(state, 3)); + + fmpz_poly_init(a); + fmpz_poly_init(b); + fmpz_poly_init(c); + fmpz_poly_init(d); + fmpz_poly_randtest(a, state, len1, bits); + fmpz_poly_randtest(b, state, len2, bits); + if (b->length == 0) + fmpz_poly_set_ui(b, 1); + + fmpz_poly_mul_SS_precache_init(pre, len1, bits + 1, b); + tr = a->length + b->length - 1; + fmpz_poly_mullow_SS_precache(c, a, pre, tr); + fmpz_poly_mul(d, a, b); + if (!fmpz_poly_equal(c, d)) + TEST_FUNCTION_FAIL("precache sd_fft: len1 = %wd, " + "len2 = %wd, bits = %wd\n", len1, len2, bits); + + fmpz_poly_mul_precache_clear(pre); + fmpz_poly_clear(a); + fmpz_poly_clear(b); + fmpz_poly_clear(c); + fmpz_poly_clear(d); + } + flint_set_num_threads(save_threads); + } + TEST_FUNCTION_END(state); } From c2fa2bbb594f66e37d86e598fb2b696b58f9eb1a Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Sat, 8 Aug 2026 11:30:31 +0200 Subject: [PATCH 3/3] More test coverage --- src/fft/convolution_precache.c | 4 +++ src/fft/fft_mfa_truncate_sqrt2_inner.c | 4 +++ src/fft/test/t-convolution_precache_sd_fft.c | 12 +++++--- src/fft/test/t-convolution_sd_fft.c | 28 ++++++++++++++++--- src/fft_small/negacyclic.c | 8 ++++-- .../test/t-sd_fft_mpn_mulmod_2expp1.c | 21 +++++++++++++- 6 files changed, 66 insertions(+), 11 deletions(-) diff --git a/src/fft/convolution_precache.c b/src/fft/convolution_precache.c index 6875dacc2d..cf7a0abd98 100644 --- a/src/fft/convolution_precache.c +++ b/src/fft/convolution_precache.c @@ -134,6 +134,10 @@ void fft_convolution_precache_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_struct slong w = (limbs*FLINT_BITS)/n; slong sqrt = (WORD(1)<<(depth/2)); + #if FLINT_HAVE_FFT_SMALL + FLINT_ASSERT(sdctx == NULL || sdctx->N == FLINT_BITS*limbs); +#endif + if (depth <= 6) { trunc = 2*((trunc + 1)/2); diff --git a/src/fft/fft_mfa_truncate_sqrt2_inner.c b/src/fft/fft_mfa_truncate_sqrt2_inner.c index aad2452767..9a691ed568 100644 --- a/src/fft/fft_mfa_truncate_sqrt2_inner.c +++ b/src/fft/fft_mfa_truncate_sqrt2_inner.c @@ -300,6 +300,10 @@ void fft_mfa_truncate_sqrt2_inner_sd_fft(struct sd_fft_mpn_mulmod_2expp1_ctx_str flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** s1, mp_size_t n1, mp_size_t trunc, mp_limb_t ** tt) { +#if FLINT_HAVE_FFT_SMALL + /* the pointwise ring is 2^(n*w) + 1 */ + FLINT_ASSERT(sdctx == NULL || sdctx->N == (slong)(n*w)); +#endif _fft_mfa_truncate_sqrt2_inner_impl(sdctx, ii, jj, n, w, t1, t2, s1, n1, trunc, tt); } diff --git a/src/fft/test/t-convolution_precache_sd_fft.c b/src/fft/test/t-convolution_precache_sd_fft.c index 1ab277b6dc..a982dbedac 100644 --- a/src/fft/test/t-convolution_precache_sd_fft.c +++ b/src/fft/test/t-convolution_precache_sd_fft.c @@ -26,11 +26,11 @@ TEST_FUNCTION_START(fft_convolution_precache_sd_fft, state) { #if FLINT_HAVE_FFT_SMALL slong depth, ti, iter; - static const int threads[] = { 1, 2, 4 }; + static const int threads[] = { 1, 2, 4, 10 }; slong save_threads = flint_get_num_threads(); for (depth = 3; depth <= 8; depth++) - for (ti = 0; ti < 3; ti++) + for (ti = 0; ti < 4; ti++) for (iter = 0; iter < 2; iter++) { mp_size_t n = (UWORD(1) << depth); @@ -126,8 +126,12 @@ TEST_FUNCTION_START(fft_convolution_precache_sd_fft, state) } jj2 = jj; - fft_convolution_precache_sd_fft(C, ii, jj, depth, limbs, trunc, - t1, t2, s1, tt); + if (n_randint(state, 4) == 0) /* NULL context = classic */ + fft_convolution_precache_sd_fft(NULL, ii, jj, depth, limbs, + trunc, t1, t2, s1, tt); + else + fft_convolution_precache_sd_fft(C, ii, jj, depth, limbs, + trunc, t1, t2, s1, tt); fft_convolution_precache(ii2, jj2, depth, limbs, trunc, t1, t2, s1, tt); diff --git a/src/fft/test/t-convolution_sd_fft.c b/src/fft/test/t-convolution_sd_fft.c index 57f4a567c7..1b0f4e8c7c 100644 --- a/src/fft/test/t-convolution_sd_fft.c +++ b/src/fft/test/t-convolution_sd_fft.c @@ -27,11 +27,11 @@ TEST_FUNCTION_START(fft_convolution_sd_fft, state) { #if FLINT_HAVE_FFT_SMALL slong depth, ti, iter; - static const int threads[] = { 1, 2, 4 }; + static const int threads[] = { 1, 2, 4, 10 }; slong save_threads = flint_get_num_threads(); for (depth = 3; depth <= 8; depth++) - for (ti = 0; ti < 3; ti++) + for (ti = 0; ti < 4; ti++) for (iter = 0; iter < 2; iter++) { mp_size_t n = (UWORD(1) << depth); @@ -137,8 +137,12 @@ TEST_FUNCTION_START(fft_convolution_sd_fft, state) } } - fft_convolution_sd_fft(C, ii, jj, depth, limbs, trunc, - t1, t2, s1, tt); + if (n_randint(state, 4) == 0) /* NULL context = classic */ + fft_convolution_sd_fft(NULL, ii, jj, depth, limbs, trunc, + t1, t2, s1, tt); + else + fft_convolution_sd_fft(C, ii, jj, depth, limbs, trunc, + t1, t2, s1, tt); fft_convolution(ii2, jj2, depth, limbs, trunc, t1, t2, s1, tt); @@ -163,6 +167,22 @@ TEST_FUNCTION_START(fft_convolution_sd_fft, state) } flint_set_num_threads(save_threads); + + /* chooser edges: a huge tight size walks the lattice through + divisibility/parity skips and capacity-infeasible digit + sizes (np < 0) before settling; the result must be a valid + ring at least as large as the request */ + { + slong mo, r; + r = fft_adjust_limbs_sd_fft(1000000, 16, &mo); + if (r < 1000000 || (mo != 0 && (64*r) % mo != 0)) + TEST_FUNCTION_FAIL("chooser edge: r = %wd, m = %wd\n", + r, mo); + r = fft_adjust_limbs_sd_fft(300, 1024, &mo); /* divisibility + skips at a large outer length */ + if (r < 300) + TEST_FUNCTION_FAIL("chooser edge 2: r = %wd\n", r); + } #endif TEST_FUNCTION_END(state); diff --git a/src/fft_small/negacyclic.c b/src/fft_small/negacyclic.c index e8b846a363..a424e08ce7 100644 --- a/src/fft_small/negacyclic.c +++ b/src/fft_small/negacyclic.c @@ -420,7 +420,7 @@ sd_fft_mpn_mulmod_2expp1(sd_fft_mpn_mulmod_2expp1_ctx_struct * C, nn_ptr z, nn_s { if (mpn_add_1(z, z, n, 1)) { - memset(z, 0, n * sizeof(ulong)); + flint_mpn_zero(z, n); z[n] = 1; } } @@ -431,11 +431,15 @@ sd_fft_mpn_mulmod_2expp1(sd_fft_mpn_mulmod_2expp1_ctx_struct * C, nn_ptr z, nn_s { if (mpn_sub_1(z, z, n, 1)) { - memset(z, 0, n * sizeof(ulong)); + flint_mpn_zero(z, n); z[n] = 1; } } } + /* result is fully reduced: the unique representative in + [0, 2^N], top limb set only for 2^N itself */ + FLINT_ASSERT(z[n] == 0 + || (z[n] == 1 && flint_mpn_zero_p(z, n))); } } } diff --git a/src/fft_small/test/t-sd_fft_mpn_mulmod_2expp1.c b/src/fft_small/test/t-sd_fft_mpn_mulmod_2expp1.c index 7564de2a85..abd0b2db8f 100644 --- a/src/fft_small/test/t-sd_fft_mpn_mulmod_2expp1.c +++ b/src/fft_small/test/t-sd_fft_mpn_mulmod_2expp1.c @@ -235,7 +235,7 @@ TEST_FUNCTION_START(sd_fft_mpn_mulmod_2expp1, state) z = flint_calloc(n + 1, sizeof(ulong)); r = flint_calloc(n + 1, sizeof(ulong)); - for (cs = 0; cs < 4; cs++) + for (cs = 0; cs < 6; cs++) { memset(x, 0, (n + 1) * sizeof(ulong)); memset(y, 0, (n + 1) * sizeof(ulong)); @@ -256,6 +256,14 @@ TEST_FUNCTION_START(sd_fft_mpn_mulmod_2expp1, state) x[n] = 1; y[n/2] = 1; break; + case 4: /* 2^N * 1 = 2^N */ + x[n] = 1; + y[0] = 1; + break; + case 5: /* 2^N * 3 = 2^N - 2 */ + x[n] = 1; + y[0] = 3; + break; } negmul_ref(r, x, y, n); sd_fft_mpn_mulmod_2expp1(C, z, x, y, S); @@ -273,5 +281,16 @@ TEST_FUNCTION_START(sd_fft_mpn_mulmod_2expp1, state) } } + /* capacity edges of the prime-count chooser: the np = 8 digit + sizes admit no larger count to grow into, so sufficiently + deep transforms are infeasible and the chooser reports so */ + if (sd_fft_mpn_mulmod_2expp1_choose_np(get_default_mpn_ctx(), + 192, 13) != -1 + || sd_fft_mpn_mulmod_2expp1_choose_np(get_default_mpn_ctx(), + 188, 21) != -1 + || sd_fft_mpn_mulmod_2expp1_choose_np(get_default_mpn_ctx(), + 64, 20) < 3) + TEST_FUNCTION_FAIL("choose_np capacity edges\n"); + TEST_FUNCTION_END(state); }