/* mpz_bc_bin_uiui - compute n over k.

Copyright 2011 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 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.

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 Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */

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

/* Algorithm:

   Plain and simply multiply things together.

   We tabulate factorials (k!/2^t)^(-1) mod B for 1= < k <= 32 (where t is
   chosen such that k!/2^t is odd).

   Improvement:

   Shift i, not c, in the loop.  Requires that we compute a starting pos in ctz
   using the low bits of n-k+1.
*/

static mp_limb_t facinv[] = {
#if GMP_NUMB_BITS == 64
  0x0000000000000001,  /*  1 */  0x0000000000000001,  /*  2 */
  0xaaaaaaaaaaaaaaab,  /*  3 */  0xaaaaaaaaaaaaaaab,  /*  4 */
  0xeeeeeeeeeeeeeeef,  /*  5 */  0x4fa4fa4fa4fa4fa5,  /*  6 */
  0x2ff2ff2ff2ff2ff3,  /*  7 */  0x2ff2ff2ff2ff2ff3,  /*  8 */
  0x938cc70553e3771b,  /*  9 */  0xb71c27cddd93e49f,  /* 10 */
  0xb38e3229fcdee63d,  /* 11 */  0xe684bb63544a4cbf,  /* 12 */
  0xc2f684917ca340fb,  /* 13 */  0xf747c9cba417526d,  /* 14 */
  0xbb26eb51d7bd49c3,  /* 15 */  0xbb26eb51d7bd49c3,  /* 16 */
  0xb0a7efb985294093,  /* 17 */  0xbe4b8c69f259eabb,  /* 18 */
  0x6854d17ed6dc4fb9,  /* 19 */  0xe1aa904c915f4325,  /* 20 */
  0x3b8206df131cead1,  /* 21 */  0x79c6009fea76fe13,  /* 22 */
  0xd8c5d381633cd365,  /* 23 */  0x4841f12b21144677,  /* 24 */
  0x4a91ff68200b0d0f,  /* 25 */  0x8f9513a58c4f9e8b,  /* 26 */
  0x2b3e690621a42251,  /* 27 */  0x4f520f00e03c04e7,  /* 28 */
  0x2edf84ee600211d3,  /* 29 */  0xadcaa2764aaacdfd,  /* 30 */
  0x161f4f9033f4fe63,  /* 31 */  0x161f4f9033f4fe63   /* 32 */
#endif
#if GMP_NUMB_BITS == 32
  0x00000001,  /*  1 */  0x00000001,  /*  2 */
  0xaaaaaaab,  /*  3 */  0xaaaaaaab,  /*  4 */
  0xeeeeeeef,  /*  5 */  0xa4fa4fa5,  /*  6 */
  0xf2ff2ff3,  /*  7 */  0xf2ff2ff3,  /*  8 */
  0x53e3771b,  /*  9 */  0xdd93e49f,  /* 10 */
  0xfcdee63d,  /* 11 */  0x544a4cbf,  /* 12 */
  0x7ca340fb,  /* 13 */  0xa417526d,  /* 14 */
  0xd7bd49c3,  /* 15 */  0xd7bd49c3   /* 16 */
#endif
};

static unsigned char ctz[] = {0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,5};

void
mpz_bc_bin_uiui (mpz_ptr r, unsigned long int n, unsigned long int k)
{
  mp_limb_t i, j, c;

  c = facinv[k - 1];
  i = n - k + 1;
  for (j = 0; j < k; j++)
    {
      c = c * i;
      c = c >> ctz[j];
      i++;
    }

  c = c & ((CNST_LIMB(1) << GMP_NUMB_BITS-4) - 1);

  PTR(r)[0] = c;
  SIZ(r) = 1;
}

#ifdef TEST
#include <stdio.h>
#include <stdlib.h>

int
main (int argc, char **argv)
{
  unsigned long n, k;
  mpz_t r, ref;

  mpz_inits (r, ref, 0);

  for (k = 1; k < GMP_NUMB_BITS / 2; k++)
    {
      for (n = k; n < GMP_LIMB_BITS; n++)
	{
	  mpz_bin_uiui (ref, n, k);
	  mpz_bc_bin_uiui (r, n, k);
	  if (mpz_cmp (ref, r) != 0)
	    {
	      printf ("\r%ld over %ld fails", n, k);
	      abort ();
	    }
	}
    }

  return 0;
}
#endif

#ifdef MAIN
#include <stdlib.h>

int
main (int argc, char **argv)
{
  unsigned long n, k;
  mpz_t r;

  n = strtoul (argv[1], 0, 0);
  k = strtoul (argv[2], 0, 0);

  mpz_init (r);

  mpz_bc_bin_uiui (r, n, k);
  gmp_printf ("%Zx\n", r);

  return 0;
}
#endif