Generic mpn_mullow and hardcoded ones for x86

.gitignore

@@ -64,3 +64,5 @@ flint.pc
 autom4te.cache/
 config.m4
 src/flint-mparam.h
+.gdb_history
+vgcore.*

Makefile.in

@@ -56,6 +56,7 @@ RM_F:=rm -f
 RM_RF:=rm -rf
 CP:=cp
 CP_A:=cp -a
+GDB:=gdb
 
 STATIC:=@STATIC@
 SHARED:=@SHARED@
@@ -498,6 +499,16 @@ else
 -include $(BUILD_DIR)/test/*.d
 -include $(BUILD_DIR)/*/test/*.d
 endif
+else ifeq ($(MAKECMDGOALS), debug)
+ifdef MOD
+$(warning Dependency tracking only set to cover the test executables of $(MOD).)
+-include $(foreach dir, $(MOD), $(BUILD_DIR)/$(dir)/test/*.d)
+else
+-include $(BUILD_DIR)/*/*.o.d
+-include $(BUILD_DIR)/*/*.lo.d
+-include $(BUILD_DIR)/test/*.d
+-include $(BUILD_DIR)/*/test/*.d
+endif
 else ifeq ($(MAKECMDGOALS), tune)
 -include $(BUILD_DIR)/*/*.o.d
 -include $(BUILD_DIR)/*/*.lo.d
@@ -798,6 +809,25 @@ check: library $(TESTS:%=%_TEST_RUN)
 	@echo 'All tests passed.'
 endif
 
+################################################################################
+# debugging
+################################################################################
+
+%_TEST_DBG_RUN_ARGS: %
+	@$(GDB) --args $< $(ARGS)
+
+ifdef MOD
+ifdef ARGS
+DEBUG:=1
+debug: library $(patsubst %,%_TEST_DBG_RUN_ARGS, $($(sort $(MOD))_TESTS))
+endif
+endif
+
+ifneq ($(DEBUG),1)
+debug:
+	$(error Can only run debugger with one module and one argument at a time)
+endif
+
 ################################################################################
 # tuning
 ################################################################################
@@ -952,11 +982,11 @@ dist:
 	dev/make_dist.sh $(FLINT_VERSION)
 
 ################################################################################
-# debugging
+# makefile debugging
 ################################################################################
 
 print-%:
 	@echo "$*=$($*)"
 
-.PHONY: all library shared static examples checkexamples profile tests check tune valgrind clean distclean install uninstall dist %_TEST_RUN %_VALGRIND_RUN print-% coverage coverage_html
-.SILENT: $(mpn_extras_S_SOURCES)
+.PHONY: all library shared static examples checkexamples profile tests check tune valgrind clean distclean install uninstall dist %_TEST_RUN %_TEST_RUN_% %_TEST_DGB_RUN_ARGS %_VALGRIND_RUN print-% coverage coverage_html debug
+.PRECIOUS: $(mpn_extras_PIC_S_SOURCES) $(mpn_extras_S_SOURCES)

src/mpn_extras.h

@@ -370,6 +370,36 @@ flint_mpn_sqr(mp_ptr r, mp_srcptr x, mp_size_t n)
         flint_mpn_mul((_z), (_y), (_yn), (_x), (_xn)); \
     }
 
+/* Low multiplication ********************************************************/
+
+#define FLINT_HAVE_MULLOW_FUNC(n) ((n) <= FLINT_MPN_MULLOW_FUNC_TAB_WIDTH)
+
+FLINT_DLL extern const flint_mpn_mul_func_t flint_mpn_mullow_func_tab[];
+
+mp_limb_t flint_mpn_mullow_basecase(mp_ptr, mp_srcptr, mp_srcptr, mp_size_t);
+
+#if FLINT_HAVE_ASSEMBLY_x86_64_adx
+# define FLINT_MPN_MULLOW_FUNC_TAB_WIDTH 8
+# define FLINT_HAVE_NATIVE_mpn_mullow_basecase 1
+#else
+# define FLINT_MPN_MULLOW_FUNC_TAB_WIDTH 0
+#endif
+
+/* TODO: Fix higher stuff */
+MPN_EXTRAS_INLINE
+mp_limb_t flint_mpn_mullow_n(mp_ptr rp, mp_srcptr xp, mp_srcptr yp, mp_size_t n)
+{
+    FLINT_ASSERT(n >= 1);
+
+    if (FLINT_HAVE_MULLOW_FUNC(n))
+    {
+        FLINT_ASSERT(rp != xp);
+        return flint_mpn_mullow_func_tab[n](rp, xp, yp);
+    }
+    else
+        return flint_mpn_mullow_basecase(rp, xp, yp, n);
+}
+
 /* High multiplication *******************************************************/
 
 #define FLINT_HAVE_MULHIGH_FUNC(n) ((n) <= FLINT_MPN_MULHIGH_FUNC_TAB_WIDTH)

