|
Packit |
5c3484 |
/* mpz_oddfac_1(RESULT, N) -- Set RESULT to the odd factor of N!.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
Contributed to the GNU project by Marco Bodrato.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.
|
|
Packit |
5c3484 |
IT IS ONLY SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.
|
|
Packit |
5c3484 |
IN FACT, IT IS ALMOST GUARANTEED THAT IT WILL CHANGE OR
|
|
Packit |
5c3484 |
DISAPPEAR IN A FUTURE GNU MP RELEASE.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
Copyright 2010-2012 Free Software Foundation, Inc.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
This file is part of the GNU MP Library.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The GNU MP Library is free software; you can redistribute it and/or modify
|
|
Packit |
5c3484 |
it under the terms of either:
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
* the GNU Lesser General Public License as published by the Free
|
|
Packit |
5c3484 |
Software Foundation; either version 3 of the License, or (at your
|
|
Packit |
5c3484 |
option) any later version.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
or
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
* the GNU General Public License as published by the Free Software
|
|
Packit |
5c3484 |
Foundation; either version 2 of the License, or (at your option) any
|
|
Packit |
5c3484 |
later version.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
or both in parallel, as here.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The GNU MP Library is distributed in the hope that it will be useful, but
|
|
Packit |
5c3484 |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
|
|
Packit |
5c3484 |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
Packit |
5c3484 |
for more details.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
You should have received copies of the GNU General Public License and the
|
|
Packit |
5c3484 |
GNU Lesser General Public License along with the GNU MP Library. If not,
|
|
Packit |
5c3484 |
see https://www.gnu.org/licenses/. */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#include "gmp.h"
|
|
Packit |
5c3484 |
#include "gmp-impl.h"
|
|
Packit |
5c3484 |
#include "longlong.h"
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* TODO:
|
|
Packit |
5c3484 |
- split this file in smaller parts with functions that can be recycled for different computations.
|
|
Packit |
5c3484 |
*/
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/**************************************************************/
|
|
Packit |
5c3484 |
/* Section macros: common macros, for mswing/fac/bin (&sieve) */
|
|
Packit |
5c3484 |
/**************************************************************/
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#define FACTOR_LIST_APPEND(PR, MAX_PR, VEC, I) \
|
|
Packit |
5c3484 |
if ((PR) > (MAX_PR)) { \
|
|
Packit |
5c3484 |
(VEC)[(I)++] = (PR); \
|
|
Packit |
5c3484 |
(PR) = 1; \
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#define FACTOR_LIST_STORE(P, PR, MAX_PR, VEC, I) \
|
|
Packit |
5c3484 |
do { \
|
|
Packit |
5c3484 |
if ((PR) > (MAX_PR)) { \
|
|
Packit |
5c3484 |
(VEC)[(I)++] = (PR); \
|
|
Packit |
5c3484 |
(PR) = (P); \
|
|
Packit |
5c3484 |
} else \
|
|
Packit |
5c3484 |
(PR) *= (P); \
|
|
Packit |
5c3484 |
} while (0)
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#define LOOP_ON_SIEVE_CONTINUE(prime,end,sieve) \
|
|
Packit |
5c3484 |
__max_i = (end); \
|
|
Packit |
5c3484 |
\
|
|
Packit |
5c3484 |
do { \
|
|
Packit |
5c3484 |
++__i; \
|
|
Packit |
5c3484 |
if (((sieve)[__index] & __mask) == 0) \
|
|
Packit |
5c3484 |
{ \
|
|
Packit |
5c3484 |
(prime) = id_to_n(__i)
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#define LOOP_ON_SIEVE_BEGIN(prime,start,end,off,sieve) \
|
|
Packit |
5c3484 |
do { \
|
|
Packit |
5c3484 |
mp_limb_t __mask, __index, __max_i, __i; \
|
|
Packit |
5c3484 |
\
|
|
Packit |
5c3484 |
__i = (start)-(off); \
|
|
Packit |
5c3484 |
__index = __i / GMP_LIMB_BITS; \
|
|
Packit |
5c3484 |
__mask = CNST_LIMB(1) << (__i % GMP_LIMB_BITS); \
|
|
Packit |
5c3484 |
__i += (off); \
|
|
Packit |
5c3484 |
\
|
|
Packit |
5c3484 |
LOOP_ON_SIEVE_CONTINUE(prime,end,sieve)
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#define LOOP_ON_SIEVE_STOP \
|
|
Packit |
5c3484 |
} \
|
|
Packit |
5c3484 |
__mask = __mask << 1 | __mask >> (GMP_LIMB_BITS-1); \
|
|
Packit |
5c3484 |
__index += __mask & 1; \
|
|
Packit |
5c3484 |
} while (__i <= __max_i) \
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#define LOOP_ON_SIEVE_END \
|
|
Packit |
5c3484 |
LOOP_ON_SIEVE_STOP; \
|
|
Packit |
5c3484 |
} while (0)
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/*********************************************************/
|
|
Packit |
5c3484 |
/* Section sieve: sieving functions and tools for primes */
|
|
Packit |
5c3484 |
/*********************************************************/
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#if WANT_ASSERT
|
|
Packit |
5c3484 |
static mp_limb_t
|
|
Packit |
5c3484 |
bit_to_n (mp_limb_t bit) { return (bit*3+4)|1; }
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* id_to_n (x) = bit_to_n (x-1) = (id*3+1)|1*/
|
|
Packit |
5c3484 |
static mp_limb_t
|
|
Packit |
5c3484 |
id_to_n (mp_limb_t id) { return id*3+1+(id&1;; }
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* n_to_bit (n) = ((n-1)&(-CNST_LIMB(2)))/3U-1 */
|
|
Packit |
5c3484 |
static mp_limb_t
|
|
Packit |
5c3484 |
n_to_bit (mp_limb_t n) { return ((n-5)|1)/3U; }
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#if WANT_ASSERT
|
|
Packit |
5c3484 |
static mp_size_t
|
|
Packit |
5c3484 |
primesieve_size (mp_limb_t n) { return n_to_bit(n) / GMP_LIMB_BITS + 1; }
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/*********************************************************/
|
|
Packit |
5c3484 |
/* Section mswing: 2-multiswing factorial */
|
|
Packit |
5c3484 |
/*********************************************************/
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Returns an approximation of the sqare root of x. *
|
|
Packit |
5c3484 |
* It gives: x <= limb_apprsqrt (x) ^ 2 < x * 9/4 */
|
|
Packit |
5c3484 |
static mp_limb_t
|
|
Packit |
5c3484 |
limb_apprsqrt (mp_limb_t x)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
int s;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
ASSERT (x > 2);
|
|
Packit |
5c3484 |
count_leading_zeros (s, x - 1);
|
|
Packit |
5c3484 |
s = GMP_LIMB_BITS - 1 - s;
|
|
Packit |
5c3484 |
return (CNST_LIMB(1) << (s >> 1)) + (CNST_LIMB(1) << ((s - 1) >> 1));
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#if 0
|
|
Packit |
5c3484 |
/* A count-then-exponentiate variant for SWING_A_PRIME */
|
|
Packit |
5c3484 |
#define SWING_A_PRIME(P, N, PR, MAX_PR, VEC, I) \
|
|
Packit |
5c3484 |
do { \
|
|
Packit |
5c3484 |
mp_limb_t __q, __prime; \
|
|
Packit |
5c3484 |
int __exp; \
|
|
Packit |
5c3484 |
__prime = (P); \
|
|
Packit |
5c3484 |
__exp = 0; \
|
|
Packit |
5c3484 |
__q = (N); \
|
|
Packit |
5c3484 |
do { \
|
|
Packit |
5c3484 |
__q /= __prime; \
|
|
Packit |
5c3484 |
__exp += __q & 1; \
|
|
Packit |
5c3484 |
} while (__q >= __prime); \
|
|
Packit |
5c3484 |
if (__exp) { /* Store $prime^{exp}$ */ \
|
|
Packit |
5c3484 |
for (__q = __prime; --__exp; __q *= __prime); \
|
|
Packit |
5c3484 |
FACTOR_LIST_STORE(__q, PR, MAX_PR, VEC, I); \
|
|
Packit |
5c3484 |
}; \
|
|
Packit |
5c3484 |
} while (0)
|
|
Packit |
5c3484 |
#else
|
|
Packit |
5c3484 |
#define SWING_A_PRIME(P, N, PR, MAX_PR, VEC, I) \
|
|
Packit |
5c3484 |
do { \
|
|
Packit |
5c3484 |
mp_limb_t __q, __prime; \
|
|
Packit |
5c3484 |
__prime = (P); \
|
|
Packit |
5c3484 |
FACTOR_LIST_APPEND(PR, MAX_PR, VEC, I); \
|
|
Packit |
5c3484 |
__q = (N); \
|
|
Packit |
5c3484 |
do { \
|
|
Packit |
5c3484 |
__q /= __prime; \
|
|
Packit |
5c3484 |
if ((__q & 1) != 0) (PR) *= __prime; \
|
|
Packit |
5c3484 |
} while (__q >= __prime); \
|
|
Packit |
5c3484 |
} while (0)
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#define SH_SWING_A_PRIME(P, N, PR, MAX_PR, VEC, I) \
|
|
Packit |
5c3484 |
do { \
|
|
Packit |
5c3484 |
mp_limb_t __prime; \
|
|
Packit |
5c3484 |
__prime = (P); \
|
|
Packit |
5c3484 |
if ((((N) / __prime) & 1) != 0) \
|
|
Packit |
5c3484 |
FACTOR_LIST_STORE(__prime, PR, MAX_PR, VEC, I); \
|
|
Packit |
5c3484 |
} while (0)
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* mpz_2multiswing_1 computes the odd part of the 2-multiswing
|
|
Packit |
5c3484 |
factorial of the parameter n. The result x is an odd positive
|
|
Packit |
5c3484 |
integer so that multiswing(n,2) = x 2^a.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
Uses the algorithm described by Peter Luschny in "Divide, Swing and
|
|
Packit |
5c3484 |
Conquer the Factorial!".
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The pointer sieve points to primesieve_size(n) limbs containing a
|
|
Packit |
5c3484 |
bit-array where primes are marked as 0.
|
|
Packit |
5c3484 |
Enough (FIXME: explain :-) limbs must be pointed by factors.
|
|
Packit |
5c3484 |
*/
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
static void
|
|
Packit |
5c3484 |
mpz_2multiswing_1 (mpz_ptr x, mp_limb_t n, mp_ptr sieve, mp_ptr factors)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mp_limb_t prod, max_prod;
|
|
Packit |
5c3484 |
mp_size_t j;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
ASSERT (n >= 26);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
j = 0;
|
|
Packit |
5c3484 |
prod = -(n & 1);
|
|
Packit |
5c3484 |
n &= ~ CNST_LIMB(1); /* n-1, if n is odd */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
prod = (prod & n) + 1; /* the original n, if it was odd, 1 otherwise */
|
|
Packit |
5c3484 |
max_prod = GMP_NUMB_MAX / (n-1);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Handle prime = 3 separately. */
|
|
Packit |
5c3484 |
SWING_A_PRIME (3, n, prod, max_prod, factors, j);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Swing primes from 5 to n/3 */
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mp_limb_t s;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mp_limb_t prime;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
s = limb_apprsqrt(n);
|
|
Packit |
5c3484 |
ASSERT (s >= 5);
|
|
Packit |
5c3484 |
s = n_to_bit (s);
|
|
Packit |
5c3484 |
LOOP_ON_SIEVE_BEGIN (prime, n_to_bit (5), s, 0,sieve);
|
|
Packit |
5c3484 |
SWING_A_PRIME (prime, n, prod, max_prod, factors, j);
|
|
Packit |
5c3484 |
LOOP_ON_SIEVE_END;
|
|
Packit |
5c3484 |
s++;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
ASSERT (max_prod <= GMP_NUMB_MAX / 3);
|
|
Packit |
5c3484 |
ASSERT (bit_to_n (s) * bit_to_n (s) > n);
|
|
Packit |
5c3484 |
ASSERT (s <= n_to_bit (n / 3));
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mp_limb_t prime;
|
|
Packit |
5c3484 |
mp_limb_t l_max_prod = max_prod * 3;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
LOOP_ON_SIEVE_BEGIN (prime, s, n_to_bit (n/3), 0, sieve);
|
|
Packit |
5c3484 |
SH_SWING_A_PRIME (prime, n, prod, l_max_prod, factors, j);
|
|
Packit |
5c3484 |
LOOP_ON_SIEVE_END;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Store primes from (n+1)/2 to n */
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mp_limb_t prime;
|
|
Packit |
5c3484 |
LOOP_ON_SIEVE_BEGIN (prime, n_to_bit (n >> 1) + 1, n_to_bit (n), 0,sieve);
|
|
Packit |
5c3484 |
FACTOR_LIST_STORE (prime, prod, max_prod, factors, j);
|
|
Packit |
5c3484 |
LOOP_ON_SIEVE_END;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (LIKELY (j != 0))
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
factors[j++] = prod;
|
|
Packit |
5c3484 |
mpz_prodlimbs (x, factors, j);
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
PTR (x)[0] = prod;
|
|
Packit |
5c3484 |
SIZ (x) = 1;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#undef SWING_A_PRIME
|
|
Packit |
5c3484 |
#undef SH_SWING_A_PRIME
|
|
Packit |
5c3484 |
#undef LOOP_ON_SIEVE_END
|
|
Packit |
5c3484 |
#undef LOOP_ON_SIEVE_STOP
|
|
Packit |
5c3484 |
#undef LOOP_ON_SIEVE_BEGIN
|
|
Packit |
5c3484 |
#undef LOOP_ON_SIEVE_CONTINUE
|
|
Packit |
5c3484 |
#undef FACTOR_LIST_APPEND
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/*********************************************************/
|
|
Packit |
5c3484 |
/* Section oddfac: odd factorial, needed also by binomial*/
|
|
Packit |
5c3484 |
/*********************************************************/
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#if TUNE_PROGRAM_BUILD
|
|
Packit |
5c3484 |
#define FACTORS_PER_LIMB (GMP_NUMB_BITS / (LOG2C(FAC_DSC_THRESHOLD_LIMIT-1)+1))
|
|
Packit |
5c3484 |
#else
|
|
Packit |
5c3484 |
#define FACTORS_PER_LIMB (GMP_NUMB_BITS / (LOG2C(FAC_DSC_THRESHOLD-1)+1))
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* mpz_oddfac_1 computes the odd part of the factorial of the
|
|
Packit |
5c3484 |
parameter n. I.e. n! = x 2^a, where x is the returned value: an
|
|
Packit |
5c3484 |
odd positive integer.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
If flag != 0 a square is skipped in the DSC part, e.g.
|
|
Packit |
5c3484 |
if n is odd, n > FAC_DSC_THRESHOLD and flag = 1, x is set to n!!.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
If n is too small, flag is ignored, and an ASSERT can be triggered.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
TODO: FAC_DSC_THRESHOLD is used here with two different roles:
|
|
Packit |
5c3484 |
- to decide when prime factorisation is needed,
|
|
Packit |
5c3484 |
- to stop the recursion, once sieving is done.
|
|
Packit |
5c3484 |
Maybe two thresholds can do a better job.
|
|
Packit |
5c3484 |
*/
|
|
Packit |
5c3484 |
void
|
|
Packit |
5c3484 |
mpz_oddfac_1 (mpz_ptr x, mp_limb_t n, unsigned flag)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
ASSERT (n <= GMP_NUMB_MAX);
|
|
Packit |
5c3484 |
ASSERT (flag == 0 || (flag == 1 && n > ODD_FACTORIAL_TABLE_LIMIT && ABOVE_THRESHOLD (n, FAC_DSC_THRESHOLD)));
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (n <= ODD_FACTORIAL_TABLE_LIMIT)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
PTR (x)[0] = __gmp_oddfac_table[n];
|
|
Packit |
5c3484 |
SIZ (x) = 1;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else if (n <= ODD_DOUBLEFACTORIAL_TABLE_LIMIT + 1)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mp_ptr px;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
px = MPZ_NEWALLOC (x, 2);
|
|
Packit |
5c3484 |
umul_ppmm (px[1], px[0], __gmp_odd2fac_table[(n - 1) >> 1], __gmp_oddfac_table[n >> 1]);
|
|
Packit |
5c3484 |
SIZ (x) = 2;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
unsigned s;
|
|
Packit |
5c3484 |
mp_ptr factors;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
s = 0;
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mp_limb_t tn;
|
|
Packit |
5c3484 |
mp_limb_t prod, max_prod, i;
|
|
Packit |
5c3484 |
mp_size_t j;
|
|
Packit |
5c3484 |
TMP_SDECL;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#if TUNE_PROGRAM_BUILD
|
|
Packit |
5c3484 |
ASSERT (FAC_DSC_THRESHOLD_LIMIT >= FAC_DSC_THRESHOLD);
|
|
Packit |
5c3484 |
ASSERT (FAC_DSC_THRESHOLD >= 2 * (ODD_DOUBLEFACTORIAL_TABLE_LIMIT + 2));
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Compute the number of recursive steps for the DSC algorithm. */
|
|
Packit |
5c3484 |
for (tn = n; ABOVE_THRESHOLD (tn, FAC_DSC_THRESHOLD); s++)
|
|
Packit |
5c3484 |
tn >>= 1;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
j = 0;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
TMP_SMARK;
|
|
Packit |
5c3484 |
factors = TMP_SALLOC_LIMBS (1 + tn / FACTORS_PER_LIMB);
|
|
Packit |
5c3484 |
ASSERT (tn >= FACTORS_PER_LIMB);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
prod = 1;
|
|
Packit |
5c3484 |
#if TUNE_PROGRAM_BUILD
|
|
Packit |
5c3484 |
max_prod = GMP_NUMB_MAX / FAC_DSC_THRESHOLD_LIMIT;
|
|
Packit |
5c3484 |
#else
|
|
Packit |
5c3484 |
max_prod = GMP_NUMB_MAX / FAC_DSC_THRESHOLD;
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
ASSERT (tn > ODD_DOUBLEFACTORIAL_TABLE_LIMIT + 1);
|
|
Packit |
5c3484 |
do {
|
|
Packit |
5c3484 |
i = ODD_DOUBLEFACTORIAL_TABLE_LIMIT + 2;
|
|
Packit |
5c3484 |
factors[j++] = ODD_DOUBLEFACTORIAL_TABLE_MAX;
|
|
Packit |
5c3484 |
do {
|
|
Packit |
5c3484 |
FACTOR_LIST_STORE (i, prod, max_prod, factors, j);
|
|
Packit |
5c3484 |
i += 2;
|
|
Packit |
5c3484 |
} while (i <= tn);
|
|
Packit |
5c3484 |
max_prod <<= 1;
|
|
Packit |
5c3484 |
tn >>= 1;
|
|
Packit |
5c3484 |
} while (tn > ODD_DOUBLEFACTORIAL_TABLE_LIMIT + 1);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
factors[j++] = prod;
|
|
Packit |
5c3484 |
factors[j++] = __gmp_odd2fac_table[(tn - 1) >> 1];
|
|
Packit |
5c3484 |
factors[j++] = __gmp_oddfac_table[tn >> 1];
|
|
Packit |
5c3484 |
mpz_prodlimbs (x, factors, j);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
TMP_SFREE;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (s != 0)
|
|
Packit |
5c3484 |
/* Use the algorithm described by Peter Luschny in "Divide,
|
|
Packit |
5c3484 |
Swing and Conquer the Factorial!".
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
Improvement: there are two temporary buffers, factors and
|
|
Packit |
5c3484 |
square, that are never used together; with a good estimate
|
|
Packit |
5c3484 |
of the maximal needed size, they could share a single
|
|
Packit |
5c3484 |
allocation.
|
|
Packit |
5c3484 |
*/
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mpz_t mswing;
|
|
Packit |
5c3484 |
mp_ptr sieve;
|
|
Packit |
5c3484 |
mp_size_t size;
|
|
Packit |
5c3484 |
TMP_DECL;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
TMP_MARK;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
flag--;
|
|
Packit |
5c3484 |
size = n / GMP_NUMB_BITS + 4;
|
|
Packit |
5c3484 |
ASSERT (primesieve_size (n - 1) <= size - (size / 2 + 1));
|
|
Packit |
5c3484 |
/* 2-multiswing(n) < 2^(n-1)*sqrt(n/pi) < 2^(n+GMP_NUMB_BITS);
|
|
Packit |
5c3484 |
one more can be overwritten by mul, another for the sieve */
|
|
Packit |
5c3484 |
MPZ_TMP_INIT (mswing, size);
|
|
Packit |
5c3484 |
/* Initialize size, so that ASSERT can check it correctly. */
|
|
Packit |
5c3484 |
ASSERT_CODE (SIZ (mswing) = 0);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Put the sieve on the second half, it will be overwritten by the last mswing. */
|
|
Packit |
5c3484 |
sieve = PTR (mswing) + size / 2 + 1;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
size = (gmp_primesieve (sieve, n - 1) + 1) / log_n_max (n) + 1;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
factors = TMP_ALLOC_LIMBS (size);
|
|
Packit |
5c3484 |
do {
|
|
Packit |
5c3484 |
mp_ptr square, px;
|
|
Packit |
5c3484 |
mp_size_t nx, ns;
|
|
Packit |
5c3484 |
mp_limb_t cy;
|
|
Packit |
5c3484 |
TMP_DECL;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
s--;
|
|
Packit |
5c3484 |
ASSERT (ABSIZ (mswing) < ALLOC (mswing) / 2); /* Check: sieve has not been overwritten */
|
|
Packit |
5c3484 |
mpz_2multiswing_1 (mswing, n >> s, sieve, factors);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
TMP_MARK;
|
|
Packit |
5c3484 |
nx = SIZ (x);
|
|
Packit |
5c3484 |
if (s == flag) {
|
|
Packit |
5c3484 |
size = nx;
|
|
Packit |
5c3484 |
square = TMP_ALLOC_LIMBS (size);
|
|
Packit |
5c3484 |
MPN_COPY (square, PTR (x), nx);
|
|
Packit |
5c3484 |
} else {
|
|
Packit |
5c3484 |
size = nx << 1;
|
|
Packit |
5c3484 |
square = TMP_ALLOC_LIMBS (size);
|
|
Packit |
5c3484 |
mpn_sqr (square, PTR (x), nx);
|
|
Packit |
5c3484 |
size -= (square[size - 1] == 0);
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
ns = SIZ (mswing);
|
|
Packit |
5c3484 |
nx = size + ns;
|
|
Packit |
5c3484 |
px = MPZ_NEWALLOC (x, nx);
|
|
Packit |
5c3484 |
ASSERT (ns <= size);
|
|
Packit |
5c3484 |
cy = mpn_mul (px, square, size, PTR(mswing), ns); /* n!= n$ * floor(n/2)!^2 */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
TMP_FREE;
|
|
Packit |
5c3484 |
SIZ(x) = nx - (cy == 0);
|
|
Packit |
5c3484 |
} while (s != 0);
|
|
Packit |
5c3484 |
TMP_FREE;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#undef FACTORS_PER_LIMB
|
|
Packit |
5c3484 |
#undef FACTOR_LIST_STORE
|