diff --git a/configure.ac b/configure.ac index 57eebc19bc..5f30147ae4 100644 --- a/configure.ac +++ b/configure.ac @@ -856,6 +856,7 @@ then AC_SEARCH_LIBS([__gmpn_addlsh1_n_ip1],[gmp],[AC_DEFINE(FLINT_HAVE_NATIVE_mpn_addlsh1_n_ip1,1,Define if GMP has mpn_addlsh1_n_ip1)]) AC_SEARCH_LIBS([__gmpn_rsh1sub_n],[gmp],[AC_DEFINE(FLINT_HAVE_NATIVE_mpn_rsh1sub_n,1,Define if GMP has mpn_rsh1sub_n)]) AC_SEARCH_LIBS([__gmpn_rsh1add_n],[gmp],[AC_DEFINE(FLINT_HAVE_NATIVE_mpn_rsh1add_n,1,Define if GMP has mpn_rsh1add_n)]) + AC_SEARCH_LIBS([__gmpn_mulmid_n],[gmp],[AC_DEFINE(FLINT_HAVE_NATIVE_mpn_mulmid_n,1,Define if GMP has mpn_mulmid_n)]) fi if test "$enable_mpfr_check" = "yes"; diff --git a/doc/source/mpn_extras.rst b/doc/source/mpn_extras.rst index f82e2d5386..2e28997fe8 100644 --- a/doc/source/mpn_extras.rst +++ b/doc/source/mpn_extras.rst @@ -217,6 +217,55 @@ More generally, we can define `n`-limb high products of `m`-limb and The low `n - 1` limbs of the output may be used as scratch space or to write the whole product when this is the fastest method. +Middle product +-------------------------------------------------------------------------------- + +The *windowed middle product* extracts a chosen limb window of a product. For +`\mathrm{an} \ge 1`, `\mathrm{bn} \ge 1` and `0 \le \mathrm{zlo} < \mathrm{zhi} +\le \mathrm{an} + \mathrm{bn}`, it writes `\mathrm{zhi} - \mathrm{zlo}` limbs to +``z`` approximating limbs `[\mathrm{zlo}, \mathrm{zhi})` of `a b`. It is a +*lower approximation*: partial products `a[p] b[q]` with `p + q < \mathrm{zlo}` +are dropped, so the computed value never exceeds the exact window, and the +deficit (a single carry from below `\mathrm{zlo}`) is bounded by +`\min(\mathrm{an}, \mathrm{bn}, \mathrm{zlo}) \cdot 2^{64}`. With +`\mathrm{zlo} = 0` the window is exact. + +.. function:: void flint_mpn_mulmid(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, mp_size_t zlo, mp_size_t zhi) + + Compute the window `[\mathrm{zlo}, \mathrm{zhi})` of `a b`, dispatching to + whichever of the routines below is expected to be fastest for the given + shape. Individual backends may return the exact window or a tighter + approximation than the classical drop; all satisfy the contract above. + +.. function:: void flint_mpn_mulmid_classical(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, mp_size_t zlo, mp_size_t zhi) + + Row-based schoolbook implementation, `O((\mathrm{zhi} - \mathrm{zlo}) \cdot + \min(\mathrm{an}, \mathrm{bn}))`. + +.. function:: void flint_mpn_mulmid_via_mul(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, mp_size_t zlo, mp_size_t zhi) + void flint_mpn_mulmid_via_mullow_n(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, mp_size_t zlo, mp_size_t zhi) + void flint_mpn_mulmid_via_mulhigh_n(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, mp_size_t zlo, mp_size_t zhi) + void flint_mpn_mulmid_via_n_padded(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, mp_size_t zlo, mp_size_t zhi) + void flint_mpn_mulmid_fft_small(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, mp_size_t zlo, mp_size_t zhi) + + Reductions of a general window to, respectively, a full product + (:func:`flint_mpn_mul`), a balanced low product (:func:`flint_mpn_mullow_n`), + a balanced high product (:func:`flint_mpn_mulhigh_n`), a balanced middle + product (:func:`flint_mpn_mulmid_n`) and the small-prime FFT. Each is valid + for arbitrary input by padding internally, but is only economical in its own + regime; :func:`flint_mpn_mulmid` chooses between them. + +.. function:: void flint_mpn_mulmid_n(mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) + + Exact balanced middle product of `\{ap, 2n-1\}` and `\{bp, n\}`, writing + `n + 2` limbs: the high `n` limbs are exact and the low two are guard limbs + (they lack the carry into the band from below). This is a wrapper around + GMP's ``mpn_mulmid_n`` and is only defined when + ``FLINT_HAVE_NATIVE_mpn_mulmid_n`` is set (that is, when the build may call + GMP internals and GMP exports the symbol; see ``configure``). When it is + unavailable, :func:`flint_mpn_mulmid_via_n_padded` is likewise unavailable + and :func:`flint_mpn_mulmid` uses its other methods. + Divisibility -------------------------------------------------------------------------------- diff --git a/src/fft_small.h b/src/fft_small.h index fbd55e1e7e..bd508a40a2 100644 --- a/src/fft_small.h +++ b/src/fft_small.h @@ -397,6 +397,7 @@ void mpn_ctx_init(mpn_ctx_t R, ulong p); void mpn_ctx_clear(mpn_ctx_t R); void* mpn_ctx_fit_buffer(mpn_ctx_t R, ulong n); void mpn_ctx_mpn_mul(mpn_ctx_t R, ulong* z, const ulong* a, ulong an, const ulong* b, ulong bn); +void _mpn_ctx_mpn_mul_range(mpn_ctx_t R, ulong* z, ulong lo, ulong hi, const ulong* a, ulong an, const ulong* b, ulong bn); void _nmod_poly_mul_mid_mpn_ctx( ulong* z, ulong zl, ulong zh, diff --git a/src/fft_small/mpn_mul.c b/src/fft_small/mpn_mul.c index bca664d324..2accab4395 100644 --- a/src/fft_small/mpn_mul.c +++ b/src/fft_small/mpn_mul.c @@ -19,6 +19,28 @@ #include "crt_helpers.h" #include "fft_small.h" +/* + If set, the truncated/wrap-around FFT optimization for partial products is + enabled: when computing only the limb window [lo, hi) of a product (a middle + or high product) we are free to fold the unused low coefficients [0, c_lo) + on top of the high ones with a cyclic (x^(2^d) - 1) transform, provided the + resulting aliasing never lands inside the window we actually read. This can + roughly halve the transform length for balanced middle products. Guarded by + an #if so its effect can be measured by toggling it to 0 (in which case the + full-length transform is used and only the output extraction is windowed). +*/ +#ifndef MPN_MUL_USE_WRAPAROUND +#define MPN_MUL_USE_WRAPAROUND 1 +#endif + +/* + Upper bound (in limbs) on the small scratch buffer used by _mpn_from_ffts to + reconstruct the "bottom band" of coefficients that straddle the limb 'lo'. + The band spans at most one BLK_SZ block plus the n+1 limb overlap of a single + coefficient; over all profiles this is < 800 limbs (checked with an assert). +*/ +#define MPN_MULMID_BOUNDBUF 1024 + /* The following profiles are hardcoded. @@ -707,12 +729,12 @@ DEFINE_IT(7, 8) #undef DEFINE_IT typedef void (*from_ffts_func)( - ulong* z, ulong zn, ulong zlen, + 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* overhang, ulong* boundbuf); /* The "n" here is the limb count Rcrts[np-1].coeff_len, which is big enough @@ -720,47 +742,65 @@ typedef void (*from_ffts_func)( intermediate dot products f[0]*x[0] + ... + f[np-1]*x[np-1]. The x[i] are single limb and the f[i] are of length "m". The number of primes is "np". - The coefficient of X^i, 0 <= i < zlen needs to be reconstructed and added - to the answer mpn (z, zn). This involves the limbs + This is the windowed reconstruction. The full integer answer is - z[floor(i*bits/64)] ... z[floor(i*bits/64)+n] + sum_{i=0}^{clen-1} coeff(X^i) * 2^(i*bits) - so is easy if floor(i*bits/64)+n < zn. + but only the limb window [lo, hi) of it is requested (a lower approximation: + carries coming from limbs below 'lo' are not recovered). The output buffer + (z, hi-lo) holds limbs [lo, hi); writing the coefficient of X^i, which spans - The the l^th fft ctx Rffts[l] is expected to have data at d + l*dstride + (z - lo)[floor(i*bits/64)] ... (z - lo)[floor(i*bits/64) + n] - if overhang = NULL + is "easy" when floor(i*bits/64) + n < hi. - handle output coefficients from [start_easy, zlen) - end_easy is still expected to be valid + Coefficients are partitioned into three groups by the caller: - if overhang != NULL + [c_lo, start_easy) bottom band: coefficients that reach into limb >= lo + but lie below the BLK_SZ-aligned easy start. Only the + worker passing boundbuf != NULL (thread 0) handles + these, reconstructing them -- including the sub-lo + limbs that carry into limb lo -- into boundbuf anchored + at limb L0 = floor(c_lo*bits/64). The driver adds the + in-range slice of boundbuf into z after all threads + join (kept serial to avoid a carry-propagation race). - overhang has space for n words + [start_easy, stop_easy) easy interior: start_easy and stop_easy are + divisible by BLK_SZ. This is the byte-identical hot + loop of the full product, just writing through the + shifted base (z - lo). - handle output coefficients from [start_easy, end_easy) where - start_easy and stop_easy are divisible by BLK_SZ + [stop_easy, clen) handled either via the overhang (parallel segments) or, + for the last segment (overhang == NULL), the hard tail + clamped at limb 'hi'. - write to output words - [start_easy*bits/64, stop_easy*bits/64) [overhang+0, overhang+n) + The l^th fft ctx Rffts[l] is expected to have data at d + l*dstride. + + if overhang == NULL + handle output coefficients from [start_easy, clen), clamped above at hi + if overhang != NULL + overhang has space for n words; handle [start_easy, stop_easy) and write + the spill of the topmost block into [overhang+0, overhang+n) */ #define DEFINE_IT(NP, N, M) \ static void CAT(_mpn_from_ffts, NP)( \ - ulong* z, ulong zn, ulong zlen, \ + 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* overhang, ulong* boundbuf) \ { \ ulong np = NP; \ ulong n = N; \ ulong m = M; \ + ulong* zbase = z - lo; /* an absolute limb 'toff' is written at z[toff - lo] */ \ ulong zn_start = start_easy*bits/64; \ - ulong zn_stop = (overhang == NULL) ? zn : stop_easy*bits/64; \ + ulong zn_stop = (overhang == NULL) ? hi : stop_easy*bits/64; \ \ FLINT_ASSERT(n == Rcrts[np-1].coeff_len); \ FLINT_ASSERT(start_easy <= stop_easy); \ + FLINT_ASSERT(zn_start >= lo); \ \ if (n == m + 1) \ { \ @@ -772,10 +812,47 @@ static void CAT(_mpn_from_ffts, NP)( \ { \ FLINT_ASSERT(n == m); \ } \ - \ - memset(z + zn_start, 0, (zn_stop - zn_start)*sizeof(ulong)); \ \ ulong Xs[BLK_SZ*NP]; \ + \ + /* bottom band: coefficients [c_lo, start_easy) -> boundbuf (thread 0 only) */ \ + if (boundbuf != NULL && c_lo < start_easy) \ + { \ + ulong L0 = (c_lo*bits)/64; \ + ulong topabs = ((start_easy - 1)*bits)/64 + n; \ + ulong tw = topabs - L0 + 1; \ + ulong* bbase = boundbuf - L0; \ + \ + FLINT_ASSERT(tw <= MPN_MULMID_BOUNDBUF); \ + for (ulong k = 0; k < tw; k++) \ + boundbuf[k] = 0; \ + \ + for (ulong blk = c_lo/BLK_SZ; blk < start_easy/BLK_SZ; blk++) \ + { \ + _convert_block(Xs, Rffts, d, dstride, np, blk); \ + ulong j0 = (blk*BLK_SZ < c_lo) ? c_lo - blk*BLK_SZ : 0; \ + for (ulong j = j0; j < BLK_SZ; j += 1) \ + { \ + ulong i = blk*BLK_SZ + j; \ + ulong r[N + 1]; \ + ulong t[N + 1]; \ + ulong l = 0; \ + \ + CAT3(_big_mul, N, M)(r, t, _crt_data_co_prime(Rcrts + np - 1, l, n), Xs[l*BLK_SZ + j]); \ + for (l++; l < np; l++) \ + CAT3(_big_addmul, N, M)(r, t, _crt_data_co_prime(Rcrts + np - 1, l, n), Xs[l*BLK_SZ + j]); \ + \ + CAT(_reduce_big_sum, N)(r, t, crt_data_prod_primes(Rcrts + np - 1)); \ + \ + ulong toff = (i*bits)/FLINT_BITS; \ + ulong tshift = (i*bits)%FLINT_BITS; \ + \ + CAT(_add_to_answer_easy, N)(bbase, r, topabs + 1, toff, tshift); \ + } \ + } \ + } \ + \ + memset(zbase + zn_start, 0, (zn_stop - zn_start)*sizeof(ulong)); \ \ if (overhang != NULL) \ { \ @@ -807,7 +884,7 @@ static void CAT(_mpn_from_ffts, NP)( \ \ FLINT_ASSERT(zn_stop > n + toff); \ \ - CAT(_add_to_answer_easy, N)(z, r, zn_stop, toff, tshift); \ + CAT(_add_to_answer_easy, N)(zbase, r, zn_stop, toff, tshift); \ } \ } \ \ @@ -833,7 +910,7 @@ static void CAT(_mpn_from_ffts, NP)( \ \ if (n + toff < zn_stop) \ { \ - CAT(_add_to_answer_easy, N)(z, r, zn_stop, toff, tshift); \ + CAT(_add_to_answer_easy, N)(zbase, r, zn_stop, toff, tshift); \ } \ else \ { \ @@ -853,7 +930,7 @@ static void CAT(_mpn_from_ffts, NP)( \ unsigned char cf = 0; \ ulong k = 0; \ for (; k < zn_stop - toff; k++) \ - cf = _addcarry_ulong(cf, z[toff + k], r[k], &z[toff + k]); \ + cf = _addcarry_ulong(cf, zbase[toff + k], r[k], &zbase[toff + k]); \ for (; k <= n; k++) \ cf = _addcarry_ulong(cf, overhang[k-(zn_stop-toff)], r[k], &overhang[k-(zn_stop-toff)]); \ } \ @@ -861,7 +938,7 @@ static void CAT(_mpn_from_ffts, NP)( \ } \ else \ { \ - for (ulong i = stop_easy; i < zlen; i++) \ + for (ulong i = stop_easy; i < clen; i++) \ { \ ulong r[N + 1]; \ ulong t[N + 1]; \ @@ -882,10 +959,10 @@ static void CAT(_mpn_from_ffts, NP)( \ ulong toff = (i*bits)/FLINT_BITS; \ ulong tshift = (i*bits)%FLINT_BITS; \ \ - if (toff >= zn) \ + if (toff >= hi) \ break; \ \ - CAT(_add_to_answer_hard, N)(z, r, zn, toff, tshift); \ + CAT(_add_to_answer_hard, N)(zbase, r, hi, toff, tshift); \ } \ } \ } @@ -1419,8 +1496,10 @@ static void mod_fft_worker_func(void* varg) typedef struct { from_ffts_func from_ffts; ulong* z; - ulong zn; - ulong zlen; + ulong lo; + ulong hi; + ulong c_lo; + ulong clen; sd_fft_ctx_struct* fctxs; double* abuf; ulong stride; @@ -1429,6 +1508,7 @@ typedef struct { ulong start_easy; ulong stop_easy; ulong* overhang; + ulong* boundbuf; ulong overhang_buffer[MPN_CTX_NCRTS]; } crt_worker_struct; @@ -1436,20 +1516,54 @@ static void crt_worker_func(void* varg) { crt_worker_struct* X = (crt_worker_struct*) varg; - X->from_ffts(X->z, X->zn, X->zlen, X->fctxs, X->abuf, X->stride, X->crts, - X->bits, X->start_easy, X->stop_easy, X->overhang); + X->from_ffts(X->z, X->lo, X->hi, X->c_lo, X->clen, X->fctxs, X->abuf, + X->stride, X->crts, X->bits, X->start_easy, X->stop_easy, + X->overhang, X->boundbuf); } -void mpn_ctx_mpn_mul(mpn_ctx_t R, ulong* z, const ulong* a, ulong an, const ulong* b, ulong bn) +/* + Compute the limb window [lo, hi) of the integer product a*b, as a lower + approximation: carries propagating up from limbs strictly below 'lo' are not + recovered, exactly as in the radix middle-product code. The low-end deficit + (true - computed) of limb 'lo' is < min(an, bn, lo) coefficients' worth, i.e. + bounded by min(an, bn, lo)*2^64; limbs sufficiently far above 'lo' are exact + (until the upper truncation at 'hi'). Writes hi - lo limbs to z. + + With lo == 0 and hi == an + bn this reduces to the full product computed by + mpn_ctx_mpn_mul (every code path below collapses to the original one). +*/ +void _mpn_ctx_mpn_mul_range(mpn_ctx_t R, ulong* z, ulong lo, ulong hi, + const ulong* a, ulong an, const ulong* b, ulong bn) { ulong zn, alen, blen, zlen, atrunc, btrunc, ztrunc, depth, stride; + ulong c_lo, c_hi, nn; double* abuf; profile_entry P; ulong sz; void* worker_struct_buffer; int squaring; + FLINT_ASSERT(an > 0); + FLINT_ASSERT(bn > 0); + + zn = an + bn; + + if (lo >= hi) + return; + + /* limbs at or above the top of the product are zero */ + if (hi > zn) + { + if (lo >= zn) + { + flint_mpn_zero(z, hi - lo); + return; + } + flint_mpn_zero(z + (zn - lo), hi - zn); + hi = zn; + } + mpn_ctx_best_profile(R, &P, an, bn); sz = sizeof(mod_worker_struct)*P.nthreads; @@ -1459,33 +1573,67 @@ void mpn_ctx_mpn_mul(mpn_ctx_t R, ulong* z, const ulong* a, ulong an, const ulon worker_struct_buffer = flint_malloc(sz); squaring = (a == b) && (an == bn); - zn = an + bn; alen = n_cdiv(FLINT_BITS*an, P.bits); blen = n_cdiv(FLINT_BITS*bn, P.bits); zlen = alen + blen - 1; atrunc = n_round_up(alen, BLK_SZ); btrunc = n_round_up(blen, BLK_SZ); + + nn = R->crts[P.np - 1].coeff_len; + + /* first coefficient whose support reaches into limb >= lo */ + c_lo = (lo > nn) ? n_min(zlen, n_cdiv((lo - nn)*FLINT_BITS, P.bits)) : 0; + /* first coefficient lying entirely at limb >= hi */ + c_hi = n_min(zlen, n_cdiv(hi*FLINT_BITS, P.bits)); + + if (c_lo >= c_hi) + { + /* no coefficient touches the window: it is all zero */ + flint_mpn_zero(z, hi - lo); + flint_free(worker_struct_buffer); + flint_give_back_threads(P.handles, P.nhandles); + return; + } + + /* + Transform-size selection. Baseline: the full-length transform, with the + output extraction windowed to [c_lo, c_hi). The wrap-around trick (guard) + may shrink the transform to a power of two 2^d when the unused low + coefficients can be aliased on top of the high ones without disturbing + [c_lo, c_hi): need max(atrunc, btrunc, c_hi) <= 2^d <= zlen and the + wrap-around of the top, zlen - 2^d, to land at or below c_lo. + */ ztrunc = n_round_up(zlen, BLK_SZ); depth = n_max(LG_BLK_SZ, n_clog2(ztrunc)); +#if MPN_MUL_USE_WRAPAROUND + { + /* + Largest power of two w = 2^d with w <= zlen. Note that this file's + n_flog2(x) returns nbits(x) = floor(log2 x) + 1 (not floor(log2 x)), + so n_pow2(n_flog2(zlen)) is the power of two *above* zlen; halve it + to land at or below zlen. + */ + ulong d = n_flog2(zlen); + ulong w = n_pow2(d); + if (w > zlen) + { + w >>= 1; + d -= 1; + } + if (d >= LG_BLK_SZ && + atrunc <= w && btrunc <= w && c_hi <= w && w <= zlen && zlen <= c_lo + w) + { + depth = d; + ztrunc = w; + } + } +#endif stride = n_round_up(sd_fft_ctx_data_size(depth), 128); - FLINT_ASSERT(an > 0); - FLINT_ASSERT(bn > 0); FLINT_ASSERT(0 <= flint_mpn_cmp_ui_2exp( crt_data_prod_primes(R->crts + P.np - 1), R->crts[P.np - 1].coeff_len, blen, 2*P.bits)); -#define TIME_THIS 0 - -#if TIME_THIS -timeit_t timer, timer_overall; -flint_printf("------------ zn = %wu, nthreads = %wu np = %wu, bits = %wu, -------------\n", zn, nthreads, np, bits); -#endif - -#if TIME_THIS -timeit_start(timer_overall); -#endif - if (P.to_ffts != NULL) { ulong bits = P.bits; @@ -1505,10 +1653,6 @@ timeit_start(timer_overall); abuf = (double*) mpn_ctx_fit_buffer(R, 2*P.np*stride*sizeof(double)); bbuf = abuf + P.np*stride; -#if TIME_THIS -timeit_start(timer); -#endif - /* some fixups for loop unrollings: round down the easy stops */ FLINT_ASSERT(bits%2 == 0); a_stop_easy &= -rounding; @@ -1550,31 +1694,6 @@ timeit_start(timer); for (slong i = P.nhandles; i > 0; i--) thread_pool_wait(global_thread_pool, P.handles[i - 1]); -#if TIME_THIS -timeit_stop(timer); -if (timer->wall > 50) -flint_printf(" mod: %wd\n", timer->wall); -#endif - -#if TIME_THIS -timeit_start(timer); -#endif - - /* - current scheduling: - np = 5, nthreads = 3: - thread0: p0, p3 - thread1: p1, p4 - thread2: p2 - - np = 3, nthreads = 5: - thread0: p0 - thread1: p1 - thread2: p2 - thread3: - - thread4: - - */ - wf = (fft_worker_struct*) worker_struct_buffer; for (ulong l = 0; l < P.np; l++) @@ -1599,12 +1718,6 @@ timeit_start(timer); for (ulong i = n_min(P.nhandles, P.np - 1); i > 0; i--) thread_pool_wait(global_thread_pool, P.handles[i - 1]); - -#if TIME_THIS -timeit_stop(timer); -if (timer->wall > 50) -flint_printf(" fft: %wd\n", timer->wall); -#endif } else { @@ -1617,9 +1730,6 @@ flint_printf(" fft: %wd\n", timer->wall); abuf = (double*) mpn_ctx_fit_buffer(R, (np+nthreads)*stride*sizeof(double)); bbuf = abuf + np*stride; -#if TIME_THIS -timeit_start(timer); -#endif for (ulong l = 0; l < np; l++) { mod_fft_worker_struct* X = w + l; @@ -1648,29 +1758,24 @@ timeit_start(timer); for (ulong i = n_min(P.nhandles, P.np - 1); i > 0; i--) thread_pool_wait(global_thread_pool, P.handles[i - 1]); - -#if TIME_THIS -timeit_stop(timer); -if (timer->wall > 50) -flint_printf("mod+fft: %wd\n", timer->wall); -#endif } -#if TIME_THIS -timeit_start(timer); -#endif - { + ulong bits = P.bits; ulong n = R->crts[P.np-1].coeff_len; crt_worker_struct* w = (crt_worker_struct*) worker_struct_buffer; ulong nthreads = P.nthreads; - ulong end_easy = (zn >= n+1 ? zn - (n+1) : UWORD(0))*64/P.bits; - - /* this is how must space was statically allocated in each struct */ - FLINT_ASSERT(n <= MPN_CTX_NCRTS); - end_easy &= -BLK_SZ; + /* BLK_SZ-aligned easy interval [E0, E1) of coefficients whose whole */ + /* span lies inside the window [lo, hi) */ + ulong E0 = n_round_up(n_cdiv(lo*FLINT_BITS, bits), BLK_SZ); + ulong E1 = ((hi >= n + 1) ? hi - (n + 1) : UWORD(0))*FLINT_BITS/bits; + E1 &= -BLK_SZ; + if (E1 < E0) + E1 = E0; + /* this is how much space was statically allocated in each struct */ + FLINT_ASSERT(n <= MPN_CTX_NCRTS); FLINT_ASSERT(4 <= P.np && P.np <= 8); static from_ffts_func tab[8-4+1] = {_mpn_from_ffts_4, _mpn_from_ffts_5, @@ -1678,71 +1783,110 @@ timeit_start(timer); _mpn_from_ffts_7, _mpn_from_ffts_8}; - for (ulong l = 0; l < nthreads; l++) + if (E1 <= E0) { - crt_worker_struct* X = w + l; - X->from_ffts = tab[P.np - 4]; - X->z = z; - X->zn = zn; - X->zlen = zlen; - X->fctxs = R->ffts; - X->abuf = abuf; - X->stride = stride; - X->crts = R->crts; - X->bits = P.bits; - X->start_easy = n_round_up((l+0)*end_easy/nthreads, BLK_SZ); - X->stop_easy = n_round_up((l+1)*end_easy/nthreads, BLK_SZ); - X->overhang = (l + 1 == nthreads) ? NULL : X->overhang_buffer; + /* + Window too small for a BLK_SZ-aligned interior: reconstruct the + whole of [c_lo, c_hi) serially into a scratch anchored at limb + L0 = floor(c_lo*bits/64) (so that no coefficient straddles its + low end), then copy out the in-range slice. + */ + ulong L0 = (c_lo*bits)/64; + ulong tlen = hi - L0; + ulong* tmp = FLINT_ARRAY_ALLOC(tlen, ulong); + + tab[P.np - 4](tmp, L0, hi, c_lo, c_hi, R->ffts, abuf, stride, + R->crts, bits, c_lo, c_lo, NULL, NULL); + + flint_mpn_copyi(z, tmp + (lo - L0), hi - lo); + flint_free(tmp); } + else + { + ulong L0 = (c_lo*bits)/64; + ulong span = E1 - E0; + ulong gap = E0*bits/64; /* absolute limb of E0 (>= lo) */ + ulong boundbuf[MPN_MULMID_BOUNDBUF]; - for (slong i = P.nhandles; i > 0; i--) - thread_pool_wake(global_thread_pool, P.handles[i - 1], 0, - crt_worker_func, w + i); - crt_worker_func(w + 0); - - for (slong i = P.nhandles; i > 0; i--) - thread_pool_wait(global_thread_pool, P.handles[i - 1]); + /* limbs [lo, gap) get only the bottom band (or stay zero) */ + if (gap > lo) + flint_mpn_zero(z, gap - lo); - unsigned char cf = 0; - for (slong i = 1; i <= P.nhandles; i++) - { - ulong start = w[i].start_easy*P.bits/64; - if (i == P.nhandles) + for (ulong l = 0; l < nthreads; l++) { - cf = flint_mpn_add_inplace_c(z + start, zn - start, - w[i - 1].overhang_buffer, n, cf); + crt_worker_struct* X = w + l; + X->from_ffts = tab[P.np - 4]; + X->z = z; + X->lo = lo; + X->hi = hi; + X->c_lo = c_lo; + X->clen = c_hi; + X->fctxs = R->ffts; + X->abuf = abuf; + X->stride = stride; + X->crts = R->crts; + X->bits = bits; + X->start_easy = E0 + n_round_up((l+0)*span/nthreads, BLK_SZ); + X->stop_easy = E0 + n_round_up((l+1)*span/nthreads, BLK_SZ); + X->overhang = (l + 1 == nthreads) ? NULL : X->overhang_buffer; + X->boundbuf = (l == 0 && c_lo < E0) ? boundbuf : NULL; } - else + + for (slong i = P.nhandles; i > 0; i--) + thread_pool_wake(global_thread_pool, P.handles[i - 1], 0, + crt_worker_func, w + i); + crt_worker_func(w + 0); + + for (slong i = P.nhandles; i > 0; i--) + thread_pool_wait(global_thread_pool, P.handles[i - 1]); + + /* stitch the per-segment overhangs (carries across boundaries) */ { - ulong stop = w[i].stop_easy*P.bits/64; - if (stop > start) + unsigned char cf = 0; + for (slong i = 1; i <= P.nhandles; i++) { - cf = flint_mpn_add_inplace_c(z + start, stop - start, - w[i - 1].overhang_buffer, n, cf); - } - else - { - for (ulong k = 0; k < n; k++) + ulong start = w[i].start_easy*bits/64; + if (i == P.nhandles) + { + cf = flint_mpn_add_inplace_c(z + (start - lo), hi - start, + w[i - 1].overhang_buffer, n, cf); + } + else { - FLINT_ASSERT(w[i].overhang_buffer[k] == 0); - w[i].overhang_buffer[k] = w[i - 1].overhang_buffer[k]; + ulong stop = w[i].stop_easy*bits/64; + if (stop > start) + { + cf = flint_mpn_add_inplace_c(z + (start - lo), stop - start, + w[i - 1].overhang_buffer, n, cf); + } + else + { + for (ulong k = 0; k < n; k++) + { + FLINT_ASSERT(w[i].overhang_buffer[k] == 0); + w[i].overhang_buffer[k] = w[i - 1].overhang_buffer[k]; + } + } } } } + + /* add the bottom band (limbs >= lo of coefficients [c_lo, E0)) */ + if (c_lo < E0) + { + ulong topabs = ((E0 - 1)*bits)/64 + n; + ulong tw = topabs - L0 + 1; + flint_mpn_add_inplace_c(z, hi - lo, + boundbuf + (lo - L0), tw - (lo - L0), 0); + } } } -#if TIME_THIS -timeit_stop(timer); -if (timer->wall > 50) -flint_printf(" crt: %wd\n", timer->wall); -timeit_stop(timer_overall); -if (timer_overall->wall > 50) -flint_printf(" +: %wd\n", timer_overall->wall); -#endif - -#undef TIME_THIS - flint_free(worker_struct_buffer); flint_give_back_threads(P.handles, P.nhandles); } + +void mpn_ctx_mpn_mul(mpn_ctx_t R, ulong* z, const ulong* a, ulong an, const ulong* b, ulong bn) +{ + _mpn_ctx_mpn_mul_range(R, z, 0, an + bn, a, an, b, bn); +} diff --git a/src/fft_small/test/main.c b/src/fft_small/test/main.c index 0a36f12735..e7b2c17ac2 100644 --- a/src/fft_small/test/main.c +++ b/src/fft_small/test/main.c @@ -13,6 +13,7 @@ #include "t-fmpz_poly_mul.c" #include "t-mpn_add_inplace_c.c" +#include "t-mpn_mulmid.c" #include "t-mul.c" #include "t-nmod_poly_divrem.c" #include "t-nmod_poly_mul.c" @@ -24,6 +25,7 @@ test_struct tests[] = { TEST_FUNCTION(_fmpz_poly_mul_mid_mpn_ctx), TEST_FUNCTION(flint_mpn_add_inplace_c), + TEST_FUNCTION(_mpn_ctx_mpn_mul_range), TEST_FUNCTION(mpn_ctx_mpn_mul), TEST_FUNCTION(_nmod_poly_divrem_mpn_ctx), TEST_FUNCTION(_nmod_poly_mul_mid_mpn_ctx), diff --git a/src/fft_small/test/t-mpn_mulmid.c b/src/fft_small/test/t-mpn_mulmid.c new file mode 100644 index 0000000000..2ea5bec541 --- /dev/null +++ b/src/fft_small/test/t-mpn_mulmid.c @@ -0,0 +1,170 @@ +/* + Copyright (C) 2024 Fredrik Johansson + + 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 "test_helpers.h" +#include "ulong_extras.h" +#include "fmpz.h" +#include "fft_small.h" +#include "machine_vectors.h" + +/* fail printer */ +static void +mulmid_fail(const char* msg, ulong an, ulong bn, ulong lo, ulong hi) +{ + flint_printf("\nFAILED: %s\n", msg); + flint_printf("an = %wu, bn = %wu, lo = %wu, hi = %wu\n", an, bn, lo, hi); + fflush(stdout); + flint_abort(); +} + +/* + Check the limb window [lo, hi) of a*b returned by _mpn_ctx_mpn_mul_range. + + 'full' must already hold _mpn_ctx_mpn_mul_range(R, full, lo, an+bn, a, an, b, bn), + i.e. the same window but with the maximal top, so that the low limbs of the + requested window can be matched against it without any high-truncation + ambiguity. +*/ +static void +check_window(mpn_ctx_t R, const ulong* a, ulong an, const ulong* b, ulong bn, + const ulong* full, ulong lo, ulong hi, ulong* scratch) +{ + ulong zn = an + bn; + + _mpn_ctx_mpn_mul_range(R, scratch, lo, hi, a, an, b, bn); + + /* limbs [lo, min(hi, zn)) must agree with the maximal-top run */ + for (ulong k = lo; k < hi; k++) + { + ulong got = scratch[k - lo]; + ulong exp = (k < zn) ? full[k - lo] : UWORD(0); + if (got != exp) + mulmid_fail("window disagrees with full-top run / zero pad", an, bn, lo, hi); + } +} + +void test_mulmid(mpn_ctx_t R, ulong maxsize, ulong nreps, flint_rand_t state) +{ + ulong *a, *b, *c, *full, *scratch; + fmpz_t fa, fb, fv, fvlo, fd, fdef, fbound; + + a = FLINT_ARRAY_ALLOC(maxsize, ulong); + b = FLINT_ARRAY_ALLOC(maxsize, ulong); + c = FLINT_ARRAY_ALLOC(maxsize, ulong); + full = FLINT_ARRAY_ALLOC(maxsize + 2, ulong); + scratch = FLINT_ARRAY_ALLOC(maxsize + 2, ulong); + + fmpz_init(fa); + fmpz_init(fb); + fmpz_init(fv); + fmpz_init(fvlo); + fmpz_init(fd); + fmpz_init(fdef); + fmpz_init(fbound); + + for (ulong rep = 0; rep < nreps; rep++) + { + ulong an = 2 + n_randint(state, maxsize - 4); + ulong bn = 1 + n_randint(state, n_min(an, maxsize - an)); + ulong zn = an + bn; + ulong lo, hi; + + for (ulong i = 0; i < maxsize; i++) + { + a[i] = n_randlimb(state); + b[i] = n_randlimb(state); + } + + /* a random low cut, biased to also hit 0 and the extremes */ + switch (n_randint(state, 4)) + { + case 0: lo = 0; break; + case 1: lo = n_randint(state, zn); break; + case 2: lo = (zn > 4) ? zn - 1 - n_randint(state, 4) : 0; break; + default: lo = n_randint(state, n_max(1, zn/2)); break; + } + + /* exact product and exact integers */ + mpn_mul(c, a, an, b, bn); + fmpz_set_ui_array(fa, a, an); + fmpz_set_ui_array(fb, b, bn); + fmpz_mul(fv, fa, fb); + + /* maximal-top window [lo, zn): a lower approximation of floor(V/2^(64 lo)) */ + _mpn_ctx_mpn_mul_range(R, full, lo, zn, a, an, b, bn); + + /* --- lower-bound + deficit-bound, as full integers (no truncation) --- */ + fmpz_fdiv_q_2exp(fvlo, fv, FLINT_BITS*lo); /* floor(V / 2^(64 lo)) */ + fmpz_set_ui_array(fd, full, zn - lo); /* our approximation */ + fmpz_sub(fdef, fvlo, fd); /* deficit = exact - ours */ + + if (fmpz_sgn(fdef) < 0) + mulmid_fail("not a lower approximation", an, bn, lo, zn); + + /* deficit must be <= min(an, bn, lo) * 2^64 */ + fmpz_set_ui(fbound, n_min(n_min(an, bn), lo)); + fmpz_mul_2exp(fbound, fbound, FLINT_BITS); + if (fmpz_cmp(fdef, fbound) > 0) + mulmid_fail("deficit exceeds min(an,bn,lo)*2^64", an, bn, lo, zn); + + /* with lo == 0 there is no dropped low part: the window is exact */ + if (lo == 0 && !fmpz_is_zero(fdef)) + mulmid_fail("low product must be exact (lo == 0)", an, bn, lo, zn); + + /* --- high-truncation consistency against the maximal-top run --- */ + for (int t = 0; t < 3; t++) + { + switch (n_randint(state, 4)) + { + case 0: hi = lo + 1; break; /* tiny window */ + case 1: hi = lo + 1 + n_randint(state, zn - lo); /* random */ + break; + case 2: hi = zn + 1 + n_randint(state, 3); break; /* past the top */ + default: hi = zn; break; + } + + check_window(R, a, an, b, bn, full, lo, hi, scratch); + } + } + + fmpz_clear(fa); + fmpz_clear(fb); + fmpz_clear(fv); + fmpz_clear(fvlo); + fmpz_clear(fd); + fmpz_clear(fdef); + fmpz_clear(fbound); + + flint_free(a); + flint_free(b); + flint_free(c); + flint_free(full); + flint_free(scratch); +} + +TEST_FUNCTION_START(_mpn_ctx_mpn_mul_range, state) +{ + { + mpn_ctx_t R; + mpn_ctx_init(R, UWORD(0x0003f00000000001)); + test_mulmid(R, 1000, 200 * flint_test_multiplier(), state); + test_mulmid(R, 10000, 40 * flint_test_multiplier(), state); + test_mulmid(R, 50000, 4 * flint_test_multiplier(), state); + + /* exercise the threading paths */ + flint_set_num_threads(1 + n_randint(state, 8)); + test_mulmid(R, 20000, 10 * flint_test_multiplier(), state); + + mpn_ctx_clear(R); + } + + TEST_FUNCTION_END(state); +} diff --git a/src/flint-config.h.in b/src/flint-config.h.in index 72740166a9..0eff278e3d 100644 --- a/src/flint-config.h.in +++ b/src/flint-config.h.in @@ -42,6 +42,9 @@ /* Define if system has mpn_modexact_1_odd */ #undef FLINT_HAVE_NATIVE_mpn_modexact_1_odd +/* Define if GMP has mpn_mulmid_n */ +#undef FLINT_HAVE_NATIVE_mpn_mulmid_n + /* Define if GMP has mpn_rsh1add_n */ #undef FLINT_HAVE_NATIVE_mpn_rsh1add_n diff --git a/src/mpn_extras.h b/src/mpn_extras.h index 94d906c4a8..aad08408c1 100644 --- a/src/mpn_extras.h +++ b/src/mpn_extras.h @@ -606,7 +606,7 @@ flint_mpn_sqr(mp_ptr r, mp_srcptr x, mp_size_t n) flint_mpn_mul((_z), (_y), (_yn), (_x), (_xn)); \ } -/* High and low multiplication *******************************************************/ +/* High, low and middle multiplication *******************************************************/ #define FLINT_HAVE_MULLOW_FUNC(n) ((n) <= FLINT_MPN_MULLOW_FUNC_TAB_WIDTH) #define FLINT_HAVE_MULHIGH_FUNC(n) ((n) <= FLINT_MPN_MULHIGH_FUNC_TAB_WIDTH) @@ -666,8 +666,41 @@ FLINT_DLL extern const flint_mpn_sqrhigh_normalised_func_t flint_mpn_sqrhigh_nor #define FLINT_MPN_MULLOW_MULDERS_CUTOFF 50 #define FLINT_MPN_MULHIGH_MULDERS_CUTOFF 40 #define FLINT_MPN_MULHIGH_MUL_CUTOFF 2000 +#define FLINT_MPN_MULHIGH_FFT_SMALL_CUTOFF 1300 #define FLINT_MPN_MULHIGH_K_TAB_SIZE 2048 +#define FLINT_MPN_SQRHIGH_MULDERS_CUTOFF 90 +#define FLINT_MPN_SQRHIGH_FFT_SMALL_CUTOFF 1500 +#define FLINT_MPN_SQRHIGH_SQR_CUTOFF 2000 +#define FLINT_MPN_SQRHIGH_K_TAB_SIZE 2048 + + +/* Windowed middle product ***************************************************/ + +void flint_mpn_mulmid(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, mp_size_t zlo, mp_size_t zhi); +void flint_mpn_mulmid_classical(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, mp_size_t zlo, mp_size_t zhi); +void flint_mpn_mulmid_via_mul(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, mp_size_t zlo, mp_size_t zhi); +void flint_mpn_mulmid_via_mullow_n(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, mp_size_t zlo, mp_size_t zhi); +void flint_mpn_mulmid_via_mulhigh_n(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, mp_size_t zlo, mp_size_t zhi); + +#if FLINT_HAVE_FFT_SMALL +void flint_mpn_mulmid_fft_small(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, mp_size_t zlo, mp_size_t zhi); +mp_limb_t _flint_mpn_mulhigh_n_fft_small(mp_ptr res, mp_srcptr u, mp_srcptr v, mp_size_t n); +mp_limb_t _flint_mpn_mullow_n_fft_small(mp_ptr res, mp_srcptr u, mp_srcptr v, mp_size_t n); +#endif + +#if FLINT_HAVE_NATIVE_mpn_mulmid_n +void __gmpn_mulmid_n(mp_ptr, mp_srcptr, mp_srcptr, mp_size_t); + +MPN_EXTRAS_INLINE +void flint_mpn_mulmid_n(mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) +{ + __gmpn_mulmid_n(rp, ap, bp, n); +} + +void flint_mpn_mulmid_via_n_padded(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, mp_size_t zlo, mp_size_t zhi); +#endif + FLINT_DLL extern const signed short flint_mpn_mulhigh_k_tab[FLINT_MPN_MULHIGH_K_TAB_SIZE]; mp_limb_t flint_mpn_mullow_basecase(mp_ptr res, mp_srcptr u, mp_srcptr v, mp_size_t n); @@ -731,10 +764,17 @@ void flint_mpn_mul_or_mullow_n(mp_ptr rp, mp_srcptr xp, mp_srcptr yp, mp_size_t if (FLINT_HAVE_MULLOW_FUNC(n)) rp[n] = flint_mpn_mullow_func_tab[n](rp, xp, yp); +#if FLINT_HAVE_FFT_SMALL + else if (n < FLINT_MPN_MULHIGH_FFT_SMALL_CUTOFF) + rp[n] = _flint_mpn_mullow_n(rp, xp, yp, n); + else + flint_mpn_mulmid_fft_small(rp, xp, n, yp, n, 0, n + 1); +#else else if (n < FLINT_MPN_MULHIGH_MUL_CUTOFF) rp[n] = _flint_mpn_mullow_n(rp, xp, yp, n); else flint_mpn_mul_n(rp, xp, yp, n); +#endif } MPN_EXTRAS_INLINE @@ -744,16 +784,19 @@ void flint_mpn_mul_or_mulhigh_n(mp_ptr rp, mp_srcptr xp, mp_srcptr yp, mp_size_t if (FLINT_HAVE_MULHIGH_FUNC(n)) rp[n - 1] = flint_mpn_mulhigh_func_tab[n](rp + n, xp, yp); +#if FLINT_HAVE_FFT_SMALL + else if (n < FLINT_MPN_MULHIGH_FFT_SMALL_CUTOFF) + rp[n - 1] = _flint_mpn_mulhigh_n(rp + n, xp, yp, n); + else + flint_mpn_mulmid_fft_small(rp + n - 2, xp, n, yp, n, n - 2, 2 * n); +#else else if (n < FLINT_MPN_MULHIGH_MUL_CUTOFF) rp[n - 1] = _flint_mpn_mulhigh_n(rp + n, xp, yp, n); else flint_mpn_mul_n(rp, xp, yp, n); +#endif } -#define FLINT_MPN_SQRHIGH_MULDERS_CUTOFF 90 -#define FLINT_MPN_SQRHIGH_SQR_CUTOFF 2000 -#define FLINT_MPN_SQRHIGH_K_TAB_SIZE 2048 - #if FLINT_HAVE_ASSEMBLY_x86_64_adx mp_limb_t _flint_mpn_sqrhigh_basecase_even(mp_ptr, mp_srcptr, mp_size_t); mp_limb_t _flint_mpn_sqrhigh_basecase_odd(mp_ptr, mp_srcptr, mp_size_t); diff --git a/src/mpn_extras/mulhigh.c b/src/mpn_extras/mulhigh.c index 725eb0ff4d..36ecc81992 100644 --- a/src/mpn_extras/mulhigh.c +++ b/src/mpn_extras/mulhigh.c @@ -204,6 +204,34 @@ _flint_mpn_mulhigh_n_mul(mp_ptr res, mp_srcptr u, mp_srcptr v, mp_size_t n) return bot; } +#if FLINT_HAVE_FFT_SMALL + +mp_limb_t +_flint_mpn_mulhigh_n_fft_small(mp_ptr res, mp_srcptr u, mp_srcptr v, mp_size_t n) +{ + mp_ptr tmp; + mp_limb_t bot; + tmp = flint_malloc(sizeof(mp_limb_t) * (n + 2)); + flint_mpn_mulmid_fft_small(tmp, u, n, v, n, n - 2, 2 * n); + memcpy(res, tmp + 2, sizeof(mp_limb_t) * n); + bot = tmp[1]; + flint_free(tmp); + return bot; +} + +mp_limb_t +_flint_mpn_mulhigh_n(mp_ptr res, mp_srcptr u, mp_srcptr v, mp_size_t n) +{ + if (n <= FLINT_MPN_MULHIGH_MULDERS_CUTOFF) + return _flint_mpn_mulhigh_n_basecase2(res, u, v, n); + else if (n <= FLINT_MPN_MULHIGH_FFT_SMALL_CUTOFF) + return _flint_mpn_mulhigh_n_mulders(res, u, v, n); + else + return _flint_mpn_mulhigh_n_fft_small(res, u, v, n); +} + +#else + mp_limb_t _flint_mpn_mulhigh_n(mp_ptr res, mp_srcptr u, mp_srcptr v, mp_size_t n) { @@ -215,6 +243,8 @@ _flint_mpn_mulhigh_n(mp_ptr res, mp_srcptr u, mp_srcptr v, mp_size_t n) return _flint_mpn_mulhigh_n_mul(res, u, v, n); } +#endif + mp_limb_pair_t _flint_mpn_mulhigh_normalised(mp_ptr rp, mp_srcptr xp, mp_srcptr yp, mp_size_t n) { FLINT_ASSERT(n >= 1); diff --git a/src/mpn_extras/mullow.c b/src/mpn_extras/mullow.c index e52d1dcaa4..d95a68fdc2 100644 --- a/src/mpn_extras/mullow.c +++ b/src/mpn_extras/mullow.c @@ -135,6 +135,34 @@ _flint_mpn_mullow_n_mul(mp_ptr res, mp_srcptr u, mp_srcptr v, mp_size_t n) return cy; } +#if FLINT_HAVE_FFT_SMALL + +mp_limb_t +_flint_mpn_mullow_n_fft_small(mp_ptr res, mp_srcptr u, mp_srcptr v, mp_size_t n) +{ + mp_ptr tmp; + mp_limb_t cy; + tmp = flint_malloc(sizeof(mp_limb_t) * (n + 1)); + flint_mpn_mulmid_fft_small(tmp, u, n, v, n, 0, n + 1); + memcpy(res, tmp, sizeof(mp_limb_t) * n); + cy = tmp[n]; + flint_free(tmp); + return cy; +} + +mp_limb_t +_flint_mpn_mullow_n(mp_ptr res, mp_srcptr u, mp_srcptr v, mp_size_t n) +{ + if (n <= FLINT_MPN_MULLOW_MULDERS_CUTOFF) + return flint_mpn_mullow_basecase(res, u, v, n); + else if (n <= FLINT_MPN_MULHIGH_FFT_SMALL_CUTOFF) + return _flint_mpn_mullow_n_mulders(res, u, v, n); + else + return _flint_mpn_mullow_n_fft_small(res, u, v, n); +} + +#else + mp_limb_t _flint_mpn_mullow_n(mp_ptr res, mp_srcptr u, mp_srcptr v, mp_size_t n) { @@ -145,3 +173,6 @@ _flint_mpn_mullow_n(mp_ptr res, mp_srcptr u, mp_srcptr v, mp_size_t n) else return _flint_mpn_mullow_n_mul(res, u, v, n); } + +#endif + diff --git a/src/mpn_extras/mulmid.c b/src/mpn_extras/mulmid.c new file mode 100644 index 0000000000..3e52eca03b --- /dev/null +++ b/src/mpn_extras/mulmid.c @@ -0,0 +1,146 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + 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 "mpn_extras.h" +#include "ulong_extras.h" + +#if FLINT_HAVE_NATIVE_mpn_mulmid_n && FLINT_HAVE_FFT_SMALL +# define MULMID_KARATSUBA_CUTOFF 36 +#else +# define MULMID_KARATSUBA_CUTOFF 20 +#endif + +#if FLINT_HAVE_FFT_SMALL +# define MULMID_FFT_SMALL_CUTOFF 250 +#endif + +#if !FLINT_HAVE_NATIVE_mpn_mulmid_n && !FLINT_HAVE_FFT_SMALL +/* bare build: above this a full product (GMP/FFT-backed flint_mpn_mul) plus a + slice beats schoolbook for a wide middle window */ +# define MULMID_BARE_MUL_CUTOFF 128 +#endif + +/* Zeros we are willing to append to a length-len operand, and the test for it + (also true when target <= len, i.e. when the operand is merely truncated). */ +#define MULMID_PAD_LIMIT(len) FLINT_MAX((mp_size_t) 4, (len) / 32) +#define MULMID_PAD_OK(target, len) ((target) - (len) <= MULMID_PAD_LIMIT(len)) + +void +flint_mpn_mulmid(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, + mp_size_t zlo, mp_size_t zhi) +{ + mp_size_t ac, bc, m, M, zn, top, lim; +#if FLINT_HAVE_NATIVE_mpn_mulmid_n + mp_size_t vac, vbc, vpa, vqb, La, Lb, zlo2, nb, vhalf; + int via_ok; +#endif + + FLINT_ASSERT(an >= 1); + FLINT_ASSERT(bn >= 1); + FLINT_ASSERT(zlo >= 0); + FLINT_ASSERT(zhi > zlo); + FLINT_ASSERT(zhi <= an + bn); + + ac = FLINT_MIN(an, zhi); + bc = FLINT_MIN(bn, zhi); + m = FLINT_MIN(ac, bc); + M = FLINT_MAX(ac, bc); + zn = zhi - zlo; + top = an + bn - zhi; /* limbs above the window */ + lim = MULMID_PAD_LIMIT(m); + + if (m < MULMID_KARATSUBA_CUTOFF || zn * zn <= m) + { + flint_mpn_mulmid_classical(z, a, an, b, bn, zlo, zhi); + return; + } + + /* (A) full / nearly-full window */ + if (zlo <= lim && top <= lim) + { + flint_mpn_mulmid_via_mul(z, a, an, b, bn, zlo, zhi); + return; + } + + /* (B) low / nearly-low window */ + if (zlo <= lim && MULMID_PAD_OK(zhi, an) && MULMID_PAD_OK(zhi, bn) +#if FLINT_HAVE_FFT_SMALL + && (m < 2 * MULMID_FFT_SMALL_CUTOFF || zn < MULMID_FFT_SMALL_CUTOFF) +#endif + ) + { + flint_mpn_mulmid_via_mullow_n(z, a, an, b, bn, zlo, zhi); + return; + } + + /* (C) high / nearly-high window */ + if (top <= lim && zlo + lim >= M) + { + flint_mpn_mulmid_via_mulhigh_n(z, a, an, b, bn, zlo, zhi); + return; + } + +#if FLINT_HAVE_FFT_SMALL + /* large operands with a large (non-full, non-low, non-high) window: FFT */ + if (m >= MULMID_FFT_SMALL_CUTOFF && zn >= MULMID_FFT_SMALL_CUTOFF) + { + flint_mpn_mulmid_fft_small(z, a, an, b, bn, zlo, zhi); + return; + } +#endif + +#if FLINT_HAVE_NATIVE_mpn_mulmid_n + /* Padding flint_mpn_mulmid_via_n_padded would incur: after clamping to zhi + and slicing the low non-contributing limbs, the longer sliced operand La + fills a 2n-1 slot and the shorter Lb an n slot. */ + vac = FLINT_MIN(an, zhi); + vbc = FLINT_MIN(bn, zhi); + vpa = (zlo > vbc - 1) ? zlo - (vbc - 1) : 0; + vqb = (zlo > vac - 1) ? zlo - (vac - 1) : 0; + La = vac - vpa; + Lb = vbc - vqb; + via_ok = 0; + if (La > 0 && Lb > 0) + { + zlo2 = zlo - vpa - vqb; + if (La < Lb) { mp_size_t t2 = La; La = Lb; Lb = t2; } + nb = zn; + if (Lb > nb) nb = Lb; + vhalf = (La + Lb - zlo2 + 1) / 2; + if (vhalf > nb) nb = vhalf; + /* via_n_padded reduces to a balanced middle product of size nb; even + with padding this is worth it whenever nb stays comfortably below the + M-sized problem the other reductions would build, so a generous + padding tolerance (max(8, len/8)) is used here. */ + via_ok = (2 * nb - 1) - La <= FLINT_MAX((mp_size_t) 8, La / 8) + && nb - Lb <= FLINT_MAX((mp_size_t) 8, Lb / 8); + } + + /* (D) balanced Karatsuba/Toom reduction, when its padding is economical */ + if (via_ok) + { + flint_mpn_mulmid_via_n_padded(z, a, an, b, bn, zlo, zhi); + return; + } +#endif + +#if FLINT_HAVE_NATIVE_mpn_mulmid_n || FLINT_HAVE_FFT_SMALL + /* medium middle window: schoolbook (large wide windows already handled) */ + flint_mpn_mulmid_classical(z, a, an, b, bn, zlo, zhi); +#else + /* bare build: a wide middle window is cheaper as a full product plus a slice + (flint_mpn_mul is GMP/FFT-backed); a narrow one stays schoolbook */ + if (m >= MULMID_BARE_MUL_CUTOFF && zn >= MULMID_BARE_MUL_CUTOFF) + flint_mpn_mulmid_via_mul(z, a, an, b, bn, zlo, zhi); + else + flint_mpn_mulmid_classical(z, a, an, b, bn, zlo, zhi); +#endif +} diff --git a/src/mpn_extras/mulmid_classical.c b/src/mpn_extras/mulmid_classical.c new file mode 100644 index 0000000000..5e0f37401d --- /dev/null +++ b/src/mpn_extras/mulmid_classical.c @@ -0,0 +1,174 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + 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 "mpn_extras.h" + +/* Propagate a single-limb carry cc into z[k..), extending the valid prefix + *top by writing fresh limbs as needed, and dropping any carry that would + land at index >= zn (i.e. at a limb position >= the top of the window). */ +static void +_flint_mpn_mulmid_carry(mp_ptr z, mp_size_t * top, mp_size_t zn, + mp_size_t k, mp_limb_t cc) +{ + while (cc != 0 && k < zn) + { + if (k < *top) + { + mp_limb_t t = z[k] + cc; + cc = (t < z[k]); + z[k] = t; + k++; + } + else + { + z[k] = cc; + *top = k + 1; + cc = 0; + } + } +} + +/* + Set (z, zhi - zlo) to the limb window [zlo, zhi) of the integer product + a * b, as a lower approximation: partial products a[p] * b[q] landing + entirely below limb zlo are dropped, so the carry they would propagate + into limb zlo is not recovered. This matches the contract of + radix_mulmid_classical with radix B = 2^64. + + The low-end deficit (exact window minus computed window) is bounded by + min(an, bn, zlo) * 2^64; limbs sufficiently above zlo are exact, up to the + truncation at zhi. With zlo == 0 the result is exact. + + Requires an >= 1, bn >= 1 and 0 <= zlo < zhi <= an + bn. + + Schoolbook: each row a * b[i] is accumulated with mpn_mul_1 / mpn_addmul_1 + over only the limbs that fall inside the window. No scratch space is used + and z is never pre-zeroed: the first write to each output limb is a + mpn_mul_1 or a carry store, and later rows accumulate on top. +*/ +void +flint_mpn_mulmid_classical(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, + mp_size_t zlo, mp_size_t zhi) +{ + mp_size_t i, o, s, e, len, zn, i0, top; + mp_limb_t cc; + + FLINT_ASSERT(an >= 1); + FLINT_ASSERT(bn >= 1); + FLINT_ASSERT(zlo >= 0); + FLINT_ASSERT(zhi > zlo); + FLINT_ASSERT(zhi <= an + bn); + + if (an < bn) + { + FLINT_SWAP(mp_srcptr, a, b); + FLINT_SWAP(mp_size_t, an, bn); + } + + /* limbs of either input at or above zhi cannot reach the window */ + an = FLINT_MIN(an, zhi); + bn = FLINT_MIN(bn, zhi); + + zn = zhi - zlo; + + /* + Phase 1: the rows i in [0, min(zlo, bn-1)] all land at output offset 0 + (their low parts a[0 .. zlo-i) are below zlo and dropped). These rows + have non-decreasing length in i, so the largest i is the longest; do it + first with mpn_mul_1 to seed the prefix, then add the shorter ones on + top (whose carries propagate within the interior). + */ + i0 = FLINT_MIN(zlo, bn - 1); + + s = zlo - i0; + e = FLINT_MIN(an, zhi - i0); + len = e - s; + if (len > 0) + { + cc = mpn_mul_1(z, a + s, len, b[i0]); + top = len; + if (len < zn) + { + z[len] = cc; + top = len + 1; + } + } + else + { + /* the longest phase-1 row is already empty, so every product lands + below zlo: the window starts out all zero */ + z[0] = 0; + top = 1; + } + + for (i = i0 - 1; i >= 0; i--) + { + s = zlo - i; + e = FLINT_MIN(an, zhi - i); + len = e - s; + if (len <= 0) + break; /* this and all smaller i are below window */ + cc = mpn_addmul_1(z, a + s, len, b[i]); + _flint_mpn_mulmid_carry(z, &top, zn, len, cc); + } + + /* + Phase 2: rows i in (zlo, bn) land at offset i - zlo, sliding right by one + limb per row. Each row's region is interior up to the current frontier + and extends it by at most one limb (written via the carry), exactly as + in a plain mpn multiplication. + */ + for (i = zlo + 1; i < bn; i++) + { + o = i - zlo; + len = FLINT_MIN(an, zhi - i); /* s == 0 */ + if (len <= 0) + break; /* zhi <= i: no further rows reach */ + + if (o + len <= top) + { + /* the whole row sits inside the current prefix */ + cc = mpn_addmul_1(z + o, a, len, b[i]); + _flint_mpn_mulmid_carry(z, &top, zn, o + len, cc); + } + else + { + /* [o, top) overlaps the prefix; [top, o+len) is fresh */ + mp_size_t ov = top - o; /* overlap length, >= 0 */ + mp_size_t right = o + len; + mp_limb_t ccf = mpn_mul_1(z + top, a + ov, len - ov, b[i]); + mp_size_t topf = right; + + if (right < zn) + { + z[right] = ccf; + topf = right + 1; + } + + if (ov > 0) + { + cc = mpn_addmul_1(z + o, a, ov, b[i]); + top = topf; + _flint_mpn_mulmid_carry(z, &top, zn, o + ov, cc); + } + else + { + top = topf; + } + } + } + + /* Any high limbs never reached by a carry are exact zeros: with zhi == + an + bn the very top limb can be one of these, so fill explicitly + rather than asserting the frontier reached the top. */ + if (top < zn) + flint_mpn_zero(z + top, zn - top); +} diff --git a/src/mpn_extras/mulmid_fft_small.c b/src/mpn_extras/mulmid_fft_small.c new file mode 100644 index 0000000000..820e6e48f9 --- /dev/null +++ b/src/mpn_extras/mulmid_fft_small.c @@ -0,0 +1,34 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + 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 "mpn_extras.h" + +#if FLINT_HAVE_FFT_SMALL + +#include "fft_small.h" + +/* + Windowed middle product through the small-prime FFT. This is a thin wrapper + over _mpn_ctx_mpn_mul_range, which returns the limb window [zlo, zhi) of a*b + as the same lower approximation as the rest of the family (partial products + landing entirely below zlo are dropped). The default (thread-local) context + is used, exactly as mpn_mul_default_mpn_ctx does for the full product. +*/ +void +flint_mpn_mulmid_fft_small(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, + mp_size_t zlo, mp_size_t zhi) +{ + _mpn_ctx_mpn_mul_range(get_default_mpn_ctx(), z, + (ulong) zlo, (ulong) zhi, + a, (ulong) an, b, (ulong) bn); +} + +#endif diff --git a/src/mpn_extras/mulmid_via_mul.c b/src/mpn_extras/mulmid_via_mul.c new file mode 100644 index 0000000000..263cb6244e --- /dev/null +++ b/src/mpn_extras/mulmid_via_mul.c @@ -0,0 +1,50 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + 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 "mpn_extras.h" + +/* + Windowed middle product via the full product. Computes the whole product + a * b with the standard (unbalanced) flint_mpn_mul and copies out the limb + window [zlo, zhi). Correct for any an >= 1, bn >= 1, 0 <= zlo < zhi <= an+bn; + the result is the *exact* window (every partial product and carry included), + which is a valid instance of the lower-approximation contract with zero + deficit. Cheapest when the window is a large fraction of the product; for a + small window it wastes work on the discarded limbs. +*/ +void +flint_mpn_mulmid_via_mul(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, + mp_size_t zlo, mp_size_t zhi) +{ + mp_size_t zn = zhi - zlo; + mp_srcptr xp, yp; + mp_size_t xn, yn; + mp_ptr t; + TMP_INIT; + + FLINT_ASSERT(an >= 1 && bn >= 1); + FLINT_ASSERT(0 <= zlo && zlo < zhi && zhi <= an + bn); + + if (an >= bn) { xp = a; xn = an; yp = b; yn = bn; } + else { xp = b; xn = bn; yp = a; yn = an; } + + if (zlo == 0 && zhi == an + bn) /* window is the entire product */ + { + flint_mpn_mul(z, xp, xn, yp, yn); + return; + } + + TMP_START; + t = TMP_ARRAY_ALLOC(an + bn, mp_limb_t); + flint_mpn_mul(t, xp, xn, yp, yn); + flint_mpn_copyi(z, t + zlo, zn); + TMP_END; +} diff --git a/src/mpn_extras/mulmid_via_mulhigh_n.c b/src/mpn_extras/mulmid_via_mulhigh_n.c new file mode 100644 index 0000000000..fa43ec5ff8 --- /dev/null +++ b/src/mpn_extras/mulmid_via_mulhigh_n.c @@ -0,0 +1,140 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + 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 "mpn_extras.h" + +/* + Windowed middle product via a balanced high product. flint_mpn_mulhigh_n + returns (a lower approximation of) the top n limbs, i.e. limbs [n, 2n), of a + balanced n x n product. + + Only the high ends of a and b reach the window: a limb a[p] can contribute to + limb >= zlo only if p + (bc - 1) >= zlo, and a limb at or above zhi never + reaches the window at all. So we first clamp to zhi and slice off the low, + non-contributing limbs + + ac = min(an, zhi), bc = min(bn, zhi), + pa = max(0, zlo - (bc - 1)), qb = max(0, zlo - (ac - 1)), + a' = a + pa (La = ac - pa limbs), b' = b + qb (Lb = bc - qb limbs), + zlo2 = zlo - pa - qb, + + which leaves a balanced-ish problem of size ~ zn rather than ~ max(an, bn). + The sliced product a'*b' still contains a few sub-zlo diagonals, but they lie + below limb zlo2 and are dropped by the high product just like the ones sliced + away -- so the net effect is the required "drop p + q < zlo". + + The window [zlo2, zlo2 + zn) of a'*b' is then exposed as a top slice: pad both + operands to n limbs and, when the window starts below M = max(La, Lb), shift + them up by L = M - zlo2 low zero limbs so that limb zlo2 lands at limb n: + + L = max(0, M - zlo2), n = M + L, off = zlo2 - M + L (>= 0), + res = mulhigh_n(0^L,a',0.. , 0^L,b',0..) ~= limbs [n,2n) of a'*b'*B^{2L}, + z = res[off .. off + zn). + + Correct for any input; the result is a lower approximation of the exact + window (mulhigh never exceeds the true high part). Cheapest when the window + sits at the top of the product (zlo2 >= M, L = 0). + + Copies are avoided where possible: when L = 0 an operand of length n is passed + verbatim and only a shorter one is padded; padded operands zero only their + genuine zero limbs; and when the window is exactly mulhigh_n's whole output + (off = 0, zn = n) it is written straight into z. +*/ +void +flint_mpn_mulmid_via_mulhigh_n(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, + mp_size_t zlo, mp_size_t zhi) +{ + mp_size_t zn = zhi - zlo; + mp_size_t ac, bc, pa, qb, La, Lb, zlo2, M, L, n, off; + mp_srcptr xp, yp; + mp_ptr res; + TMP_INIT; + + FLINT_ASSERT(an >= 1 && bn >= 1); + FLINT_ASSERT(0 <= zlo && zlo < zhi && zhi <= an + bn); + + ac = FLINT_MIN(an, zhi); + bc = FLINT_MIN(bn, zhi); + pa = (zlo > bc - 1) ? zlo - (bc - 1) : 0; + qb = (zlo > ac - 1) ? zlo - (ac - 1) : 0; + La = ac - pa; + Lb = bc - qb; + + if (La <= 0 || Lb <= 0) + { + /* every partial product lands below zlo: the lower-approximation window + is all zero */ + flint_mpn_zero(z, zn); + return; + } + + a += pa; + b += qb; + zlo2 = zlo - pa - qb; + + M = FLINT_MAX(La, Lb); + L = (zlo2 < M) ? (M - zlo2) : 0; + n = M + L; + off = (zlo2 - M) + L; + xp = a; + yp = b; + + FLINT_ASSERT(off >= 0 && off + zn <= n); + + TMP_START; + + if (L == 0) + { + if (La < n) + { + mp_ptr X = TMP_ARRAY_ALLOC(n, mp_limb_t); + flint_mpn_copyi(X, a, La); + flint_mpn_zero(X + La, n - La); + xp = X; + } + if (Lb < n) + { + mp_ptr Y = TMP_ARRAY_ALLOC(n, mp_limb_t); + flint_mpn_copyi(Y, b, Lb); + flint_mpn_zero(Y + Lb, n - Lb); + yp = Y; + } + } + else + { + mp_ptr X = TMP_ARRAY_ALLOC(n, mp_limb_t); + mp_ptr Y = TMP_ARRAY_ALLOC(n, mp_limb_t); + + flint_mpn_zero(X, L); + flint_mpn_copyi(X + L, a, La); + flint_mpn_zero(X + L + La, n - L - La); + + flint_mpn_zero(Y, L); + flint_mpn_copyi(Y + L, b, Lb); + flint_mpn_zero(Y + L + Lb, n - L - Lb); + + xp = X; + yp = Y; + } + + if (off == 0 && zn == n) + { + flint_mpn_mulhigh_n(z, xp, yp, n); + } + else + { + res = TMP_ARRAY_ALLOC(n, mp_limb_t); + flint_mpn_mulhigh_n(res, xp, yp, n); + flint_mpn_copyi(z, res + off, zn); + } + + TMP_END; +} diff --git a/src/mpn_extras/mulmid_via_mullow_n.c b/src/mpn_extras/mulmid_via_mullow_n.c new file mode 100644 index 0000000000..3d1d08b3d3 --- /dev/null +++ b/src/mpn_extras/mulmid_via_mullow_n.c @@ -0,0 +1,65 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + 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 "mpn_extras.h" + +/* + Windowed middle product via a balanced low product. The low zhi limbs of + a * b are (a mod B^zhi) * (b mod B^zhi) mod B^zhi, computed exactly by + flint_mpn_mullow_n on operands zero-padded / truncated to n = zhi limbs; the + window [zlo, zhi) is then the top zn of those. Correct for any input; the + result is the *exact* window (with carry-in from below zlo), a zero-deficit + instance of the contract. Cheapest when zhi is small (a low or nearly-low + window); for large zhi it computes an almost full n x n product. +*/ +void +flint_mpn_mulmid_via_mullow_n(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, + mp_size_t zlo, mp_size_t zhi) +{ + mp_size_t n = zhi; + mp_size_t zn = zhi - zlo; + mp_srcptr xp = a, yp = b; + mp_ptr rp; + TMP_INIT; + + FLINT_ASSERT(an >= 1 && bn >= 1); + FLINT_ASSERT(0 <= zlo && zlo < zhi && zhi <= an + bn); + + TMP_START; + + /* present both operands as n limbs: truncate if longer, zero-extend if shorter */ + if (an < n) + { + mp_ptr X = TMP_ARRAY_ALLOC(n, mp_limb_t); + flint_mpn_copyi(X, a, an); + flint_mpn_zero(X + an, n - an); + xp = X; + } + if (bn < n) + { + mp_ptr Y = TMP_ARRAY_ALLOC(n, mp_limb_t); + flint_mpn_copyi(Y, b, bn); + flint_mpn_zero(Y + bn, n - bn); + yp = Y; + } + + if (zlo == 0) /* window is exactly the low n limbs */ + { + flint_mpn_mullow_n(z, xp, yp, n); + TMP_END; + return; + } + + rp = TMP_ARRAY_ALLOC(n, mp_limb_t); + flint_mpn_mullow_n(rp, xp, yp, n); + flint_mpn_copyi(z, rp + zlo, zn); + TMP_END; +} diff --git a/src/mpn_extras/mulmid_via_n_padded.c b/src/mpn_extras/mulmid_via_n_padded.c new file mode 100644 index 0000000000..336619c961 --- /dev/null +++ b/src/mpn_extras/mulmid_via_n_padded.c @@ -0,0 +1,147 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + 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 "mpn_extras.h" + +#if FLINT_HAVE_NATIVE_mpn_mulmid_n + +/* + Set (z, zhi - zlo) to the limb window [zlo, zhi) of the integer product + a * b, as a lower approximation: partial products a[p] * b[q] whose low + limb lies below limb zlo (i.e. p + q < zlo) are dropped, so the carry they + would propagate into limb zlo is not recovered. This matches the contract + of radix_mulmid_classical with radix B = 2^64 (and the schoolbook + flint_mpn_mulmid_classical). Equivalently + + z = ( sum_{p+q >= zlo} a[p]*b[q] * B^(p+q-zlo) ) mod B^(zhi-zlo). + + Requires an >= 1, bn >= 1 and 0 <= zlo < zhi <= an + bn. + + Implementation: reduce the window to a single *balanced* middle product and + hand it to flint_mpn_mulmid_n (Karatsuba/Toom-42). flint_mpn_mulmid_n + computes, for inputs {A, 2n-1} and {B, n}, the band n-1 <= p'+q' < 2n-1 of + A*B shifted down by n-1, dropping all products with p'+q' < n-1 -- which is + exactly the same "drop everything below the floor" rule as our window's + lower-approximation, so the deficits coincide and the result is the *same* + value the schoolbook would produce (it is not merely another approximation). + + To map the window onto that shape we: + + (1) clamp an, bn to zhi (limbs at/above zhi cannot reach the window); + + (2) slice off the low limbs that can never contribute: a[p] with + p + (bn-1) < zlo, and b[q] with q + (an-1) < zlo. Every product + removed here has p+q < zlo (already dropped) or low limb >= zhi + (outside the window), so slicing changes nothing. This also keeps + the balanced size from being tied to the full an; + + (3) place the (longer) sliced operand of length La in the 2n-1 slot at + offset alpha and the shorter one of length Lb in the n slot at the + top (offset beta = n - Lb), with alpha + beta = n - 1 - zlo2 so that + output limb k lands exactly at product position zlo + k. Here zlo2 + is the window floor in sliced coordinates; one shows zlo2 <= Lb - 1, + hence alpha = (Lb-1) - zlo2 >= 0. + + We choose the smallest n that fits both operands and covers the window: + + n = max( zn, Lb, ceil((La + Lb - zlo2) / 2) ). + + NOTE ON THE ALTERNATIVE. This is the "generous n + zero-pad" route: one + balanced call, no fix-up. It is a good fit when the window is wide (zn + comparable to Lb), i.e. exactly the regime where Karatsuba/Toom pays off. + It is *not* a good fit for a narrow window sitting in a tall, skewed region: + there n is forced up to ~Lb (B must hold all of b), so the balanced call + costs ~O(Lb^2) while the window only depends on ~zn diagonals and plain + schoolbook (flint_mpn_mulmid_classical) costs ~O(zn * Lb). A + production dispatcher should therefore fall back to schoolbook for narrow + windows, or chunk the region the way the GMP-derived general middle product + (flint_mpn_mulmid_unbalanced, mulmid_gmp.c) does with MULMID_CHUNK, and only enter this + balanced reduction once a chunk is square enough to benefit. +*/ +void +flint_mpn_mulmid_via_n_padded(mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, mp_size_t bn, + mp_size_t zlo, mp_size_t zhi) +{ + mp_size_t zn = zhi - zlo; + mp_size_t pa, qb, La, Lb, zlo2, n, half, alpha, beta; + mp_srcptr ap, bp; + mp_ptr A, B, rp; + TMP_INIT; + + FLINT_ASSERT(an >= 1); + FLINT_ASSERT(bn >= 1); + FLINT_ASSERT(zlo >= 0); + FLINT_ASSERT(zhi > zlo); + FLINT_ASSERT(zhi <= an + bn); + + /* (1) limbs of either input at or above zhi cannot reach the window */ + an = FLINT_MIN(an, zhi); + bn = FLINT_MIN(bn, zhi); + + /* (2) drop low limbs that cannot contribute to any in-window diagonal */ + pa = (zlo > bn - 1) ? zlo - (bn - 1) : 0; + qb = (zlo > an - 1) ? zlo - (an - 1) : 0; + La = an - pa; + Lb = bn - qb; + + /* nothing survives: every contributing product was dropped or out of range */ + if (La <= 0 || Lb <= 0) + { + flint_mpn_zero(z, zn); + return; + } + + ap = a + pa; + bp = b + qb; + zlo2 = zlo - pa - qb; /* window floor in sliced coordinates, >= 0 */ + + /* longer sliced operand goes in the 2n-1 slot */ + if (La < Lb) + { + FLINT_SWAP(mp_srcptr, ap, bp); + FLINT_SWAP(mp_size_t, La, Lb); + } + + /* (3) smallest balanced size covering the window and holding both operands */ + n = zn; + if (Lb > n) + n = Lb; + half = (La + Lb - zlo2 + 1) / 2; /* ceil((La + Lb - zlo2) / 2) */ + if (half > n) + n = half; + + alpha = (Lb - 1) - zlo2; /* offset of the long operand in A */ + beta = n - Lb; /* short operand sits atop B */ + + FLINT_ASSERT(alpha >= 0 && beta >= 0); + FLINT_ASSERT(zn <= n); + FLINT_ASSERT(alpha + La <= 2 * n - 1); + FLINT_ASSERT(beta + Lb <= n); + + TMP_START; + A = TMP_ARRAY_ALLOC(2 * n - 1, mp_limb_t); + B = TMP_ARRAY_ALLOC(n, mp_limb_t); + rp = TMP_ARRAY_ALLOC(n + 2, mp_limb_t); + + flint_mpn_zero(A, 2 * n - 1); + flint_mpn_copyi(A + alpha, ap, La); + flint_mpn_zero(B, n); + flint_mpn_copyi(B + beta, bp, Lb); + + flint_mpn_mulmid_n(rp, A, B, n); + + /* output limb k is product position zlo + k; copy the window verbatim */ + flint_mpn_copyi(z, rp, zn); + + TMP_END; +} + +#endif diff --git a/src/mpn_extras/profile/p-mulmid.c b/src/mpn_extras/profile/p-mulmid.c new file mode 100644 index 0000000000..98eb5f3a41 --- /dev/null +++ b/src/mpn_extras/profile/p-mulmid.c @@ -0,0 +1,274 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + 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 "flint.h" +#include "mpn_extras.h" +#include "ulong_extras.h" +#include "profiler.h" + +typedef void (*mulmid_fn)(mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t, + mp_size_t, mp_size_t); + +#define IDX_DISPATCH 0 +#define NALG 7 + +static const char * const alg_name[NALG] = +{ + "dispatch", "classical", "via_mul", "via_mullow", + "via_mulhigh", "via_n_pad", "fft_small" +}; + +static const mulmid_fn alg_fn[NALG] = +{ + flint_mpn_mulmid, + flint_mpn_mulmid_classical, + flint_mpn_mulmid_via_mul, + flint_mpn_mulmid_via_mullow_n, + flint_mpn_mulmid_via_mulhigh_n, + flint_mpn_mulmid_via_n_padded, + flint_mpn_mulmid_fft_small +}; + +/* seconds per call, adaptively repeated until the cpu clock is meaningful */ +static double +time_call(mulmid_fn f, mp_ptr z, mp_srcptr a, mp_size_t an, mp_srcptr b, + mp_size_t bn, mp_size_t zlo, mp_size_t zhi) +{ + timeit_t timer; + slong reps; + + f(z, a, an, b, bn, zlo, zhi); /* warm up / prime caches */ + + TIMEIT_REPEAT(timer, reps) + f(z, a, an, b, bn, zlo, zhi); + TIMEIT_END_REPEAT(timer, reps); + + return (double) timer->wall * 0.001 / reps; +} + +/* shape categories, chosen to bracket the algorithm crossovers */ +enum { + SH_LOW, SH_NEAR_LOW, + SH_HIGH, SH_NEAR_HIGH, SH_BAL_MID, SH_NEAR_BAL_MID, + SH_RANDOM, SH_NEAR_FULL, + SH_NCAT +}; + +static const char * const cat_name[SH_NCAT] = +{ + "low", "near-low", + "high", "near-high", "bal-mid", "near-bal-mid", + "near-full", + "random", +}; + +/* small signed jitter in [-j, j] */ +static mp_size_t +jitter(flint_rand_t state, mp_size_t j) +{ + return (mp_size_t) (n_randint(state, 2 * j + 1)) - j; +} + +static mp_size_t +clampsz(mp_size_t x, mp_size_t lo, mp_size_t hi) +{ + if (x < lo) x = lo; + if (x > hi) x = hi; + return x; +} + +/* fill (an, bn, zlo, zhi) for a category at a given scale */ +static void +gen_shape(flint_rand_t state, int cat, mp_size_t scale, + mp_size_t * an, mp_size_t * bn, mp_size_t * zlo, mp_size_t * zhi) +{ + mp_size_t a, b, n, M, tot, lo, hi; + + switch (cat) + { + case SH_RANDOM: + a = 1 + n_randint(state, scale); + b = 1 + n_randint(state, scale); + tot = a + b; + lo = n_randint(state, tot); + hi = lo + 1 + n_randint(state, tot - lo); + break; + + case SH_NEAR_FULL: /* whole product but a few limbs */ + a = scale / 2 + 1 + n_randint(state, scale / 2 + 1); + b = scale / 2 + 1 + n_randint(state, scale / 2 + 1); + tot = a + b; + lo = n_randint(state, 5); + hi = tot - n_randint(state, 5); + break; + + case SH_LOW: /* exact low product, mullow shape */ + a = scale / 2 + 1 + n_randint(state, scale / 2 + 1); + b = scale / 2 + 1 + n_randint(state, scale / 2 + 1); + tot = a + b; + lo = 0; + hi = clampsz(FLINT_MIN(a, b) + jitter(state, 4), 1, tot); + break; + + case SH_NEAR_LOW: /* low but with a small nonzero zlo */ + a = scale / 2 + 1 + n_randint(state, scale / 2 + 1); + b = scale / 2 + 1 + n_randint(state, scale / 2 + 1); + tot = a + b; + lo = 1 + n_randint(state, 4); + hi = clampsz(FLINT_MIN(a, b) + jitter(state, 4), lo + 1, tot); + break; + + case SH_HIGH: /* top slice, mulhigh shape */ + a = scale / 2 + 1 + n_randint(state, scale / 2 + 1); + b = scale / 2 + 1 + n_randint(state, scale / 2 + 1); + tot = a + b; + M = FLINT_MAX(a, b); + lo = M + n_randint(state, FLINT_MIN(a, b)); + hi = tot; + break; + + case SH_NEAR_HIGH: /* near the mulhigh boundary */ + a = scale / 2 + 1 + n_randint(state, scale / 2 + 1); + b = scale / 2 + 1 + n_randint(state, scale / 2 + 1); + tot = a + b; + M = FLINT_MAX(a, b); + lo = clampsz(M + jitter(state, 4), 0, tot - 2); + hi = tot - n_randint(state, 5); + break; + + case SH_BAL_MID: /* balanced middle: mulmid_n shape */ + n = FLINT_MAX((mp_size_t) 2, scale); + a = 2 * n - 1; + b = n; + lo = n - 1; + hi = 2 * n - 1; + break; + + case SH_NEAR_BAL_MID: /* perturbed balanced middle */ + default: + n = FLINT_MAX((mp_size_t) 8, scale); + a = 2 * n - 1 + jitter(state, 5); + b = n + jitter(state, 5); + tot = a + b; + lo = clampsz(n - 1 + jitter(state, 5), 0, tot - 2); + hi = clampsz(2 * n - 1 + jitter(state, 5), lo + 1, tot); + break; + } + + if (hi > a + b) hi = a + b; + if (lo >= hi) lo = hi - 1; + if (lo < 0) lo = 0; + + *an = a; *bn = b; *zlo = lo; *zhi = hi; +} + +int main(int argc, char ** argv) +{ + flint_rand_t state; + mp_size_t scales[] = { 4, 8, 16, 32, 64, 160, 400, 1000, 2500, 10000, }; + slong nscales = sizeof(scales) / sizeof(scales[0]); + slong reps_per_cell = (argc > 1) ? atol(argv[1]) : 2; + mp_size_t maxn = 2 * scales[nscales - 1] + 16; + mp_ptr a, b, z, zref, d; + double worst_ratio = 1.0; + char worst_desc[128]; + int cat; + slong s, r; + + flint_rand_init(state); + + a = flint_malloc(sizeof(mp_limb_t) * maxn); + b = flint_malloc(sizeof(mp_limb_t) * maxn); + z = flint_malloc(sizeof(mp_limb_t) * 2 * maxn); + zref = flint_malloc(sizeof(mp_limb_t) * 2 * maxn); + d = flint_malloc(sizeof(mp_limb_t) * 2 * maxn); + strcpy(worst_desc, "(none)"); + + for (cat = 0; cat < SH_NCAT; cat++) + { + flint_printf("%-14s %6s %6s %8s %8s |", "shape", "an", "bn", "zlo", "zn"); + for (int j = 0; j < NALG; j++) + flint_printf(" %10s", alg_name[j]); + flint_printf(" | %8s %s\n", "disp/best", "best"); + + for (s = 0; s < nscales; s++) + for (r = 0; r < reps_per_cell; r++) + { + mp_size_t an, bn, zlo, zhi, zn; + double t[NALG]; + double best = 0.0; + int best_j = 0, j; + + do { + gen_shape(state, cat, scales[s], &an, &bn, &zlo, &zhi); + } while (zhi == 0); + + zn = zhi - zlo; + + flint_mpn_rrandom(a, state, an); + flint_mpn_rrandom(b, state, bn); + + /* exact reference (high limbs) to catch gross errors while tuning */ + flint_mpn_mulmid_via_mul(zref, a, an, b, bn, zlo, zhi); + + for (j = 0; j < NALG; j++) + { + t[j] = time_call(alg_fn[j], z, a, an, b, bn, zlo, zhi); + + /* deficit = exact - computed must be non-negative and confined to + the low two window limbs (borrow through zero runs may still make + individual higher limbs differ, so compare the deficit, not the + limbs) */ + { + int bad = 0; + mp_size_t i; + mpn_sub_n(d, zref, z, zn); /* deficit mod 2^(64*zn) */ + for (i = 2; i < zn && !bad; i++) + if (d[i] != 0) + bad = 1; + if (zn >= 2 && d[1] > (mp_limb_t) (an + bn)) + bad = 1; + if (bad) + flint_printf("MISMATCH %s an=%wd bn=%wd zlo=%wd zhi=%wd\n", + alg_name[j], an, bn, zlo, zhi); + } + + if (j == 0 || t[j] < best) { best = t[j]; best_j = j; } + } + + { + double ratio = (best > 0.0) ? t[IDX_DISPATCH] / best : 1.0; + + flint_printf("%-14s %6wd %6wd %8wd %8wd |", + cat_name[cat], an, bn, zlo, zn); + for (j = 0; j < NALG; j++) + flint_printf(" %10.2e", t[j]); + flint_printf(" | %8.2f %s\n", ratio, alg_name[best_j]); + + if (ratio > worst_ratio && best_j != IDX_DISPATCH) + { + worst_ratio = ratio; + flint_sprintf(worst_desc, "%s an=%wd bn=%wd zlo=%wd zn=%wd (best %s)", + cat_name[cat], an, bn, zlo, zn, alg_name[best_j]); + } + } + } + } + + flint_printf("\nworst dispatch slowdown: %.2fx at %s\n", worst_ratio, worst_desc); + + flint_free(a); flint_free(b); flint_free(z); flint_free(zref); flint_free(d); + flint_rand_clear(state); + flint_cleanup_master(); + return 0; +} diff --git a/src/mpn_extras/sqrhigh.c b/src/mpn_extras/sqrhigh.c index a6c0d4eec9..43b6ffa7d2 100644 --- a/src/mpn_extras/sqrhigh.c +++ b/src/mpn_extras/sqrhigh.c @@ -176,6 +176,21 @@ _flint_mpn_sqrhigh_sqr(mp_ptr res, mp_srcptr u, mp_size_t n) return bot; } +#if FLINT_HAVE_FFT_SMALL + +mp_limb_t +_flint_mpn_sqrhigh(mp_ptr res, mp_srcptr u, mp_size_t n) +{ + if (n <= FLINT_MPN_SQRHIGH_MULDERS_CUTOFF) + return _flint_mpn_sqrhigh_basecase(res, u, n); + else if (n <= FLINT_MPN_SQRHIGH_FFT_SMALL_CUTOFF) + return _flint_mpn_sqrhigh_mulders(res, u, n); + else + return _flint_mpn_mulhigh_n_fft_small(res, u, u, n); +} + +#else + mp_limb_t _flint_mpn_sqrhigh(mp_ptr res, mp_srcptr u, mp_size_t n) { @@ -187,6 +202,9 @@ _flint_mpn_sqrhigh(mp_ptr res, mp_srcptr u, mp_size_t n) return _flint_mpn_sqrhigh_sqr(res, u, n); } +#endif + + mp_limb_pair_t _flint_mpn_sqrhigh_normalised(mp_ptr rp, mp_srcptr xp, mp_size_t n) { mp_limb_pair_t ret; diff --git a/src/mpn_extras/test/main.c b/src/mpn_extras/test/main.c index 2c1c54c050..138a47398e 100644 --- a/src/mpn_extras/test/main.c +++ b/src/mpn_extras/test/main.c @@ -29,6 +29,8 @@ #include "t-mulhigh_n_tab.c" #include "t-mulhigh_n_recursive.c" #include "t-mulhigh_normalised.c" +#include "t-mulmid.c" +#include "t-mulmid_n.c" #include "t-mulmod_2expp1.c" #include "t-mulmod_precond_matrix.c" #include "t-mulmod_precond_shoup.c" @@ -62,6 +64,8 @@ test_struct tests[] = TEST_FUNCTION(flint_mpn_mulhigh_n), TEST_FUNCTION(flint_mpn_mulhigh_n_recursive), TEST_FUNCTION(flint_mpn_mulhigh_normalised), + TEST_FUNCTION(flint_mpn_mulmid), + TEST_FUNCTION(flint_mpn_mulmid_n), TEST_FUNCTION(flint_mpn_mulmod_2expp1), TEST_FUNCTION(flint_mpn_mulmod_precond_matrix), TEST_FUNCTION(flint_mpn_mulmod_precond_shoup), diff --git a/src/mpn_extras/test/t-mulmid.c b/src/mpn_extras/test/t-mulmid.c new file mode 100644 index 0000000000..c7a199833e --- /dev/null +++ b/src/mpn_extras/test/t-mulmid.c @@ -0,0 +1,108 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + 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 "test_helpers.h" +#include "mpn_extras.h" + +/* + flint_mpn_mulmid returns the window [zlo, zhi) of a*b as a lower + approximation: computed = exact - D, where the deficit D is a single carry + from the products dropped below zlo, bounded by min(an,bn,zlo)*2^64 (a couple + of backends, e.g. via_mulhigh, use a different but comparably small D). Since + D is subtracted modulo 2^(64*zn), a narrow window can wrap so that the result + exceeds the exact window numerically -- that is allowed. What we can check is + that the deficit d = (exact - computed) mod 2^(64*zn) is confined to the low + two limbs (it fits there whenever the window has room), and that with + zlo == 0, where nothing is dropped, the window is exact. +*/ +TEST_FUNCTION_START(flint_mpn_mulmid, state) +{ + slong ix; + + for (ix = 0; ix < 30000 * flint_test_multiplier(); ix++) + { + mp_size_t an, bn, zn, zlo, zhi, k; + mp_ptr a, b, z, full, d; + + an = 1 + n_randint(state, 30); + bn = 1 + n_randint(state, 30); + if (n_randint(state, 1000) == 0) + { + an = 1 + n_randint(state, 600); + bn = 1 + n_randint(state, 600); + } + + zlo = n_randint(state, an + bn); + zhi = zlo + 1 + n_randint(state, an + bn - zlo); + zn = zhi - zlo; + + a = flint_malloc(sizeof(mp_limb_t) * an); + b = flint_malloc(sizeof(mp_limb_t) * bn); + z = flint_malloc(sizeof(mp_limb_t) * zn); + d = flint_malloc(sizeof(mp_limb_t) * zn); + full = flint_malloc(sizeof(mp_limb_t) * (an + bn)); + + flint_mpn_rrandom(a, state, an); + flint_mpn_rrandom(b, state, bn); + flint_mpn_rrandom(z, state, zn); /* poison */ + + flint_mpn_mulmid(z, a, an, b, bn, zlo, zhi); + + if (an >= bn) + mpn_mul(full, a, an, b, bn); + else + mpn_mul(full, b, bn, a, an); + + /* d = (exact - computed) mod 2^(64*zn) = deficit (mod 2^(64*zn)) */ + mpn_sub_n(d, full + zlo, z, zn); + + if (zn >= 2) + { + for (k = 2; k < zn; k++) + if (d[k] != 0) + TEST_FUNCTION_FAIL( + "deficit exceeds the low two window limbs\n" + "ix = %wd, (an, bn) = (%wd, %wd), (zlo, zhi) = (%wd, %wd)\n" + "a = %{ulong*}\nb = %{ulong*}\n" + "exact = %{ulong*}\ngot = %{ulong*}\n", + ix, an, bn, zlo, zhi, a, an, b, bn, + full + zlo, zn, z, zn); + + if (d[1] > (mp_limb_t) (an + bn)) + TEST_FUNCTION_FAIL( + "deficit too large (second limb %wu > an + bn)\n" + "ix = %wd, (an, bn) = (%wd, %wd), (zlo, zhi) = (%wd, %wd)\n" + "a = %{ulong*}\nb = %{ulong*}\n" + "exact = %{ulong*}\ngot = %{ulong*}\n", + d[1], ix, an, bn, zlo, zhi, a, an, b, bn, + full + zlo, zn, z, zn); + } + + if (zlo == 0) + for (k = 0; k < zn; k++) + if (d[k] != 0) + TEST_FUNCTION_FAIL( + "zlo == 0 window is not exact\n" + "ix = %wd, (an, bn) = (%wd, %wd), zhi = %wd\n" + "a = %{ulong*}\nb = %{ulong*}\n" + "exact = %{ulong*}\ngot = %{ulong*}\n", + ix, an, bn, zhi, a, an, b, bn, + full, zhi, z, zhi); + + flint_free(a); + flint_free(b); + flint_free(z); + flint_free(d); + flint_free(full); + } + + TEST_FUNCTION_END(state); +} diff --git a/src/mpn_extras/test/t-mulmid_n.c b/src/mpn_extras/test/t-mulmid_n.c new file mode 100644 index 0000000000..7de2710c5c --- /dev/null +++ b/src/mpn_extras/test/t-mulmid_n.c @@ -0,0 +1,82 @@ +/* + Copyright (C) 2026 Fredrik Johansson + + 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 "test_helpers.h" +#include "longlong.h" +#include "mpn_extras.h" + +/* + flint_mpn_mulmid_n(rp, a, b, n) is the balanced exact middle product of + {a, 2n-1} and {b, n}, writing n + 2 limbs. It is the sum of the diagonals + n-1 <= p+q < 2n-1 shifted down by n-1; the high n limbs (rp[2..n+2)) are + exact and the low two are guard limbs. Reference: accumulate exactly that + band and compare the high n limbs. +*/ +TEST_FUNCTION_START(flint_mpn_mulmid_n, state) +{ + slong ix; + +#if !FLINT_HAVE_NATIVE_mpn_mulmid_n + TEST_FUNCTION_END_SKIPPED(state); +#else + for (ix = 0; ix < 10000 * flint_test_multiplier(); ix++) + { + mp_size_t n, an, p, q; + mp_ptr a, b, rp, acc; + + n = 2 + n_randint(state, 40); + if (n_randint(state, 500) == 0) + n = 2 + n_randint(state, 400); + an = 2 * n - 1; + + a = flint_malloc(sizeof(mp_limb_t) * an); + b = flint_malloc(sizeof(mp_limb_t) * n); + rp = flint_malloc(sizeof(mp_limb_t) * (n + 2)); + acc = flint_calloc(n + 3, sizeof(mp_limb_t)); + + flint_mpn_rrandom(a, state, an); + flint_mpn_rrandom(b, state, n); + flint_mpn_rrandom(rp, state, n + 2); /* poison */ + + flint_mpn_mulmid_n(rp, a, b, n); + + for (p = 0; p < an; p++) + for (q = 0; q < n; q++) + { + mp_size_t s = p + q, off; + mp_limb_t hi, lo; + + if (s < n - 1 || s >= 2 * n - 1) + continue; + + umul_ppmm(hi, lo, a[p], b[q]); + off = s - (n - 1); + mpn_add_1(acc + off, acc + off, (n + 3) - off, lo); + mpn_add_1(acc + off + 1, acc + off + 1, (n + 3) - off - 1, hi); + } + + if (mpn_cmp(rp + 2, acc + 2, n) != 0) + TEST_FUNCTION_FAIL( + "ix = %wd, n = %wd\n" + "a = %{ulong*}\nb = %{ulong*}\n" + "exact high limbs = %{ulong*}\n" + "got = %{ulong*}\n", + ix, n, a, an, b, n, acc + 2, n, rp + 2, n); + + flint_free(a); + flint_free(b); + flint_free(rp); + flint_free(acc); + } + + TEST_FUNCTION_END(state); +#endif +}