ghc/rts/gmp/mpz/bin_uiui.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 /* Avoid reallocs by rounding up any new size */
  28 #define ROUNDUP_MASK  15
  29
  30 /* Enhancement: use mpn_divexact_1 when it exists */
  31 #define MULDIV()                                                \
  32   MPZ_REALLOC (r, (SIZ(r)+1)|ROUNDUP_MASK);                     \
  33   PTR(r)[SIZ(r)] = mpn_mul_1 (PTR(r), PTR(r), SIZ(r), nacc);    \
  34   ASSERT_NOCARRY (mpn_divrem_1 (PTR(r), (mp_size_t) 0,          \
  35                                 PTR(r), SIZ(r)+1, kacc));       \
  36   SIZ(r) += (PTR(r)[SIZ(r)] != 0);
  37
  38 void
  39 #if __STDC__
  40 mpz_bin_uiui (mpz_ptr r, unsigned long int n, unsigned long int k)
  41 #else
  42 mpz_bin_uiui (r, n, k)
  43      mpz_ptr r;
  44      unsigned long int n;
  45      unsigned long int k;
  46 #endif
  47 {
  48   unsigned long int  i, j;
  49   mp_limb_t          nacc, kacc;
  50   unsigned long int  cnt;
  51
  52   /* bin(n,k) = 0 if k>n. */
  53   if (n < k)
  54     {
  55       mpz_set_ui (r, 0);
  56       return;
  57     }
  58
  59   /* Rewrite bin(n,k) as bin(n,n-k) if that is smaller. */
  60   k = MIN (k, n-k);
  61
  62   /* bin(n,0) = 1 */
  63   if (k == 0)
  64     {
  65       mpz_set_ui (r, 1);
  66       return;
  67     }
  68
  69   j = n - k + 1;
  70   mpz_set_ui (r, j);
  71
  72   /* Initialize accumulators.  */
  73   nacc = 1;
  74   kacc = 1;
  75
  76   cnt = 0;
  77   for (i = 2; i <= k; i++)
  78     {
  79       mp_limb_t n1, n0, k1, k0;
  80
  81       j++;
  82 #if 0
  83       /* Remove common multiples of 2.  This will allow us to accumulate
  84          more in nacc and kacc before we need a bignum step.  It would make
  85          sense to cancel factors of 3, 5, etc too, but this would be best
  86          handled by sieving out factors.  Alternatively, we could perform a
  87          gcd of the accumulators just as they have overflown, and keep
  88          accumulating until the gcd doesn't remove a significant factor.  */
  89       while (((nacc | kacc) & 1) == 0)
  90         {
  91           nacc >>= 1;
  92           kacc >>= 1;
  93         }
  94 #else
  95       cnt = ((nacc | kacc) & 1) ^ 1;
  96       nacc >>= cnt;
  97       kacc >>= cnt;
  98 #endif
  99       /* Accumulate next multiples.  */
 100       umul_ppmm (n1, n0, nacc, j);
 101       umul_ppmm (k1, k0, kacc, i);
 102       if (n1 != 0)
 103         {
 104           /* Accumulator overflow.  Perform bignum step.  */
 105           MULDIV ();
 106           nacc = j;
 107           kacc = i;
 108         }
 109       else
 110         {
 111           if (k1 != 0) abort ();
 112           /* Save new products in accumulators to keep accumulating.  */
 113           nacc = n0;
 114           kacc = k0;
 115         }
 116     }
 117
 118   /* Take care of whatever is left in accumulators.  */
 119   MULDIV ();
 120 }