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.2

1.1       maekawa     1: /* mpz_bin_uiui - compute n over k.
                      2:
1.1.1.2 ! ohara       3: Copyright 1998, 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
1.1       maekawa     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 */
1.1.1.2 ! ohara      36: #define DIVIDE()                                                \
        !            37:   do {                                                          \
        !            38:     ASSERT (SIZ(r) > 0);                                        \
        !            39:     MPN_DIVREM_OR_DIVEXACT_1 (PTR(r), PTR(r), SIZ(r), kacc);    \
        !            40:     SIZ(r) -= (PTR(r)[SIZ(r)-1] == 0);                          \
        !            41:   } while (0)
1.1       maekawa    42:
                     43: void
                     44: mpz_bin_ui (mpz_ptr r, mpz_srcptr n, unsigned long int k)
                     45: {
                     46:   mpz_t      ni;
                     47:   mp_limb_t  i;
                     48:   mpz_t      nacc;
                     49:   mp_limb_t  kacc;
                     50:   mp_size_t  negate;
1.1.1.2 ! ohara      51:
1.1       maekawa    52:   if (mpz_sgn (n) < 0)
                     53:     {
                     54:       /* bin(n,k) = (-1)^k * bin(-n+k-1,k), and set ni = -n+k-1 - k = -n-1 */
                     55:       mpz_init (ni);
                     56:       mpz_neg (ni, n);
                     57:       mpz_sub_ui (ni, ni, 1L);
                     58:       negate = (k & 1);   /* (-1)^k */
                     59:     }
                     60:   else
                     61:     {
                     62:       /* bin(n,k) == 0 if k>n
                     63:          (no test for this under the n<0 case, since -n+k-1 >= k there) */
                     64:       if (mpz_cmp_ui (n, k) < 0)
                     65:         {
                     66:           mpz_set_ui (r, 0L);
                     67:           return;
                     68:         }
                     69:
                     70:       /* set ni = n-k */
                     71:       mpz_init (ni);
                     72:       mpz_sub_ui (ni, n, k);
                     73:       negate = 0;
                     74:     }
                     75:
                     76:   /* Now wanting bin(ni+k,k), with ni positive, and "negate" is the sign (0
                     77:      for positive, 1 for negative). */
                     78:   mpz_set_ui (r, 1L);
                     79:
                     80:   /* Rewrite bin(n,k) as bin(n,n-k) if that is smaller.  In this case it's
                     81:      whether ni+k-k < k meaning ni<k, and if so change to denominator ni+k-k
                     82:      = ni, and new ni of ni+k-ni = k.  */
                     83:   if (mpz_cmp_ui (ni, k) < 0)
                     84:     {
                     85:       unsigned long  tmp;
                     86:       tmp = k;
                     87:       k = mpz_get_ui (ni);
                     88:       mpz_set_ui (ni, tmp);
                     89:     }
                     90:
                     91:   kacc = 1;
                     92:   mpz_init_set_ui (nacc, 1);
                     93:
                     94:   for (i = 1; i <= k; i++)
                     95:     {
                     96:       mp_limb_t k1, k0;
                     97:
                     98: #if 0
                     99:       mp_limb_t nacclow;
                    100:       int c;
                    101:
                    102:       nacclow = PTR(nacc)[0];
                    103:       for (c = 0; (((kacc | nacclow) & 1) == 0); c++)
                    104:        {
                    105:          kacc >>= 1;
                    106:          nacclow >>= 1;
                    107:        }
                    108:       mpz_div_2exp (nacc, nacc, c);
                    109: #endif
                    110:
                    111:       mpz_add_ui (ni, ni, 1);
                    112:       mpz_mul (nacc, nacc, ni);
1.1.1.2 ! ohara     113:       umul_ppmm (k1, k0, kacc, i << GMP_NAIL_BITS);
        !           114:       k0 >>= GMP_NAIL_BITS;
1.1       maekawa   115:       if (k1 != 0)
                    116:        {
                    117:          /* Accumulator overflow.  Perform bignum step.  */
                    118:          mpz_mul (r, r, nacc);
                    119:          mpz_set_ui (nacc, 1);
                    120:           DIVIDE ();
                    121:          kacc = i;
                    122:        }
                    123:       else
                    124:        {
                    125:          /* Save new products in accumulators to keep accumulating.  */
                    126:          kacc = k0;
                    127:        }
                    128:     }
                    129:
                    130:   mpz_mul (r, r, nacc);
                    131:   DIVIDE ();
                    132:   SIZ(r) = (SIZ(r) ^ -negate) + negate;
                    133:
                    134:   mpz_clear (nacc);
                    135:   mpz_clear (ni);
                    136: }