Return to bin_uiui.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpz

Annotation of OpenXM_contrib/gmp/mpz/bin_uiui.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:
1.1.1.2 ! ohara      27: /* Enhancement: It ought to be possible to calculate the size of the final
        !            28:    result in advance, to a rough approximation at least, and use it to do
        !            29:    just one realloc.  Stirling's approximation n! ~= sqrt(2*pi*n)*(n/e)^n
        !            30:    (Knuth section 1.2.5) might be of use.  */
        !            31:
        !            32: /* "inc" in the main loop allocates a chunk more space if not already
        !            33:    enough, so as to avoid repeated reallocs.  The final step on the other
        !            34:    hand requires only one more limb.  */
        !            35: #define MULDIV(inc)                                                     \
        !            36:   do {                                                                  \
        !            37:     ASSERT (rsize <= ralloc);                                           \
        !            38:                                                                         \
        !            39:     if (rsize == ralloc)                                                \
        !            40:       {                                                                 \
        !            41:         mp_size_t  new_ralloc = ralloc + (inc);                         \
        !            42:         rp = __GMP_REALLOCATE_FUNC_LIMBS (rp, ralloc, new_ralloc);      \
        !            43:         ralloc = new_ralloc;                                            \
        !            44:       }                                                                 \
        !            45:                                                                         \
        !            46:     rp[rsize] = mpn_mul_1 (rp, rp, rsize, nacc);                        \
        !            47:     MPN_DIVREM_OR_DIVEXACT_1 (rp, rp, rsize+1, kacc);                   \
        !            48:     rsize += (rp[rsize] != 0);                                          \
        !            49:                                                                         \
        !            50: } while (0)
1.1       maekawa    51:
                     52: void
                     53: mpz_bin_uiui (mpz_ptr r, unsigned long int n, unsigned long int k)
                     54: {
                     55:   unsigned long int  i, j;
                     56:   mp_limb_t          nacc, kacc;
                     57:   unsigned long int  cnt;
1.1.1.2 ! ohara      58:   mp_size_t          rsize, ralloc;
        !            59:   mp_ptr             rp;
1.1       maekawa    60:
                     61:   /* bin(n,k) = 0 if k>n. */
                     62:   if (n < k)
                     63:     {
1.1.1.2 ! ohara      64:       SIZ(r) = 0;
1.1       maekawa    65:       return;
                     66:     }
                     67:
1.1.1.2 ! ohara      68:   rp = PTR(r);
        !            69:
1.1       maekawa    70:   /* Rewrite bin(n,k) as bin(n,n-k) if that is smaller. */
                     71:   k = MIN (k, n-k);
                     72:
                     73:   /* bin(n,0) = 1 */
                     74:   if (k == 0)
                     75:     {
1.1.1.2 ! ohara      76:       SIZ(r) = 1;
        !            77:       rp[0] = 1;
1.1       maekawa    78:       return;
                     79:     }
                     80:
                     81:   j = n - k + 1;
1.1.1.2 ! ohara      82:   rp[0] = j;
        !            83:   rsize = 1;
        !            84:   ralloc = ALLOC(r);
1.1       maekawa    85:
                     86:   /* Initialize accumulators.  */
                     87:   nacc = 1;
                     88:   kacc = 1;
                     89:
                     90:   cnt = 0;
                     91:   for (i = 2; i <= k; i++)
                     92:     {
                     93:       mp_limb_t n1, n0, k1, k0;
                     94:
                     95:       j++;
                     96: #if 0
                     97:       /* Remove common multiples of 2.  This will allow us to accumulate
                     98:          more in nacc and kacc before we need a bignum step.  It would make
                     99:          sense to cancel factors of 3, 5, etc too, but this would be best
                    100:          handled by sieving out factors.  Alternatively, we could perform a
                    101:          gcd of the accumulators just as they have overflown, and keep
                    102:          accumulating until the gcd doesn't remove a significant factor.  */
                    103:       while (((nacc | kacc) & 1) == 0)
                    104:         {
                    105:           nacc >>= 1;
                    106:           kacc >>= 1;
                    107:         }
                    108: #else
                    109:       cnt = ((nacc | kacc) & 1) ^ 1;
                    110:       nacc >>= cnt;
                    111:       kacc >>= cnt;
                    112: #endif
                    113:       /* Accumulate next multiples.  */
1.1.1.2 ! ohara     114:       umul_ppmm (n1, n0, nacc, j << GMP_NAIL_BITS);
        !           115:       umul_ppmm (k1, k0, kacc, i << GMP_NAIL_BITS);
        !           116:       n0 >>= GMP_NAIL_BITS;
        !           117:       k0 >>= GMP_NAIL_BITS;
1.1       maekawa   118:       if (n1 != 0)
                    119:         {
1.1.1.2 ! ohara     120:           /* Accumulator overflow.  Perform bignum step. */
        !           121:           MULDIV (32);
1.1       maekawa   122:           nacc = j;
                    123:           kacc = i;
                    124:         }
                    125:       else
                    126:         {
1.1.1.2 ! ohara     127:           ASSERT (k1 == 0); /* n>=k, so high k zero when high n zero */
        !           128:
1.1       maekawa   129:           /* Save new products in accumulators to keep accumulating.  */
                    130:           nacc = n0;
                    131:           kacc = k0;
                    132:         }
                    133:     }
                    134:
                    135:   /* Take care of whatever is left in accumulators.  */
1.1.1.2 ! ohara     136:   MULDIV (1);
        !           137:
        !           138:   ALLOC(r) = ralloc;
        !           139:   SIZ(r) = rsize;
        !           140:   PTR(r) = rp;
1.1       maekawa   141: }