gcc/gmp/mpz/fac_ui.c

/* mpz_fac_ui(result, n) -- Set RESULT to N!.

Copyright 1991, 1993, 1994, 1995, 2000, 2001, 2002, 2003 Free Software
Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */

#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"

#include "fac_ui.h"


static void odd_product __GMP_PROTO ((unsigned long, unsigned long, mpz_t *));
static void ap_product_small __GMP_PROTO ((mpz_t, mp_limb_t, mp_limb_t, unsigned long, unsigned long));


/* must be >=2	*/
#define APCONST	5

/* 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 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)


#if BITS_PER_ULONG == GMP_LIMB_BITS
#define BSWAP_ULONG(x,y)	BSWAP_LIMB(x,y)
#endif

/* We used to have a case here for limb==2*long, doing a BSWAP_LIMB followed
   by a shift down to get the high part.  But it provoked incorrect code
   from "HP aC++/ANSI C B3910B A.05.52 [Sep 05 2003]" in ILP32 mode.  This
   case would have been nice for gcc ia64 where BSWAP_LIMB is a mux1, but we
   can get that directly muxing a 4-byte ulong if it matters enough.  */

#if ! defined (BSWAP_ULONG)
#define BSWAP_ULONG(dst, src)						\
  do {									\
    unsigned long  __bswapl_src = (src);				\
    unsigned long  __bswapl_dst = 0;					\
    int	       __i;							\
    for (__i = 0; __i < sizeof(unsigned long); __i++)			\
      {									\
	__bswapl_dst = (__bswapl_dst << 8) | (__bswapl_src & 0xFF);	\
	__bswapl_src >>= 8;						\
      }									\
    (dst) = __bswapl_dst;						\
  } while (0)
#endif

/* x is bit reverse of y */
/* Note the divides below are all exact */
#define BITREV_ULONG(x,y)						   \
  do {									   \
   unsigned long __dst;							   \
   BSWAP_ULONG(__dst,y);						   \
   __dst = ((__dst>>4)&(ULONG_MAX/17)) | ((__dst<<4)&((ULONG_MAX/17)*16)); \
   __dst = ((__dst>>2)&(ULONG_MAX/5) ) | ((__dst<<2)&((ULONG_MAX/5)*4)  ); \
   __dst = ((__dst>>1)&(ULONG_MAX/3) ) | ((__dst<<1)&((ULONG_MAX/3)*2)  ); \
   (x) = __dst;								   \
  } while(0)
/* above could be improved if cpu has a nibble/bit swap/muxing instruction */
/* above code is serialized, possible to write as a big parallel expression */



