Imported gcc-4.4.3

[msp430-gcc.git] / gmp / mpz / bin_ui.c

diff --git a/gmp/mpz/bin_ui.c b/gmp/mpz/bin_ui.c
new file mode 100644 (file)
index 0000000..c97ed7c
--- /dev/null
+++ b/gmp/mpz/bin_ui.c
@@ -0,0 +1,133 @@
+/* mpz_bin_ui - compute n over k.
+
+Copyright 1998, 1999, 2000, 2001, 2002 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"
+
+
+/* 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. */
+
+
+#define DIVIDE()                                                              \
+  do {                                                                        \
+    ASSERT (SIZ(r) > 0);                                                      \
+    MPN_DIVREM_OR_DIVEXACT_1 (PTR(r), PTR(r), (mp_size_t) SIZ(r), kacc);      \
+    SIZ(r) -= (PTR(r)[SIZ(r)-1] == 0);                                        \
+  } while (0)
+
+void
+mpz_bin_ui (mpz_ptr r, mpz_srcptr n, unsigned long int k)
+{
+  mpz_t      ni;
+  mp_limb_t  i;
+  mpz_t      nacc;
+  mp_limb_t  kacc;
+  mp_size_t  negate;
+
+  if (mpz_sgn (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_neg (ni, n);
+      mpz_sub_ui (ni, ni, 1L);
+      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)
+        {
+          mpz_set_ui (r, 0L);
+          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). */
+  mpz_set_ui (r, 1L);
+
+  /* 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);
+    }
+
+  kacc = 1;
+  mpz_init_set_ui (nacc, 1L);
+
+  for (i = 1; i <= k; i++)
+    {
+      mp_limb_t k1, k0;
+
+#if 0
+      mp_limb_t nacclow;
+      int c;
+
+      nacclow = PTR(nacc)[0];
+      for (c = 0; (((kacc | nacclow) & 1) == 0); c++)
+       {
+         kacc >>= 1;
+         nacclow >>= 1;
+       }
+      mpz_div_2exp (nacc, nacc, c);
+#endif
+
+      mpz_add_ui (ni, ni, 1L);
+      mpz_mul (nacc, nacc, ni);
+      umul_ppmm (k1, k0, kacc, i << GMP_NAIL_BITS);
+      k0 >>= GMP_NAIL_BITS;
+      if (k1 != 0)
+       {
+         /* Accumulator overflow.  Perform bignum step.  */
+         mpz_mul (r, r, nacc);
+         mpz_set_ui (nacc, 1L);
+          DIVIDE ();
+         kacc = i;
+       }
+      else
+       {
+         /* Save new products in accumulators to keep accumulating.  */
+         kacc = k0;
+       }
+    }
+
+  mpz_mul (r, r, nacc);
+  DIVIDE ();
+  SIZ(r) = (SIZ(r) ^ -negate) + negate;
+
+  mpz_clear (nacc);
+  mpz_clear (ni);
+}