Home | History | Annotate | Line # | Download | only in mpz 
      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, 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 /* TODO: Remove duplicated constants / macros / static functions...
     36  */
     37 
     38 /*************************************************************/
     39 /* Section macros: common macros, for swing/fac/bin (&sieve) */
     40 /*************************************************************/
     41 
     42 #define FACTOR_LIST_STORE(P, PR, MAX_PR, VEC, I)		\
     43   do {								\
     44     if ((PR) > (MAX_PR)) {					\
     45       (VEC)[(I)++] = (PR);					\
     46       (PR) = (P);						\
     47     } else							\
     48       (PR) *= (P);						\
     49   } while (0)
     50 
     51 #define LOOP_ON_SIEVE_CONTINUE(prime,end,sieve)			\
     52     __max_i = (end);						\
     53 								\
     54     do {							\
     55       ++__i;							\
     56       if (((sieve)[__index] & __mask) == 0)			\
     57 	{							\
     58 	  mp_limb_t prime;					\
     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   ASSERT (n <= GMP_NUMB_MAX);
    112 
    113   if (n < 5)
    114     {
    115       MPZ_NEWALLOC (x, 1)[0] = (066211 >> (n*3)) & 7;
    116       SIZ (x) = 1;
    117     }
    118   else
    119     {
    120       mp_limb_t *sieve, *factors;
    121       mp_size_t size, j;
    122       mp_limb_t prod;
    123       TMP_DECL;
    124 
    125       size = n / GMP_NUMB_BITS;
    126       size = size + (size >> 1) + 1;
    127       ASSERT (size >= primesieve_size (n));
    128       sieve = MPZ_NEWALLOC (x, size);
    129       size = (gmp_primesieve (sieve, n) + 1) / log_n_max (n) + 1;
    130 
    131       TMP_MARK;
    132       factors = TMP_ALLOC_LIMBS (size);
    133 
    134       j = 0;
    135 
    136       prod = 6;
    137 
    138       /* Store primes from 5 to n */
    139       {
    140 	mp_limb_t max_prod;
    141 
    142 	max_prod = GMP_NUMB_MAX / n;
    143 
    144 	LOOP_ON_SIEVE_BEGIN (prime, n_to_bit(5), n_to_bit (n), 0, sieve);
    145 	FACTOR_LIST_STORE (prime, prod, max_prod, factors, j);
    146 	LOOP_ON_SIEVE_END;
    147       }
    148 
    149       if (j != 0)
    150 	{
    151 	  factors[j++] = prod;
    152 	  mpz_prodlimbs (x, factors, j);
    153 	}
    154       else
    155 	{
    156 	  PTR (x)[0] = prod;
    157 	  SIZ (x) = 1;
    158 	}
    159 
    160       TMP_FREE;
    161     }
    162 }
    163