Quelle  t-div.c

/* Copyright 2006, 2007, 2009, 2010, 2013-2015, 2018 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>  /* for strtol */
#include <stdio.h>  /* for printf */

#include "gmp-impl.h"
#include "longlong.h"
#include "tests/tests.h"


static void
dumpy (mp_srcptr p, mp_size_t n)
{
  mp_size_t i;
  if (n > 20)
    {
      for (i = n - 1; i >= n - 4; i--)
 {
   printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
   printf (" ");
 }
      printf ("... ");
      for (i = 3; i >= 0; i--)
 {
   printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
   printf (i == 0 ? "" : " ");
 }
    }
  else
    {
      for (i = n - 1; i >= 0; i--)
 {
   printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]);
   printf (i == 0 ? "" : " ");
 }
    }
  puts ("");
}

static signed long test;

static void
check_one (mp_ptr qp, mp_srcptr rp,
    mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn,
    const char *fname, mp_limb_t q_allowed_err)
{
  mp_size_t qn = nn - dn + 1;
  mp_ptr tp;
  const char *msg;
  const char *tvalue;
  mp_limb_t i;
  TMP_DECL;
  TMP_MARK;

  tp = TMP_ALLOC_LIMBS (nn + 1);
  if (dn >= qn)
    refmpn_mul (tp, dp, dn, qp, qn);
  else
    refmpn_mul (tp, qp, qn, dp, dn);

  for (i = 0; i < q_allowed_err && (tp[nn] > 0 || mpn_cmp (tp, np, nn) > 0); i++)
    ASSERT_NOCARRY (refmpn_sub (tp, tp, nn+1, dp, dn));

  if (tp[nn] > 0 || mpn_cmp (tp, np, nn) > 0)
    {
      msg = "q too large";
      tvalue = "Q*D";
    error:
      printf ("\r*******************************************************************************\n");
      printf ("%s failed test %ld: %s\n", fname, test, msg);
      printf ("N=    "); dumpy (np, nn);
      printf ("D=    "); dumpy (dp, dn);
      printf ("Q=    "); dumpy (qp, qn);
      if (rp)
 { printf ("R=    "); dumpy (rp, dn); }
      printf ("%5s=", tvalue); dumpy (tp, nn+1);
      printf ("nn = %ld, dn = %ld, qn = %ld\n", nn, dn, qn);
      abort ();
    }

  ASSERT_NOCARRY (refmpn_sub_n (tp, np, tp, nn));
  tvalue = "N-Q*D";
  if (!(nn == dn || mpn_zero_p (tp + dn, nn - dn)) || mpn_cmp (tp, dp, dn) >= 0)
    {
      msg = "q too small";
      goto error;
    }

  if (rp && mpn_cmp (rp, tp, dn) != 0)
    {
      msg = "r incorrect";
      goto error;
    }

  TMP_FREE;
}


/* These are *bit* sizes. */
#ifndef SIZE_LOG
#define SIZE_LOG 17
#endif
#define MAX_DN (1L << SIZE_LOG)
#define MAX_NN (1L << (SIZE_LOG + 1))

#define COUNT 200

