1 /* mpz_bin_uiui - compute n over k.
3 Copyright (C) 1998, 1999, 2000 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
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.
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.
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. */
27 /* Avoid reallocs by rounding up any new size */
28 #define ROUNDUP_MASK 15
30 /* Enhancement: use mpn_divexact_1 when it exists */
32 MPZ_REALLOC (r, (SIZ(r)+1)|ROUNDUP_MASK); \
33 PTR(r)[SIZ(r)] = mpn_mul_1 (PTR(r), PTR(r), SIZ(r), nacc); \
34 ASSERT_NOCARRY (mpn_divrem_1 (PTR(r), (mp_size_t) 0, \
35 PTR(r), SIZ(r)+1, kacc)); \
36 SIZ(r) += (PTR(r)[SIZ(r)] != 0);
40 mpz_bin_uiui (mpz_ptr r, unsigned long int n, unsigned long int k)
42 mpz_bin_uiui (r, n, k)
48 unsigned long int i, j;
50 unsigned long int cnt;
52 /* bin(n,k) = 0 if k>n. */
59 /* Rewrite bin(n,k) as bin(n,n-k) if that is smaller. */
72 /* Initialize accumulators. */
77 for (i = 2; i <= k; i++)
79 mp_limb_t n1, n0, k1, k0;
83 /* Remove common multiples of 2. This will allow us to accumulate
84 more in nacc and kacc before we need a bignum step. It would make
85 sense to cancel factors of 3, 5, etc too, but this would be best
86 handled by sieving out factors. Alternatively, we could perform a
87 gcd of the accumulators just as they have overflown, and keep
88 accumulating until the gcd doesn't remove a significant factor. */
89 while (((nacc | kacc) & 1) == 0)
95 cnt = ((nacc | kacc) & 1) ^ 1;
99 /* Accumulate next multiples. */
100 umul_ppmm (n1, n0, nacc, j);
101 umul_ppmm (k1, k0, kacc, i);
104 /* Accumulator overflow. Perform bignum step. */
111 if (k1 != 0) abort ();
112 /* Save new products in accumulators to keep accumulating. */
118 /* Take care of whatever is left in accumulators. */