new mpz_fac_ui code
Kevin Ryde
user42@zip.com.au
Thu, 19 Jun 2003 10:04:50 +1000
--=-=-=
David Newman <david.newman@jesus.ox.ac.uk> writes:
>
> Here's some factorial code which uses the factorization of n! and repeated
> squaring.
Jason Moxham has contributed a new mpz_fac_ui (below) based on
stripping factors of two, and splitting the remaining products.
> It runs about twice as fast as the original code for large enough n.
> Unfortunately, 'large enough' is about 15000 on my Athlon 1.2G so a bit of
> threshold handling is needed.
That's rather big, perhaps it can be brought down.
I'd wondered if a full prime factorization would be necessary, perhaps
just stripping factors of for instance 2, 3, 5, 7 and 11 would result
in most of the possible squaring.
For reference, mpz_fac_ui is not a very important function, but good
algorithms well implemented are always welcome.
Factorial
---------
Factorials are calculated by a combination of removal of twos,
powering, and binary splitting. The procedure can be best illustrated
with an example,
23! = 1.2.3.4.5.6.7.8.9.10.11.12.13.14.15.16.17.18.19.20.21.22.23
has factors of two removed,
23! = 2^19.1.1.3.1.5.3.7.1.9.5.11.3.13.7.15.1.17.9.19.5.21.11.23
and the resulting terms collected up according to their multiplicity,
23! = 2^19.(3.5)^3.(7.9.11)^2.(13.15.17.19.21.23)
Each sequence such as 13.15.17.19.21.23 is evaluated by splitting
into every second term, as for instance (13.17.21).(15.19.23), and the
same recursively on each half. This is implemented iteratively using
some bit twiddling.
Such splitting is more efficient than repeated Nx1 multiplies since
it forms big multiplies, allowing Karatsuba and higher algorithms to be
used. And even below the Karatsuba threshold a big block of work can
be more efficient for the basecase algorithm.
Splitting into subsequences of every second term keeps the resulting
products more nearly equal in size than would the simpler approach of
say taking the first half and second half of the sequence. Nearly
equal products are more efficient for the current multiply
implementation.
--=-=-=
Content-Type: text/x-csrc
Content-Disposition: attachment; filename=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 2.1 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; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
#include "fac_ui.h"
static void odd_product _PROTO ((unsigned long low, unsigned long high, mpz_t * st));
static void ap_product_small _PROTO ((mpz_t ret, mp_limb_t start, mp_limb_t step, unsigned long count, unsigned long nm));
/* 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==BITS_PER_MP_LIMB
#define BSWAP_ULONG(x,y) BSWAP_LIMB(x,y)
#endif
#if BITS_PER_ULONG*2==BITS_PER_MP_LIMB
#define BSWAP_ULONG(x,y) \
do{mp_limb_t __src=(y); \
BSWAP_LIMB(__src,__src); \
__src>>=BITS_PER_ULONG; \
(x)=__src; \
}while(0)
#endif
#if ! defined (BSWAP_ULONG)
#define BSWAP_ULONG(dst, src) \
do { \
unsigned long __bswapl_src = (src); \
unsigned long __dst = 0; \
int __i; \
for (__i = 0; __i < sizeof(unsigned long); __i++) \
{ \
__dst = (__dst << 8) | (__bswapl_src & 0xFF); \
__bswapl_src >>= 8; \
} \
(dst) = __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 formulars are only approxiamtions , also we ignore
the few "side" calculuations , corect allocation seems to speed up the
small sizes better , having very little effect on the large sizes */
/* estimate space for stack entrys 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, 2, 1, n - 1, 4);
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, 2, 1, n / 2 - 1, 4);
ap_product_small (t1, 3, 2, (n - 1) / 2, 4);
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 entrys */
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, 2, 1, n / 4 - 1, 4);
/* 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], low, 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 entrys */
for (z = mask; z >= 0; z--)
{
BITREV_ULONG (y, z);
y >>= (BITS_PER_ULONG - n);
ap_product_small (st[stn], low + 2 * ((~y) & mask), 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;
}
--=-=-=
Content-Type: text/x-chdr
Content-Disposition: attachment; filename=fac_ui.h
/* This file is automatically generated by gen-fac_ui.c */
#if GMP_NUMB_BITS != 32
Error , error this data is for 32 GMP_NUMB_BITS only
#endif
#if GMP_LIMB_BITS != 32
Error , error this data is for 32 GMP_LIMB_BITS only
#endif
/* This table is 0!,1!,2!,3!,...,n! where n! has <= GMP_NUMB_BITS bits */
#define ONE_LIMB_FACTORIAL_TABLE CNST_LIMB(0x1),CNST_LIMB(0x1),CNST_LIMB(0x2),CNST_LIMB(0x6),CNST_LIMB(0x18),CNST_LIMB(0x78),CNST_LIMB(0x2d0),CNST_LIMB(0x13b0),CNST_LIMB(0x9d80),CNST_LIMB(0x58980),CNST_LIMB(0x375f00),CNST_LIMB(0x2611500),CNST_LIMB(0x1c8cfc00)
/* is 2^(GMP_LIMB_BITS+1)/exp(1) */
#define FAC2OVERE CNST_LIMB(0xbc5c254b)
/* FACMULn is largest odd x such that x*(x+2)*...*(x+2(n-1))<=2^GMP_NUMB_BITS-1 */
#define FACMUL2 CNST_LIMB(0xffff)
#define FACMUL3 CNST_LIMB(0x657)
#define FACMUL4 CNST_LIMB(0xfd)
--=-=-=--