third_party/gmp/mpz/bin_ui.c

/* mpz_bin_ui(RESULT, N, K) -- Set RESULT to N over K.

Copyright 1998-2002, 2012, 2013, 2015, 2017-2018 Free Software Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of either:

  * the GNU Lesser General Public License as published by the Free
    Software Foundation; either version 3 of the License, or (at your
    option) any later version.

or

  * the GNU General Public License as published by the Free Software
    Foundation; either version 2 of the License, or (at your option) any
    later version.

or both in parallel, as here.

The GNU MP Library 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 copies of the GNU General Public License and the
GNU Lesser General Public License along with the GNU MP Library.  If not,
see https://www.gnu.org/licenses/.  */

#include "gmp-impl.h"

/* How many special cases? Minimum is 2: 0 and 1;
 * also 3 {0,1,2} and 5 {0,1,2,3,4} are implemented.
 */
#define APARTAJ_KALKULOJ 2

/* Whether to use (1) or not (0) the function mpz_bin_uiui whenever
 * the operands fit.
 */
#define UZU_BIN_UIUI 0

/* Whether to use a shortcut to precompute the product of four
 * elements (1), or precompute only the product of a couple (0).
 *
 * In both cases the precomputed product is then updated with some
 * linear operations to obtain the product of the next four (1)
 * [or two (0)] operands.
 */
#define KVAROPE 1

static void
posmpz_init (mpz_ptr r)
{
  mp_ptr rp;
  ASSERT (SIZ (r) > 0);
  rp = SIZ (r) + MPZ_REALLOC (r, SIZ (r) + 2);
  *rp = 0;
  *++rp = 0;
}

/* Equivalent to mpz_add_ui (r, r, in), but faster when
   0 < SIZ (r) < ALLOC (r) and limbs above SIZ (r) contain 0. */
static void
posmpz_inc_ui (mpz_ptr r, unsigned long in)
{
#if BITS_PER_ULONG > GMP_NUMB_BITS
  mpz_add_ui (r, r, in);
#else
  ASSERT (SIZ (r) > 0);
  MPN_INCR_U (PTR (r), SIZ (r) + 1, in);
  SIZ (r) += (PTR (r)[SIZ (r)] != 0);
#endif
}

/* Equivalent to mpz_sub_ui (r, r, in), but faster when
   0 < SIZ (r) and we know in advance that the result is positive. */
static void
posmpz_dec_ui (mpz_ptr r, unsigned long in)
{
#if BITS_PER_ULONG > GMP_NUMB_BITS
  mpz_sub_ui (r, r, in);
#else
  ASSERT (mpz_cmp_ui (r, in) >= 0);
  MPN_DECR_U (PTR (r), SIZ (r), in);
  SIZ (r) -= (PTR (r)[SIZ (r)-1] == 0);
#endif
}

/* Equivalent to mpz_tdiv_q_2exp (r, r, 1), but faster when
   0 < SIZ (r) and we know in advance that the result is positive. */
static void
posmpz_rsh1 (mpz_ptr r)
{
  mp_ptr rp;
  mp_size_t rn;

  rn = SIZ (r);
  rp = PTR (r);
  ASSERT (rn > 0);
  mpn_rshift (rp, rp, rn, 1);
  SIZ (r) -= rp[rn - 1] == 0;
}

/* Computes r = n(n+(2*k-1))/2
   It uses a sqare instead of a product, computing
   r = ((n+k-1)^2 + n - (k-1)^2)/2
   As a side effect, sets t = n+k-1
 */
static void
mpz_hmul_nbnpk (mpz_ptr r, mpz_srcptr n, unsigned long int k, mpz_ptr t)
{
  ASSERT (k > 0 && SIZ(n) > 0);
  --k;
  mpz_add_ui (t, n, k);
  mpz_mul (r, t, t);
  mpz_add (r, r, n);
  posmpz_rsh1 (r);
  if (LIKELY (k <= (1UL << (BITS_PER_ULONG / 2))))
    posmpz_dec_ui (r, (k + (k & 1))*(k >> 1));
  else
    {
      mpz_t tmp;
      mpz_init_set_ui (tmp, (k + (k & 1)));
      mpz_mul_ui (tmp, tmp, k >> 1);
      mpz_sub (r, r, tmp);
      mpz_clear (tmp);
    }
}

