1  /* mpz_mfac_uiui(RESULT, N, M) -- Set RESULT to N!^(M) = N(N-M)(N-2M)...
       2  
       3  Contributed to the GNU project by Marco Bodrato.
       4  
       5  Copyright 2012, 2013, 2015, 2016 Free Software Foundation, Inc.
       6  
       7  This file is part of the GNU MP Library.
       8  
       9  The GNU MP Library is free software; you can redistribute it and/or modify
      10  it under the terms of either:
      11  
      12    * the GNU Lesser General Public License as published by the Free
      13      Software Foundation; either version 3 of the License, or (at your
      14      option) any later version.
      15  
      16  or
      17  
      18    * the GNU General Public License as published by the Free Software
      19      Foundation; either version 2 of the License, or (at your option) any
      20      later version.
      21  
      22  or both in parallel, as here.
      23  
      24  The GNU MP Library is distributed in the hope that it will be useful, but
      25  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
      26  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
      27  for more details.
      28  
      29  You should have received copies of the GNU General Public License and the
      30  GNU Lesser General Public License along with the GNU MP Library.  If not,
      31  see https://www.gnu.org/licenses/.  */
      32  
      33  #include "gmp-impl.h"
      34  
      35  /*************************************************************/
      36  /* Section macros: common macros, for swing/fac/bin (&sieve) */
      37  /*************************************************************/
      38  
      39  #define FACTOR_LIST_STORE(P, PR, MAX_PR, VEC, I)		\
      40    do {								\
      41      if ((PR) > (MAX_PR)) {					\
      42        (VEC)[(I)++] = (PR);					\
      43        (PR) = (P);						\
      44      } else							\
      45        (PR) *= (P);						\
      46    } while (0)
      47  
      48  /*********************************************************/
      49  /* Section oder factorials:                              */
      50  /*********************************************************/
      51  
      52  /* mpz_mfac_uiui (x, n, m) computes x = n!^(m) = n*(n-m)*(n-2m)*...   */
      53  
      54  void
      55  mpz_mfac_uiui (mpz_ptr x, unsigned long n, unsigned long m)
      56  {
      57    ASSERT (n <= GMP_NUMB_MAX);
      58    ASSERT (m != 0);
      59  
      60    if ((n < 3) | (n - 3 < m - 1)) { /* (n < 3 || n - 1 <= m || m == 0) */
      61      MPZ_NEWALLOC (x, 1)[0] = n + (n == 0);
      62      SIZ (x) = 1;
      63    } else { /* 0 < m < n - 1 < GMP_NUMB_MAX */
      64      mp_limb_t g, sn;
      65      mpz_t     t;
      66  
      67      sn = n;
      68      g = mpn_gcd_1 (&sn, 1, m);
      69      if (g > 1) { n/=g; m/=g; }
      70  
      71      if (m <= 2) { /* fac or 2fac */
      72        if (m == 1) {
      73  	if (g > 2) {
      74  	  mpz_init (t);
      75  	  mpz_fac_ui (t, n);
      76  	  sn = n;
      77  	} else {
      78  	  if (g == 2)
      79  	    mpz_2fac_ui (x, n << 1);
      80  	  else
      81  	    mpz_fac_ui (x, n);
      82  	  return;
      83  	}
      84        } else { /* m == 2 */
      85  	if (g > 1) {
      86  	  mpz_init (t);
      87  	  mpz_2fac_ui (t, n);
      88  	  sn = n / 2 + 1;
      89  	} else {
      90  	  mpz_2fac_ui (x, n);
      91  	  return;
      92  	}
      93        }
      94      } else { /* m >= 3, gcd(n,m) = 1 */
      95        mp_limb_t *factors;
      96        mp_limb_t prod, max_prod;
      97        mp_size_t j;
      98        TMP_DECL;
      99  
     100        sn = n / m + 1;
     101  
     102        j = 0;
     103        prod = n;
     104        n -= m;
     105        max_prod = GMP_NUMB_MAX / n;
     106  
     107        if (g > 1)
     108  	factors = MPZ_NEWALLOC (x, sn / log_n_max (n) + 2);
     109        else {
     110  	TMP_MARK;
     111  	factors = TMP_ALLOC_LIMBS (sn / log_n_max (n) + 2);
     112        }
     113  
     114        for (; n > m; n -= m)
     115  	FACTOR_LIST_STORE (n, prod, max_prod, factors, j);
     116  
     117        factors[j++] = n;
     118        factors[j++] = prod;
     119  
     120        if (g > 1) {
     121  	mpz_init (t);
     122  	mpz_prodlimbs (t, factors, j);
     123        } else {
     124  	mpz_prodlimbs (x, factors, j);
     125  	TMP_FREE;
     126  	return;
     127        }
     128      }
     129  
     130      {
     131        mpz_t p;
     132  
     133        mpz_init (p);
     134        mpz_ui_pow_ui (p, g, sn); /* g^sn */
     135        mpz_mul (x, p, t);
     136        mpz_clear (p);
     137        mpz_clear (t);
     138      }
     139    }
     140  }