source/libs/gmp/gmp-src/mpz/bin_ui.c

   1 /* mpz_bin_ui - compute n over k.
   2
   3 Copyright 1998-2002, 2012, 2013 Free Software Foundation, Inc.
   4
   5 This file is part of the GNU MP Library.
   6
   7 The GNU MP Library is free software; you can redistribute it and/or modify
   8 it under the terms of either:
   9
  10   * the GNU Lesser General Public License as published by the Free
  11     Software Foundation; either version 3 of the License, or (at your
  12     option) any later version.
  13
  14 or
  15
  16   * the GNU General Public License as published by the Free Software
  17     Foundation; either version 2 of the License, or (at your option) any
  18     later version.
  19
  20 or both in parallel, as here.
  21
  22 The GNU MP Library is distributed in the hope that it will be useful, but
  23 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  24 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  25 for more details.
  26
  27 You should have received copies of the GNU General Public License and the
  28 GNU Lesser General Public License along with the GNU MP Library.  If not,
  29 see https://www.gnu.org/licenses/.  */
  30
  31 #include "gmp.h"
  32 #include "gmp-impl.h"
  33 #include "longlong.h"
  34
  35
  36 /* This is a poor implementation.  Look at bin_uiui.c for improvement ideas.
  37    In fact consider calling mpz_bin_uiui() when the arguments fit, leaving
  38    the code here only for big n.
  39
  40    The identity bin(n,k) = (-1)^k * bin(-n+k-1,k) can be found in Knuth vol
  41    1 section 1.2.6 part G. */
  42
  43
  44 #define DIVIDE()                                                              \
  45   do {                                                                        \
  46     ASSERT (SIZ(r) > 0);                                                      \
  47     MPN_DIVREM_OR_DIVEXACT_1 (PTR(r), PTR(r), (mp_size_t) SIZ(r), kacc);      \
  48     SIZ(r) -= (PTR(r)[SIZ(r)-1] == 0);                                        \
  49   } while (0)
  50
  51 void
  52 mpz_bin_ui (mpz_ptr r, mpz_srcptr n, unsigned long int k)
  53 {
  54   mpz_t      ni;
  55   mp_limb_t  i;
  56   mpz_t      nacc;
  57   mp_limb_t  kacc;
  58   mp_size_t  negate;
  59
  60   if (SIZ (n) < 0)
  61     {
  62       /* bin(n,k) = (-1)^k * bin(-n+k-1,k), and set ni = -n+k-1 - k = -n-1 */
  63       mpz_init (ni);
  64       mpz_add_ui (ni, n, 1L);
  65       mpz_neg (ni, ni);
  66       negate = (k & 1);   /* (-1)^k */
  67     }
  68   else
  69     {
  70       /* bin(n,k) == 0 if k>n
  71          (no test for this under the n<0 case, since -n+k-1 >= k there) */
  72       if (mpz_cmp_ui (n, k) < 0)
  73         {
  74           SIZ (r) = 0;
  75           return;
  76         }
  77
  78       /* set ni = n-k */
  79       mpz_init (ni);
  80       mpz_sub_ui (ni, n, k);
  81       negate = 0;
  82     }
  83
  84   /* Now wanting bin(ni+k,k), with ni positive, and "negate" is the sign (0
  85      for positive, 1 for negative). */
  86   SIZ (r) = 1; PTR (r)[0] = 1;
  87
  88   /* Rewrite bin(n,k) as bin(n,n-k) if that is smaller.  In this case it's
  89      whether ni+k-k < k meaning ni<k, and if so change to denominator ni+k-k
  90      = ni, and new ni of ni+k-ni = k.  */
  91   if (mpz_cmp_ui (ni, k) < 0)
  92     {
  93       unsigned long  tmp;
  94       tmp = k;
  95       k = mpz_get_ui (ni);
  96       mpz_set_ui (ni, tmp);
  97     }
  98
  99   kacc = 1;
 100   mpz_init_set_ui (nacc, 1L);
 101
 102   for (i = 1; i <= k; i++)
 103     {
 104       mp_limb_t k1, k0;
 105
 106 #if 0
 107       mp_limb_t nacclow;
 108       int c;
 109
 110       nacclow = PTR(nacc)[0];
 111       for (c = 0; (((kacc | nacclow) & 1) == 0); c++)
 112         {
 113           kacc >>= 1;
 114           nacclow >>= 1;
 115         }
 116       mpz_div_2exp (nacc, nacc, c);
 117 #endif
 118
 119       mpz_add_ui (ni, ni, 1L);
 120       mpz_mul (nacc, nacc, ni);
 121       umul_ppmm (k1, k0, kacc, i << GMP_NAIL_BITS);
 122       if (k1 != 0)
 123         {
 124           /* Accumulator overflow.  Perform bignum step.  */
 125           mpz_mul (r, r, nacc);
 126           SIZ (nacc) = 1; PTR (nacc)[0] = 1;
 127           DIVIDE ();
 128           kacc = i;
 129         }
 130       else
 131         {
 132           /* Save new products in accumulators to keep accumulating.  */
 133           kacc = k0 >> GMP_NAIL_BITS;
 134         }
 135     }
 136
 137   mpz_mul (r, r, nacc);
 138   DIVIDE ();
 139   SIZ(r) = (SIZ(r) ^ -negate) + negate;
 140
 141   mpz_clear (nacc);
 142   mpz_clear (ni);
 143 }