github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/mpz/mfac_uiui.c (about)

     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 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.h"
    34  #include "gmp-impl.h"
    35  
    36  /*************************************************************/
    37  /* Section macros: common macros, for swing/fac/bin (&sieve) */
    38  /*************************************************************/
    39  
    40  #define FACTOR_LIST_STORE(P, PR, MAX_PR, VEC, I)		\
    41    do {								\
    42      if ((PR) > (MAX_PR)) {					\
    43        (VEC)[(I)++] = (PR);					\
    44        (PR) = (P);						\
    45      } else							\
    46        (PR) *= (P);						\
    47    } while (0)
    48  
    49  /*********************************************************/
    50  /* Section oder factorials:                              */
    51  /*********************************************************/
    52  
    53  /* mpz_mfac_uiui (x, n, m) computes x = n!^(m) = n*(n-m)*(n-2m)*...   */
    54  
    55  void
    56  mpz_mfac_uiui (mpz_ptr x, unsigned long n, unsigned long m)
    57  {
    58    ASSERT (n <= GMP_NUMB_MAX);
    59    ASSERT (m != 0);
    60  
    61    if ((n < 3) | (n - 3 < m - 1)) { /* (n < 3 || n - 1 <= m || m == 0) */
    62      PTR (x)[0] = n + (n == 0);
    63      SIZ (x) = 1;
    64    } else { /* m < n - 1 < GMP_NUMB_MAX */
    65      mp_limb_t g, sn;
    66      mpz_t     t;
    67  
    68      sn = n;
    69      g = mpn_gcd_1 (&sn, 1, m);
    70      if (g > 1) { n/=g; m/=g; }
    71  
    72      if (m <= 2) { /* fac or 2fac */
    73        if (m == 1) {
    74  	if (g > 2) {
    75  	  mpz_init (t);
    76  	  mpz_fac_ui (t, n);
    77  	  sn = n;
    78  	} else {
    79  	  if (g == 2)
    80  	    mpz_2fac_ui (x, n << 1);
    81  	  else
    82  	    mpz_fac_ui (x, n);
    83  	  return;
    84  	}
    85        } else { /* m == 2 */
    86  	if (g > 1) {
    87  	  mpz_init (t);
    88  	  mpz_2fac_ui (t, n);
    89  	  sn = n / 2 + 1;
    90  	} else {
    91  	  mpz_2fac_ui (x, n);
    92  	  return;
    93  	}
    94        }
    95      } else { /* m >= 3, gcd(n,m) = 1 */
    96        mp_limb_t *factors;
    97        mp_limb_t prod, max_prod, 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  }