Blame mpn/generic/mod_1.c

Packit 5c3484
/* mpn_mod_1(dividend_ptr, dividend_size, divisor_limb) --
Packit 5c3484
   Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
Packit 5c3484
   Return the single-limb remainder.
Packit 5c3484
   There are no constraints on the value of the divisor.
Packit 5c3484
Packit 5c3484
Copyright 1991, 1993, 1994, 1999, 2000, 2002, 2007-2009, 2012 Free Software
Packit 5c3484
Foundation, Inc.
Packit 5c3484
Packit 5c3484
This file is part of the GNU MP Library.
Packit 5c3484
Packit 5c3484
The GNU MP Library is free software; you can redistribute it and/or modify
Packit 5c3484
it under the terms of either:
Packit 5c3484
Packit 5c3484
  * the GNU Lesser General Public License as published by the Free
Packit 5c3484
    Software Foundation; either version 3 of the License, or (at your
Packit 5c3484
    option) any later version.
Packit 5c3484
Packit 5c3484
or
Packit 5c3484
Packit 5c3484
  * the GNU General Public License as published by the Free Software
Packit 5c3484
    Foundation; either version 2 of the License, or (at your option) any
Packit 5c3484
    later version.
Packit 5c3484
Packit 5c3484
or both in parallel, as here.
Packit 5c3484
Packit 5c3484
The GNU MP Library is distributed in the hope that it will be useful, but
Packit 5c3484
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
Packit 5c3484
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
Packit 5c3484
for more details.
Packit 5c3484
Packit 5c3484
You should have received copies of the GNU General Public License and the
Packit 5c3484
GNU Lesser General Public License along with the GNU MP Library.  If not,
Packit 5c3484
see https://www.gnu.org/licenses/.  */
Packit 5c3484
Packit 5c3484
#include "gmp.h"
Packit 5c3484
#include "gmp-impl.h"
Packit 5c3484
#include "longlong.h"
Packit 5c3484
Packit 5c3484
Packit 5c3484
/* The size where udiv_qrnnd_preinv should be used rather than udiv_qrnnd,
Packit 5c3484
   meaning the quotient size where that should happen, the quotient size
Packit 5c3484
   being how many udiv divisions will be done.
Packit 5c3484
Packit 5c3484
   The default is to use preinv always, CPUs where this doesn't suit have
Packit 5c3484
   tuned thresholds.  Note in particular that preinv should certainly be
Packit 5c3484
   used if that's the only division available (USE_PREINV_ALWAYS).  */
Packit 5c3484
Packit 5c3484
#ifndef MOD_1_NORM_THRESHOLD
Packit 5c3484
#define MOD_1_NORM_THRESHOLD  0
Packit 5c3484
#endif
Packit 5c3484
Packit 5c3484
#ifndef MOD_1_UNNORM_THRESHOLD
Packit 5c3484
#define MOD_1_UNNORM_THRESHOLD  0
Packit 5c3484
#endif
Packit 5c3484
Packit 5c3484
#ifndef MOD_1U_TO_MOD_1_1_THRESHOLD
Packit 5c3484
#define MOD_1U_TO_MOD_1_1_THRESHOLD  MP_SIZE_T_MAX /* default is not to use mpn_mod_1s */
Packit 5c3484
#endif
Packit 5c3484
Packit 5c3484
#ifndef MOD_1N_TO_MOD_1_1_THRESHOLD
Packit 5c3484
#define MOD_1N_TO_MOD_1_1_THRESHOLD  MP_SIZE_T_MAX /* default is not to use mpn_mod_1s */
Packit 5c3484
#endif
Packit 5c3484
Packit 5c3484
#ifndef MOD_1_1_TO_MOD_1_2_THRESHOLD
Packit 5c3484
#define MOD_1_1_TO_MOD_1_2_THRESHOLD  10
Packit 5c3484
#endif
Packit 5c3484
Packit 5c3484
#ifndef MOD_1_2_TO_MOD_1_4_THRESHOLD
Packit 5c3484
#define MOD_1_2_TO_MOD_1_4_THRESHOLD  20
Packit 5c3484
#endif
Packit 5c3484
Packit 5c3484
#if TUNE_PROGRAM_BUILD && !HAVE_NATIVE_mpn_mod_1_1p
Packit 5c3484
/* Duplicates declarations in tune/speed.h */
Packit 5c3484
mp_limb_t mpn_mod_1_1p_1 (mp_srcptr, mp_size_t, mp_limb_t, mp_limb_t [4]);
Packit 5c3484
mp_limb_t mpn_mod_1_1p_2 (mp_srcptr, mp_size_t, mp_limb_t, mp_limb_t [4]);
Packit 5c3484
Packit 5c3484
void mpn_mod_1_1p_cps_1 (mp_limb_t [4], mp_limb_t);
Packit 5c3484
void mpn_mod_1_1p_cps_2 (mp_limb_t [4], mp_limb_t);
Packit 5c3484
Packit 5c3484
#undef mpn_mod_1_1p
Packit 5c3484
#define mpn_mod_1_1p(ap, n, b, pre)			     \
Packit 5c3484
  (mod_1_1p_method == 1 ? mpn_mod_1_1p_1 (ap, n, b, pre)     \
Packit 5c3484
   : (mod_1_1p_method == 2 ? mpn_mod_1_1p_2 (ap, n, b, pre)  \
Packit 5c3484
      : __gmpn_mod_1_1p (ap, n, b, pre)))
