1*86d7f5d3SJohn Marino /* mpz_bin_ui - compute n over k.
2*86d7f5d3SJohn Marino
3*86d7f5d3SJohn Marino Copyright 1998, 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
4*86d7f5d3SJohn Marino
5*86d7f5d3SJohn Marino This file is part of the GNU MP Library.
6*86d7f5d3SJohn Marino
7*86d7f5d3SJohn Marino The GNU MP Library is free software; you can redistribute it and/or modify
8*86d7f5d3SJohn Marino it under the terms of the GNU Lesser General Public License as published by
9*86d7f5d3SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
10*86d7f5d3SJohn Marino option) any later version.
11*86d7f5d3SJohn Marino
12*86d7f5d3SJohn Marino The GNU MP Library is distributed in the hope that it will be useful, but
13*86d7f5d3SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14*86d7f5d3SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15*86d7f5d3SJohn Marino License for more details.
16*86d7f5d3SJohn Marino
17*86d7f5d3SJohn Marino You should have received a copy of the GNU Lesser General Public License
18*86d7f5d3SJohn Marino along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
19*86d7f5d3SJohn Marino
20*86d7f5d3SJohn Marino #include "gmp.h"
21*86d7f5d3SJohn Marino #include "gmp-impl.h"
22*86d7f5d3SJohn Marino #include "longlong.h"
23*86d7f5d3SJohn Marino
24*86d7f5d3SJohn Marino
25*86d7f5d3SJohn Marino /* This is a poor implementation. Look at bin_uiui.c for improvement ideas.
26*86d7f5d3SJohn Marino In fact consider calling mpz_bin_uiui() when the arguments fit, leaving
27*86d7f5d3SJohn Marino the code here only for big n.
28*86d7f5d3SJohn Marino
29*86d7f5d3SJohn Marino The identity bin(n,k) = (-1)^k * bin(-n+k-1,k) can be found in Knuth vol
30*86d7f5d3SJohn Marino 1 section 1.2.6 part G. */
31*86d7f5d3SJohn Marino
32*86d7f5d3SJohn Marino
33*86d7f5d3SJohn Marino #define DIVIDE() \
34*86d7f5d3SJohn Marino do { \
35*86d7f5d3SJohn Marino ASSERT (SIZ(r) > 0); \
36*86d7f5d3SJohn Marino MPN_DIVREM_OR_DIVEXACT_1 (PTR(r), PTR(r), (mp_size_t) SIZ(r), kacc); \
37*86d7f5d3SJohn Marino SIZ(r) -= (PTR(r)[SIZ(r)-1] == 0); \
38*86d7f5d3SJohn Marino } while (0)
39*86d7f5d3SJohn Marino
40*86d7f5d3SJohn Marino void
mpz_bin_ui(mpz_ptr r,mpz_srcptr n,unsigned long int k)41*86d7f5d3SJohn Marino mpz_bin_ui (mpz_ptr r, mpz_srcptr n, unsigned long int k)
42*86d7f5d3SJohn Marino {
43*86d7f5d3SJohn Marino mpz_t ni;
44*86d7f5d3SJohn Marino mp_limb_t i;
45*86d7f5d3SJohn Marino mpz_t nacc;
46*86d7f5d3SJohn Marino mp_limb_t kacc;
47*86d7f5d3SJohn Marino mp_size_t negate;
48*86d7f5d3SJohn Marino
49*86d7f5d3SJohn Marino if (mpz_sgn (n) < 0)
50*86d7f5d3SJohn Marino {
51*86d7f5d3SJohn Marino /* bin(n,k) = (-1)^k * bin(-n+k-1,k), and set ni = -n+k-1 - k = -n-1 */
52*86d7f5d3SJohn Marino mpz_init (ni);
53*86d7f5d3SJohn Marino mpz_neg (ni, n);
54*86d7f5d3SJohn Marino mpz_sub_ui (ni, ni, 1L);
55*86d7f5d3SJohn Marino negate = (k & 1); /* (-1)^k */
56*86d7f5d3SJohn Marino }
57*86d7f5d3SJohn Marino else
58*86d7f5d3SJohn Marino {
59*86d7f5d3SJohn Marino /* bin(n,k) == 0 if k>n
60*86d7f5d3SJohn Marino (no test for this under the n<0 case, since -n+k-1 >= k there) */
61*86d7f5d3SJohn Marino if (mpz_cmp_ui (n, k) < 0)
62*86d7f5d3SJohn Marino {
63*86d7f5d3SJohn Marino mpz_set_ui (r, 0L);
64*86d7f5d3SJohn Marino return;
65*86d7f5d3SJohn Marino }
66*86d7f5d3SJohn Marino
67*86d7f5d3SJohn Marino /* set ni = n-k */
68*86d7f5d3SJohn Marino mpz_init (ni);
69*86d7f5d3SJohn Marino mpz_sub_ui (ni, n, k);
70*86d7f5d3SJohn Marino negate = 0;
71*86d7f5d3SJohn Marino }
72*86d7f5d3SJohn Marino
73*86d7f5d3SJohn Marino /* Now wanting bin(ni+k,k), with ni positive, and "negate" is the sign (0
74*86d7f5d3SJohn Marino for positive, 1 for negative). */
75*86d7f5d3SJohn Marino mpz_set_ui (r, 1L);
76*86d7f5d3SJohn Marino
77*86d7f5d3SJohn Marino /* Rewrite bin(n,k) as bin(n,n-k) if that is smaller. In this case it's
78*86d7f5d3SJohn Marino whether ni+k-k < k meaning ni<k, and if so change to denominator ni+k-k
79*86d7f5d3SJohn Marino = ni, and new ni of ni+k-ni = k. */
80*86d7f5d3SJohn Marino if (mpz_cmp_ui (ni, k) < 0)
81*86d7f5d3SJohn Marino {
82*86d7f5d3SJohn Marino unsigned long tmp;
83*86d7f5d3SJohn Marino tmp = k;
84*86d7f5d3SJohn Marino k = mpz_get_ui (ni);
85*86d7f5d3SJohn Marino mpz_set_ui (ni, tmp);
86*86d7f5d3SJohn Marino }
87*86d7f5d3SJohn Marino
88*86d7f5d3SJohn Marino kacc = 1;
89*86d7f5d3SJohn Marino mpz_init_set_ui (nacc, 1L);
90*86d7f5d3SJohn Marino
91*86d7f5d3SJohn Marino for (i = 1; i <= k; i++)
92*86d7f5d3SJohn Marino {
93*86d7f5d3SJohn Marino mp_limb_t k1, k0;
94*86d7f5d3SJohn Marino
95*86d7f5d3SJohn Marino #if 0
96*86d7f5d3SJohn Marino mp_limb_t nacclow;
97*86d7f5d3SJohn Marino int c;
98*86d7f5d3SJohn Marino
99*86d7f5d3SJohn Marino nacclow = PTR(nacc)[0];
100*86d7f5d3SJohn Marino for (c = 0; (((kacc | nacclow) & 1) == 0); c++)
101*86d7f5d3SJohn Marino {
102*86d7f5d3SJohn Marino kacc >>= 1;
103*86d7f5d3SJohn Marino nacclow >>= 1;
104*86d7f5d3SJohn Marino }
105*86d7f5d3SJohn Marino mpz_div_2exp (nacc, nacc, c);
106*86d7f5d3SJohn Marino #endif
107*86d7f5d3SJohn Marino
108*86d7f5d3SJohn Marino mpz_add_ui (ni, ni, 1L);
109*86d7f5d3SJohn Marino mpz_mul (nacc, nacc, ni);
110*86d7f5d3SJohn Marino umul_ppmm (k1, k0, kacc, i << GMP_NAIL_BITS);
111*86d7f5d3SJohn Marino k0 >>= GMP_NAIL_BITS;
112*86d7f5d3SJohn Marino if (k1 != 0)
113*86d7f5d3SJohn Marino {
114*86d7f5d3SJohn Marino /* Accumulator overflow. Perform bignum step. */
115*86d7f5d3SJohn Marino mpz_mul (r, r, nacc);
116*86d7f5d3SJohn Marino mpz_set_ui (nacc, 1L);
117*86d7f5d3SJohn Marino DIVIDE ();
118*86d7f5d3SJohn Marino kacc = i;
119*86d7f5d3SJohn Marino }
120*86d7f5d3SJohn Marino else
121*86d7f5d3SJohn Marino {
122*86d7f5d3SJohn Marino /* Save new products in accumulators to keep accumulating. */
123*86d7f5d3SJohn Marino kacc = k0;
124*86d7f5d3SJohn Marino }
125*86d7f5d3SJohn Marino }
126*86d7f5d3SJohn Marino
127*86d7f5d3SJohn Marino mpz_mul (r, r, nacc);
128*86d7f5d3SJohn Marino DIVIDE ();
129*86d7f5d3SJohn Marino SIZ(r) = (SIZ(r) ^ -negate) + negate;
130*86d7f5d3SJohn Marino
131*86d7f5d3SJohn Marino mpz_clear (nacc);
132*86d7f5d3SJohn Marino mpz_clear (ni);
133*86d7f5d3SJohn Marino }
134