1 /* mpz_fac_ui(result, n) -- Set RESULT to N!.
3 Copyright (C) 1991, 1993, 1994, 1995 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. */
32 mpz_fac_ui (mpz_ptr result, unsigned long int n)
34 mpz_fac_ui (result, n)
41 /* Be silly. Just multiply the numbers in ascending order. O(n**2). */
45 mpz_set_ui (result, 1L);
47 for (k = 2; k <= n; k++)
48 mpz_mul_ui (result, result, k);
51 /* Be smarter. Multiply groups of numbers in ascending order until the
52 product doesn't fit in a limb. Multiply these partial product in a
53 balanced binary tree fashion, to make the operand have as equal sizes
54 as possible. When the operands have about the same size, mpn_mul
57 unsigned long int p, k;
60 /* Stack of partial products, used to make the computation balanced
61 (i.e. make the sizes of the multiplication operands equal). The
62 topmost position of MP_STACK will contain a one-limb partial product,
63 the second topmost will contain a two-limb partial product, and so
64 on. MP_STACK[0] will contain a partial product with 2**t limbs.
65 To compute n! MP_STACK needs to be less than
66 log(n)**2/log(BITS_PER_MP_LIMB), so 30 is surely enough. */
67 #define MP_STACK_SIZE 30
68 mpz_t mp_stack[MP_STACK_SIZE];
70 /* TOP is an index into MP_STACK, giving the topmost element.
71 TOP_LIMIT_SO_FAR is the largets value it has taken so far. */
72 int top, top_limit_so_far;
74 /* Count of the total number of limbs put on MP_STACK so far. This
75 variable plays an essential role in making the compututation balanced.
77 unsigned int tree_cnt;
79 top = top_limit_so_far = -1;
82 for (k = 2; k <= n; k++)
84 /* Multiply the partial product in P with K. */
85 umul_ppmm (p1, p0, (mp_limb_t) p, (mp_limb_t) k);
87 /* Did we get overflow into the high limb, i.e. is the partial
88 product now more than one limb? */
93 if (tree_cnt % 2 == 0)
97 /* TREE_CNT is even (i.e. we have generated an even number of
98 one-limb partial products), which means that we have a
99 single-limb product on the top of MP_STACK. */
101 mpz_mul_ui (mp_stack[top], mp_stack[top], p);
103 /* If TREE_CNT is divisable by 4, 8,..., we have two
104 similar-sized partial products with 2, 4,... limbs at
105 the topmost two positions of MP_STACK. Multiply them
106 to form a new partial product with 4, 8,... limbs. */
107 for (i = 4; (tree_cnt & (i - 1)) == 0; i <<= 1)
109 mpz_mul (mp_stack[top - 1],
110 mp_stack[top], mp_stack[top - 1]);
116 /* Put the single-limb partial product in P on the stack.
117 (The next time we get a single-limb product, we will
118 multiply the two together.) */
120 if (top > top_limit_so_far)
122 if (top > MP_STACK_SIZE)
124 /* The stack is now bigger than ever, initialize the top
126 mpz_init_set_ui (mp_stack[top], p);
130 mpz_set_ui (mp_stack[top], p);
133 /* We ignored the last result from umul_ppmm. Put K in P as the
134 first component of the next single-limb partial product. */
138 /* We didn't get overflow in umul_ppmm. Put p0 in P and try
139 with one more value of K. */
140 p = p0; /* bogus if long != mp_limb_t */
143 /* We have partial products in mp_stack[0..top], in descending order.
144 We also have a small partial product in p.
145 Their product is the final result. */
147 mpz_set_ui (result, p);
149 mpz_mul_ui (result, mp_stack[top--], p);
151 mpz_mul (result, result, mp_stack[top--]);
153 /* Free the storage allocated for MP_STACK. */
154 for (top = top_limit_so_far; top >= 0; top--)
155 mpz_clear (mp_stack[top]);