Blame stdlib/mod_1.c

Packit Service 82fcde
/* mpn_mod_1(dividend_ptr, dividend_size, divisor_limb) --
Packit Service 82fcde
   Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
Packit Service 82fcde
   Return the single-limb remainder.
Packit Service 82fcde
   There are no constraints on the value of the divisor.
Packit Service 82fcde
Packit Service 82fcde
Copyright (C) 1991-2018 Free Software Foundation, Inc.
Packit Service 82fcde
Packit Service 82fcde
This file is part of the GNU MP Library.
Packit Service 82fcde
Packit Service 82fcde
The GNU MP Library is free software; you can redistribute it and/or modify
Packit Service 82fcde
it under the terms of the GNU Lesser General Public License as published by
Packit Service 82fcde
the Free Software Foundation; either version 2.1 of the License, or (at your
Packit Service 82fcde
option) any later version.
Packit Service 82fcde
Packit Service 82fcde
The GNU MP Library is distributed in the hope that it will be useful, but
Packit Service 82fcde
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
Packit Service 82fcde
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
Packit Service 82fcde
License for more details.
Packit Service 82fcde
Packit Service 82fcde
You should have received a copy of the GNU Lesser General Public License
Packit Service 82fcde
along with the GNU MP Library; see the file COPYING.LIB.  If not, see
Packit Service 82fcde
<http://www.gnu.org/licenses/>.  */
Packit Service 82fcde
Packit Service 82fcde
#include <gmp.h>
Packit Service 82fcde
#include "gmp-impl.h"
Packit Service 82fcde
#include "longlong.h"
Packit Service 82fcde
Packit Service 82fcde
#ifndef UMUL_TIME
Packit Service 82fcde
#define UMUL_TIME 1
Packit Service 82fcde
#endif
Packit Service 82fcde
Packit Service 82fcde
#ifndef UDIV_TIME
Packit Service 82fcde
#define UDIV_TIME UMUL_TIME
Packit Service 82fcde
#endif
Packit Service 82fcde
Packit Service 82fcde
/* FIXME: We should be using invert_limb (or invert_normalized_limb)
Packit Service 82fcde
   here (not udiv_qrnnd).  */
Packit Service 82fcde
Packit Service 82fcde
mp_limb_t
Packit Service 82fcde
mpn_mod_1 (mp_srcptr dividend_ptr, mp_size_t dividend_size,
Packit Service 82fcde
	   mp_limb_t divisor_limb)
