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

     1  /* mpz_primorial_ui(RESULT, N) -- Set RESULT to N# the product of primes <= N.
     2  
     3  Contributed to the GNU project by Marco Bodrato.
     4  
     5  Copyright 2012 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  /* TODO: Remove duplicated constants / macros / static functions...
    37   */
    38  
    39  /*************************************************************/
    40  /* Section macros: common macros, for swing/fac/bin (&sieve) */
    41  /*************************************************************/
    42  
    43  #define FACTOR_LIST_STORE(P, PR, MAX_PR, VEC, I)		\
    44    do {								\
    45      if ((PR) > (MAX_PR)) {					\
    46        (VEC)[(I)++] = (PR);					\
    47        (PR) = (P);						\
    48      } else							\
    49        (PR) *= (P);						\
    50    } while (0)
    51  
    52  #define LOOP_ON_SIEVE_CONTINUE(prime,end,sieve)			\
    53      __max_i = (end);						\
    54  								\
    55      do {							\
    56        ++__i;							\
    57        if (((sieve)[__index] & __mask) == 0)			\
    58  	{							\
    59  	  (prime) = id_to_n(__i)
    60  
    61  #define LOOP_ON_SIEVE_BEGIN(prime,start,end,off,sieve)		\
    62    do {								\
    63      mp_limb_t __mask, __index, __max_i, __i;			\
    64  								\
    65      __i = (start)-(off);					\
    66      __index = __i / GMP_LIMB_BITS;				\
    67      __mask = CNST_LIMB(1) << (__i % GMP_LIMB_BITS);		\
    68      __i += (off);						\
    69  								\
    70      LOOP_ON_SIEVE_CONTINUE(prime,end,sieve)
    71  
    72  #define LOOP_ON_SIEVE_STOP					\
    73  	}							\
    74        __mask = __mask << 1 | __mask >> (GMP_LIMB_BITS-1);	\
    75        __index += __mask & 1;					\
    76      }  while (__i <= __max_i)					\
    77  
    78  #define LOOP_ON_SIEVE_END					\
    79      LOOP_ON_SIEVE_STOP;						\
    80    } while (0)
    81  
    82  /*********************************************************/
    83  /* Section sieve: sieving functions and tools for primes */
    84  /*********************************************************/
    85  
    86  #if 0
    87  static mp_limb_t
    88  bit_to_n (mp_limb_t bit) { return (bit*3+4)|1; }
    89  #endif
    90  
    91  /* id_to_n (x) = bit_to_n (x-1) = (id*3+1)|1*/
    92  static mp_limb_t
    93  id_to_n  (mp_limb_t id)  { return id*3+1+(id&1); }
    94  
    95  /* n_to_bit (n) = ((n-1)&(-CNST_LIMB(2)))/3U-1 */
    96  static mp_limb_t
    97  n_to_bit (mp_limb_t n) { return ((n-5)|1)/3U; }
    98  
    99  #if WANT_ASSERT
   100  static mp_size_t
   101  primesieve_size (mp_limb_t n) { return n_to_bit(n) / GMP_LIMB_BITS + 1; }
   102  #endif
   103  
   104  /*********************************************************/
   105  /* Section primorial: implementation                     */
   106  /*********************************************************/
   107  
   108  void
   109  mpz_primorial_ui (mpz_ptr x, unsigned long n)
   110  {
   111    static const mp_limb_t table[] = { 1, 1, 2, 6, 6 };
   112  
   113    ASSERT (n <= GMP_NUMB_MAX);
   114  
   115    if (n < numberof (table))
   116      {
   117        PTR (x)[0] = table[n];
   118        SIZ (x) = 1;
   119      }
   120    else
   121      {
   122        mp_limb_t *sieve, *factors;
   123        mp_size_t size;
   124        mp_limb_t prod;
   125        mp_limb_t j;
   126        TMP_DECL;
   127  
   128        size = 1 + n / GMP_NUMB_BITS + n / (2*GMP_NUMB_BITS);
   129        ASSERT (size >= primesieve_size (n));
   130        sieve = MPZ_NEWALLOC (x, size);
   131        size = (gmp_primesieve (sieve, n) + 1) / log_n_max (n) + 1;
   132  
   133        TMP_MARK;
   134        factors = TMP_ALLOC_LIMBS (size);
   135  
   136        j = 0;
   137  
   138        prod = table[numberof (table)-1];
   139  
   140        /* Store primes from 5 to n */
   141        {
   142  	mp_limb_t prime, max_prod;
   143  
   144  	max_prod = GMP_NUMB_MAX / n;
   145  
   146  	LOOP_ON_SIEVE_BEGIN (prime, n_to_bit(numberof (table)), n_to_bit (n), 0, sieve);
   147  	FACTOR_LIST_STORE (prime, prod, max_prod, factors, j);
   148  	LOOP_ON_SIEVE_END;
   149        }
   150  
   151        if (LIKELY (j != 0))
   152  	{
   153  	  factors[j++] = prod;
   154  	  mpz_prodlimbs (x, factors, j);
   155  	}
   156        else
   157  	{
   158  	  PTR (x)[0] = prod;
   159  	  SIZ (x) = 1;
   160  	}
   161  
   162        TMP_FREE;
   163      }
   164  }