Quelle  t-mulmod_bknp1.c   Sprache: C

/* 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 <stdlib.h>
#include <stdio.h>

#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 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);

      if (test == 0)
 {
   mpn_random2 (ap, n);
   mpn_add_1 (ap + n, ap, n, 1); /* {ap,an} = -1 mod B+1 */
   MPN_ZERO (ap + 2 * n, rn - 2 * n + 1);
 }
      else
 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;
}