[BACK]Return to fac_ui.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gmp / mpz

Diff for /OpenXM_contrib/gmp/mpz/Attic/fac_ui.c between version 1.1.1.2 and 1.1.1.3

version 1.1.1.2, 2000/09/09 14:12:50 version 1.1.1.3, 2003/08/25 16:06:33
Line 1 
Line 1 
 /* mpz_fac_ui(result, n) -- Set RESULT to N!.  /* mpz_fac_ui(result, n) -- Set RESULT to N!.
   
 Copyright (C) 1991, 1993, 1994, 1995 Free Software Foundation, Inc.  Copyright 1991, 1993, 1994, 1995, 2000, 2001, 2002 Free Software Foundation,
   Inc.
   
 This file is part of the GNU MP Library.  This file is part of the GNU MP Library.
   
Line 19  along with the GNU MP Library; see the file COPYING.LI
Line 20  along with the GNU MP Library; see the file COPYING.LI
 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,  the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 MA 02111-1307, USA. */  MA 02111-1307, USA. */
   
 #ifdef DBG  
 #include <stdio.h>  
 #endif  
   
 #include "gmp.h"  #include "gmp.h"
 #include "gmp-impl.h"  #include "gmp-impl.h"
 #include "longlong.h"  #include "longlong.h"
   
   
   /* Enhancements:
   
      Data tables could be used for results up to 3 or 4 limbs to avoid
      fiddling around with small quantities.
   
      The product accumulation might be worth splitting out into something that
      could be used elsewhere too.
   
      The tree of partial products should be done with TMP_ALLOC, not mpz_init.
      It should be possible to know a maximum size at each level.
   
      Factors of two could be stripped from k to save some multiplying (but not
      very much).  The same could be done with factors of 3, perhaps.
   
      The prime factorization of n! is easy to determine, it might be worth
      using that rather than a simple 1 to n.  The powering of primes could do
      some squaring instead of multiplying.  There's probably other ways to use
      some squaring too.  */
   
   
   /* for single non-zero limb */
   #define MPZ_SET_1_NZ(z,n)       \
     do {                          \
       mpz_ptr  __z = (z);         \
       ASSERT ((n) != 0);          \
       PTR(__z)[0] = (n);          \
       SIZ(__z) = 1;               \
     } while (0)
   
   /* for single non-zero limb */
   #define MPZ_INIT_SET_1_NZ(z,n)                  \
     do {                                          \
       mpz_ptr  __iz = (z);                        \
       ALLOC(__iz) = 1;                            \
       PTR(__iz) = __GMP_ALLOCATE_FUNC_LIMBS (1);  \
       MPZ_SET_1_NZ (__iz, n);                     \
     } while (0)
   
   /* for src>0 and n>0 */
   #define MPZ_MUL_1_POS(dst,src,n)                        \
     do {                                                  \
       mpz_ptr    __dst = (dst);                           \
       mpz_srcptr __src = (src);                           \
       mp_size_t  __size = SIZ(__src);                     \
       mp_ptr     __dst_p;                                 \
       mp_limb_t  __c;                                     \
                                                           \
       ASSERT (__size > 0);                                \
       ASSERT ((n) != 0);                                  \
                                                           \
       MPZ_REALLOC (__dst, __size+1);                      \
       __dst_p = PTR(__dst);                               \
                                                           \
       __c = mpn_mul_1 (__dst_p, PTR(__src), __size, n);   \
       __dst_p[__size] = __c;                              \
       SIZ(__dst) = __size + (__c != 0);                   \
                                                           \
     } while (0)
   
   
 void  void
 #if __STDC__  
 mpz_fac_ui (mpz_ptr result, unsigned long int n)  mpz_fac_ui (mpz_ptr result, unsigned long int n)
 #else  
 mpz_fac_ui (result, n)  
      mpz_ptr result;  
      unsigned long int n;  
 #endif  
 {  {
 #if SIMPLE_FAC  #if SIMPLE_FAC
   
   /* Be silly.  Just multiply the numbers in ascending order.  O(n**2).  */    /* Be silly.  Just multiply the numbers in ascending order.  O(n**2).  */
   
   unsigned long int k;    unsigned long int k;
   
   mpz_set_ui (result, 1L);    mpz_set_ui (result, 1L);
   
   for (k = 2; k <= n; k++)    for (k = 2; k <= n; k++)
     mpz_mul_ui (result, result, k);      mpz_mul_ui (result, result, k);
 #else  #else
Line 54  mpz_fac_ui (result, n)
Line 102  mpz_fac_ui (result, n)
      as possible.  When the operands have about the same size, mpn_mul       as possible.  When the operands have about the same size, mpn_mul
      becomes faster.  */       becomes faster.  */
   
   unsigned long int p, k;    unsigned long  k;
   mp_limb_t p1, p0;    mp_limb_t      p, p1, p0;
   
   /* Stack of partial products, used to make the computation balanced    /* Stack of partial products, used to make the computation balanced
      (i.e. make the sizes of the multiplication operands equal).  The       (i.e. make the sizes of the multiplication operands equal).  The
Line 76  mpz_fac_ui (result, n)
Line 124  mpz_fac_ui (result, n)
      See below.  */       See below.  */
   unsigned int tree_cnt;    unsigned int tree_cnt;
   
     /* for testing with small limbs */
     if (MP_LIMB_T_MAX < ULONG_MAX)
       ASSERT_ALWAYS (n <= MP_LIMB_T_MAX);
   
   top = top_limit_so_far = -1;    top = top_limit_so_far = -1;
   tree_cnt = 0;    tree_cnt = 0;
   p = 1;    p = 1;
   for (k = 2; k <= n; k++)    for (k = 2; k <= n; k++)
     {      {
       /* Multiply the partial product in P with K.  */        /* Multiply the partial product in P with K.  */
       umul_ppmm (p1, p0, (mp_limb_t) p, (mp_limb_t) k);        umul_ppmm (p1, p0, p, (mp_limb_t) k);
   
   #if GMP_NAIL_BITS == 0
   #define OVERFLOW (p1 != 0)
   #else
   #define OVERFLOW ((p1 | (p0 >> GMP_NUMB_BITS)) != 0)
   #endif
       /* Did we get overflow into the high limb, i.e. is the partial        /* Did we get overflow into the high limb, i.e. is the partial
          product now more than one limb?  */           product now more than one limb?  */
       if (p1 != 0)        if OVERFLOW
         {          {
           tree_cnt++;            tree_cnt++;
   
Line 98  mpz_fac_ui (result, n)
Line 155  mpz_fac_ui (result, n)
                  one-limb partial products), which means that we have a                   one-limb partial products), which means that we have a
                  single-limb product on the top of MP_STACK.  */                   single-limb product on the top of MP_STACK.  */
   
               mpz_mul_ui (mp_stack[top], mp_stack[top], p);                MPZ_MUL_1_POS (mp_stack[top], mp_stack[top], p);
   
               /* If TREE_CNT is divisable by 4, 8,..., we have two                /* If TREE_CNT is divisable by 4, 8,..., we have two
                  similar-sized partial products with 2, 4,... limbs at                   similar-sized partial products with 2, 4,... limbs at
Line 123  mpz_fac_ui (result, n)
Line 180  mpz_fac_ui (result, n)
                     abort();                      abort();
                   /* The stack is now bigger than ever, initialize the top                    /* The stack is now bigger than ever, initialize the top
                      element.  */                       element.  */
                   mpz_init_set_ui (mp_stack[top], p);                    MPZ_INIT_SET_1_NZ (mp_stack[top], p);
                   top_limit_so_far++;                    top_limit_so_far++;
                 }                  }
               else                else
                 mpz_set_ui (mp_stack[top], p);                  MPZ_SET_1_NZ (mp_stack[top], p);
             }              }
   
           /* We ignored the last result from umul_ppmm.  Put K in P as the            /* We ignored the last result from umul_ppmm.  Put K in P as the
Line 137  mpz_fac_ui (result, n)
Line 194  mpz_fac_ui (result, n)
       else        else
         /* We didn't get overflow in umul_ppmm.  Put p0 in P and try          /* We didn't get overflow in umul_ppmm.  Put p0 in P and try
            with one more value of K.  */             with one more value of K.  */
         p = p0;                 /* bogus if long != mp_limb_t */          p = p0;
     }      }
   
   /* We have partial products in mp_stack[0..top], in descending order.    /* We have partial products in mp_stack[0..top], in descending order.
      We also have a small partial product in p.       We also have a small partial product in p.
      Their product is the final result.  */       Their product is the final result.  */
   if (top < 0)    if (top < 0)
     mpz_set_ui (result, p);      MPZ_SET_1_NZ (result, p);
   else    else
     mpz_mul_ui (result, mp_stack[top--], p);      MPZ_MUL_1_POS (result, mp_stack[top--], p);
   while (top >= 0)    while (top >= 0)
     mpz_mul (result, result, mp_stack[top--]);      mpz_mul (result, result, mp_stack[top--]);
   
Legend:
Removed from v.1.1.1.2  
changed lines
  Added in v.1.1.1.3
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>