#if KVAROPE
static void
rek_raising_fac4 (mpz_ptr r, mpz_ptr p, mpz_ptr P, unsigned long int k, unsigned long int lk, mpz_ptr t)
{
  if (k - lk < 5)
    {
      do {
	posmpz_inc_ui (p, 4*k+2);
	mpz_addmul_ui (P, p, 4*k);
	posmpz_dec_ui (P, k);
	mpz_mul (r, r, P);
      } while (--k > lk);
    }
  else
    {
      mpz_t lt;
      unsigned long int m;

      m = ((k + lk) >> 1) + 1;
      rek_raising_fac4 (r, p, P, k, m, t);

      posmpz_inc_ui (p, 4*m+2);
      mpz_addmul_ui (P, p, 4*m);
      posmpz_dec_ui (P, m);
      if (t == NULL)
	{
	  mpz_init_set (lt, P);
	  t = lt;
	}
      else
	{
	  ALLOC (lt) = 0;
	  mpz_set (t, P);
	}
      rek_raising_fac4 (t, p, P, m - 1, lk, NULL);

      mpz_mul (r, r, t);
      mpz_clear (lt);
    }
}

/* Computes (n+1)(n+2)...(n+k)/2^(k/2 +k/4) using the helper function
   rek_raising_fac4, and exploiting an idea inspired by a piece of
   code that Fredrik Johansson wrote and by a comment by Niels Möller.

   Assume k = 4i then compute:
     p  = (n+1)(n+4i)/2 - i
	  (n+1+1)(n+4i)/2 = p + i + (n+4i)/2
	  (n+1+1)(n+4i-1)/2 = p + i + ((n+4i)-(n+1+1))/2 = p + i + (n-n+4i-2)/2 = p + 3i-1
     P  = (p + i)*(p+3i-1)/2 = (n+1)(n+2)(n+4i-1)(n+4i)/8
     n' = n + 2
     i' = i - 1
	  (n'-1)(n')(n'+4i'+1)(n'+4i'+2)/8 = P
	  (n'-1)(n'+4i'+2)/2 - i' - 1 = p
	  (n'-1+2)(n'+4i'+2)/2 - i' - 1 = p + (n'+4i'+2)
	  (n'-1+2)(n'+4i'+2-2)/2 - i' - 1 = p + (n'+4i'+2) - (n'-1+2) =  p + 4i' + 1
	  (n'-1+2)(n'+4i'+2-2)/2 - i' = p + 4i' + 2
     p' = p + 4i' + 2 = (n'+1)(n'+4i')/2 - i'
	  p' - 4i' - 2 = p
	  (p' - 4i' - 2 + i)*(p' - 4i' - 2+3i-1)/2 = P
	  (p' - 4i' - 2 + i' + 1)*(p' - 4i' - 2 + 3i' + 3 - 1)/2 = P
	  (p' - 3i' - 1)*(p' - i')/2 = P
	  (p' - 3i' - 1 + 4i' + 1)*(p' - i' + 4i' - 1)/2 = P + (4i' + 1)*(p' - i')/2 + (p' - 3i' - 1 + 4i' + 1)*(4i' - 1)/2
	  (p' + i')*(p' + 3i' - 1)/2 = P + (4i')*(p' + p')/2 + (p' - i' - (p' + i'))/2
	  (p' + i')*(p' + 3i' - 1)/2 = P + 4i'p' + (p' - i' - p' - i')/2
	  (p' + i')*(p' + 3i' - 1)/2 = P + 4i'p' - i'
     P' = P + 4i'p' - i'

   And compute the product P * P' * P" ...
 */

