1 #ifndef CADO_UTILS_GMP_AUX_H_
2 #define CADO_UTILS_GMP_AUX_H_
3
4 #include "cado_config.h" // for ULONG_BITS
5 #include <gmp.h>
6 #include <stddef.h> // for size_t
7 #include <stdint.h>
8 #include "macros.h"
9
10 /* the following function are missing in GMP */
11 #ifndef mpz_add_si
12 static inline void
mpz_add_si(mpz_ptr a,mpz_srcptr b,const long c)13 mpz_add_si (mpz_ptr a, mpz_srcptr b, const long c) {
14 if (c >= 0)
15 mpz_add_ui (a, b, (unsigned long) c);
16 else
17 mpz_sub_ui (a, b, -(unsigned long) c);
18 }
19 #endif
20
21 #ifndef mpz_sub_si
22 static inline void
mpz_sub_si(mpz_ptr a,mpz_srcptr b,const long c)23 mpz_sub_si (mpz_ptr a, mpz_srcptr b, const long c) {
24 if (c >= 0)
25 mpz_sub_ui (a, b, (unsigned long) c);
26 else
27 mpz_add_ui (a, b, -(unsigned long) c);
28 }
29 #endif
30
31 #ifndef mpz_addmul_si
32 static inline void
mpz_addmul_si(mpz_ptr a,mpz_srcptr b,const long c)33 mpz_addmul_si (mpz_ptr a, mpz_srcptr b, const long c) {
34 if (c >= 0)
35 mpz_addmul_ui (a, b, (unsigned long) c);
36 else
37 mpz_submul_ui (a, b, -(unsigned long) c);
38 }
39 #endif
40
41 #ifndef mpz_submul_si
42 static inline void
mpz_submul_si(mpz_ptr a,mpz_srcptr b,const long c)43 mpz_submul_si (mpz_ptr a, mpz_srcptr b, const long c) {
44 if (c >= 0)
45 mpz_submul_ui (a, b, (unsigned long) c);
46 else mpz_addmul_ui (a, b, -(unsigned long) c);
47 }
48 #endif
49
50 #ifdef __cplusplus
51 extern "C" {
52 #endif
53
54 /* gmp_aux */
55 extern unsigned long ulong_nextprime (unsigned long);
56
57 #if ULONG_BITS < 64
58 extern void mpz_init_set_uint64 (mpz_ptr, uint64_t);
59 extern void mpz_init_set_int64 (mpz_ptr, int64_t);
60 extern void mpz_set_uint64 (mpz_ptr, uint64_t);
61 extern void mpz_set_int64 (mpz_ptr, int64_t);
62 extern uint64_t mpz_get_uint64 (mpz_srcptr);
63 extern int64_t mpz_get_int64 (mpz_srcptr);
64 extern void mpz_mul_uint64 (mpz_ptr a, mpz_srcptr b, uint64_t c);
65 extern int mpz_cmp_uint64 (mpz_srcptr a, uint64_t c);
66 extern int mpz_cmp_int64 (mpz_srcptr a, int64_t c);
67 extern void mpz_add_uint64 (mpz_ptr a, mpz_srcptr b, uint64_t c);
68 extern void mpz_add_int64 (mpz_ptr a, mpz_srcptr b, int64_t c);
69 extern void mpz_sub_uint64 (mpz_ptr a, mpz_srcptr b, uint64_t c);
70 extern void mpz_sub_int64 (mpz_ptr a, mpz_srcptr b, int64_t c);
71 extern void mpz_addmul_uint64 (mpz_ptr a, mpz_srcptr b, uint64_t c);
72 extern void mpz_submul_uint64 (mpz_ptr a, mpz_srcptr b, uint64_t c);
73 extern void mpz_divexact_uint64 (mpz_ptr a, mpz_srcptr b, uint64_t c);
74 extern void mpz_mul_int64 (mpz_ptr a, mpz_srcptr b, int64_t c);
75 extern int mpz_fits_uint64_p(mpz_srcptr);
76 extern int mpz_fits_sint64_p(mpz_srcptr);
77 extern int mpz_divisible_uint64_p (mpz_ptr a, uint64_t c);
78 extern uint64_t uint64_nextprime (uint64_t);
mpz_uint64_sub(mpz_ptr a,const uint64_t b,mpz_srcptr c)79 static inline void mpz_uint64_sub (mpz_ptr a, const uint64_t b, mpz_srcptr c) {
80 mpz_sub_uint64 (a, c, b);
81 mpz_neg (a, a);
82 }
83 extern uint64_t mpz_tdiv_qr_uint64 (mpz_ptr Q, mpz_ptr R, mpz_srcptr N, uint64_t d);
84 extern uint64_t mpz_tdiv_q_uint64 (mpz_ptr Q, mpz_srcptr N, uint64_t d);
85 extern uint64_t mpz_tdiv_r_uint64 (mpz_ptr R, mpz_srcptr N, uint64_t d);
86 extern uint64_t mpz_tdiv_uint64 (mpz_srcptr N, uint64_t d);
87 #else
mpz_init_set_uint64(mpz_ptr a,const uint64_t b)88 static inline void mpz_init_set_uint64 (mpz_ptr a, const uint64_t b) {mpz_init_set_ui(a, b);}
mpz_init_set_int64(mpz_ptr a,const int64_t b)89 static inline void mpz_init_set_int64 (mpz_ptr a, const int64_t b) {mpz_init_set_si(a, b);}
mpz_set_uint64(mpz_ptr a,const uint64_t b)90 static inline void mpz_set_uint64 (mpz_ptr a, const uint64_t b) {mpz_set_ui(a, b);}
mpz_set_int64(mpz_ptr a,const int64_t b)91 static inline void mpz_set_int64 (mpz_ptr a, const int64_t b) {mpz_set_si(a, b);}
mpz_get_uint64(mpz_srcptr a)92 static inline uint64_t mpz_get_uint64 (mpz_srcptr a) {return (uint64_t) mpz_get_ui(a);}
mpz_get_int64(mpz_srcptr a)93 static inline int64_t mpz_get_int64 (mpz_srcptr a) {return (int64_t) mpz_get_si(a);}
mpz_mul_uint64(mpz_ptr a,mpz_srcptr b,const uint64_t c)94 static inline void mpz_mul_uint64 (mpz_ptr a, mpz_srcptr b, const uint64_t c) {mpz_mul_ui(a, b, c);}
mpz_cmp_uint64(mpz_srcptr a,const uint64_t c)95 static inline int mpz_cmp_uint64 (mpz_srcptr a, const uint64_t c) {return mpz_cmp_ui(a, c);}
mpz_cmp_int64(mpz_srcptr a,const int64_t c)96 static inline int mpz_cmp_int64 (mpz_srcptr a, const int64_t c) {return mpz_cmp_si(a, c);}
mpz_add_uint64(mpz_ptr a,mpz_srcptr b,const uint64_t c)97 static inline void mpz_add_uint64 (mpz_ptr a, mpz_srcptr b, const uint64_t c) {mpz_add_ui(a, b, c);}
mpz_add_int64(mpz_ptr a,mpz_srcptr b,const int64_t c)98 static inline void mpz_add_int64 (mpz_ptr a, mpz_srcptr b, const int64_t c) {mpz_add_si(a, b, c);}
mpz_sub_uint64(mpz_ptr a,mpz_srcptr b,const uint64_t c)99 static inline void mpz_sub_uint64 (mpz_ptr a, mpz_srcptr b, const uint64_t c) {mpz_sub_ui(a, b, c);}
mpz_sub_int64(mpz_ptr a,mpz_srcptr b,const int64_t c)100 static inline void mpz_sub_int64 (mpz_ptr a, mpz_srcptr b, const int64_t c) {mpz_sub_si(a, b, c);}
mpz_addmul_uint64(mpz_ptr a,mpz_srcptr b,const uint64_t c)101 static inline void mpz_addmul_uint64 (mpz_ptr a, mpz_srcptr b, const uint64_t c) {mpz_addmul_ui(a, b, c);}
mpz_submul_uint64(mpz_ptr a,mpz_srcptr b,const uint64_t c)102 static inline void mpz_submul_uint64 (mpz_ptr a, mpz_srcptr b, const uint64_t c) {mpz_submul_ui(a, b, c);}
mpz_divexact_uint64(mpz_ptr a,mpz_srcptr b,const uint64_t c)103 static inline void mpz_divexact_uint64 (mpz_ptr a, mpz_srcptr b, const uint64_t c) {mpz_divexact_ui(a, b, c);}
mpz_mul_int64(mpz_ptr a,mpz_srcptr b,const int64_t c)104 static inline void mpz_mul_int64 (mpz_ptr a, mpz_srcptr b, const int64_t c) {mpz_mul_si(a, b, c);}
mpz_fits_uint64_p(mpz_srcptr a)105 static inline int mpz_fits_uint64_p(mpz_srcptr a) {return mpz_fits_ulong_p(a);}
mpz_fits_sint64_p(mpz_srcptr a)106 static inline int mpz_fits_sint64_p(mpz_srcptr a) {return mpz_fits_slong_p(a);}
mpz_divisible_uint64_p(mpz_ptr a,const uint64_t c)107 static inline int mpz_divisible_uint64_p (mpz_ptr a, const uint64_t c) {return mpz_divisible_ui_p(a, c);}
uint64_nextprime(uint64_t a)108 static inline uint64_t uint64_nextprime (uint64_t a) {return (uint64_t) ulong_nextprime(a);}
mpz_uint64_sub(mpz_ptr a,const uint64_t b,mpz_srcptr c)109 static inline void mpz_uint64_sub(mpz_ptr a, const uint64_t b, mpz_srcptr c) {mpz_ui_sub(a, b, c);}
110
mpz_tdiv_qr_uint64(mpz_ptr Q,mpz_ptr R,mpz_srcptr N,uint64_t d)111 static inline uint64_t mpz_tdiv_qr_uint64 (mpz_ptr Q, mpz_ptr R, mpz_srcptr N, uint64_t d) {return mpz_tdiv_qr_ui(Q, R, N, d);}
mpz_tdiv_q_uint64(mpz_ptr Q,mpz_srcptr N,uint64_t d)112 static inline uint64_t mpz_tdiv_q_uint64 (mpz_ptr Q, mpz_srcptr N, uint64_t d) {return mpz_tdiv_q_ui(Q, N, d);}
mpz_tdiv_r_uint64(mpz_ptr R,mpz_srcptr N,uint64_t d)113 static inline uint64_t mpz_tdiv_r_uint64 (mpz_ptr R, mpz_srcptr N, uint64_t d) {return mpz_tdiv_r_ui(R, N, d);}
mpz_tdiv_uint64(mpz_srcptr N,uint64_t d)114 static inline uint64_t mpz_tdiv_uint64 (mpz_srcptr N, uint64_t d) {return mpz_tdiv_ui(N, d);}
115 #endif
116
117 extern void mpz_submul_int64 (mpz_ptr a, mpz_srcptr b, int64_t c);
118 extern void mpz_addmul_int64 (mpz_ptr a, mpz_srcptr b, int64_t c);
119 extern int ulong_isprime (unsigned long);
120 extern unsigned long ulong_nextcomposite (unsigned long, unsigned long);
121 extern void mpz_ndiv_qr (mpz_ptr q, mpz_ptr r, mpz_srcptr n, mpz_srcptr d);
122 extern void mpz_ndiv_r (mpz_ptr r, mpz_srcptr n, mpz_srcptr d);
123 extern void mpz_ndiv_qr_ui (mpz_ptr q, mpz_ptr r, mpz_srcptr n, unsigned long int d);
124 extern void mpz_ndiv_q (mpz_ptr q, mpz_srcptr n, mpz_srcptr d);
125 extern void mpz_ndiv_q_ui (mpz_ptr q, mpz_srcptr n, unsigned long int d);
126 extern int mpz_coprime_p (mpz_srcptr a, mpz_srcptr b);
127
128 extern int mpz_rdiv_q(mpz_ptr q, mpz_srcptr a, mpz_srcptr b);
129
130 /* Put in r the smallest legitimate value that it at least s + diff (note
131 that if s+diff is already legitimate, then r = s+diff will result.
132
133 Here, legitimate means prime or squarefree composite, with the constraint
134 that all the prime factors must be in [pmin, pmax[ .
135
136 The prime factors of r are put in factor_r, and the number of them is
137 returned. The caller must have allocated factor_r with enough space.
138 */
139 int
140 next_mpz_with_factor_constraints(mpz_ptr r,
141 unsigned long factor_r[],
142 mpz_srcptr s,
143 unsigned long diff,
144 unsigned long pmin,
145 unsigned long pmax);
146
147 /* return the number of bits of p, counting from the least significant end */
148 extern int nbits (uintmax_t p);
149 extern long double mpz_get_ld (mpz_srcptr z);
150
151 extern int mpz_p_valuation(mpz_srcptr a, mpz_srcptr p);
152 extern int mpz_p_valuation_ui(mpz_srcptr a, unsigned long p);
153
154 #if !GMP_VERSION_ATLEAST(5,0,0)
155 mp_limb_t mpn_neg (mp_limb_t *rp, const mp_limb_t *sp, mp_size_t n);
156 void mpn_xor_n (mp_limb_t *rp, const mp_limb_t *s1p, const mp_limb_t *s2p,
157 mp_size_t n);
158 void mpn_zero(mp_limb_t *rp, mp_size_t n);
159 void mpn_copyi(mp_limb_t *rp, const mp_limb_t * up, mp_size_t n);
160 void mpn_copyd(mp_limb_t *rp, const mp_limb_t * up, mp_size_t n);
161 #endif
162
163 #if !GMP_VERSION_ATLEAST(6,1,0)
164 int mpn_zero_p(const mp_limb_t *rp, mp_size_t n);
165 #endif
166
167 #ifndef HAVE_MPIR
168
169 /* Yes, it's a bit ugly of course. */
mpn_rrandom(mp_limb_t * rp,gmp_randstate_t rstate,mp_size_t N)170 static inline void mpn_rrandom (mp_limb_t *rp, gmp_randstate_t rstate, mp_size_t N)
171 {
172 mpz_t dummy;
173 dummy->_mp_d = rp;
174 dummy->_mp_alloc = N;
175 dummy->_mp_size = N;
176 mpz_rrandomb(dummy, rstate, (mp_bitcnt_t) N * GMP_LIMB_BITS);
177 }
178
179
mpn_randomb(mp_limb_t * rp,gmp_randstate_t rstate,mp_size_t N)180 static inline void mpn_randomb (mp_limb_t *rp, gmp_randstate_t rstate, mp_size_t N)
181 {
182 mpz_t dummy;
183 dummy->_mp_d = rp;
184 dummy->_mp_alloc = N;
185 dummy->_mp_size = N;
186 mpz_urandomb(dummy, rstate, (mp_bitcnt_t) N * GMP_LIMB_BITS);
187 }
188
189 #endif
190
191 void memfill_random(void *p, size_t s, gmp_randstate_t rstate);
192
193
194 /* The lack of these has become really annoying. There's a patch floating
195 * on the gmp list (may 2018, with rare reminders and no apparent action
196 * taken in the following months. Maybe it will be merged for gmp 17.42
197 * eventually...)
198 */
199 #if !GMP_VERSION_ATLEAST(17,42,0)
200 typedef __gmp_randstate_struct * gmp_randstate_ptr;
201 typedef const __gmp_randstate_struct * gmp_randstate_srcptr;
202 #endif
203
204 #ifdef __cplusplus
205 }
206 #endif
207
208 #ifdef __cplusplus
209 /* Same idea as for cxx_mpz and friends */
210 struct cxx_gmp_randstate {
211 gmp_randstate_t x;
cxx_gmp_randstatecxx_gmp_randstate212 cxx_gmp_randstate() { gmp_randinit_default(x); }
cxx_gmp_randstatecxx_gmp_randstate213 cxx_gmp_randstate(gmp_randstate_srcptr f) { gmp_randinit_set(x, f); }
~cxx_gmp_randstatecxx_gmp_randstate214 ~cxx_gmp_randstate() { gmp_randclear(x); }
cxx_gmp_randstatecxx_gmp_randstate215 cxx_gmp_randstate(cxx_gmp_randstate const & o) {
216 gmp_randinit_set(x, o.x);
217 }
218 cxx_gmp_randstate & operator=(cxx_gmp_randstate const & o) {
219 gmp_randclear(x);
220 gmp_randinit_set(x, o.x);
221 return *this;
222 }
gmp_randstate_ptrcxx_gmp_randstate223 operator gmp_randstate_ptr() { return x; }
gmp_randstate_srcptrcxx_gmp_randstate224 operator gmp_randstate_srcptr() const { return x; }
225 gmp_randstate_ptr operator->() { return x; }
226 gmp_randstate_srcptr operator->() const { return x; }
227 };
228
229 #if GNUC_VERSION_ATLEAST(4,3,0)
230 extern void gmp_randinit_default(cxx_gmp_randstate & pl) __attribute__((error("gmp_randinit_default must not be called on a cxx_gmp_randstate reference -- it is the caller's business (via a ctor)")));
231 extern void gmp_randclear(cxx_gmp_randstate & pl) __attribute__((error("gmp_randclear must not be called on a cxx_gmp_randstate reference -- it is the caller's business (via a dtor)")));
232 #endif
233
234 #endif
235
236
237 #endif /* CADO_UTILS_GMP_AUX_H_ */
238