Hybrid Schonhage-Strassen-NTT convolution

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.
+

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.
+

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,29 @@ 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_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,
+                          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 +252,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
 }