Packit 5c3484
Packit 5c3484
#undef mpn_mod_1_1p_cps
Packit 5c3484
#define mpn_mod_1_1p_cps(pre, b)				\
Packit 5c3484
  (mod_1_1p_method == 1 ? mpn_mod_1_1p_cps_1 (pre, b)		\
Packit 5c3484
   : (mod_1_1p_method == 2 ? mpn_mod_1_1p_cps_2 (pre, b)	\
Packit 5c3484
      : __gmpn_mod_1_1p_cps (pre, b)))
Packit 5c3484
#endif /* TUNE_PROGRAM_BUILD && !HAVE_NATIVE_mpn_mod_1_1p */
Packit 5c3484
Packit 5c3484
Packit 5c3484
/* The comments in mpn/generic/divrem_1.c apply here too.
Packit 5c3484
Packit 5c3484
   As noted in the algorithms section of the manual, the shifts in the loop
Packit 5c3484
   for the unnorm case can be avoided by calculating r = a%(d*2^n), followed
Packit 5c3484
   by a final (r*2^n)%(d*2^n).  In fact if it happens that a%(d*2^n) can
Packit 5c3484
   skip a division where (a*2^n)%(d*2^n) can't then there's the same number
Packit 5c3484
   of divide steps, though how often that happens depends on the assumed
Packit 5c3484
   distributions of dividend and divisor.  In any case this idea is left to
Packit 5c3484
   CPU specific implementations to consider.  */
