/* 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;

}