src/mpn_extras/mullow.c

@@ -0,0 +1,57 @@
+/*
+    Copyright (C) 2024 Albin Ahlbäck
+
+    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 <https://www.gnu.org/licenses/>.
+*/
+
+#include "mpn_extras.h"
+
+#if FLINT_HAVE_ASSEMBLY_x86_64_adx
+mp_limb_t flint_mpn_mullow_1(mp_ptr, mp_srcptr, mp_srcptr);
+mp_limb_t flint_mpn_mullow_2(mp_ptr, mp_srcptr, mp_srcptr);
+mp_limb_t flint_mpn_mullow_3(mp_ptr, mp_srcptr, mp_srcptr);
+mp_limb_t flint_mpn_mullow_4(mp_ptr, mp_srcptr, mp_srcptr);
+mp_limb_t flint_mpn_mullow_5(mp_ptr, mp_srcptr, mp_srcptr);
+mp_limb_t flint_mpn_mullow_6(mp_ptr, mp_srcptr, mp_srcptr);
+mp_limb_t flint_mpn_mullow_7(mp_ptr, mp_srcptr, mp_srcptr);
+mp_limb_t flint_mpn_mullow_8(mp_ptr, mp_srcptr, mp_srcptr);
+
+const flint_mpn_mul_func_t flint_mpn_mullow_func_tab[] =
+{
+    NULL,
+    flint_mpn_mullow_1,
+    flint_mpn_mullow_2,
+    flint_mpn_mullow_3,
+    flint_mpn_mullow_4,
+    flint_mpn_mullow_5,
+    flint_mpn_mullow_6,
+    flint_mpn_mullow_7,
+    flint_mpn_mullow_8
+};
+#else
+const flint_mpn_mul_func_t flint_mpn_mullow_func_tab[] = { NULL };
+#endif
+
+#if !FLINT_HAVE_NATIVE_mpn_mullow_basecase
+mp_limb_t
+flint_mpn_mullow_basecase(mp_ptr rp, mp_srcptr xp, mp_srcptr yp, mp_size_t n)
+{
+    mp_limb_t ret;
+    mp_size_t ix;
+
+    ret = mpn_mul_1(rp, xp, n, yp[0]);
+
+    for (ix = 1; ix < n; ix++)
+    {
+        ret += mpn_addmul_1(rp + ix, xp, n - ix, yp[ix]);
+        ret += xp[n - ix] * yp[ix];
+    }
+
+    return ret;
+}
+#endif

src/mpn_extras/profile/p-mullow.c