int
main (int argc, char **argv)
{
  gmp_randstate_ptr rands;
  unsigned long maxnbits, maxdbits, nbits, dbits;
  mpz_t n, d, q, r, tz, junk;
  mp_size_t maxnn, maxdn, nn, dn, clearn, i;
  mp_ptr np, dup, dnp, qp, rp, junkp;
  mp_limb_t t;
  gmp_pi1_t dinv;
  long count = COUNT;
  mp_ptr scratch;
  mp_limb_t ran;
  mp_size_t alloc, itch;
  mp_limb_t rran0, rran1, qran0, qran1;
  TMP_DECL;

  TESTS_REPS (count, argv, argc);

  maxdbits = MAX_DN;
  maxnbits = MAX_NN;

  tests_start ();
  rands = RANDS;

  mpz_init (n);
  mpz_init (d);
  mpz_init (q);
  mpz_init (r);
  mpz_init (tz);
  mpz_init (junk);

  maxnn = maxnbits / GMP_NUMB_BITS + 1;
  maxdn = maxdbits / GMP_NUMB_BITS + 1;

  TMP_MARK;

  qp = TMP_ALLOC_LIMBS (maxnn + 2) + 1;
  rp = TMP_ALLOC_LIMBS (maxnn + 2) + 1;
  dnp = TMP_ALLOC_LIMBS (maxdn);

  alloc = 1;
  scratch = __GMP_ALLOCATE_FUNC_LIMBS (alloc);

  for (test = -300; test < count; test++)
    {
      nbits = urandom () % (maxnbits - GMP_NUMB_BITS) + 2 * GMP_NUMB_BITS;

      if (test < 0)
 dbits = (test + 300) % (nbits - 1) + 1;
      else
 dbits = urandom () % (nbits - 1) % maxdbits + 1;

#if RAND_UNIFORM
#define RANDFUNC mpz_urandomb
#else
#define RANDFUNC mpz_rrandomb
#endif

      do
 RANDFUNC (d, rands, dbits);
      while (mpz_sgn (d) == 0);
      dn = SIZ (d);
      dup = PTR (d);
      MPN_COPY (dnp, dup, dn);
      dnp[dn - 1] |= GMP_NUMB_HIGHBIT;

      if (test % 2 == 0)
 {
   RANDFUNC (n, rands, nbits);
   nn = SIZ (n);
   ASSERT_ALWAYS (nn >= dn);
 }
      else
 {
   do
     {
       RANDFUNC (q, rands, urandom () % (nbits - dbits + 1));
       RANDFUNC (r, rands, urandom () % mpz_sizeinbase (d, 2));
       mpz_mul (n, q, d);
       mpz_add (n, n, r);
       nn = SIZ (n);
     }
   while (nn > maxnn || nn < dn);
 }

      ASSERT_ALWAYS (nn <= maxnn);
      ASSERT_ALWAYS (dn <= maxdn);

      mpz_urandomb (junk, rands, nbits);
      junkp = PTR (junk);

      np = PTR (n);

      mpz_urandomb (tz, rands, 32);
      t = mpz_get_ui (tz);

      if (t % 17 == 0)
 {
   dnp[dn - 1] = GMP_NUMB_MAX;
   dup[dn - 1] = GMP_NUMB_MAX;
 }

      switch ((int) t % 16)
 {
 case 0:
   clearn = urandom () % nn;
   for (i = clearn; i < nn; i++)
     np[i] = 0;
   break;
 case 1:
   mpn_sub_1 (np + nn - dn, dnp, dn, urandom ());
   break;
 case 2:
   mpn_add_1 (np + nn - dn, dnp, dn, urandom ());
   break;
 }

      if (dn >= 2)
 invert_pi1 (dinv, dnp[dn - 1], dnp[dn - 2]);

      rran0 = urandom ();
      rran1 = urandom ();
      qran0 = urandom ();
      qran1 = urandom ();

      qp[-1] = qran0;
      qp[nn - dn + 1] = qran1;
      rp[-1] = rran0;

      ran = urandom ();

      if ((double) (nn - dn) * dn < 1e5)
 {
   /* Test mpn_sbpi1_div_qr */
   if (dn > 2)
     {
       MPN_COPY (rp, np, nn);
       if (nn > dn)
  MPN_COPY (qp, junkp, nn - dn);
       qp[nn - dn] = mpn_sbpi1_div_qr (qp, rp, nn, dnp, dn, dinv.inv32);
       check_one (qp, rp, np, nn, dnp, dn, "mpn_sbpi1_div_qr", 0);
     }

   /* Test mpn_sbpi1_divappr_q */
   if (dn > 2)
     {
       MPN_COPY (rp, np, nn);
       if (nn > dn)
  MPN_COPY (qp, junkp, nn - dn);
       qp[nn - dn] = mpn_sbpi1_divappr_q (qp, rp, nn, dnp, dn, dinv.inv32);
       check_one (qp, NULL, np, nn, dnp, dn, "mpn_sbpi1_divappr_q", 1);
     }

   /* Test mpn_sbpi1_div_q */
   if (dn > 2)
     {
       MPN_COPY (rp, np, nn);
       if (nn > dn)
  MPN_COPY (qp, junkp, nn - dn);
       qp[nn - dn] = mpn_sbpi1_div_q (qp, rp, nn, dnp, dn, dinv.inv32);
       check_one (qp, NULL, np, nn, dnp, dn, "mpn_sbpi1_div_q", 0);
     }

   /* Test mpn_sec_div_qr */
   itch = mpn_sec_div_qr_itch (nn, dn);
   if (itch + 1 > alloc)
     {
       scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
       alloc = itch + 1;
     }
   scratch[itch] = ran;
   MPN_COPY (rp, np, nn);
   if (nn >= dn)
     MPN_COPY (qp, junkp, nn - dn + 1);
   qp[nn - dn] = mpn_sec_div_qr (qp, rp, nn, dup, dn, scratch);
   ASSERT_ALWAYS (ran == scratch[itch]);
   check_one (qp, rp, np, nn, dup, dn, "mpn_sec_div_qr (unnorm)", 0);

   /* Test mpn_sec_div_r */
   itch = mpn_sec_div_r_itch (nn, dn);
   if (itch + 1 > alloc)
     {
       scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
       alloc = itch + 1;
     }
   scratch[itch] = ran;
   MPN_COPY (rp, np, nn);
   mpn_sec_div_r (rp, nn, dup, dn, scratch);
   ASSERT_ALWAYS (ran == scratch[itch]);
   /* Note: Since check_one cannot cope with remainder-only functions, we
     pass qp[] from the previous function, mpn_sec_div_qr.  */
   check_one (qp, rp, np, nn, dup, dn, "mpn_sec_div_r (unnorm)", 0);

   /* Normalised case, mpn_sec_div_qr */
   itch = mpn_sec_div_qr_itch (nn, dn);
   scratch[itch] = ran;

   MPN_COPY (rp, np, nn);
   if (nn >= dn)
     MPN_COPY (qp, junkp, nn - dn + 1);
   qp[nn - dn] = mpn_sec_div_qr (qp, rp, nn, dnp, dn, scratch);
   ASSERT_ALWAYS (ran == scratch[itch]);
   check_one (qp, rp, np, nn, dnp, dn, "mpn_sec_div_qr (norm)", 0);

   /* Normalised case, mpn_sec_div_r */
   itch = mpn_sec_div_r_itch (nn, dn);
   scratch[itch] = ran;
   MPN_COPY (rp, np, nn);
   mpn_sec_div_r (rp, nn, dnp, dn, scratch);
   ASSERT_ALWAYS (ran == scratch[itch]);
   /* Note: Since check_one cannot cope with remainder-only functions, we
     pass qp[] from the previous function, mpn_sec_div_qr.  */
   check_one (qp, rp, np, nn, dnp, dn, "mpn_sec_div_r (norm)", 0);
 }

      /* Test mpn_dcpi1_div_qr */
      if (dn >= 6 && nn - dn >= 3)
 {
   MPN_COPY (rp, np, nn);
   if (nn > dn)
     MPN_COPY (qp, junkp, nn - dn);
   qp[nn - dn] = mpn_dcpi1_div_qr (qp, rp, nn, dnp, dn, &dinv);
   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
   ASSERT_ALWAYS (rp[-1] == rran0);
   check_one (qp, rp, np, nn, dnp, dn, "mpn_dcpi1_div_qr", 0);
 }

      /* Test mpn_dcpi1_divappr_q */
      if (dn >= 6 && nn - dn >= 3)
 {
   MPN_COPY (rp, np, nn);
   if (nn > dn)
     MPN_COPY (qp, junkp, nn - dn);
   qp[nn - dn] = mpn_dcpi1_divappr_q (qp, rp, nn, dnp, dn, &dinv);
   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
   ASSERT_ALWAYS (rp[-1] == rran0);
   check_one (qp, NULL, np, nn, dnp, dn, "mpn_dcpi1_divappr_q", 1);
 }

      /* Test mpn_dcpi1_div_q */
      if (dn >= 6 && nn - dn >= 3)
 {
   MPN_COPY (rp, np, nn);
   if (nn > dn)
     MPN_COPY (qp, junkp, nn - dn);
   qp[nn - dn] = mpn_dcpi1_div_q (qp, rp, nn, dnp, dn, &dinv);
   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
   ASSERT_ALWAYS (rp[-1] == rran0);
   check_one (qp, NULL, np, nn, dnp, dn, "mpn_dcpi1_div_q", 0);
 }

     /* Test mpn_mu_div_qr */
      if (nn - dn > 2 && dn >= 2)
 {
   itch = mpn_mu_div_qr_itch (nn, dn, 0);
   if (itch + 1 > alloc)
     {
       scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
       alloc = itch + 1;
     }
   scratch[itch] = ran;
   MPN_COPY (qp, junkp, nn - dn);
   MPN_ZERO (rp, dn);
   rp[dn] = rran1;
   qp[nn - dn] = mpn_mu_div_qr (qp, rp, np, nn, dnp, dn, scratch);
   ASSERT_ALWAYS (ran == scratch[itch]);
   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
   ASSERT_ALWAYS (rp[-1] == rran0);  ASSERT_ALWAYS (rp[dn] == rran1);
   check_one (qp, rp, np, nn, dnp, dn, "mpn_mu_div_qr", 0);
 }

      /* Test mpn_mu_divappr_q */
      if (nn - dn > 2 && dn >= 2)
 {
   itch = mpn_mu_divappr_q_itch (nn, dn, 0);
   if (itch + 1 > alloc)
     {
       scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
       alloc = itch + 1;
     }
   scratch[itch] = ran;
   MPN_COPY (qp, junkp, nn - dn);
   qp[nn - dn] = mpn_mu_divappr_q (qp, np, nn, dnp, dn, scratch);
   ASSERT_ALWAYS (ran == scratch[itch]);
   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
   check_one (qp, NULL, np, nn, dnp, dn, "mpn_mu_divappr_q", 4);
 }

      /* Test mpn_mu_div_q */
      if (nn - dn > 2 && dn >= 2)
 {
   itch = mpn_mu_div_q_itch (nn, dn, 0);
   if (itch + 1> alloc)
     {
       scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
       alloc = itch + 1;
     }
   scratch[itch] = ran;
   MPN_COPY (qp, junkp, nn - dn);
   qp[nn - dn] = mpn_mu_div_q (qp, np, nn, dnp, dn, scratch);
   ASSERT_ALWAYS (ran == scratch[itch]);
   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
   check_one (qp, NULL, np, nn, dnp, dn, "mpn_mu_div_q", 0);
 }

      if (1)
 {
   itch = nn + 1;
   if (itch + 1> alloc)
     {
       scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1);
       alloc = itch + 1;
     }
   scratch[itch] = ran;
   MPN_COPY (qp, junkp, nn - dn + 1);
   mpn_div_q (qp, np, nn, dup, dn, scratch);
   ASSERT_ALWAYS (ran == scratch[itch]);
   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
   check_one (qp, NULL, np, nn, dup, dn, "mpn_div_q", 0);
 }

      if (dn >= 2 && nn >= 2)
 {
   mp_limb_t qh;

   /* mpn_divrem_2 */
   MPN_COPY (rp, np, nn);
   qp[nn - 2] = qp[nn-1] = qran1;

   qh = mpn_divrem_2 (qp, 0, rp, nn, dnp + dn - 2);
   ASSERT_ALWAYS (qp[nn - 2] == qran1);
   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - 1] == qran1);
   qp[nn - 2] = qh;
   check_one (qp, rp, np, nn, dnp + dn - 2, 2, "mpn_divrem_2", 0);

   /* Missing: divrem_2 with fraction limbs. */

   /* mpn_div_qr_2 */
   qp[nn - 2] = qran1;

   qh = mpn_div_qr_2 (qp, rp, np, nn, dup + dn - 2);
   ASSERT_ALWAYS (qp[nn - 2] == qran1);
   ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - 1] == qran1);
   qp[nn - 2] = qh;
   check_one (qp, rp, np, nn, dup + dn - 2, 2, "mpn_div_qr_2", 0);
 }
      if (dn >= 1 && nn >= 1)
 {
   /* mpn_div_qr_1 */
   mp_limb_t qh;
   qp[nn-1] = qran1;
   rp[0] = mpn_div_qr_1 (qp, &qh, np, nn, dnp[dn - 1]);
   ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - 1] == qran1);
   qp[nn - 1] = qh;
   check_one (qp, rp, np, nn,  dnp + dn - 1, 1, "mpn_div_qr_1", 0);
 }
    }

  __GMP_FREE_FUNC_LIMBS (scratch, alloc);

  TMP_FREE;

  mpz_clear (n);
  mpz_clear (d);
  mpz_clear (q);
  mpz_clear (r);
  mpz_clear (tz);
  mpz_clear (junk);

  tests_end ();
  return 0;
}