← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/ecl/intrepid

« back to all changes in this revision

Viewing changes to src/gmp/mpz/bin_uiui.c

  • Committer: Bazaar Package Importer
  • Author(s): Peter Van Eynde
  • Date: 2006-05-17 02:46:26 UTC
  • Revision ID: james.westby@ubuntu.com-20060517024626-lljr08ftv9g9vefl
Tags: upstream-0.9h-20060510
Import upstream version 0.9h-20060510

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/gmp/mpz/bin_uiui.c
 
1
/* mpz_bin_uiui - compute n over k.
 
2
 
 
3
Copyright 1998, 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
 
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
 
 
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)
 
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;
 
58
  mp_size_t          rsize, ralloc;
 
59
  mp_ptr             rp;
 
60
 
 
61
  /* bin(n,k) = 0 if k>n. */
 
62
  if (n < k)
 
63
    {
 
64
      SIZ(r) = 0;
 
65
      return;
 
66
    }
 
67
 
 
68
  rp = PTR(r);
 
69
 
 
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
    {
 
76
      SIZ(r) = 1;
 
77
      rp[0] = 1;
 
78
      return;
 
79
    }
 
80
 
 
81
  j = n - k + 1;
 
82
  rp[0] = j;
 
83
  rsize = 1;
 
84
  ralloc = ALLOC(r);
 
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.  */
 
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;
 
118
      if (n1 != 0)
 
119
        {
 
120
          /* Accumulator overflow.  Perform bignum step. */
 
121
          MULDIV (32);
 
122
          nacc = j;
 
123
          kacc = i;
 
124
        }
 
125
      else
 
126
        {
 
127
          ASSERT (k1 == 0); /* n>=k, so high k zero when high n zero */
 
128
 
 
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.  */
 
136
  MULDIV (1);
 
137
 
 
138
  ALLOC(r) = ralloc;
 
139
  SIZ(r) = rsize;
 
140
  PTR(r) = rp;
 
141
}