← Back to branch summary

~ubuntu-branches/ubuntu/quantal/gclcvs/quantal

« back to all changes in this revision

Viewing changes to gmp3/mpz/fac_ui.c

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2004-06-24 15:13:46 UTC
  • Revision ID: james.westby@ubuntu.com-20040624151346-xh0xaaktyyp7aorc
Tags: 2.7.0-26
C_GC_OFFSET is 2 on m68k-linux

Show diffs side-by-side

added added

removed removed

Lines of Context:
gmp3/mpz/fac_ui.c
 
1
/* mpz_fac_ui(result, n) -- Set RESULT to N!.
 
2
 
 
3
Copyright 1991, 1993, 1994, 1995, 2000, 2001 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
/* Enhancements:
 
28
 
 
29
   Data tables could be used for results up to 3 or 4 limbs to avoid
 
30
   fiddling around with small quantities.
 
31
 
 
32
   The product accumulation might be worth splitting out into something that
 
33
   could be used elsewhere too.
 
34
 
 
35
   The tree of partial products should be done with TMP_ALLOC, not mpz_init.
 
36
   It should be possible to know a maximum size at each level.
 
37
 
 
38
   Factors of two could be stripped from k to save some multiplying (but not
 
39
   very much).  The same could be done with factors of 3, perhaps.
 
40
 
 
41
   The prime factorization of n! is easy to determine, it might be worth
 
42
   using that rather than a simple 1 to n.  The powering of primes could do
 
43
   some squaring instead of multiplying.  There's probably other ways to use
 
44
   some squaring too.  */
 
45
 
 
46
 
 
47
/* for single non-zero limb */
 
48
#define MPZ_SET_1_NZ(z,n)       \
 
49
  do {                          \
 
50
    mpz_ptr  __z = (z);         \
 
51
    ASSERT ((n) != 0);          \
 
52
    PTR(__z)[0] = (n);          \
 
53
    SIZ(__z) = 1;               \
 
54
  } while (0)
 
55
 
 
56
/* for single non-zero limb */
 
57
#define MPZ_INIT_SET_1_NZ(z,n)                  \
 
58
  do {                                          \
 
59
    mpz_ptr  __iz = (z);                        \
 
60
    ALLOC(__iz) = 1;                            \
 
61
    PTR(__iz) = __GMP_ALLOCATE_FUNC_LIMBS (1);  \
 
62
    MPZ_SET_1_NZ (__iz, n);                     \
 
63
  } while (0)
 
64
 
 
65
/* for src>0 and n>0 */
 
66
#define MPZ_MUL_1_POS(dst,src,n)                        \
 
67
  do {                                                  \
 
68
    mpz_ptr    __dst = (dst);                           \
 
69
    mpz_srcptr __src = (src);                           \
 
70
    mp_size_t  __size = SIZ(__src);                     \
 
71
    mp_ptr     __dst_p;                                 \
 
72
    mp_limb_t  __c;                                     \
 
73
                                                        \
 
74
    ASSERT (__size > 0);                                \
 
75
    ASSERT ((n) != 0);                                  \
 
76
                                                        \
 
77
    MPZ_REALLOC (__dst, __size+1);                      \
 
78
    __dst_p = PTR(__dst);                               \
 
79
                                                        \
 
80
    __c = mpn_mul_1 (__dst_p, PTR(__src), __size, n);   \
 
81
    __dst_p[__size] = __c;                              \
 
82
    SIZ(__dst) = __size + (__c != 0);                   \
 
83
                                                        \
 
84
  } while (0)
 
85
 
 
86
 
 
87
void
 
88
mpz_fac_ui (mpz_ptr result, unsigned long int n)
 