@@ -0,0 +1,74 @@
+/*
+    Copyright (C) 2024 Albin Ahlbäck
+
+    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 <https://www.gnu.org/licenses/>.
+*/
+
+#include "mpn_extras.h"
+#include "profiler.h"
+
+#define mpn_mullo_basecase __gmpn_mullo_basecase
+void mpn_mullo_basecase(mp_ptr, mp_srcptr, mp_srcptr, mp_size_t);
+
+#define N_MIN 1
+#define N_MAX 30
+
+#if N_MAX < 9
+# define P 1
+#elif N_MAX < 99
+# define P 2
+#elif N_MAX < 999
+# define P 3
+#elif N_MAX < 9999
+# define P 4
+#else
+# define P 5
+#endif
+
+#define _STR(x) #x
+#define STR(x) _STR(x)
+
+int
+main(void)
+{
+    mp_limb_t rp[2 * N_MAX];
+    mp_limb_t xp[N_MAX];
+    mp_limb_t yp[N_MAX];
+    mp_size_t n;
+
+    flint_printf("%.*s      mullo_basecase / FLINT || mul / mullow\n", P, "                              ");
+    for (n = N_MIN; n <= N_MAX; n++)
+    {
+        double t1, t2, t3, FLINT_SET_BUT_UNUSED(__);
+        flint_printf("n = %" STR(P) "wd:", n);
+
+        mpn_random2(xp, n);
+        mpn_random2(yp, n);
+
+        TIMEIT_START
+        mpn_mullo_basecase(rp, xp, yp, n);
+        TIMEIT_STOP_VALUES(__, t1)
+
+        TIMEIT_START
+        flint_mpn_mul_n(rp, xp, yp, n);
+        TIMEIT_STOP_VALUES(__, t2)
+
+        TIMEIT_START
+        flint_mpn_mullow_n(rp, xp, yp, n);
+        TIMEIT_STOP_VALUES(__, t3)
+
+        flint_printf("         %7.2fx       ||  %7.2f\n", t1 / t3, t2 / t3);
+    }
+
+    flint_cleanup_master();
+
+    return 0;
+}
+
+#undef N_MIN
+#undef N_MAX

src/mpn_extras/test/main.c

@@ -24,6 +24,7 @@
 #include "t-mul.c"
 #include "t-mul_n.c"
 #include "t-mul_toom22.c"
+#include "t-mullow_n.c"
 #include "t-mulhigh.c"
 #include "t-mulhigh_normalised.c"
 #include "t-mulmod_2expp1.c"
@@ -48,6 +49,7 @@ test_struct tests[] =
     TEST_FUNCTION(flint_mpn_mul),
     TEST_FUNCTION(flint_mpn_mul_n),
     TEST_FUNCTION(flint_mpn_mul_toom22),
+    TEST_FUNCTION(flint_mpn_mullow_n),
     TEST_FUNCTION(flint_mpn_mulhigh),
     TEST_FUNCTION(flint_mpn_mulhigh_normalised),
     TEST_FUNCTION(flint_mpn_mulmod_2expp1),

src/mpn_extras/test/t-mullow_n.c

@@ -0,0 +1,69 @@
+/*
+    Copyright (C) 2024 Albin Ahlbäck
+
+    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 <https://www.gnu.org/licenses/>.
+*/
+
+#include "test_helpers.h"
+#include "mpn_extras.h"
+
+#define N_MIN 1
+#define N_MAX 30
+
+TEST_FUNCTION_START(flint_mpn_mullow_n, state)
+{
+    slong ix;
+    int result;
+
+    mp_ptr rp, rpf, xp, yp;
+
+    for (ix = 0; ix < 100000 * flint_test_multiplier(); ix++)
+    {
+        mp_limb_t ret;
+        mp_size_t n;
+
+        n = N_MIN + n_randint(state, N_MAX - N_MIN + 1);
+
+        rp = flint_malloc(sizeof(mp_limb_t) * n);
+        rpf = flint_malloc(2 * sizeof(mp_limb_t) * n);
+        xp = flint_malloc(sizeof(mp_limb_t) * n);
+        yp = flint_malloc(sizeof(mp_limb_t) * n);
+
+        mpn_random2(xp, n);
+        mpn_random2(yp, n);
+
+        if (ix < 1)
+            continue;
+
+        ret = flint_mpn_mullow_n(rp, xp, yp, n);
+        flint_mpn_mul_n(rpf, xp, yp, n);
+
+        result = (mpn_cmp(rp, rpf, n) == 0 && ret == rpf[n]);
+        if (!result)
+            TEST_FUNCTION_FAIL(
+                    "ix = %wd\n"
+                    "n = %wd\n"
+                    "xp = %{ulong*}\n"
+                    "yp = %{ulong*}\n"
+                    "Exp ret: %{ulong}\n"
+                    "Got ret: %{ulong}\n"
+                    "Expected: %{ulong*}\n"
+                    "Got:      %{ulong*}\n",
+                    ix, n, xp, n, yp, n, rpf[n], ret, rpf, n, rp, n);
+
+        flint_free(rp);
+        flint_free(rpf);
+        flint_free(xp);
+        flint_free(yp);
+    }
+
+    TEST_FUNCTION_END(state);
+}
+
+#undef N_MIN
+#undef N_MAX