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

Diff for /OpenXM_contrib/gmp/mpz/Attic/bin_uiui.c between version 1.1.1.1 and 1.1.1.2

version 1.1.1.1, 2000/09/09 14:12:47 version 1.1.1.2, 2003/08/25 16:06:32
Line 1 
Line 1 
 /* mpz_bin_uiui - compute n over k.  /* mpz_bin_uiui - compute n over k.
   
 Copyright (C) 1998, 1999, 2000 Free Software Foundation, Inc.  Copyright 1998, 1999, 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 24  MA 02111-1307, USA. */
Line 24  MA 02111-1307, USA. */
 #include "longlong.h"  #include "longlong.h"
   
   
 /* Avoid reallocs by rounding up any new size */  /* Enhancement: It ought to be possible to calculate the size of the final
 #define ROUNDUP_MASK  15     result in advance, to a rough approximation at least, and use it to do
      just one realloc.  Stirling's approximation n! ~= sqrt(2*pi*n)*(n/e)^n
      (Knuth section 1.2.5) might be of use.  */
   
 /* Enhancement: use mpn_divexact_1 when it exists */  /* "inc" in the main loop allocates a chunk more space if not already
 #define MULDIV()                                                \     enough, so as to avoid repeated reallocs.  The final step on the other
   MPZ_REALLOC (r, (SIZ(r)+1)|ROUNDUP_MASK);                     \     hand requires only one more limb.  */
   PTR(r)[SIZ(r)] = mpn_mul_1 (PTR(r), PTR(r), SIZ(r), nacc);    \  #define MULDIV(inc)                                                     \
   ASSERT_NOCARRY (mpn_divrem_1 (PTR(r), (mp_size_t) 0,          \    do {                                                                  \
                                 PTR(r), SIZ(r)+1, kacc));       \      ASSERT (rsize <= ralloc);                                           \
   SIZ(r) += (PTR(r)[SIZ(r)] != 0);                                                                          \
       if (rsize == ralloc)                                                \
         {                                                                 \
           mp_size_t  new_ralloc = ralloc + (inc);                         \
           rp = __GMP_REALLOCATE_FUNC_LIMBS (rp, ralloc, new_ralloc);      \
           ralloc = new_ralloc;                                            \
         }                                                                 \
                                                                           \
       rp[rsize] = mpn_mul_1 (rp, rp, rsize, nacc);                        \
       MPN_DIVREM_OR_DIVEXACT_1 (rp, rp, rsize+1, kacc);                   \
       rsize += (rp[rsize] != 0);                                          \
                                                                           \
   } while (0)
   
 void  void
 #if __STDC__  
 mpz_bin_uiui (mpz_ptr r, unsigned long int n, unsigned long int k)  mpz_bin_uiui (mpz_ptr r, unsigned long int n, unsigned long int k)
 #else  
 mpz_bin_uiui (r, n, k)  
      mpz_ptr r;  
      unsigned long int n;  
      unsigned long int k;  
 #endif  
 {  {
   unsigned long int  i, j;    unsigned long int  i, j;
   mp_limb_t          nacc, kacc;    mp_limb_t          nacc, kacc;
   unsigned long int  cnt;    unsigned long int  cnt;
     mp_size_t          rsize, ralloc;
     mp_ptr             rp;
   
   /* bin(n,k) = 0 if k>n. */    /* bin(n,k) = 0 if k>n. */
   if (n < k)    if (n < k)
     {      {
       mpz_set_ui (r, 0);        SIZ(r) = 0;
       return;        return;
     }      }
   
     rp = PTR(r);
   
   /* Rewrite bin(n,k) as bin(n,n-k) if that is smaller. */    /* Rewrite bin(n,k) as bin(n,n-k) if that is smaller. */
   k = MIN (k, n-k);    k = MIN (k, n-k);
   
   /* bin(n,0) = 1 */    /* bin(n,0) = 1 */
   if (k == 0)    if (k == 0)
     {      {
       mpz_set_ui (r, 1);        SIZ(r) = 1;
         rp[0] = 1;
       return;        return;
     }      }
   
   j = n - k + 1;    j = n - k + 1;
   mpz_set_ui (r, j);    rp[0] = j;
     rsize = 1;
     ralloc = ALLOC(r);
   
   /* Initialize accumulators.  */    /* Initialize accumulators.  */
   nacc = 1;    nacc = 1;
Line 97  mpz_bin_uiui (r, n, k)
Line 111  mpz_bin_uiui (r, n, k)
       kacc >>= cnt;        kacc >>= cnt;
 #endif  #endif
       /* Accumulate next multiples.  */        /* Accumulate next multiples.  */
       umul_ppmm (n1, n0, nacc, j);        umul_ppmm (n1, n0, nacc, j << GMP_NAIL_BITS);
       umul_ppmm (k1, k0, kacc, i);        umul_ppmm (k1, k0, kacc, i << GMP_NAIL_BITS);
         n0 >>= GMP_NAIL_BITS;
         k0 >>= GMP_NAIL_BITS;
       if (n1 != 0)        if (n1 != 0)
         {          {
           /* Accumulator overflow.  Perform bignum step.  */            /* Accumulator overflow.  Perform bignum step. */
           MULDIV ();            MULDIV (32);
           nacc = j;            nacc = j;
           kacc = i;            kacc = i;
         }          }
       else        else
         {          {
           if (k1 != 0) abort ();            ASSERT (k1 == 0); /* n>=k, so high k zero when high n zero */
   
           /* Save new products in accumulators to keep accumulating.  */            /* Save new products in accumulators to keep accumulating.  */
           nacc = n0;            nacc = n0;
           kacc = k0;            kacc = k0;
Line 116  mpz_bin_uiui (r, n, k)
Line 133  mpz_bin_uiui (r, n, k)
     }      }
   
   /* Take care of whatever is left in accumulators.  */    /* Take care of whatever is left in accumulators.  */
   MULDIV ();    MULDIV (1);
   
     ALLOC(r) = ralloc;
     SIZ(r) = rsize;
     PTR(r) = rp;
 }  }
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.1.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>