Return to bin_ui.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpz

Annotation of OpenXM_contrib/gmp/mpz/bin_ui.c, Revision 1.1.1.1

1.1       maekawa     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: }