1 /* mpz_bin_uiui - compute n over k. 2 3 Copyright 1998, 1999, 2000, 2001, 2002, 2003, 2006 Free Software Foundation, 4 Inc. 5 6 This file is part of the GNU MP Library. 7 8 The GNU MP Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MP Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 20 21 #include "gmp.h" 22 #include "gmp-impl.h" 23 #include "longlong.h" 24 25 26 /* Enhancement: It ought to be possible to calculate the size of the final 27 result in advance, to a rough approximation at least, and use it to do 28 just one realloc. Stirling's approximation n! ~= sqrt(2*pi*n)*(n/e)^n 29 (Knuth section 1.2.5) might be of use. */ 30 31 /* "inc" in the main loop allocates a chunk more space if not already 32 enough, so as to avoid repeated reallocs. The final step on the other 33 hand requires only one more limb. */ 34 #define MULDIV(inc) \ 35 do { \ 36 ASSERT (rsize <= ralloc); \ 37 \ 38 if (rsize == ralloc) \ 39 { \ 40 mp_size_t new_ralloc = ralloc + (inc); \ 41 rp = __GMP_REALLOCATE_FUNC_LIMBS (rp, ralloc, new_ralloc); \ 42 ralloc = new_ralloc; \ 43 } \ 44 \ 45 rp[rsize] = mpn_mul_1 (rp, rp, rsize, nacc); \ 46 MPN_DIVREM_OR_DIVEXACT_1 (rp, rp, rsize+1, kacc); \ 47 rsize += (rp[rsize] != 0); \ 48 \ 49 } while (0) 50 51 void 52 mpz_bin_uiui (mpz_ptr r, unsigned long int n, unsigned long int k) 53 { 54 unsigned long int i, j; 55 mp_limb_t nacc, kacc; 56 unsigned long int cnt; 57 mp_size_t rsize, ralloc; 58 mp_ptr rp; 59 60 /* bin(n,k) = 0 if k>n. */ 61 if (n < k) 62 { 63 SIZ(r) = 0; 64 return; 65 } 66 67 rp = PTR(r); 68 69 /* Rewrite bin(n,k) as bin(n,n-k) if that is smaller. */ 70 k = MIN (k, n-k); 71 72 /* bin(n,0) = 1 */ 73 if (k == 0) 74 { 75 SIZ(r) = 1; 76 rp[0] = 1; 77 return; 78 } 79 80 j = n - k + 1; 81 rp[0] = j; 82 rsize = 1; 83 ralloc = ALLOC(r); 84 85 /* Initialize accumulators. */ 86 nacc = 1; 87 kacc = 1; 88 89 for (i = 2; i <= k; i++) 90 { 91 mp_limb_t n1, n0; 92 93 /* Remove common 2 factors. */ 94 cnt = ((nacc | kacc) & 1) ^ 1; 95 nacc >>= cnt; 96 kacc >>= cnt; 97 98 j++; 99 /* Accumulate next multiples. */ 100 umul_ppmm (n1, n0, nacc, (mp_limb_t) j << GMP_NAIL_BITS); 101 n0 >>= GMP_NAIL_BITS; 102 if (n1 == 0) 103 { 104 /* Save new products in accumulators to keep accumulating. */ 105 nacc = n0; 106 kacc = kacc * i; 107 } 108 else 109 { 110 /* Accumulator overflow. Perform bignum step. */ 111 MULDIV (32); 112 nacc = j; 113 kacc = i; 114 } 115 } 116 117 /* Take care of whatever is left in accumulators. */ 118 MULDIV (1); 119 120 ALLOC(r) = ralloc; 121 SIZ(r) = rsize; 122 PTR(r) = rp; 123 } 124