static void
mpz_raising_fac4 (mpz_ptr r, mpz_ptr n, unsigned long int k, mpz_ptr t, mpz_ptr p)
{
  ASSERT ((k >= APARTAJ_KALKULOJ) && (APARTAJ_KALKULOJ > 0));
  posmpz_init (n);
  posmpz_inc_ui (n, 1);
  SIZ (r) = 0;
  if (k & 1)
    {
      mpz_set (r, n);
      posmpz_inc_ui (n, 1);
    }
  k >>= 1;
  if (APARTAJ_KALKULOJ < 2 && k == 0)
    return;

  mpz_hmul_nbnpk (p, n, k, t);
  posmpz_init (p);

  if (k & 1)
    {
      if (SIZ (r))
	mpz_mul (r, r, p);
      else
	mpz_set (r, p);
      posmpz_inc_ui (p, k - 1);
    }
  k >>= 1;
  if (APARTAJ_KALKULOJ < 4 && k == 0)
    return;

  mpz_hmul_nbnpk (t, p, k, n);
  if (SIZ (r))
    mpz_mul (r, r, t);
  else
    mpz_set (r, t);

  if (APARTAJ_KALKULOJ > 8 || k > 1)
    {
      posmpz_dec_ui (p, k);
      rek_raising_fac4 (r, p, t, k - 1, 0, n);
    }
}

#else /* KVAROPE */

static void
rek_raising_fac (mpz_ptr r, mpz_ptr n, unsigned long int k, unsigned long int lk, mpz_ptr t1, mpz_ptr t2)
{
  /* Should the threshold depend on SIZ (n) ? */
  if (k - lk < 10)
    {
      do {
	posmpz_inc_ui (n, k);
	mpz_mul (r, r, n);
	--k;
      } while (k > lk);
    }
  else
    {
      mpz_t t3;
      unsigned long int m;

      m = ((k + lk) >> 1) + 1;
      rek_raising_fac (r, n, k, m, t1, t2);

      posmpz_inc_ui (n, m);
      if (t1 == NULL)
	{
	  mpz_init_set (t3, n);
	  t1 = t3;
	}
      else
	{
	  ALLOC (t3) = 0;
	  mpz_set (t1, n);
	}
      rek_raising_fac (t1, n, m - 1, lk, t2, NULL);

      mpz_mul (r, r, t1);
      mpz_clear (t3);
    }
}

/* Computes (n+1)(n+2)...(n+k)/2^(k/2) using the helper function
   rek_raising_fac, and exploiting an idea inspired by a piece of
   code that Fredrik Johansson wrote.

   Force an even k = 2i then compute:
     p  = (n+1)(n+2i)/2
     i' = i - 1
     p == (n+1)(n+2i'+2)/2
     p' = p + i' == (n+2)(n+2i'+1)/2
     n' = n + 1
     p'== (n'+1)(n'+2i')/2 == (n+1 +1)(n+2i -1)/2

   And compute the product p * p' * p" ...
*/

static void
mpz_raising_fac (mpz_ptr r, mpz_ptr n, unsigned long int k, mpz_ptr t, mpz_ptr p)
{
  unsigned long int hk;
  ASSERT ((k >= APARTAJ_KALKULOJ) && (APARTAJ_KALKULOJ > 1));
  mpz_add_ui (n, n, 1);
  hk = k >> 1;
  mpz_hmul_nbnpk (p, n, hk, t);

  if ((k & 1) != 0)
    {
      mpz_add_ui (t, t, hk + 1);
      mpz_mul (r, t, p);
    }
  else
    {
      mpz_set (r, p);
    }

  if ((APARTAJ_KALKULOJ > 3) || (hk > 1))
    {
      posmpz_init (p);
      rek_raising_fac (r, p, hk - 1, 0, t, n);
    }
}
#endif /* KVAROPE */

/* This is a poor implementation.  Look at bin_uiui.c for improvement ideas.
   In fact consider calling mpz_bin_uiui() when the arguments fit, leaving
   the code here only for big n.

   The identity bin(n,k) = (-1)^k * bin(-n+k-1,k) can be found in Knuth vol
   1 section 1.2.6 part G. */