void
mpz_fac_ui (mpz_ptr x, unsigned long n)
{
  unsigned long z, stt;
  int i, j;
  mpz_t t1, st[8 * sizeof (unsigned long) + 1 - APCONST];
  mp_limb_t d[4];

  static const mp_limb_t table[] = { ONE_LIMB_FACTORIAL_TABLE };

  if (n < numberof (table))
    {
      MPZ_SET_1_NZ (x, table[n]);
      return;
    }

  /*  NOTE : MUST have n>=3 here */
  ASSERT (n >= 3);
  /* for estimating the alloc sizes the calculation of these formula's is not
     exact and also the formulas are only approximations, also we ignore
     the few "side" calculations, correct allocation seems to speed up the
     small sizes better, having very little effect on the large sizes */

  /* estimate space for stack entries see below
     number of bits for n! is
     (1+log_2(2*pi)/2)-n*log_2(exp(1))+(n+1/2)*log_2(n)=
     2.325748065-n*1.442695041+(n+0.5)*log_2(n)  */
  umul_ppmm (d[1], d[0], (mp_limb_t) n, (mp_limb_t) FAC2OVERE);
  /* d[1] is 2n/e, d[0] ignored        */
  count_leading_zeros (z, d[1]);
  z = GMP_LIMB_BITS - z - 1;	/* z=floor(log_2(2n/e))   */
  umul_ppmm (d[1], d[0], (mp_limb_t) n, (mp_limb_t) z);
  /* d=n*floor(log_2(2n/e))   */
  d[0] = (d[0] >> 2) | (d[1] << (GMP_LIMB_BITS - 2));
  d[1] >>= 2;
  /* d=n*floor(log_2(2n/e))/4   */
  z = d[0] + 1;			/* have to ignore any overflow */
  /* so z is the number of bits wanted for st[0]    */


  if (n <= ((unsigned long) 1) << (APCONST))
    {
      mpz_realloc2 (x, 4 * z);
      ap_product_small (x, CNST_LIMB(2), CNST_LIMB(1), n - 1, 4L);
      return;
    }
  if (n <= ((unsigned long) 1) << (APCONST + 1))
    {				/*  use n!=odd(1,n)*(n/2)!*2^(n/2)         */
      mpz_init2 (t1, 2 * z);
      mpz_realloc2 (x, 4 * z);
      ap_product_small (x, CNST_LIMB(2), CNST_LIMB(1), n / 2 - 1, 4L);
      ap_product_small (t1, CNST_LIMB(3), CNST_LIMB(2), (n - 1) / 2, 4L);
      mpz_mul (x, x, t1);
      mpz_clear (t1);
      mpz_mul_2exp (x, x, n / 2);
      return;
    }
  if (n <= ((unsigned long) 1) << (APCONST + 2))
    {
      /* use n!=C_2(1,n/2)^2*C_2(n/2,n)*(n/4)!*2^(n/2+n/4) all int divs
	 so need (BITS_IN_N-APCONST+1)=(APCONST+3-APCONST+1)=4 stack entries */
      mpz_init2 (t1, 2 * z);
      mpz_realloc2 (x, 4 * z);
      for (i = 0; i < 4; i++)
	{
	  mpz_init2 (st[i], z);
	  z >>= 1;
	}
      odd_product (1, n / 2, st);
      mpz_set (x, st[0]);
      odd_product (n / 2, n, st);
      mpz_mul (x, x, x);
      ASSERT (n / 4 <= FACMUL4 + 6);
      ap_product_small (t1, CNST_LIMB(2), CNST_LIMB(1), n / 4 - 1, 4L);
      /* must have 2^APCONST odd numbers max */
      mpz_mul (t1, t1, st[0]);
      for (i = 0; i < 4; i++)
	mpz_clear (st[i]);
      mpz_mul (x, x, t1);
      mpz_clear (t1);
      mpz_mul_2exp (x, x, n / 2 + n / 4);
      return;
    }

  count_leading_zeros (stt, (mp_limb_t) n);
  stt = GMP_LIMB_BITS - stt + 1 - APCONST;

  for (i = 0; i < (signed long) stt; i++)
    {
      mpz_init2 (st[i], z);
      z >>= 1;
    }

  count_leading_zeros (z, (mp_limb_t) (n / 3));
  /* find z st 2^z>n/3 range for z is 1 <= z <= 8 * sizeof(unsigned long)-1 */
  z = GMP_LIMB_BITS - z;

  /*
     n! = 2^e * PRODUCT_{i=0}^{i=z-1} C_2( n/2^{i+1}, n/2^i )^{i+1}
     where 2^e || n!   3.2^z>n   C_2(a,b)=PRODUCT of odd z such that a<z<=b
   */


  mpz_init_set_ui (t1, 1);
  for (j = 8 * sizeof (unsigned long) / 2; j != 0; j >>= 1)
    {
      MPZ_SET_1_NZ (x, 1);
      for (i = 8 * sizeof (unsigned long) - j; i >= j; i -= 2 * j)
	if ((signed long) z >= i)
	  {
	    odd_product (n >> i, n >> (i - 1), st);
	    /* largest odd product when j=i=1 then we have
	       odd_product(n/2,n,st) which is approx (2n/e)^(n/4)
	       so log_base2(largest oddproduct)=n*log_base2(2n/e)/4
	       number of bits is n*log_base2(2n/e)/4+1  */
	    if (i != j)
	      mpz_pow_ui (st[0], st[0], i / j);
	    mpz_mul (x, x, st[0]);
	  }
      if ((signed long) z >= j && j != 1)
	{
	  mpz_mul (t1, t1, x);
	  mpz_mul (t1, t1, t1);
	}
    }
  for (i = 0; i < (signed long) stt; i++)
    mpz_clear (st[i]);
  mpz_mul (x, x, t1);
  mpz_clear (t1);
  popc_limb (i, (mp_limb_t) n);
  mpz_mul_2exp (x, x, n - i);
  return;
}

