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