89
{
 
90
#if SIMPLE_FAC
 
91
  /* Be silly.  Just multiply the numbers in ascending order.  O(n**2).  */
 
92
  unsigned long int k;
 
93
  mpz_set_ui (result, 1L);
 
94
  for (k = 2; k <= n; k++)
 
95
    mpz_mul_ui (result, result, k);
 
96
#else
 
97
 
 
98
  /* Be smarter.  Multiply groups of numbers in ascending order until the
 
99
     product doesn't fit in a limb.  Multiply these partial product in a
 
100
     balanced binary tree fashion, to make the operand have as equal sizes
 
101
     as possible.  When the operands have about the same size, mpn_mul
 
102
     becomes faster.  */
 
103
 
 
104
  unsigned long  k;
 
105
  mp_limb_t      p, p1, p0;
 
106
 
 
107
  /* Stack of partial products, used to make the computation balanced
 
108
     (i.e. make the sizes of the multiplication operands equal).  The
 
109
     topmost position of MP_STACK will contain a one-limb partial product,
 
110
     the second topmost will contain a two-limb partial product, and so
 
111
     on.  MP_STACK[0] will contain a partial product with 2**t limbs.
 
112
     To compute n! MP_STACK needs to be less than
 
113
     log(n)**2/log(BITS_PER_MP_LIMB), so 30 is surely enough.  */
 
114
#define MP_STACK_SIZE 30
 
115
  mpz_t mp_stack[MP_STACK_SIZE];
 
116
 
 
117
  /* TOP is an index into MP_STACK, giving the topmost element.
 
118
     TOP_LIMIT_SO_FAR is the largets value it has taken so far.  */
 
119
  int top, top_limit_so_far;
 
120
 
 
121
  /* Count of the total number of limbs put on MP_STACK so far.  This
 
122
     variable plays an essential role in making the compututation balanced.
 
123
     See below.  */
 
124
  unsigned int tree_cnt;
 
125
 
 
126
  /* for testing with small limbs */
 
127
  if (MP_LIMB_T_MAX < ULONG_MAX)
 
128
    ASSERT_ALWAYS (n <= MP_LIMB_T_MAX);
 
129
 
 
130
  top = top_limit_so_far = -1;
 
131
  tree_cnt = 0;
 
132
  p = 1;
 
133
  for (k = 2; k <= n; k++)
 
134
    {
 
135
      /* Multiply the partial product in P with K.  */
 
136
      umul_ppmm (p1, p0, p, (mp_limb_t) k);
 
137
 
 
138
      /* Did we get overflow into the high limb, i.e. is the partial
 
139
         product now more than one limb?  */
 
140
      if (p1 != 0)
 
141
        {
 
142
          tree_cnt++;
 
143
 
 
144
          if (tree_cnt % 2 == 0)
 
145
            {
 
146
              mp_size_t i;
 
147
 
 
148
              /* TREE_CNT is even (i.e. we have generated an even number of
 
149
                 one-limb partial products), which means that we have a
 
150
                 single-limb product on the top of MP_STACK.  */
 
151
 
 
152
              MPZ_MUL_1_POS (mp_stack[top], mp_stack[top], p);
 
153
 
 
154
              /* If TREE_CNT is divisable by 4, 8,..., we have two
 
155
                 similar-sized partial products with 2, 4,... limbs at
 
156
                 the topmost two positions of MP_STACK.  Multiply them
 
157
                 to form a new partial product with 4, 8,... limbs.  */
 
158
              for (i = 4; (tree_cnt & (i - 1)) == 0; i <<= 1)
 
159
                {
 
160
                  mpz_mul (mp_stack[top - 1],
 
161
                           mp_stack[top], mp_stack[top - 1]);
 
162
                  top--;
 
163
                }
 
164
            }
 
165
          else
 
166
            {
 
167
              /* Put the single-limb partial product in P on the stack.
 
168
                 (The next time we get a single-limb product, we will
 
169
                 multiply the two together.)  */
 
170
              top++;
 
171
              if (top > top_limit_so_far)
 
172
                {
 
173
                  if (top > MP_STACK_SIZE)
 
174
                    abort();
 
175
                  /* The stack is now bigger than ever, initialize the top
 
176
                     element.  */
 
177
                  MPZ_INIT_SET_1_NZ (mp_stack[top], p);
 
178
                  top_limit_so_far++;
 
179
                }
 
180
              else
 
181
                MPZ_SET_1_NZ (mp_stack[top], p);
 
182
            }
 
183
 
 
184
          /* We ignored the last result from umul_ppmm.  Put K in P as the
 
185
             first component of the next single-limb partial product.  */
 
186
          p = k;
 
187
        }
 
188
      else
 
189
        /* We didn't get overflow in umul_ppmm.  Put p0 in P and try
 
190
           with one more value of K.  */
 
191
        p = p0;
 
192
    }
 
193
 
 
194
  /* We have partial products in mp_stack[0..top], in descending order.
 
195
     We also have a small partial product in p.
 
196
     Their product is the final result.  */
 
197
  if (top < 0)
 
198
    MPZ_SET_1_NZ (result, p);
 
199
  else
 
200
    MPZ_MUL_1_POS (result, mp_stack[top--], p);
 
201
  while (top >= 0)
 
202
    mpz_mul (result, result, mp_stack[top--]);
 
203
 
 
204
  /* Free the storage allocated for MP_STACK.  */
 
205
  for (top = top_limit_so_far; top >= 0; top--)
 
206
    mpz_clear (mp_stack[top]);
 
207
#endif
 
208
}