/* start,step are mp_limb_t although they will fit in unsigned long	*/
static void
ap_product_small (mpz_t ret, mp_limb_t start, mp_limb_t step,
		  unsigned long count, unsigned long nm)
{
  unsigned long a;
  mp_limb_t b;

  ASSERT (count <= (((unsigned long) 1) << APCONST));
/* count can never be zero ? check this and remove test below */
  if (count == 0)
    {
      MPZ_SET_1_NZ (ret, 1);
      return;
    }
  if (count == 1)
    {
      MPZ_SET_1_NZ (ret, start);
      return;
    }
  switch (nm)
    {
    case 1:
      MPZ_SET_1_NZ (ret, start);
      b = start + step;
      for (a = 0; a < count - 1; b += step, a++)
	MPZ_MUL_1_POS (ret, ret, b);
      return;
    case 2:
      MPZ_SET_1_NZ (ret, start * (start + step));
      if (count == 2)
	return;
      for (b = start + 2 * step, a = count / 2 - 1; a != 0;
	   a--, b += 2 * step)
	MPZ_MUL_1_POS (ret, ret, b * (b + step));
      if (count % 2 == 1)
	MPZ_MUL_1_POS (ret, ret, b);
      return;
    case 3:
      if (count == 2)
	{
	  MPZ_SET_1_NZ (ret, start * (start + step));
	  return;
	}
      MPZ_SET_1_NZ (ret, start * (start + step) * (start + 2 * step));
      if (count == 3)
	return;
      for (b = start + 3 * step, a = count / 3 - 1; a != 0;
	   a--, b += 3 * step)
	MPZ_MUL_1_POS (ret, ret, b * (b + step) * (b + 2 * step));
      if (count % 3 == 2)
	b = b * (b + step);
      if (count % 3 != 0)
	MPZ_MUL_1_POS (ret, ret, b);
      return;
    default:			/* ie nm=4      */
      if (count == 2)
	{
	  MPZ_SET_1_NZ (ret, start * (start + step));
	  return;
	}
      if (count == 3)
	{
	  MPZ_SET_1_NZ (ret, start * (start + step) * (start + 2 * step));
	  return;
	}
      MPZ_SET_1_NZ (ret,
		    start * (start + step) * (start + 2 * step) * (start +
								   3 * step));
      if (count == 4)
	return;
      for (b = start + 4 * step, a = count / 4 - 1; a != 0;
	   a--, b += 4 * step)
	MPZ_MUL_1_POS (ret, ret,
		       b * (b + step) * (b + 2 * step) * (b + 3 * step));
      if (count % 4 == 2)
	b = b * (b + step);
      if (count % 4 == 3)
	b = b * (b + step) * (b + 2 * step);
      if (count % 4 != 0)
	MPZ_MUL_1_POS (ret, ret, b);
      return;
    }
}

/* return value in st[0]
   odd_product(l,h)=sqrt((h/e)^h/(l/e)^l) using Stirling approx and e=exp(1)
   so st[0] needs enough bits for above, st[1] needs half these bits and
   st[2] needs 1/4 of these bits etc	*/
static void
odd_product (unsigned long low, unsigned long high, mpz_t * st)
{
  unsigned long stc = 1, stn = 0, n, y, mask, a, nm = 1;
  signed long z;

  low++;
  if (low % 2 == 0)
    low++;
  if (high == 0)
    high = 1;
  if (high % 2 == 0)
    high--;
/* must have high>=low ? check this and remove test below */
  if (high < low)
    {
      MPZ_SET_1_NZ (st[0], 1);
      return;
    }
  if (high == low)
    {
      MPZ_SET_1_NZ (st[0], low);
      return;
    }
  if (high <= FACMUL2 + 2)
    {
      nm = 2;
      if (high <= FACMUL3 + 4)
	{
	  nm = 3;
	  if (high <= FACMUL4 + 6)
	    nm = 4;
	}
    }
  high = (high - low) / 2 + 1;	/* high is now count,high<=2^(BITS_PER_ULONG-1) */
  if (high <= (((unsigned long) 1) << APCONST))
    {
      ap_product_small (st[0], (mp_limb_t) low, CNST_LIMB(2), high, nm);
      return;
    }
  count_leading_zeros (n, (mp_limb_t) high);
/* assumes clz above is LIMB based not NUMB based */
  n = GMP_LIMB_BITS - n - APCONST;
  mask = (((unsigned long) 1) << n);
  a = mask << 1;
  mask--;
/* have 2^(BITS_IN_N-APCONST) iterations so need
   (BITS_IN_N-APCONST+1) stack entries	*/
  for (z = mask; z >= 0; z--)
    {
      BITREV_ULONG (y, z);
      y >>= (BITS_PER_ULONG - n);
      ap_product_small (st[stn],
			(mp_limb_t) (low + 2 * ((~y) & mask)), (mp_limb_t) a,
			(high + y) >> n, nm);
      ASSERT (((high + y) >> n) <= (((unsigned long) 1) << APCONST));
      stn++;
      y = stc++;
      while ((y & 1) == 0)
	{
	  mpz_mul (st[stn - 2], st[stn - 2], st[stn - 1]);
	  stn--;
	  y >>= 1;
	}
    }
  ASSERT (stn == 1);
  return;
}