void
mpz_bin_ui (mpz_ptr r, mpz_srcptr n, unsigned long int k)
{
  mpz_t      ni;
  mp_size_t  negate;

  if (SIZ (n) < 0)
    {
      /* bin(n,k) = (-1)^k * bin(-n+k-1,k), and set ni = -n+k-1 - k = -n-1 */
      mpz_init (ni);
      mpz_add_ui (ni, n, 1L);
      mpz_neg (ni, ni);
      negate = (k & 1);   /* (-1)^k */
    }
  else
    {
      /* bin(n,k) == 0 if k>n
	 (no test for this under the n<0 case, since -n+k-1 >= k there) */
      if (mpz_cmp_ui (n, k) < 0)
	{
	  SIZ (r) = 0;
	  return;
	}

      /* set ni = n-k */
      mpz_init (ni);
      mpz_sub_ui (ni, n, k);
      negate = 0;
    }

  /* Now wanting bin(ni+k,k), with ni positive, and "negate" is the sign (0
     for positive, 1 for negative). */

  /* Rewrite bin(n,k) as bin(n,n-k) if that is smaller.  In this case it's
     whether ni+k-k < k meaning ni<k, and if so change to denominator ni+k-k
     = ni, and new ni of ni+k-ni = k.  */
  if (mpz_cmp_ui (ni, k) < 0)
    {
      unsigned long  tmp;
      tmp = k;
      k = mpz_get_ui (ni);
      mpz_set_ui (ni, tmp);
    }

  if (k < APARTAJ_KALKULOJ)
    {
      if (k == 0)
	{
	  SIZ (r) = 1;
	  MPZ_NEWALLOC (r, 1)[0] = 1;
	}
#if APARTAJ_KALKULOJ > 2
      else if (k == 2)
	{
	  mpz_add_ui (ni, ni, 1);
	  mpz_mul (r, ni, ni);
	  mpz_add (r, r, ni);
	  posmpz_rsh1 (r);
	}
#endif
#if APARTAJ_KALKULOJ > 3
      else if (k > 2)
	{ /* k = 3, 4 */
	  mpz_add_ui (ni, ni, 2); /* n+1 */
	  mpz_mul (r, ni, ni); /* (n+1)^2 */
	  mpz_sub_ui (r, r, 1); /* (n+1)^2-1 */
	  if (k == 3)
	    {
	      mpz_mul (r, r, ni); /* ((n+1)^2-1)(n+1) = n(n+1)(n+2) */
	      /* mpz_divexact_ui (r, r, 6); /\* 6=3<<1; div_by3 ? *\/ */
	      mpn_pi1_bdiv_q_1 (PTR(r), PTR(r), SIZ(r), 3, GMP_NUMB_MASK/3*2+1, 1);
	      MPN_NORMALIZE_NOT_ZERO (PTR(r), SIZ(r));
	    }
	  else /* k = 4 */
	    {
	      mpz_add (ni, ni, r); /* (n+1)^2+n */
	      mpz_mul (r, ni, ni); /* ((n+1)^2+n)^2 */
	      mpz_sub_ui (r, r, 1); /* ((n+1)^2+n)^2-1 = n(n+1)(n+2)(n+3) */
	      /* mpz_divexact_ui (r, r, 24); /\* 24=3<<3; div_by3 ? *\/ */
	      mpn_pi1_bdiv_q_1 (PTR(r), PTR(r), SIZ(r), 3, GMP_NUMB_MASK/3*2+1, 3);
	      MPN_NORMALIZE_NOT_ZERO (PTR(r), SIZ(r));
	    }
	}
#endif
      else
	{ /* k = 1 */
	  mpz_add_ui (r, ni, 1);
	}
    }
#if UZU_BIN_UIUI
  else if (mpz_cmp_ui (ni, ULONG_MAX - k) <= 0)
    {
      mpz_bin_uiui (r, mpz_get_ui (ni) + k, k);
    }
#endif
  else
    {
      mp_limb_t count;
      mpz_t num, den;

      mpz_init (num);
      mpz_init (den);

#if KVAROPE
      mpz_raising_fac4 (num, ni, k, den, r);
      popc_limb (count, k);
      ASSERT (k - (k >> 1) - (k >> 2) - count >= 0);
      mpz_tdiv_q_2exp (num, num, k - (k >> 1) - (k >> 2) - count);
#else
      mpz_raising_fac (num, ni, k, den, r);
      popc_limb (count, k);
      ASSERT (k - (k >> 1) - count >= 0);
      mpz_tdiv_q_2exp (num, num, k - (k >> 1) - count);
#endif

      mpz_oddfac_1(den, k, 0);

      mpz_divexact(r, num, den);
      mpz_clear (num);
      mpz_clear (den);
    }
  mpz_clear (ni);

  SIZ(r) = (SIZ(r) ^ -negate) + negate;
}