version 1.1.1.1, 2000/09/09 14:12:47 |
version 1.1.1.2, 2003/08/25 16:06:32 |
|
|
/* mpz_bin_uiui - compute n over k. |
/* mpz_bin_uiui - compute n over k. |
|
|
Copyright (C) 1998, 1999, 2000 Free Software Foundation, Inc. |
Copyright 1998, 1999, 2000, 2001, 2002 Free Software Foundation, Inc. |
|
|
This file is part of the GNU MP Library. |
This file is part of the GNU MP Library. |
|
|
Line 24 MA 02111-1307, USA. */ |
|
Line 24 MA 02111-1307, USA. */ |
|
#include "longlong.h" |
#include "longlong.h" |
|
|
|
|
/* Avoid reallocs by rounding up any new size */ |
/* Enhancement: It ought to be possible to calculate the size of the final |
#define ROUNDUP_MASK 15 |
result in advance, to a rough approximation at least, and use it to do |
|
just one realloc. Stirling's approximation n! ~= sqrt(2*pi*n)*(n/e)^n |
|
(Knuth section 1.2.5) might be of use. */ |
|
|
/* Enhancement: use mpn_divexact_1 when it exists */ |
/* "inc" in the main loop allocates a chunk more space if not already |
#define MULDIV() \ |
enough, so as to avoid repeated reallocs. The final step on the other |
MPZ_REALLOC (r, (SIZ(r)+1)|ROUNDUP_MASK); \ |
hand requires only one more limb. */ |
PTR(r)[SIZ(r)] = mpn_mul_1 (PTR(r), PTR(r), SIZ(r), nacc); \ |
#define MULDIV(inc) \ |
ASSERT_NOCARRY (mpn_divrem_1 (PTR(r), (mp_size_t) 0, \ |
do { \ |
PTR(r), SIZ(r)+1, kacc)); \ |
ASSERT (rsize <= ralloc); \ |
SIZ(r) += (PTR(r)[SIZ(r)] != 0); |
\ |
|
if (rsize == ralloc) \ |
|
{ \ |
|
mp_size_t new_ralloc = ralloc + (inc); \ |
|
rp = __GMP_REALLOCATE_FUNC_LIMBS (rp, ralloc, new_ralloc); \ |
|
ralloc = new_ralloc; \ |
|
} \ |
|
\ |
|
rp[rsize] = mpn_mul_1 (rp, rp, rsize, nacc); \ |
|
MPN_DIVREM_OR_DIVEXACT_1 (rp, rp, rsize+1, kacc); \ |
|
rsize += (rp[rsize] != 0); \ |
|
\ |
|
} while (0) |
|
|
void |
void |
#if __STDC__ |
|
mpz_bin_uiui (mpz_ptr r, unsigned long int n, unsigned long int k) |
mpz_bin_uiui (mpz_ptr r, unsigned long int n, unsigned long int k) |
#else |
|
mpz_bin_uiui (r, n, k) |
|
mpz_ptr r; |
|
unsigned long int n; |
|
unsigned long int k; |
|
#endif |
|
{ |
{ |
unsigned long int i, j; |
unsigned long int i, j; |
mp_limb_t nacc, kacc; |
mp_limb_t nacc, kacc; |
unsigned long int cnt; |
unsigned long int cnt; |
|
mp_size_t rsize, ralloc; |
|
mp_ptr rp; |
|
|
/* bin(n,k) = 0 if k>n. */ |
/* bin(n,k) = 0 if k>n. */ |
if (n < k) |
if (n < k) |
{ |
{ |
mpz_set_ui (r, 0); |
SIZ(r) = 0; |
return; |
return; |
} |
} |
|
|
|
rp = PTR(r); |
|
|
/* Rewrite bin(n,k) as bin(n,n-k) if that is smaller. */ |
/* Rewrite bin(n,k) as bin(n,n-k) if that is smaller. */ |
k = MIN (k, n-k); |
k = MIN (k, n-k); |
|
|
/* bin(n,0) = 1 */ |
/* bin(n,0) = 1 */ |
if (k == 0) |
if (k == 0) |
{ |
{ |
mpz_set_ui (r, 1); |
SIZ(r) = 1; |
|
rp[0] = 1; |
return; |
return; |
} |
} |
|
|
j = n - k + 1; |
j = n - k + 1; |
mpz_set_ui (r, j); |
rp[0] = j; |
|
rsize = 1; |
|
ralloc = ALLOC(r); |
|
|
/* Initialize accumulators. */ |
/* Initialize accumulators. */ |
nacc = 1; |
nacc = 1; |
Line 97 mpz_bin_uiui (r, n, k) |
|
Line 111 mpz_bin_uiui (r, n, k) |
|
kacc >>= cnt; |
kacc >>= cnt; |
#endif |
#endif |
/* Accumulate next multiples. */ |
/* Accumulate next multiples. */ |
umul_ppmm (n1, n0, nacc, j); |
umul_ppmm (n1, n0, nacc, j << GMP_NAIL_BITS); |
umul_ppmm (k1, k0, kacc, i); |
umul_ppmm (k1, k0, kacc, i << GMP_NAIL_BITS); |
|
n0 >>= GMP_NAIL_BITS; |
|
k0 >>= GMP_NAIL_BITS; |
if (n1 != 0) |
if (n1 != 0) |
{ |
{ |
/* Accumulator overflow. Perform bignum step. */ |
/* Accumulator overflow. Perform bignum step. */ |
MULDIV (); |
MULDIV (32); |
nacc = j; |
nacc = j; |
kacc = i; |
kacc = i; |
} |
} |
else |
else |
{ |
{ |
if (k1 != 0) abort (); |
ASSERT (k1 == 0); /* n>=k, so high k zero when high n zero */ |
|
|
/* Save new products in accumulators to keep accumulating. */ |
/* Save new products in accumulators to keep accumulating. */ |
nacc = n0; |
nacc = n0; |
kacc = k0; |
kacc = k0; |
Line 116 mpz_bin_uiui (r, n, k) |
|
Line 133 mpz_bin_uiui (r, n, k) |
|
} |
} |
|
|
/* Take care of whatever is left in accumulators. */ |
/* Take care of whatever is left in accumulators. */ |
MULDIV (); |
MULDIV (1); |
|
|
|
ALLOC(r) = ralloc; |
|
SIZ(r) = rsize; |
|
PTR(r) = rp; |
} |
} |