1 /*
2 Copyright (C) 2010 Fredrik Johansson
3 Copyright (C) 2010 William Hart
4
5 This file is part of FLINT.
6
7 FLINT is free software: you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License (LGPL) as published
9 by the Free Software Foundation; either version 2.1 of the License, or
10 (at your option) any later version. See <http://www.gnu.org/licenses/>.
11 */
12
13 #include "arith.h"
14
15 #if FLINT64
16 #define LARGEST_ULONG_PRIMORIAL 52
17 #else
18 #define LARGEST_ULONG_PRIMORIAL 28
19 #endif
20
21 /* Only those with odd index */
22 const ulong ULONG_PRIMORIALS[] =
23 {
24 UWORD(6),
25 UWORD(30),
26 UWORD(210),
27 UWORD(210),
28 UWORD(2310),
29 UWORD(30030),
30 UWORD(30030),
31 UWORD(510510),
32 UWORD(9699690),
33 UWORD(9699690),
34 UWORD(223092870),
35 UWORD(223092870),
36 UWORD(223092870),
37 #if FLINT64
38 UWORD(6469693230),
39 UWORD(200560490130),
40 UWORD(200560490130),
41 UWORD(200560490130),
42 UWORD(7420738134810),
43 UWORD(7420738134810),
44 UWORD(304250263527210),
45 UWORD(13082761331670030),
46 UWORD(13082761331670030),
47 UWORD(614889782588491410),
48 UWORD(614889782588491410),
49 UWORD(614889782588491410)
50 #endif
51 };
52
53
54 #define PROD_LIMBS_DIRECT_CUTOFF 50
55
mpn_prod_limbs_direct(mp_limb_t * result,const mp_limb_t * factors,mp_size_t n)56 mp_size_t mpn_prod_limbs_direct(mp_limb_t * result, const mp_limb_t * factors,
57 mp_size_t n)
58 {
59 mp_size_t k, len;
60 mp_limb_t top;
61 if (n < 1)
62 {
63 result[0] = UWORD(1);
64 return 1;
65 }
66 result[0] = factors[0];
67 len = 1;
68 for (k=1; k<n; k++)
69 {
70 top = mpn_mul_1(result, result, len, factors[k]);
71 if (top)
72 {
73 result[len] = top;
74 len++;
75 }
76 }
77 return len;
78 }
79
mpn_prod_limbs_balanced(mp_limb_t * result,mp_limb_t * scratch,const mp_limb_t * factors,mp_size_t n,ulong bits)80 mp_size_t mpn_prod_limbs_balanced(mp_limb_t * result, mp_limb_t * scratch,
81 const mp_limb_t * factors, mp_size_t n, ulong bits)
82 {
83 mp_size_t an, bn, alen, blen, len;
84 mp_limb_t top;
85
86 if (n < PROD_LIMBS_DIRECT_CUTOFF)
87 return mpn_prod_limbs_direct(result, factors, n);
88
89 an = n/2;
90 bn = n - an;
91
92 alen = mpn_prod_limbs_balanced(scratch, result, factors, an, bits);
93 blen = mpn_prod_limbs_balanced(scratch + alen, result, factors + an, bn, bits);
94 len = alen + blen;
95
96 if (alen <= blen)
97 top = mpn_mul(result, scratch + alen, blen, scratch, alen);
98 else
99 top = mpn_mul(result, scratch, alen, scratch + alen, blen);
100
101 if (!top)
102 len--;
103
104 return len;
105 }
106
107 /*
108 Set result to the product of the given factors, return the
109 length of the result. It is assumed that no factors are zero.
110 bits must be set to some bound on the bit size of the entries
111 in factors. If no bound is known, simply use FLINT_BITS.
112 */
mpn_prod_limbs(mp_limb_t * result,const mp_limb_t * factors,mp_size_t n,ulong bits)113 mp_size_t mpn_prod_limbs(mp_limb_t * result, const mp_limb_t * factors,
114 mp_size_t n, ulong bits)
115 {
116 mp_size_t len, limbs;
117 mp_limb_t * scratch;
118
119 if (n < PROD_LIMBS_DIRECT_CUTOFF)
120 return mpn_prod_limbs_direct(result, factors, n);
121
122 limbs = (n * bits - 1)/FLINT_BITS + 2;
123
124 scratch = flint_malloc(sizeof(mp_limb_t) * limbs);
125 len = mpn_prod_limbs_balanced(result, scratch, factors, n, bits);
126 flint_free(scratch);
127
128 return len;
129 }
130
131 void
fmpz_primorial(fmpz_t res,ulong n)132 fmpz_primorial(fmpz_t res, ulong n)
133 {
134 mp_size_t len, pi;
135 ulong bits;
136 __mpz_struct * mpz_ptr;
137 const mp_limb_t * primes;
138
139 if (n <= LARGEST_ULONG_PRIMORIAL)
140 {
141 if (n <= 2)
142 fmpz_set_ui(res, 1 + (n==2));
143 else
144 fmpz_set_ui(res, ULONG_PRIMORIALS[(n-1)/2-1]);
145 return;
146 }
147
148 pi = n_prime_pi(n);
149
150 primes = n_primes_arr_readonly(pi);
151 bits = FLINT_BIT_COUNT(primes[pi - 1]);
152
153 mpz_ptr = _fmpz_promote(res);
154 mpz_realloc2(mpz_ptr, pi*bits);
155
156 len = mpn_prod_limbs(mpz_ptr->_mp_d, primes, pi, bits);
157 mpz_ptr->_mp_size = len;
158 }
159
160