rts/gmp/mpz/fac_ui.c

   1 /* mpz_fac_ui(result, n) -- Set RESULT to N!.
   2
   3 Copyright (C) 1991, 1993, 1994, 1995 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 #ifdef DBG
  23 #include <stdio.h>
  24 #endif
  25
  26 #include "gmp.h"
  27 #include "gmp-impl.h"
  28 #include "longlong.h"
  29
  30 void
  31 #if __STDC__
  32 mpz_fac_ui (mpz_ptr result, unsigned long int n)
  33 #else
  34 mpz_fac_ui (result, n)
  35      mpz_ptr result;
  36      unsigned long int n;
  37 #endif
  38 {
  39 #if SIMPLE_FAC
  40
  41   /* Be silly.  Just multiply the numbers in ascending order.  O(n**2).  */
  42
  43   unsigned long int k;
  44
  45   mpz_set_ui (result, 1L);
  46
  47   for (k = 2; k <= n; k++)
  48     mpz_mul_ui (result, result, k);
  49 #else
  50
  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
  55      becomes faster.  */
  56
  57   unsigned long int p, k;
  58   mp_limb_t p1, p0;
  59
  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];
  69
  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;
  73
  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.
  76      See below.  */
  77   unsigned int tree_cnt;
  78
  79   top = top_limit_so_far = -1;
  80   tree_cnt = 0;
  81   p = 1;
  82   for (k = 2; k <= n; k++)
  83     {
  84       /* Multiply the partial product in P with K.  */
  85       umul_ppmm (p1, p0, (mp_limb_t) p, (mp_limb_t) k);
  86
  87       /* Did we get overflow into the high limb, i.e. is the partial
  88          product now more than one limb?  */
  89       if (p1 != 0)
  90         {
  91           tree_cnt++;
  92
  93           if (tree_cnt % 2 == 0)
  94             {
  95               mp_size_t i;
  96
  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.  */
 100
 101               mpz_mul_ui (mp_stack[top], mp_stack[top], p);
 102
 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)
 108                 {
 109                   mpz_mul (mp_stack[top - 1],
 110                            mp_stack[top], mp_stack[top - 1]);
 111                   top--;
 112                 }
 113             }
 114           else
 115             {
 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.)  */
 119               top++;
 120               if (top > top_limit_so_far)
 121                 {
 122                   if (top > MP_STACK_SIZE)
 123                     abort();
 124                   /* The stack is now bigger than ever, initialize the top
 125                      element.  */
 126                   mpz_init_set_ui (mp_stack[top], p);
 127                   top_limit_so_far++;
 128                 }
 129               else
 130                 mpz_set_ui (mp_stack[top], p);
 131             }
 132
 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.  */
 135           p = k;
 136         }
 137       else
 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 */
 141     }
 142
 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.  */
 146   if (top < 0)
 147     mpz_set_ui (result, p);
 148   else
 149     mpz_mul_ui (result, mp_stack[top--], p);
 150   while (top >= 0)
 151     mpz_mul (result, result, mp_stack[top--]);
 152
 153   /* Free the storage allocated for MP_STACK.  */
 154   for (top = top_limit_so_far; top >= 0; top--)
 155     mpz_clear (mp_stack[top]);
 156 #endif
 157 }