From 7ec32571f3a3dbeafa591d9e558b62fcd01ee3ff Mon Sep 17 00:00:00 2001 From: Marco Bodrato Date: Tue, 15 Feb 2022 08:57:56 +0100 Subject: tests/mpn/t-{mul,sqr}mod_bknp1.c: New tests for mpn_{mul,sqr}mod_bknp1 tests/mpn/Makefile.ami (check_PROGRAMS): Compile and run new tests --- tests/mpn/Makefile.am | 1 + tests/mpn/t-mulmod_bknp1.c | 195 +++++++++++++++++++++++++++++++++++++ tests/mpn/t-sqrmod_bknp1.c | 235 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 431 insertions(+) create mode 100644 tests/mpn/t-mulmod_bknp1.c create mode 100644 tests/mpn/t-sqrmod_bknp1.c (limited to 'tests') diff --git a/tests/mpn/Makefile.am b/tests/mpn/Makefile.am index 7e370b68c..eb7993335 100644 --- a/tests/mpn/Makefile.am +++ b/tests/mpn/Makefile.am @@ -28,6 +28,7 @@ check_PROGRAMS = t-asmtype t-aors_1 t-divrem_1 t-mod_1 t-fat t-get_d \ t-toom52 t-toom53 t-toom54 t-toom62 t-toom63 t-toom6h t-toom8h \ t-toom2-sqr t-toom3-sqr t-toom4-sqr t-toom6-sqr t-toom8-sqr \ t-div t-mul t-mullo t-sqrlo t-mulmod_bnm1 t-sqrmod_bnm1 t-mulmid \ + t-mulmod_bknp1 t-sqrmod_bknp1 \ t-addaddmul t-hgcd t-hgcd_appr t-matrix22 t-invert t-bdiv t-fib2m \ t-broot t-brootinv t-minvert t-sizeinbase t-gcd_11 t-gcd_22 t-gcdext_1 diff --git a/tests/mpn/t-mulmod_bknp1.c b/tests/mpn/t-mulmod_bknp1.c new file mode 100644 index 000000000..f07b85489 --- /dev/null +++ b/tests/mpn/t-mulmod_bknp1.c @@ -0,0 +1,195 @@ +/* Test for mulmod_bknp1 function. + + Contributed to the GNU project by Marco Bodrato. + +Copyright 2009, 2020-2022 Free Software Foundation, Inc. + +This file is part of the GNU MP Library test suite. + +The GNU MP Library test suite is free software; you can redistribute it +and/or modify it under the terms of the GNU General Public License as +published by the Free Software Foundation; either version 3 of the License, +or (at your option) any later version. + +The GNU MP Library test suite is distributed in the hope that it will be +useful, but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +Public License for more details. + +You should have received a copy of the GNU General Public License along with +the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */ + + +#include +#include + +#include "gmp-impl.h" +#include "tests.h" + +#if MOD_BKNP1_USE11 +#define USE11 11, +#else +#define USE11 +#endif + + +#if GMP_NUMB_BITS % 32 == 0 +#define MAX_K 17 +#define SUPPORTED_K {3, 5, 7, 13, USE11 MAX_K} +#else +#if GMP_NUMB_BITS % 16 == 0 +#define MAX_K 13 +#define SUPPORTED_K {3, 5, 7, USE11 MAX_K} +#else +#if GMP_NUMB_BITS % 8 == 0 +#define MAX_K 7 +#define SUPPORTED_K {3, USE11 MAX_K} +#else +#define SUPPORTED_K {USE11} /* Supported ? */ +#endif /* GMP_NUMB_BITS % 8 == 0 */ +#endif /* GMP_NUMB_BITS % 16 == 0 */ +#endif /* GMP_NUMB_BITS % 32 == 0 */ + +#if MOD_BKNP1_ONLY3 +#undef SUPPORTED_K +#undef MAX_K +#define MAX_K 3 +#define SUPPORTED_K {3} +#endif + +/* Sizes are up to MAX_K * 2^SIZE_LOG limbs */ +#ifndef SIZE_LOG +#define SIZE_LOG 7 +#endif + +#ifndef COUNT +#define COUNT 5000 +#endif + +#define MAX_N (MAX_K << SIZE_LOG) +#define MIN_N 1 + +/* + Reference function for multiplication modulo B^{k*rn}+1. +*/ + +static void +ref_mulmod_bnp1 (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t rn) +{ + mp_limb_t cy; + + mpn_mul_n (rp, ap, bp, rn + 1); + cy = rp[2 * rn]; + MPN_INCR_U (rp, 2 * rn + 1, rp[2 * rn]); + cy = rp[2 * rn] - cy + mpn_sub_n (rp, rp, rp + rn, rn); + rp[rn] = 0; + MPN_INCR_U (rp, rn + 1, cy); +} + +/* + Compare the result of the mpn_mulmod_bnp1 function in the library + with the reference function above. +*/ +unsigned supported_k[] = SUPPORTED_K; + +int +main (int argc, char **argv) +{ + mp_ptr ap, bp, refp, pp, scratch; + int count = COUNT; + int test; + gmp_randstate_ptr rands; + TMP_DECL; + TMP_MARK; + + TESTS_REPS (count, argv, argc); + + tests_start (); + rands = RANDS; + + ap = TMP_ALLOC_LIMBS (MAX_N + 1); + bp = TMP_ALLOC_LIMBS (MAX_N + 1); + refp = TMP_ALLOC_LIMBS (MAX_N * 2 + 2); + pp = 1 + TMP_ALLOC_LIMBS (MAX_N + 3); + scratch + = 1 + TMP_ALLOC_LIMBS (mpn_mulmod_bknp1_itch (MAX_N) + 2); + + for (test = 0; test < count; test++) + { + unsigned size_min; + unsigned size_range; + unsigned k; + mp_size_t an,bn,rn,n; + mp_size_t itch; + mp_limb_t p_before, p_after, s_before, s_after; + + for (size_min = 1; (1L << size_min) < MIN_N; size_min++) + ; + + /* We generate rn in the MIN_N <= n <= (1 << size_range). */ + size_range = size_min + + gmp_urandomm_ui (rands, SIZE_LOG + 1 - size_min); + + k = supported_k[test % numberof (supported_k)]; + n = MIN_N + + gmp_urandomm_ui (rands, (1L << size_range) + 1 - MIN_N); + rn = k * n; + if ((GMP_NUMB_MAX % k != 0) && (rn % 3 == 0)) + n = rn / (k = 3); + + mpn_random2 (ap, rn + 1); + mpn_random2 (bp, rn + 1); + + bp [rn] &= 1; + ap [rn] &= 1; + + mpn_random2 (pp-1, rn + 3); + p_before = pp[-1]; + p_after = pp[rn + 1]; + + itch = mpn_mulmod_bknp1_itch (rn); + ASSERT_ALWAYS (itch <= mpn_mulmod_bknp1_itch (MAX_N)); + mpn_random2 (scratch - 1, itch + 2); + s_before = scratch[-1]; + s_after = scratch[itch]; + + mpn_mulmod_bknp1 ( pp, ap, bp, n, k, scratch); + ref_mulmod_bnp1 (refp, ap, bp, rn); + if (pp[-1] != p_before || pp[rn + 1] != p_after + || scratch[-1] != s_before || scratch[itch] != s_after + || mpn_cmp (refp, pp, rn + 1) != 0) + { + printf ("ERROR in test %d, rn = %d, n = %d, k = %d\n", + test, (int) rn, (int) n, (int) k); + if (pp[-1] != p_before) + { + printf ("before pp:"); mpn_dump (pp - 1, 1); + printf ("keep: "); mpn_dump (&p_before, 1); + } + if (pp[rn + 1] != p_after) + { + printf ("after pp:"); mpn_dump (pp + rn + 1, 1); + printf ("keep: "); mpn_dump (&p_after, 1); + } + if (scratch[-1] != s_before) + { + printf ("before scratch:"); mpn_dump (scratch - 1, 1); + printf ("keep: "); mpn_dump (&s_before, 1); + } + if (scratch[itch] != s_after) + { + printf ("after scratch:"); mpn_dump (scratch + itch, 1); + printf ("keep: "); mpn_dump (&s_after, 1); + } + mpn_dump (ap, rn + 1); + mpn_dump (bp, rn + 1); + mpn_dump (pp, rn + 1); + mpn_dump (refp, rn + 1); + + abort(); + } + } + TMP_FREE; + tests_end (); + return 0; +} diff --git a/tests/mpn/t-sqrmod_bknp1.c b/tests/mpn/t-sqrmod_bknp1.c new file mode 100644 index 000000000..b28a55a05 --- /dev/null +++ b/tests/mpn/t-sqrmod_bknp1.c @@ -0,0 +1,235 @@ +/* Test for mulmod_bknp1 function. + + Contributed to the GNU project by Marco Bodrato. + +Copyright 2009, 2020-2022 Free Software Foundation, Inc. + +This file is part of the GNU MP Library test suite. + +The GNU MP Library test suite is free software; you can redistribute it +and/or modify it under the terms of the GNU General Public License as +published by the Free Software Foundation; either version 3 of the License, +or (at your option) any later version. + +The GNU MP Library test suite is distributed in the hope that it will be +useful, but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +Public License for more details. + +You should have received a copy of the GNU General Public License along with +the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */ + + +#include +#include + +#include "gmp-impl.h" +#include "tests.h" + + +#if MOD_BKNP1_USE11 +#define USE11 11, +#else +#define USE11 +#endif + +#if GMP_NUMB_BITS % 32 == 0 +#define MAX_K 17 +#define SUPPORTED_K {3, 5, 7, 13, USE11 MAX_K} +#else +#if GMP_NUMB_BITS % 16 == 0 +#define MAX_K 13 +#define SUPPORTED_K {3, 5, 7, USE11 MAX_K} +#else +#if GMP_NUMB_BITS % 8 == 0 +#define MAX_K 7 +#define SUPPORTED_K {3, USE11 MAX_K} +#else +#define SUPPORTED_K {USE11} /* Supported ? */ +#endif /* GMP_NUMB_BITS % 8 == 0 */ +#endif /* GMP_NUMB_BITS % 16 == 0 */ +#endif /* GMP_NUMB_BITS % 32 == 0 */ + +#if MOD_BKNP1_ONLY3 +#undef SUPPORTED_K +#undef MAX_K +#define MAX_K 3 +#define SUPPORTED_K {3} +#endif + +/* Sizes are up to MAX_K * 2^SIZE_LOG limbs */ +#ifndef SIZE_LOG +#define SIZE_LOG 7 +#endif + +#ifndef COUNT +#define COUNT 5000 +#endif + +#define MAX_N (MAX_K << SIZE_LOG) +#define MIN_N 1 + +/* + Reference function for multiplication modulo B^{k*rn}+1. +*/ + +static void +ref_sqrmod_bnp1 (mp_ptr rp, mp_srcptr ap, mp_size_t rn) +{ + mp_limb_t cy; + + mpn_sqr (rp, ap, rn + 1); + cy = rp[2 * rn]; + MPN_INCR_U (rp, 2 * rn + 1, rp[2 * rn]); + cy = rp[2 * rn] - cy + mpn_sub_n (rp, rp, rp + rn, rn); + rp[rn] = 0; + MPN_INCR_U (rp, rn + 1, cy); +} + +/* + Compare the result of the mpn_mulmod_bnp1 function in the library + with the reference function above. +*/ +unsigned supported_k[] = SUPPORTED_K; + +int +main (int argc, char **argv) +{ + mp_ptr ap, refp, pp, scratch; + int count = COUNT; + int test; + gmp_randstate_ptr rands; + TMP_DECL; + TMP_MARK; + + TESTS_REPS (count, argv, argc); + + tests_start (); + rands = RANDS; + + ap = TMP_ALLOC_LIMBS (MAX_N + 1); + refp = TMP_ALLOC_LIMBS (MAX_N * 2 + 2); + pp = 1 + TMP_ALLOC_LIMBS (MAX_N + 3); + scratch + = 1 + TMP_ALLOC_LIMBS (mpn_mulmod_bknp1_itch (MAX_N) + 2); + + for (test = 0; test < count; test++) + { + unsigned size_min; + unsigned size_range; + unsigned k; + mp_size_t an,rn,n; + mp_size_t itch; + mp_limb_t p_before, p_after, s_before, s_after; + + for (size_min = 1; (1L << size_min) < MIN_N; size_min++) + ; + + /* We generate rn in the MIN_N <= n <= (1 << size_range). */ + size_range = size_min + + gmp_urandomm_ui (rands, SIZE_LOG + 1 - size_min); + + k = supported_k[test % numberof (supported_k)]; + n = MIN_N + + gmp_urandomm_ui (rands, (1L << size_range) + 1 - MIN_N); + rn = k * n; + if ((GMP_NUMB_MAX % k != 0) && (rn % 3 == 0)) + n = rn / (k = 3); + + mpn_random2 (ap, rn + 1); + + ap [rn] &= 1; + + mpn_random2 (pp-1, rn + 3); + p_before = pp[-1]; + p_after = pp[rn + 1]; + + itch = mpn_sqrmod_bknp1_itch (rn); + ASSERT_ALWAYS (itch <= mpn_mulmod_bknp1_itch (MAX_N)); + mpn_random2 (scratch - 1, itch + 2); + s_before = scratch[-1]; + s_after = scratch[itch]; + + mpn_sqrmod_bknp1 ( pp, ap, n, k, scratch); + ref_sqrmod_bnp1 (refp, ap, rn); + if (pp[-1] != p_before || pp[rn + 1] != p_after + || scratch[-1] != s_before || scratch[itch] != s_after + || mpn_cmp (refp, pp, rn + 1) != 0) + { + printf ("ERROR in test %d(sqr), rn = %d, n = %d, k = %d\n", + test, (int) rn, (int) n, (int) k); + if (pp[-1] != p_before) + { + printf ("before pp:"); mpn_dump (pp - 1, 1); + printf ("keep: "); mpn_dump (&p_before, 1); + } + if (pp[rn + 1] != p_after) + { + printf ("after pp:"); mpn_dump (pp + rn + 1, 1); + printf ("keep: "); mpn_dump (&p_after, 1); + } + if (scratch[-1] != s_before) + { + printf ("before scratch:"); mpn_dump (scratch - 1, 1); + printf ("keep: "); mpn_dump (&s_before, 1); + } + if (scratch[itch] != s_after) + { + printf ("after scratch:"); mpn_dump (scratch + itch, 1); + printf ("keep: "); mpn_dump (&s_after, 1); + } + mpn_dump (ap, rn + 1); + mpn_dump (pp, rn + 1); + mpn_dump (refp, rn + 1); + + abort(); + } + + mpn_random2 (pp-1, rn + 3); + p_before = pp[-1]; + p_after = pp[rn + 1]; + + itch = mpn_mulmod_bknp1_itch (rn); + ASSERT_ALWAYS (itch <= mpn_mulmod_bknp1_itch (MAX_N)); + mpn_random2 (scratch - 1, itch + 2); + s_before = scratch[-1]; + s_after = scratch[itch]; + + mpn_mulmod_bknp1 ( pp, ap, ap, n, k, scratch); + if (pp[-1] != p_before || pp[rn + 1] != p_after + || scratch[-1] != s_before || scratch[itch] != s_after + || mpn_cmp (refp, pp, rn + 1) != 0) + { + printf ("ERROR in test %d(mul), rn = %d, n = %d, k = %d\n", + test, (int) rn, (int) n, (int) k); + if (pp[-1] != p_before) + { + printf ("before pp:"); mpn_dump (pp - 1, 1); + printf ("keep: "); mpn_dump (&p_before, 1); + } + if (pp[rn + 1] != p_after) + { + printf ("after pp:"); mpn_dump (pp + rn + 1, 1); + printf ("keep: "); mpn_dump (&p_after, 1); + } + if (scratch[-1] != s_before) + { + printf ("before scratch:"); mpn_dump (scratch - 1, 1); + printf ("keep: "); mpn_dump (&s_before, 1); + } + if (scratch[itch] != s_after) + { + printf ("after scratch:"); mpn_dump (scratch + itch, 1); + printf ("keep: "); mpn_dump (&s_after, 1); + } + mpn_dump (ap, rn + 1); + mpn_dump (pp, rn + 1); + mpn_dump (refp, rn + 1); + + abort(); + } + } + TMP_FREE; + tests_end (); + return 0; +} -- cgit v1.2.1