Packit 5c3484
Packit 5c3484
static mp_limb_t
Packit 5c3484
mpn_mod_1_unnorm (mp_srcptr up, mp_size_t un, mp_limb_t d)
Packit 5c3484
{
Packit 5c3484
  mp_size_t  i;
Packit 5c3484
  mp_limb_t  n1, n0, r;
Packit 5c3484
  mp_limb_t  dummy;
Packit 5c3484
  int cnt;
Packit 5c3484
Packit 5c3484
  ASSERT (un > 0);
Packit 5c3484
  ASSERT (d != 0);
Packit 5c3484
Packit 5c3484
  d <<= GMP_NAIL_BITS;
Packit 5c3484
Packit 5c3484
  /* Skip a division if high < divisor.  Having the test here before
Packit 5c3484
     normalizing will still skip as often as possible.  */
Packit 5c3484
  r = up[un - 1] << GMP_NAIL_BITS;
Packit 5c3484
  if (r < d)
Packit 5c3484
    {
Packit 5c3484
      r >>= GMP_NAIL_BITS;
Packit 5c3484
      un--;
Packit 5c3484
      if (un == 0)
Packit 5c3484
	return r;
Packit 5c3484
    }
Packit 5c3484
  else
Packit 5c3484
    r = 0;
Packit 5c3484
Packit 5c3484
  /* If udiv_qrnnd doesn't need a normalized divisor, can use the simple
Packit 5c3484
     code above. */
Packit 5c3484
  if (! UDIV_NEEDS_NORMALIZATION
Packit 5c3484
      && BELOW_THRESHOLD (un, MOD_1_UNNORM_THRESHOLD))
Packit 5c3484
    {
Packit 5c3484
      for (i = un - 1; i >= 0; i--)
Packit 5c3484
	{
Packit 5c3484
	  n0 = up[i] << GMP_NAIL_BITS;
Packit 5c3484
	  udiv_qrnnd (dummy, r, r, n0, d);
Packit 5c3484
	  r >>= GMP_NAIL_BITS;
Packit 5c3484
	}
Packit 5c3484
      return r;
Packit 5c3484
    }
Packit 5c3484
Packit 5c3484
  count_leading_zeros (cnt, d);
Packit 5c3484
  d <<= cnt;
Packit 5c3484
Packit 5c3484
  n1 = up[un - 1] << GMP_NAIL_BITS;
Packit 5c3484
  r = (r << cnt) | (n1 >> (GMP_LIMB_BITS - cnt));
Packit 5c3484
Packit 5c3484
  if (UDIV_NEEDS_NORMALIZATION
Packit 5c3484
      && BELOW_THRESHOLD (un, MOD_1_UNNORM_THRESHOLD))
Packit 5c3484
    {
Packit 5c3484
      mp_limb_t nshift;
Packit 5c3484
      for (i = un - 2; i >= 0; i--)
Packit 5c3484
	{
Packit 5c3484
	  n0 = up[i] << GMP_NAIL_BITS;
Packit 5c3484
	  nshift = (n1 << cnt) | (n0 >> (GMP_NUMB_BITS - cnt));
Packit 5c3484
	  udiv_qrnnd (dummy, r, r, nshift, d);
Packit 5c3484
	  r >>= GMP_NAIL_BITS;
Packit 5c3484
	  n1 = n0;
Packit 5c3484
	}
Packit 5c3484
      udiv_qrnnd (dummy, r, r, n1 << cnt, d);
Packit 5c3484
      r >>= GMP_NAIL_BITS;
Packit 5c3484
      return r >> cnt;
Packit 5c3484
    }
Packit 5c3484
  else
Packit 5c3484
    {
Packit 5c3484
      mp_limb_t inv, nshift;
Packit 5c3484
      invert_limb (inv, d);
Packit 5c3484
Packit 5c3484
      for (i = un - 2; i >= 0; i--)
Packit 5c3484
	{
Packit 5c3484
	  n0 = up[i] << GMP_NAIL_BITS;
Packit 5c3484
	  nshift = (n1 << cnt) | (n0 >> (GMP_NUMB_BITS - cnt));
Packit 5c3484
	  udiv_rnnd_preinv (r, r, nshift, d, inv);
Packit 5c3484
	  r >>= GMP_NAIL_BITS;
Packit 5c3484
	  n1 = n0;
Packit 5c3484
	}
Packit 5c3484
      udiv_rnnd_preinv (r, r, n1 << cnt, d, inv);
Packit 5c3484
      r >>= GMP_NAIL_BITS;
Packit 5c3484
      return r >> cnt;
Packit 5c3484
    }
Packit 5c3484
}
Packit 5c3484
Packit 5c3484
static mp_limb_t
Packit 5c3484
mpn_mod_1_norm (mp_srcptr up, mp_size_t un, mp_limb_t d)
Packit 5c3484
{
Packit 5c3484
  mp_size_t  i;
Packit 5c3484
  mp_limb_t  n0, r;
Packit 5c3484
  mp_limb_t  dummy;
Packit 5c3484
Packit 5c3484
  ASSERT (un > 0);
Packit 5c3484
Packit 5c3484
  d <<= GMP_NAIL_BITS;
Packit 5c3484
Packit 5c3484
  ASSERT (d & GMP_LIMB_HIGHBIT);
Packit 5c3484
Packit 5c3484
  /* High limb is initial remainder, possibly with one subtract of
Packit 5c3484
     d to get r
Packit 5c3484
  r = up[un - 1] << GMP_NAIL_BITS;
Packit 5c3484
  if (r >= d)
Packit 5c3484
    r -= d;
Packit 5c3484
  r >>= GMP_NAIL_BITS;
Packit 5c3484
  un--;
Packit 5c3484
  if (un == 0)
Packit 5c3484
    return r;
Packit 5c3484
Packit 5c3484
  if (BELOW_THRESHOLD (un, MOD_1_NORM_THRESHOLD))
Packit 5c3484
    {
Packit 5c3484
      for (i = un - 1; i >= 0; i--)
Packit 5c3484
	{
Packit 5c3484
	  n0 = up[i] << GMP_NAIL_BITS;
Packit 5c3484
	  udiv_qrnnd (dummy, r, r, n0, d);
Packit 5c3484
	  r >>= GMP_NAIL_BITS;
Packit 5c3484
	}
Packit 5c3484
      return r;
Packit 5c3484
    }
Packit 5c3484
  else
Packit 5c3484
    {
Packit 5c3484
      mp_limb_t  inv;
Packit 5c3484
      invert_limb (inv, d);
Packit 5c3484
      for (i = un - 1; i >= 0; i--)
Packit 5c3484
	{
Packit 5c3484
	  n0 = up[i] << GMP_NAIL_BITS;
Packit 5c3484
	  udiv_rnnd_preinv (r, r, n0, d, inv);
Packit 5c3484
	  r >>= GMP_NAIL_BITS;
Packit 5c3484
	}
Packit 5c3484
      return r;
Packit 5c3484
    }
Packit 5c3484
}
Packit 5c3484
Packit 5c3484
mp_limb_t
Packit 5c3484
mpn_mod_1 (mp_srcptr ap, mp_size_t n, mp_limb_t b)
Packit 5c3484
{
Packit 5c3484
  ASSERT (n >= 0);
Packit 5c3484
  ASSERT (b != 0);
Packit 5c3484
Packit 5c3484
  /* Should this be handled at all?  Rely on callers?  Note un==0 is currently
Packit 5c3484
     required by mpz/fdiv_r_ui.c and possibly other places.  */
Packit 5c3484
  if (n == 0)
Packit 5c3484
    return 0;
Packit 5c3484
Packit 5c3484
  if (UNLIKELY ((b & GMP_NUMB_HIGHBIT) != 0))
Packit 5c3484
    {
Packit 5c3484
      if (BELOW_THRESHOLD (n, MOD_1N_TO_MOD_1_1_THRESHOLD))
Packit 5c3484
	{
Packit 5c3484
	  return mpn_mod_1_norm (ap, n, b);
Packit 5c3484
	}
Packit 5c3484
      else
Packit 5c3484
	{
Packit 5c3484
	  mp_limb_t pre[4];
Packit 5c3484
	  mpn_mod_1_1p_cps (pre, b);
Packit 5c3484
	  return mpn_mod_1_1p (ap, n, b, pre);
Packit 5c3484
	}
Packit 5c3484
    }
Packit 5c3484
  else
Packit 5c3484
    {
Packit 5c3484
      if (BELOW_THRESHOLD (n, MOD_1U_TO_MOD_1_1_THRESHOLD))
Packit 5c3484
	{
Packit 5c3484
	  return mpn_mod_1_unnorm (ap, n, b);
Packit 5c3484
	}
Packit 5c3484
      else if (BELOW_THRESHOLD (n, MOD_1_1_TO_MOD_1_2_THRESHOLD))
Packit 5c3484
	{
Packit 5c3484
	  mp_limb_t pre[4];
Packit 5c3484
	  mpn_mod_1_1p_cps (pre, b);
Packit 5c3484
	  return mpn_mod_1_1p (ap, n, b << pre[1], pre);
Packit 5c3484
	}
Packit 5c3484
      else if (BELOW_THRESHOLD (n, MOD_1_2_TO_MOD_1_4_THRESHOLD) || UNLIKELY (b > GMP_NUMB_MASK / 4))
Packit 5c3484
	{
Packit 5c3484
	  mp_limb_t pre[5];
Packit 5c3484
	  mpn_mod_1s_2p_cps (pre, b);
Packit 5c3484
	  return mpn_mod_1s_2p (ap, n, b << pre[1], pre);
Packit 5c3484
	}
Packit 5c3484
      else
Packit 5c3484
	{
Packit 5c3484
	  mp_limb_t pre[7];
Packit 5c3484
	  mpn_mod_1s_4p_cps (pre, b);
Packit 5c3484
	  return mpn_mod_1s_4p (ap, n, b << pre[1], pre);
Packit 5c3484
	}
Packit 5c3484
    }
Packit 5c3484
}