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     ! 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: }