Packit Service 82fcde
{
Packit Service 82fcde
  mp_size_t i;
Packit Service 82fcde
  mp_limb_t n1, n0, r;
Packit Service 82fcde
  mp_limb_t dummy __attribute__ ((unused));
Packit Service 82fcde
Packit Service 82fcde
  /* Botch: Should this be handled at all?  Rely on callers?  */
Packit Service 82fcde
  if (dividend_size == 0)
Packit Service 82fcde
    return 0;
Packit Service 82fcde
Packit Service 82fcde
  /* If multiplication is much faster than division, and the
Packit Service 82fcde
     dividend is large, pre-invert the divisor, and use
Packit Service 82fcde
     only multiplications in the inner loop.  */
Packit Service 82fcde
Packit Service 82fcde
  /* This test should be read:
Packit Service 82fcde
       Does it ever help to use udiv_qrnnd_preinv?
Packit Service 82fcde
	 && Does what we save compensate for the inversion overhead?  */
Packit Service 82fcde
  if (UDIV_TIME > (2 * UMUL_TIME + 6)
Packit Service 82fcde
      && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME)
Packit Service 82fcde
    {
Packit Service 82fcde
      int normalization_steps;
Packit Service 82fcde
Packit Service 82fcde
      count_leading_zeros (normalization_steps, divisor_limb);
Packit Service 82fcde
      if (normalization_steps != 0)
Packit Service 82fcde
	{
Packit Service 82fcde
	  mp_limb_t divisor_limb_inverted;
Packit Service 82fcde
Packit Service 82fcde
	  divisor_limb <<= normalization_steps;
Packit Service 82fcde
Packit Service 82fcde
	  /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB.  The
Packit Service 82fcde
	     result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
Packit Service 82fcde
	     most significant bit (with weight 2**N) implicit.  */
Packit Service 82fcde
Packit Service 82fcde
	  /* Special case for DIVISOR_LIMB == 100...000.  */
Packit Service 82fcde
	  if (divisor_limb << 1 == 0)
Packit Service 82fcde
	    divisor_limb_inverted = ~(mp_limb_t) 0;
Packit Service 82fcde
	  else
Packit Service 82fcde
	    udiv_qrnnd (divisor_limb_inverted, dummy,
Packit Service 82fcde
			-divisor_limb, 0, divisor_limb);
Packit Service 82fcde
Packit Service 82fcde
	  n1 = dividend_ptr[dividend_size - 1];
Packit Service 82fcde
	  r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
Packit Service 82fcde
Packit Service 82fcde
	  /* Possible optimization:
Packit Service 82fcde
	     if (r == 0
Packit Service 82fcde
	     && divisor_limb > ((n1 << normalization_steps)
Packit Service 82fcde
			     | (dividend_ptr[dividend_size - 2] >> ...)))
Packit Service 82fcde
	     ...one division less... */
Packit Service 82fcde
Packit Service 82fcde
	  for (i = dividend_size - 2; i >= 0; i--)
Packit Service 82fcde
	    {
Packit Service 82fcde
	      n0 = dividend_ptr[i];
Packit Service 82fcde
	      udiv_qrnnd_preinv (dummy, r, r,
Packit Service 82fcde
				 ((n1 << normalization_steps)
Packit Service 82fcde
				  | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
Packit Service 82fcde
				 divisor_limb, divisor_limb_inverted);
Packit Service 82fcde
	      n1 = n0;
Packit Service 82fcde
	    }
Packit Service 82fcde
	  udiv_qrnnd_preinv (dummy, r, r,
Packit Service 82fcde
			     n1 << normalization_steps,
Packit Service 82fcde
			     divisor_limb, divisor_limb_inverted);
Packit Service 82fcde
	  return r >> normalization_steps;
Packit Service 82fcde
	}
Packit Service 82fcde
      else
Packit Service 82fcde
	{
Packit Service 82fcde
	  mp_limb_t divisor_limb_inverted;
Packit Service 82fcde
Packit Service 82fcde
	  /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB.  The
Packit Service 82fcde
	     result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
Packit Service 82fcde
	     most significant bit (with weight 2**N) implicit.  */
Packit Service 82fcde
Packit Service 82fcde
	  /* Special case for DIVISOR_LIMB == 100...000.  */
Packit Service 82fcde
	  if (divisor_limb << 1 == 0)
Packit Service 82fcde
	    divisor_limb_inverted = ~(mp_limb_t) 0;
Packit Service 82fcde
	  else
Packit Service 82fcde
	    udiv_qrnnd (divisor_limb_inverted, dummy,
Packit Service 82fcde
			-divisor_limb, 0, divisor_limb);
Packit Service 82fcde
Packit Service 82fcde
	  i = dividend_size - 1;
Packit Service 82fcde
	  r = dividend_ptr[i];
Packit Service 82fcde
Packit Service 82fcde
	  if (r >= divisor_limb)
Packit Service 82fcde
	    r = 0;
Packit Service 82fcde
	  else
Packit Service 82fcde
	    i--;
Packit Service 82fcde
Packit Service 82fcde
	  for (; i >= 0; i--)
Packit Service 82fcde
	    {
Packit Service 82fcde
	      n0 = dividend_ptr[i];
Packit Service 82fcde
	      udiv_qrnnd_preinv (dummy, r, r,
Packit Service 82fcde
				 n0, divisor_limb, divisor_limb_inverted);
Packit Service 82fcde
	    }
Packit Service 82fcde
	  return r;
Packit Service 82fcde
	}
Packit Service 82fcde
    }
Packit Service 82fcde
  else
Packit Service 82fcde
    {
Packit Service 82fcde
      if (UDIV_NEEDS_NORMALIZATION)
Packit Service 82fcde
	{
Packit Service 82fcde
	  int normalization_steps;
Packit Service 82fcde
Packit Service 82fcde
	  count_leading_zeros (normalization_steps, divisor_limb);
Packit Service 82fcde
	  if (normalization_steps != 0)
Packit Service 82fcde
	    {
Packit Service 82fcde
	      divisor_limb <<= normalization_steps;
Packit Service 82fcde
Packit Service 82fcde
	      n1 = dividend_ptr[dividend_size - 1];
Packit Service 82fcde
	      r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
Packit Service 82fcde
Packit Service 82fcde
	      /* Possible optimization:
Packit Service 82fcde
		 if (r == 0
Packit Service 82fcde
		 && divisor_limb > ((n1 << normalization_steps)
Packit Service 82fcde
				 | (dividend_ptr[dividend_size - 2] >> ...)))
Packit Service 82fcde
		 ...one division less... */
Packit Service 82fcde
Packit Service 82fcde
	      for (i = dividend_size - 2; i >= 0; i--)
Packit Service 82fcde
		{
Packit Service 82fcde
		  n0 = dividend_ptr[i];
Packit Service 82fcde
		  udiv_qrnnd (dummy, r, r,
Packit Service 82fcde
			      ((n1 << normalization_steps)
Packit Service 82fcde
			       | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
Packit Service 82fcde
			      divisor_limb);
Packit Service 82fcde
		  n1 = n0;
Packit Service 82fcde
		}
Packit Service 82fcde
	      udiv_qrnnd (dummy, r, r,
Packit Service 82fcde
			  n1 << normalization_steps,
Packit Service 82fcde
			  divisor_limb);
Packit Service 82fcde
	      return r >> normalization_steps;
Packit Service 82fcde
	    }
Packit Service 82fcde
	}
Packit Service 82fcde
      /* No normalization needed, either because udiv_qrnnd doesn't require
Packit Service 82fcde
	 it, or because DIVISOR_LIMB is already normalized.  */
Packit Service 82fcde
Packit Service 82fcde
      i = dividend_size - 1;
Packit Service 82fcde
      r = dividend_ptr[i];
Packit Service 82fcde
Packit Service 82fcde
      if (r >= divisor_limb)
Packit Service 82fcde
	r = 0;
Packit Service 82fcde
      else
Packit Service 82fcde
	i--;
Packit Service 82fcde
Packit Service 82fcde
      for (; i >= 0; i--)
Packit Service 82fcde
	{
Packit Service 82fcde
	  n0 = dividend_ptr[i];
Packit Service 82fcde
	  udiv_qrnnd (dummy, r, r, n0, divisor_limb);
Packit Service 82fcde
	}
Packit Service 82fcde
      return r;
Packit Service 82fcde
    }
Packit Service 82fcde
}