rts/gmp/mpz/bin_ui.c

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