home *** CD-ROM | disk | FTP | other *** search
/ Fresh Fish 4 / FreshFish_May-June1994.bin / bbs / gnu / gmp-1.3.2-src.lha / src / amiga / gmp-1.3.2 / mpz_fac_ui.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-03  |  4.6 KB  |  157 lines
  1. /* mpz_fac_ui(result, n) -- Set RESULT to N!.
  2.  
  3. Copyright (C) 1991 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 General Public License as published by
  9. the Free Software Foundation; either version 2, or (at your option)
  10. any later version.
  11.  
  12. The GNU MP Library is distributed in the hope that it will be useful,
  13. but WITHOUT ANY WARRANTY; without even the implied warranty of
  14. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15. GNU General Public License for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with the GNU MP Library; see the file COPYING.  If not, write to
  19. the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  20.  
  21. #ifdef DBG
  22. #include <stdio.h>
  23. #endif
  24.  
  25. #include "gmp.h"
  26. #include "gmp-impl.h"
  27. #include "longlong.h"
  28.  
  29. void
  30. #ifdef __STDC__
  31. mpz_fac_ui (MP_INT *result, unsigned long int n)
  32. #else
  33. mpz_fac_ui (result, n)
  34.      MP_INT *result;
  35.      unsigned long int n;
  36. #endif
  37. {
  38. #if SIMPLE_FAC
  39.  
  40.   /* Be silly.  Just multiply the numbers in ascending order.  O(n**2).  */
  41.  
  42.   mp_limb k;
  43.  
  44.   mpz_set_ui (result, (mp_limb) 1);
  45.  
  46.   for (k = 2; k <= n; k++)
  47.     mpz_mul_ui (result, result, k);
  48. #else
  49.  
  50.   /* Be smarter.  Multiply groups of numbers in ascending order until the
  51.      product doesn't fit in a limb.  Multiply these partial products in a
  52.      balanced binary tree fashion, to make the operand have as equal sizes
  53.      as possible.  (When the operands have about the same size, mpn_mul
  54.      becomes faster.)  */
  55.  
  56.   mp_limb k;
  57.   mp_limb p1, p0, p;
  58.  
  59.   /* Stack of partial products, used to make the computation balanced
  60.      (i.e. make the sizes of the multiplication operands equal).  The
  61.      topmost position of MP_STACK will contain a one-limb partial product,
  62.      the second topmost will contain a two-limb partial product, and so
  63.      on.  MP_STACK[0] will contain a partial product with 2**t limbs.
  64.      To compute n! MP_STACK needs to be less than
  65.      log(n)**2/log(BITS_PER_MP_LIMB), so 30 is surely enough.  */
  66. #define MP_STACK_SIZE 30
  67.   MP_INT mp_stack[MP_STACK_SIZE];
  68.  
  69.   /* TOP is an index into MP_STACK, giving the topmost element.
  70.      TOP_LIMIT_SO_FAR is the largets value it has taken so far.  */
  71.   int top, top_limit_so_far;
  72.  
  73.   /* Count of the total number of limbs put on MP_STACK so far.  This
  74.      variable plays an essential role in making the compututation balanced.
  75.      See below.  */
  76.   unsigned int tree_cnt;
  77.  
  78.   top = top_limit_so_far = -1;
  79.   tree_cnt = 0;
  80.   p = 1;
  81.   for (k = 2; k <= n; k++)
  82.     {
  83.       /* Multiply the partial product in P with K.  */
  84.       umul_ppmm (p1, p0, p, k);
  85.  
  86.       /* Did we get overflow into the high limb, i.e. is the partial
  87.      product now more than one limb?  */
  88.       if (p1 != 0)
  89.     {
  90.       tree_cnt++;
  91.  
  92.       if (tree_cnt % 2 == 0)
  93.         {
  94.           mp_size i;
  95.  
  96.           /* TREE_CNT is even (i.e. we have generated an even number of
  97.          one-limb partial products), which means that we have a
  98.          single-limb product on the top of MP_STACK.  */
  99.  
  100.           mpz_mul_ui (&mp_stack[top], &mp_stack[top], p);
  101.  
  102.           /* If TREE_CNT is divisable by 4, 8,..., we have two
  103.          similar-sized partial products with 2, 4,... limbs at
  104.          the topmost two positions of MP_STACK.  Multiply them
  105.          to form a new partial product with 4, 8,... limbs.  */
  106.           for (i = 4; (tree_cnt & (i - 1)) == 0; i <<= 1)
  107.         {
  108.           mpz_mul (&mp_stack[top - 1],
  109.                &mp_stack[top], &mp_stack[top - 1]);
  110.           top--;
  111.         }
  112.         }
  113.       else
  114.         {
  115.           /* Put the single-limb partial product in P on the stack.
  116.          (The next time we get a single-limb product, we will
  117.          multiply the two together.)  */
  118.           top++;
  119.           if (top > top_limit_so_far)
  120.         {
  121.           if (top > MP_STACK_SIZE)
  122.             abort();
  123.           /* The stack is now bigger than ever, initialize the top
  124.              element.  */
  125.           mpz_init_set_ui (&mp_stack[top], p);
  126.           top_limit_so_far++;
  127.         }
  128.           else
  129.         mpz_set_ui (&mp_stack[top], p);
  130.         }
  131.  
  132.       /* We ignored the last result from umul_ppmm.  Put K in P as the
  133.          first component of the next single-limb partial product.  */
  134.       p = k;
  135.     }
  136.       else
  137.     /* We didn't get overflow in umul_ppmm.  Put p0 in P and try
  138.        with one more value of K.  */
  139.     p = p0;
  140.     }
  141.  
  142.   /* We have partial products in mp_stack[0..top], in descending order.
  143.      We also have a small partial product in p.
  144.      Their product is the final result.  */
  145.   if (top < 0)
  146.     mpz_set_ui (result, p);
  147.   else
  148.     mpz_mul_ui (result, &mp_stack[top--], p);
  149.   while (top >= 0)
  150.     mpz_mul (result, result, &mp_stack[top--]);
  151.  
  152.   /* Free the storage allocated for MP_STACK.  */
  153.   for (top = top_limit_so_far; top >= 0; top--)
  154.     mpz_clear (&mp_stack[top]);
  155. #endif
  156. }
  157.