mp_class.h
1 #ifndef SYMENGINE_INTEGER_CLASS_H
2 #define SYMENGINE_INTEGER_CLASS_H
3 
4 #include <symengine/symengine_config.h>
5 #include <symengine/symengine_casts.h>
6 #if SYMENGINE_INTEGER_CLASS != SYMENGINE_BOOSTMP
7 #include <symengine/mp_wrapper.h>
8 #endif
9 
10 #if SYMENGINE_INTEGER_CLASS == SYMENGINE_BOOSTMP
11 #include <boost/multiprecision/cpp_int.hpp>
12 #include <boost/random/mersenne_twister.hpp>
13 #include <boost/random/uniform_int.hpp>
14 #include <symengine/symengine_rcp.h>
15 #elif SYMENGINE_INTEGER_CLASS == SYMENGINE_PIRANHA
16 #include <piranha/mp_integer.hpp>
17 #include <piranha/mp_rational.hpp>
18 #elif SYMENGINE_INTEGER_CLASS == SYMENGINE_GMPXX
19 #define __GMPXX_USE_CXX11 1
20 #include <gmpxx.h>
21 #endif
22 
23 namespace SymEngine
24 {
25 /*
26  * integer_class and rational_class are the storage classes for
27  * SymEngine::Integer
28  * and SymEngine::Rational respectively. There are 4 choices for the
29  * integer_class.
30  * mpz_class from libgmpxx, piranha::integer from piranha, mpz_wrapper which
31  * wraps
32  * mpz_t from libgmp and fmpz_wrapper which wraps fmpz_t from libflint. This
33  * choice
34  * is made at compile time with SYMENGINE_INTEGER_CLASS.
35  *
36  * Each of these classes has to have all the arithmetic operators overloaded
37  * with
38  * operands of these classes and C++ integer types except long long. Also, shift
39  * operators, move operators, string, integer and mpz_t constructors are also
40  * required.
41  *
42  * To add a new type, several non-member functions need to be defined for the
43  * new
44  * type. See mpz_wrapper implementation for all the non-member functions that
45  * need
46  * to be defined.
47  */
48 
49 #if SYMENGINE_INTEGER_CLASS == SYMENGINE_BOOSTMP
50 typedef boost::multiprecision::number<boost::multiprecision::cpp_int_backend<>,
51  boost::multiprecision::et_off>
52  integer_class;
53 typedef boost::multiprecision::number<
54  boost::multiprecision::cpp_rational_backend, boost::multiprecision::et_off>
55  rational_class;
56 #elif SYMENGINE_INTEGER_CLASS == SYMENGINE_PIRANHA
57 typedef piranha::integer integer_class;
58 typedef piranha::rational rational_class;
59 #elif SYMENGINE_INTEGER_CLASS == SYMENGINE_FLINT
60 typedef fmpz_wrapper integer_class;
61 typedef fmpq_wrapper rational_class;
62 #elif SYMENGINE_INTEGER_CLASS == SYMENGINE_GMP
63 typedef mpz_wrapper integer_class;
64 typedef mpq_wrapper rational_class;
65 #elif SYMENGINE_INTEGER_CLASS == SYMENGINE_GMPXX
66 typedef mpz_class integer_class;
67 typedef mpq_class rational_class;
68 #endif
69 
70 // needs to be in a separate namespace to import the literals.
71 // eg: using namespace SymEngine::literals;
72 inline namespace literals
73 {
75 inline integer_class operator"" _z(const char *str)
76 {
77  return integer_class(str);
78 }
79 
80 inline rational_class operator"" _q(const char *str)
81 {
82  return rational_class(integer_class(str));
83 }
84 } // namespace literals
85 
86 #if SYMENGINE_INTEGER_CLASS == SYMENGINE_GMPXX \
87  || SYMENGINE_INTEGER_CLASS == SYMENGINE_GMP
88 // Helper functions for mpz
89 inline integer_class mp_abs(const integer_class &i)
90 {
91  integer_class res;
92  mpz_abs(res.get_mpz_t(), i.get_mpz_t());
93  return res;
94 }
95 
96 inline int mp_sign(const integer_class &i)
97 {
98  return mpz_sgn(i.get_mpz_t());
99 }
100 
101 inline integer_class mp_sqrt(const integer_class &i)
102 {
103  integer_class res;
104  mpz_sqrt(res.get_mpz_t(), i.get_mpz_t());
105  return res;
106 }
107 
108 inline double mp_get_d(const integer_class &i)
109 {
110  return static_cast<double>(i.get_d());
111 }
112 
113 inline void mp_set_d(integer_class &i, double a)
114 {
115  mpz_set_d(i.get_mpz_t(), a);
116 }
117 
118 inline void mp_demote(integer_class &i) {}
119 
120 inline bool mp_fits_ulong_p(const integer_class &i)
121 {
122  return i.fits_ulong_p() != 0;
123 }
124 
125 inline bool mp_fits_slong_p(const integer_class &i)
126 {
127  return i.fits_slong_p() != 0;
128 }
129 
130 inline unsigned long mp_get_ui(const integer_class &i)
131 {
132  return i.get_ui();
133 }
134 
135 inline long mp_get_si(const integer_class &i)
136 {
137  return i.get_si();
138 }
139 
140 inline mpz_srcptr get_mpz_t(const integer_class &i)
141 {
142  return i.get_mpz_t();
143 }
144 
145 inline mpz_ptr get_mpz_t(integer_class &i)
146 {
147  return i.get_mpz_t();
148 }
149 
150 inline void mp_pow_ui(integer_class &res, const integer_class &i,
151  unsigned long n)
152 {
153  mpz_pow_ui(res.get_mpz_t(), i.get_mpz_t(), n);
154 }
155 
156 inline void mp_pow_ui(rational_class &res, const rational_class &i,
157  unsigned long n)
158 {
159  mpz_pow_ui(res.get_num().get_mpz_t(), i.get_num().get_mpz_t(), n);
160  mpz_pow_ui(res.get_den().get_mpz_t(), i.get_den().get_mpz_t(), n);
161 }
162 
163 inline void mp_powm(integer_class &res, const integer_class &a,
164  const integer_class &b, const integer_class &m)
165 {
166  mpz_powm(res.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t(), m.get_mpz_t());
167 }
168 
169 inline bool mp_invert(integer_class &res, const integer_class &a,
170  const integer_class &m)
171 {
172  return mpz_invert(res.get_mpz_t(), a.get_mpz_t(), m.get_mpz_t()) != 0;
173 }
174 
175 inline void mp_gcd(integer_class &res, const integer_class &a,
176  const integer_class &b)
177 {
178  mpz_gcd(res.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
179 }
180 
181 inline void mp_gcdext(integer_class &res, integer_class &r, integer_class &s,
182  const integer_class &a, const integer_class &b)
183 {
184  mpz_gcdext(res.get_mpz_t(), r.get_mpz_t(), s.get_mpz_t(), a.get_mpz_t(),
185  b.get_mpz_t());
186 }
187 
188 inline void mp_and(integer_class &res, const integer_class &a,
189  const integer_class &b)
190 {
191  mpz_and(res.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
192 }
193 
194 inline void mp_fdiv_r(integer_class &res, const integer_class &a,
195  const integer_class &b)
196 {
197  mpz_fdiv_r(res.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
198 }
199 
200 inline void mp_fdiv_q(integer_class &res, const integer_class &a,
201  const integer_class &b)
202 {
203  mpz_fdiv_q(res.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
204 }
205 
206 inline void mp_cdiv_q(integer_class &res, const integer_class &a,
207  const integer_class &b)
208 {
209  mpz_cdiv_q(res.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
210 }
211 
212 inline void mp_tdiv_q(integer_class &res, const integer_class &a,
213  const integer_class &b)
214 {
215  mpz_tdiv_q(res.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
216 }
217 
218 inline void mp_fdiv_qr(integer_class &q, integer_class &r,
219  const integer_class &a, const integer_class &b)
220 {
221  mpz_fdiv_qr(q.get_mpz_t(), r.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
222 }
223 
224 inline void mp_divexact(integer_class &q, const integer_class &a,
225  const integer_class &b)
226 {
227  mpz_divexact(q.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
228 }
229 
230 inline void mp_lcm(integer_class &q, const integer_class &a,
231  const integer_class &b)
232 {
233  mpz_lcm(q.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
234 }
235 
236 inline void mp_tdiv_qr(integer_class &q, integer_class &r,
237  const integer_class &a, const integer_class &b)
238 {
239  mpz_tdiv_qr(q.get_mpz_t(), r.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
240 }
241 
242 inline void mp_addmul(integer_class &r, const integer_class &a,
243  const integer_class &b)
244 {
245  mpz_addmul(r.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
246 }
247 
248 // Helper functions for mpq
249 inline const integer_class &get_den(const rational_class &i)
250 {
251  return i.get_den();
252 }
253 
254 inline const integer_class &get_num(const rational_class &i)
255 {
256  return i.get_num();
257 }
258 
259 inline integer_class &get_den(rational_class &i)
260 {
261  return i.get_den();
262 }
263 
264 inline integer_class &get_num(rational_class &i)
265 {
266  return i.get_num();
267 }
268 
269 inline mpq_srcptr get_mpq_t(const rational_class &i)
270 {
271  return i.get_mpq_t();
272 }
273 
274 inline void canonicalize(rational_class &i)
275 {
276  i.canonicalize();
277 }
278 
279 inline double mp_get_d(const rational_class &i)
280 {
281  return i.get_d();
282 }
283 
284 inline int mp_sign(const rational_class &i)
285 {
286  return mpq_sgn(i.get_mpq_t());
287 }
288 
289 inline rational_class mp_abs(const rational_class &i)
290 {
291  rational_class res;
292  mpq_abs(res.get_mpq_t(), i.get_mpq_t());
293  return res;
294 }
295 
296 inline integer_class mp_primorial(unsigned long n)
297 {
298  integer_class res;
299  mpz_primorial_ui(res.get_mpz_t(), n);
300  return res;
301 }
302 
303 #elif SYMENGINE_INTEGER_CLASS == SYMENGINE_PIRANHA
304 // Helper functions for piranha::integer
305 inline piranha::integer mp_abs(const piranha::integer &i)
306 {
307  return i.abs();
308 }
309 
310 inline piranha::integer mp_sqrt(const piranha::integer &i)
311 {
312  return i.sqrt();
313 }
314 
315 inline void mp_demote(piranha::integer &i) {}
316 
317 inline mpz_ptr get_mpz_t(piranha::integer &i)
318 {
319  return i._get_mpz_ptr();
320 }
321 
322 inline auto get_mpz_t(const piranha::integer &i) -> decltype(i.get_mpz_view())
323 {
324  return i.get_mpz_view();
325 }
326 
327 inline void mp_pow_ui(piranha::integer &res, const piranha::integer &i,
328  unsigned long n)
329 {
330  res = i.pow(n);
331 }
332 
333 inline void mp_pow_ui(piranha::rational &res, const piranha::rational &i,
334  unsigned long n)
335 {
336  res = i.pow(n);
337 }
338 
339 inline void mp_powm(piranha::integer &res, const piranha::integer &a,
340  const piranha::integer &b, const piranha::integer &m)
341 {
342  auto _res = get_mpz_t(res);
343  mpz_powm(_res, get_mpz_t(a), get_mpz_t(b), get_mpz_t(m));
344 }
345 
346 inline bool mp_invert(piranha::integer &res, const piranha::integer &a,
347  const piranha::integer &m)
348 {
349  auto _res = get_mpz_t(res);
350  return mpz_invert(_res, get_mpz_t(a), get_mpz_t(m)) != 0;
351 }
352 
353 inline void mp_gcd(piranha::integer &g, const piranha::integer &a,
354  const piranha::integer &b)
355 {
356  piranha::integer::gcd(g, a, b);
357 }
358 
359 inline void mp_gcdext(piranha::integer &g, piranha::integer &r,
360  piranha::integer &s, const piranha::integer &a,
361  const piranha::integer &b)
362 {
363  auto _g = get_mpz_t(g);
364  auto _r = get_mpz_t(r);
365  auto _s = get_mpz_t(s);
366  mpz_gcdext(_g, _r, _s, get_mpz_t(a), get_mpz_t(b));
367 }
368 
369 inline void mp_and(piranha::integer &res, const piranha::integer &a,
370  const piranha::integer &b)
371 {
372  auto _res = get_mpz_t(res);
373  mpz_and(_res, get_mpz_t(a), get_mpz_t(b));
374 }
375 
376 inline void mp_fdiv_r(piranha::integer &res, const piranha::integer &a,
377  const piranha::integer &b)
378 {
379  auto _res = get_mpz_t(res);
380  mpz_fdiv_r(_res, get_mpz_t(a), get_mpz_t(b));
381 }
382 
383 inline void mp_fdiv_q(piranha::integer &res, const piranha::integer &a,
384  const piranha::integer &b)
385 {
386  auto _res = get_mpz_t(res);
387  mpz_fdiv_q(_res, get_mpz_t(a), get_mpz_t(b));
388 }
389 
390 inline void mp_cdiv_q(piranha::integer &res, const piranha::integer &a,
391  const piranha::integer &b)
392 {
393  auto _res = get_mpz_t(res);
394  mpz_cdiv_q(_res, get_mpz_t(a), get_mpz_t(b));
395 }
396 
397 inline void mp_tdiv_q(piranha::integer &res, const piranha::integer &a,
398  const piranha::integer &b)
399 {
400  auto _res = get_mpz_t(res);
401  mpz_tdiv_q(_res, get_mpz_t(a), get_mpz_t(b));
402 }
403 
404 inline void mp_fdiv_qr(piranha::integer &q, piranha::integer &r,
405  const piranha::integer &a, const piranha::integer &b)
406 {
407  auto _q = get_mpz_t(q);
408  mpz_fdiv_qr(_q, get_mpz_t(r), get_mpz_t(a), get_mpz_t(b));
409 }
410 
411 inline void mp_divexact(piranha::integer &q, const piranha::integer &a,
412  const piranha::integer &b)
413 {
414  piranha::integer::_divexact(q, a, b);
415 }
416 
417 inline void mp_lcm(piranha::integer &q, const piranha::integer &a,
418  const piranha::integer &b)
419 {
420  auto _q = get_mpz_t(q);
421  mpz_lcm(_q, get_mpz_t(a), get_mpz_t(b));
422 }
423 
424 inline void mp_tdiv_qr(piranha::integer &q, piranha::integer &r,
425  const piranha::integer &a, const piranha::integer &b)
426 {
427  piranha::integer::divrem(q, r, a, b);
428 }
429 
430 inline int mp_sign(const piranha::integer &i)
431 {
432  return i.sign();
433 }
434 
435 inline long mp_get_si(const piranha::integer &i)
436 {
437  return mpz_get_si(i.get_mpz_view());
438 }
439 
440 inline unsigned long mp_get_ui(const piranha::integer &i)
441 {
442  return mpz_get_ui(i.get_mpz_view());
443 }
444 
445 inline double mp_get_d(const piranha::integer &i)
446 {
447  return mpz_get_d(i.get_mpz_view());
448 }
449 
450 inline void mp_set_d(piranha::integer &i, double a)
451 {
452  i = a;
453 }
454 
455 inline bool mp_fits_ulong_p(const piranha::integer &i)
456 {
457  return mpz_fits_ulong_p(i.get_mpz_view()) != 0;
458 }
459 
460 inline bool mp_fits_slong_p(const piranha::integer &i)
461 {
462  return mpz_fits_slong_p(i.get_mpz_view()) != 0;
463 }
464 
465 inline void mp_addmul(integer_class &r, const integer_class &a,
466  const integer_class &b)
467 {
468  piranha::math::multiply_accumulate(r, a, b);
469 }
470 
471 inline integer_class mp_primorial(unsigned long n)
472 {
473  integer_class res;
474  mpz_primorial_ui(get_mpz_t(res), n);
475  return res;
476 }
477 
478 // Helper functions for piranha::rational
479 
480 inline piranha::rational mp_abs(const piranha::rational &i)
481 {
482  return i.abs();
483 }
484 
485 inline const piranha::integer &get_den(const piranha::rational &i)
486 {
487  return i.den();
488 }
489 
490 inline const piranha::integer &get_num(const piranha::rational &i)
491 {
492  return i.num();
493 }
494 
495 inline piranha::integer &get_den(piranha::rational &i)
496 {
497  return i._den();
498 }
499 
500 inline piranha::integer &get_num(piranha::rational &i)
501 {
502  return i._num();
503 }
504 
505 inline void canonicalize(piranha::rational &i)
506 {
507  i.canonicalise();
508 }
509 
510 inline double mp_get_d(const piranha::rational &i)
511 {
512  return mpq_get_d(i.get_mpq_view().get());
513 }
514 
515 inline auto get_mpq_t(const piranha::rational &i) -> decltype(i.get_mpq_view())
516 {
517  return i.get_mpq_view();
518 }
519 
520 inline int mp_sign(const piranha::rational &i)
521 {
522  return i.num().sign();
523 }
524 #elif SYMENGINE_INTEGER_CLASS == SYMENGINE_FLINT
525 
526 // helper functions for fmpz
527 
528 inline mpz_view_flint get_mpz_t(const fmpz_wrapper &i)
529 {
530  return mpz_view_flint(i);
531 }
532 
533 inline mpz_ptr get_mpz_t(fmpz_wrapper &i)
534 {
535  return _fmpz_promote_val(i.get_fmpz_t());
536 }
537 
538 inline void mp_demote(fmpz_wrapper &i)
539 {
540  _fmpz_demote_val(i.get_fmpz_t());
541 }
542 
543 inline int mp_sign(const fmpz_wrapper &i)
544 {
545  return fmpz_sgn(i.get_fmpz_t());
546 }
547 
548 inline long mp_get_si(const fmpz_wrapper &i)
549 {
550  return fmpz_get_si(i.get_fmpz_t());
551 }
552 
553 inline unsigned long mp_get_ui(const fmpz_wrapper &i)
554 {
555  return fmpz_get_ui(i.get_fmpz_t());
556 }
557 
558 inline bool mp_fits_slong_p(const fmpz_wrapper &i)
559 {
560  return fmpz_fits_si(i.get_fmpz_t());
561 }
562 
563 inline bool mp_fits_ulong_p(const fmpz_wrapper &i)
564 {
565  return fmpz_sgn(i.get_fmpz_t()) >= 0 && fmpz_abs_fits_ui(i.get_fmpz_t());
566 }
567 
568 inline double mp_get_d(const fmpz_wrapper &i)
569 {
570  return fmpz_get_d(i.get_fmpz_t());
571 }
572 
573 inline void mp_set_d(fmpz_wrapper &i, double a)
574 {
575  return fmpz_set_d(i.get_fmpz_t(), a);
576 }
577 
578 inline fmpz_wrapper mp_abs(const fmpz_wrapper &i)
579 {
580  fmpz_wrapper res;
581  fmpz_abs(res.get_fmpz_t(), i.get_fmpz_t());
582  return res;
583 }
584 
585 inline fmpz_wrapper mp_sqrt(const fmpz_wrapper &i)
586 {
587  fmpz_wrapper res;
588  fmpz_sqrt(res.get_fmpz_t(), i.get_fmpz_t());
589  return res;
590 }
591 
592 inline void mp_pow_ui(fmpz_wrapper &res, const fmpz_wrapper &i, unsigned long n)
593 {
594  fmpz_pow_ui(res.get_fmpz_t(), i.get_fmpz_t(), n);
595 }
596 
597 inline void mp_pow_ui(fmpq_wrapper &res, const fmpq_wrapper &i, unsigned long n)
598 {
599  fmpq_pow_si(res.get_fmpq_t(), i.get_fmpq_t(), n);
600 }
601 
602 inline void mp_powm(fmpz_wrapper &res, const fmpz_wrapper &a,
603  const fmpz_wrapper &b, const fmpz_wrapper &m)
604 {
605  if (b >= 0) {
606  fmpz_powm(res.get_fmpz_t(), a.get_fmpz_t(), b.get_fmpz_t(),
607  m.get_fmpz_t());
608  } else {
609  fmpz_neg(res.get_fmpz_t(), b.get_fmpz_t());
610  fmpz_powm(res.get_fmpz_t(), a.get_fmpz_t(), res.get_fmpz_t(),
611  m.get_fmpz_t());
612  fmpz_invmod(res.get_fmpz_t(), res.get_fmpz_t(), m.get_fmpz_t());
613  }
614 }
615 
616 inline bool mp_invert(fmpz_wrapper &res, const fmpz_wrapper &a,
617  const fmpz_wrapper &m)
618 {
619  return fmpz_invmod(res.get_fmpz_t(), a.get_fmpz_t(), m.get_fmpz_t()) != 0;
620 }
621 
622 inline void mp_gcd(fmpz_wrapper &res, const fmpz_wrapper &a,
623  const fmpz_wrapper &b)
624 {
625  fmpz_gcd(res.get_fmpz_t(), a.get_fmpz_t(), b.get_fmpz_t());
626 }
627 
628 inline void mp_gcdext(fmpz_wrapper &g, fmpz_wrapper &r, fmpz_wrapper &s,
629  const fmpz_wrapper &a, const fmpz_wrapper &b)
630 {
631  fmpz_xgcd(g.get_fmpz_t(), r.get_fmpz_t(), s.get_fmpz_t(), a.get_fmpz_t(),
632  b.get_fmpz_t());
633 }
634 
635 inline void mp_and(fmpz_wrapper &res, const fmpz_wrapper &a,
636  const fmpz_wrapper &b)
637 {
638  fmpz_and(res.get_fmpz_t(), a.get_fmpz_t(), b.get_fmpz_t());
639 }
640 
641 inline void mp_fdiv_r(fmpz_wrapper &res, const fmpz_wrapper &a,
642  const fmpz_wrapper &b)
643 {
644  fmpz_fdiv_r(res.get_fmpz_t(), a.get_fmpz_t(), b.get_fmpz_t());
645 }
646 
647 inline void mp_fdiv_q(fmpz_wrapper &res, const fmpz_wrapper &a,
648  const fmpz_wrapper &b)
649 {
650  fmpz_fdiv_q(res.get_fmpz_t(), a.get_fmpz_t(), b.get_fmpz_t());
651 }
652 
653 inline void mp_cdiv_q(fmpz_wrapper &res, const fmpz_wrapper &a,
654  const fmpz_wrapper &b)
655 {
656  fmpz_cdiv_q(res.get_fmpz_t(), a.get_fmpz_t(), b.get_fmpz_t());
657 }
658 
659 inline void mp_tdiv_q(fmpz_wrapper &res, const fmpz_wrapper &a,
660  const fmpz_wrapper &b)
661 {
662  fmpz_tdiv_q(res.get_fmpz_t(), a.get_fmpz_t(), b.get_fmpz_t());
663 }
664 
665 inline void mp_fdiv_qr(fmpz_wrapper &q, fmpz_wrapper &r, const fmpz_wrapper &a,
666  const fmpz_wrapper &b)
667 {
668  fmpz_fdiv_qr(q.get_fmpz_t(), r.get_fmpz_t(), a.get_fmpz_t(),
669  b.get_fmpz_t());
670 }
671 
672 inline void mp_divexact(fmpz_wrapper &q, const fmpz_wrapper &a,
673  const fmpz_wrapper &b)
674 {
675  fmpz_divexact(q.get_fmpz_t(), a.get_fmpz_t(), b.get_fmpz_t());
676 }
677 
678 inline void mp_lcm(fmpz_wrapper &q, const fmpz_wrapper &a,
679  const fmpz_wrapper &b)
680 {
681  fmpz_lcm(q.get_fmpz_t(), a.get_fmpz_t(), b.get_fmpz_t());
682 }
683 
684 inline void mp_tdiv_qr(fmpz_wrapper &q, fmpz_wrapper &r, const fmpz_wrapper &a,
685  const fmpz_wrapper &b)
686 {
687  fmpz_tdiv_qr(q.get_fmpz_t(), r.get_fmpz_t(), a.get_fmpz_t(),
688  b.get_fmpz_t());
689 }
690 
691 inline void mp_addmul(fmpz_wrapper &r, const fmpz_wrapper &a,
692  const fmpz_wrapper &b)
693 {
694  fmpz_addmul(r.get_fmpz_t(), a.get_fmpz_t(), b.get_fmpz_t());
695 }
696 
697 inline integer_class mp_primorial(unsigned long n)
698 {
699  fmpz_wrapper res;
700  fmpz_primorial(res.get_fmpz_t(), n);
701  return res;
702 }
703 
704 // helper functions for fmpq
705 
706 inline const fmpz_wrapper &get_den(const fmpq_wrapper &i)
707 {
708  return i.get_den();
709 }
710 
711 inline const fmpz_wrapper &get_num(const fmpq_wrapper &i)
712 {
713  return i.get_num();
714 }
715 
716 inline fmpz_wrapper &get_den(fmpq_wrapper &i)
717 {
718  return i.get_den();
719 }
720 
721 inline fmpz_wrapper &get_num(fmpq_wrapper &i)
722 {
723  return i.get_num();
724 }
725 
726 inline mpq_view_flint get_mpq_t(const fmpq_wrapper &i)
727 {
728  return mpq_view_flint(i);
729 }
730 
731 inline void canonicalize(fmpq_wrapper &i)
732 {
733  fmpq_canonicalise(i.get_fmpq_t());
734 }
735 
736 inline double mp_get_d(const fmpq_wrapper &i)
737 {
738  return mp_get_d(i.get_num()) / mp_get_d(i.get_den());
739 }
740 
741 inline int mp_sign(const fmpq_wrapper &i)
742 {
743  return fmpq_sgn(i.get_fmpq_t());
744 }
745 
746 inline fmpq_wrapper mp_abs(const fmpq_wrapper &i)
747 {
748  fmpq_wrapper res;
749  fmpq_abs(res.get_fmpq_t(), i.get_fmpq_t());
750  return res;
751 }
752 
753 #elif SYMENGINE_INTEGER_CLASS == SYMENGINE_BOOSTMP
754 
755 inline integer_class mp_abs(const integer_class &i)
756 {
757  // boost::multiprecision::abs(i) returns
758  // an expression template, not a cpp_int
759  // but it's ok: cpp_int is constructible from an expression template
760  return boost::multiprecision::abs(i);
761 }
762 
763 inline int mp_cmpabs(const integer_class &a, const integer_class &b)
764 {
765  if (mp_abs(a) > mp_abs(b)) {
766  return 1;
767  }
768  if (mp_abs(a) == mp_abs(b)) {
769  return 0;
770  }
771  return -1;
772 }
773 
774 inline int mp_sign(const integer_class &i)
775 {
776  return boost::math::sign(i);
777 }
778 
779 inline double mp_get_d(const integer_class &i)
780 {
781  return i.convert_to<double>();
782 }
783 
784 inline void mp_set_d(integer_class &i, double a)
785 {
786  i.assign(a);
787 }
788 
789 inline unsigned long mp_get_ui(const integer_class &i)
790 {
791  return mp_abs(i).convert_to<unsigned long>();
792 }
793 
794 inline long mp_get_si(const integer_class &i)
795 {
796  return i.convert_to<long>();
797 }
798 
799 inline bool mp_fits_ulong_p(const integer_class &i)
800 {
801  return (i >= 0) && (i <= ULONG_MAX);
802 }
803 
804 inline bool mp_fits_slong_p(const integer_class &i)
805 {
806  return (i >= LONG_MIN) && (i <= LONG_MAX);
807 }
808 
809 // bitwise and
810 inline void mp_and(integer_class &res, const integer_class &a,
811  const integer_class &b)
812 {
813  res = boost::multiprecision::operator&(a, b);
814 }
815 
816 inline void mp_pow_ui(integer_class &res, const integer_class &i,
817  unsigned long n)
818 {
819  res = boost::multiprecision::pow(i, numeric_cast<unsigned>(n));
820 }
821 
822 inline void mp_gcd(integer_class &res, const integer_class &a,
823  const integer_class &b)
824 {
825  res = boost::multiprecision::gcd(a, b);
826 }
827 
828 void mp_fdiv_qr(integer_class &q, integer_class &r, const integer_class &a,
829  const integer_class &b);
830 
831 void mp_cdiv_qr(integer_class &q, integer_class &r, const integer_class &a,
832  const integer_class &b);
833 
834 inline void mp_fdiv_r(integer_class &res, const integer_class &a,
835  const integer_class &b)
836 {
837  // TODO: benchmark this speed
838  integer_class quo;
839  mp_fdiv_qr(quo, res, a, b);
840 }
841 
842 inline void mp_fdiv_q(integer_class &res, const integer_class &a,
843  const integer_class &b)
844 {
845  // TODO: benchmark this speed
846  integer_class rem;
847  mp_fdiv_qr(res, rem, a, b);
848 }
849 
850 inline void mp_cdiv_q(integer_class &res, const integer_class &a,
851  const integer_class &b)
852 {
853  // TODO: benchmark this speed
854  integer_class rem;
855  mp_cdiv_qr(res, rem, a, b);
856 }
857 
858 inline void mp_tdiv_qr(integer_class &q, integer_class &r,
859  const integer_class &a, const integer_class &b)
860 {
861  boost::multiprecision::divide_qr(a, b, q, r);
862 }
863 
864 inline void mp_tdiv_q(integer_class &res, const integer_class &a,
865  const integer_class &b)
866 {
867  integer_class rem;
868  mp_tdiv_qr(res, rem, a, b);
869 }
870 
871 inline void mp_divexact(integer_class &q, const integer_class &a,
872  const integer_class &b)
873 {
874  // TODO: make faster
875  q = a / b;
876 }
877 
878 inline void mp_lcm(integer_class &q, const integer_class &a,
879  const integer_class &b)
880 {
881  q = boost::multiprecision::lcm(a, b);
882 }
883 
884 inline void mp_addmul(integer_class &r, const integer_class &a,
885  const integer_class &b)
886 {
887  // boost::multiprecision::default_ops::eval_multiply_add(r,a,b);
888  // //segfaults.
889  r += a * b;
890 }
891 
892 // Helper functions for cpp_rational
893 inline const integer_class get_den(const rational_class &i)
894 {
895  return boost::multiprecision::denominator(i);
896 }
897 
898 inline const integer_class get_num(const rational_class &i)
899 {
900  return boost::multiprecision::numerator(i);
901 }
902 
903 inline void canonicalize(rational_class &i)
904 {
905  // do nothing; boost::multiprecision::cpp_int
906  // is always stored in canonical form
907  // numerator and denominator share no common factors
908  // denominator is positive.
909 }
910 
911 inline double mp_get_d(const rational_class &i)
912 {
913  return i.convert_to<double>();
914 }
915 
916 inline int mp_sign(const rational_class &i)
917 {
918  return i.sign();
919 }
920 
921 inline rational_class mp_abs(const rational_class &i)
922 {
923  return boost::multiprecision::abs(i);
924 }
925 
926 inline bool mp_divisible_p(const integer_class &a, const integer_class &b)
927 {
928  return a % b == 0;
929 }
930 
931 void mp_pow_ui(rational_class &res, const rational_class &i, unsigned long n);
932 
933 void mp_powm(integer_class &res, const integer_class &a, const integer_class &b,
934  const integer_class &m);
935 
936 /* Extended Euclidean algorithm in Z
937  * inargs: integers a, b
938  * outargs: gcd, the greatest common divisor of a and b
939  * s, t such that sa + tb = gcd
940  */
941 void mp_gcdext(integer_class &gcd, integer_class &s, integer_class &t,
942  const integer_class &a, const integer_class &b);
943 
944 bool mp_invert(integer_class &res, const integer_class &a,
945  const integer_class &m);
946 
947 bool mp_root(integer_class &res, const integer_class &i, unsigned long n);
948 
949 integer_class mp_sqrt(const integer_class &i);
950 
951 void mp_rootrem(integer_class &a, integer_class &b, const integer_class &i,
952  unsigned long n);
953 
954 void mp_sqrtrem(integer_class &a, integer_class &b, const integer_class &i);
955 
956 int mp_probab_prime_p(const integer_class &i, unsigned retries);
957 
958 void mp_nextprime(integer_class &res, const integer_class &i);
959 
960 unsigned long mp_scan1(const integer_class &i);
961 
962 void mp_fib_ui(integer_class &res, unsigned long n);
963 
964 void mp_fib2_ui(integer_class &a, integer_class &b, unsigned long n);
965 
966 void mp_lucnum_ui(integer_class &res, unsigned long n);
967 
968 void mp_lucnum2_ui(integer_class &a, integer_class &b, unsigned long n);
969 
970 void mp_fac_ui(integer_class &res, unsigned long n);
971 
972 void mp_bin_ui(integer_class &res, const integer_class &n, unsigned long r);
973 
974 bool mp_perfect_power_p(const integer_class &i);
975 
976 bool mp_perfect_square_p(const integer_class &i);
977 
978 int mp_legendre(const integer_class &a, const integer_class &n);
979 
980 int mp_jacobi(const integer_class &a, const integer_class &n);
981 
982 int mp_kronecker(const integer_class &a, const integer_class &n);
983 
984 integer_class mp_primorial(unsigned long n);
985 
986 class mp_randstate
987 {
988 public:
989  // returns a uniformly distributed random integer between 0 and a-1,
990  // inclusive
991  void urandomint(integer_class &res, const integer_class &a)
992  {
993  boost::random::uniform_int_distribution<integer_class> ui(0, a);
994  res = ui(_twister);
995  }
996 
997  void seed(const uint32_t &i)
998  {
999  _twister.seed(i);
1000  }
1001 
1002 private:
1003  boost::random::mt19937 _twister;
1004 };
1005 
1006 #endif // SYMENGINE_INTEGER_CLASS == SYMENGINE_BOOSTMP
1007 
1008 // The implementation of each of the following
1009 // requires only get_mpz_t(integer_class), mp_demote
1010 // and functions from GMP. They don't depend on
1011 // any members of the backend class mpz, fmpz, or piranha::integer.
1012 // These functions will all need to be implemented separately
1013 // for a GMP-free build (with boost::multiprecision, for instance)
1014 #if SYMENGINE_INTEGER_CLASS == SYMENGINE_PIRANHA \
1015  || SYMENGINE_INTEGER_CLASS == SYMENGINE_FLINT \
1016  || SYMENGINE_INTEGER_CLASS == SYMENGINE_GMP \
1017  || SYMENGINE_INTEGER_CLASS == SYMENGINE_GMPXX
1018 
1019 inline bool mp_root(integer_class &res, const integer_class &i, unsigned long n)
1020 {
1021  auto _res = get_mpz_t(res);
1022  int ret = mpz_root(_res, get_mpz_t(i), n);
1023  mp_demote(res);
1024  return ret != 0;
1025 }
1026 
1027 inline void mp_nextprime(integer_class &res, const integer_class &i)
1028 {
1029  auto _res = get_mpz_t(res);
1030  mpz_nextprime(_res, get_mpz_t(i));
1031  mp_demote(res);
1032 }
1033 
1034 inline void mp_sqrtrem(integer_class &a, integer_class &b,
1035  const integer_class &i)
1036 {
1037  auto _a = get_mpz_t(a);
1038  auto _b = get_mpz_t(b);
1039  mpz_sqrtrem(_a, _b, get_mpz_t(i));
1040  mp_demote(a);
1041  mp_demote(b);
1042 }
1043 
1044 inline void mp_rootrem(integer_class &a, integer_class &b,
1045  const integer_class &i, unsigned long n)
1046 {
1047  auto _a = get_mpz_t(a);
1048  auto _b = get_mpz_t(b);
1049  mpz_rootrem(_a, _b, get_mpz_t(i), n);
1050  mp_demote(a);
1051  mp_demote(b);
1052 }
1053 
1054 inline unsigned long mp_scan1(const integer_class &i)
1055 {
1056  return mpz_scan1(get_mpz_t(i), 0);
1057 }
1058 
1059 inline void mp_fib_ui(integer_class &res, unsigned long n)
1060 {
1061  mpz_fib_ui(get_mpz_t(res), n);
1062  mp_demote(res);
1063 }
1064 
1065 inline void mp_fib2_ui(integer_class &a, integer_class &b, unsigned long n)
1066 {
1067  mpz_fib2_ui(get_mpz_t(a), get_mpz_t(b), n);
1068  mp_demote(a);
1069  mp_demote(b);
1070 }
1071 
1072 inline void mp_lucnum_ui(integer_class &res, unsigned long n)
1073 {
1074  mpz_lucnum_ui(get_mpz_t(res), n);
1075  mp_demote(res);
1076 }
1077 
1078 inline void mp_lucnum2_ui(integer_class &a, integer_class &b, unsigned long n)
1079 {
1080  mpz_lucnum2_ui(get_mpz_t(a), get_mpz_t(b), n);
1081  mp_demote(a);
1082  mp_demote(b);
1083 }
1084 
1085 inline void mp_bin_ui(integer_class &res, const integer_class &n,
1086  unsigned long r)
1087 {
1088  auto _res = get_mpz_t(res);
1089  mpz_bin_ui(_res, get_mpz_t(n), r);
1090  mp_demote(res);
1091 }
1092 
1093 inline void mp_fac_ui(integer_class &res, unsigned long n)
1094 {
1095  mpz_fac_ui(get_mpz_t(res), n);
1096  mp_demote(res);
1097 }
1098 
1099 inline int mp_legendre(const integer_class &a, const integer_class &n)
1100 {
1101  return mpz_legendre(get_mpz_t(a), get_mpz_t(n));
1102 }
1103 
1104 inline int mp_kronecker(const integer_class &a, const integer_class &n)
1105 {
1106  return mpz_kronecker(get_mpz_t(a), get_mpz_t(n));
1107 }
1108 
1109 inline int mp_jacobi(const integer_class &a, const integer_class &n)
1110 {
1111  return mpz_jacobi(get_mpz_t(a), get_mpz_t(n));
1112 }
1113 
1114 inline bool mp_perfect_power_p(const integer_class &i)
1115 {
1116  return mpz_perfect_power_p(get_mpz_t(i)) != 0;
1117 }
1118 
1119 inline bool mp_perfect_square_p(const integer_class &i)
1120 {
1121  return mpz_perfect_square_p(get_mpz_t(i)) != 0;
1122 }
1123 
1124 inline int mp_probab_prime_p(const integer_class &i, unsigned retries)
1125 {
1126  return mpz_probab_prime_p(get_mpz_t(i), retries);
1127 }
1128 
1129 inline bool mp_divisible_p(const integer_class &a, const integer_class &b)
1130 {
1131  return mpz_divisible_p(get_mpz_t(a), get_mpz_t(b)) != 0;
1132 }
1133 
1134 inline void mp_urandomm(integer_class &a, gmp_randstate_t &t,
1135  const integer_class &b)
1136 {
1137  auto _a = get_mpz_t(a);
1138  mpz_urandomm(_a, t, get_mpz_t(b));
1139  mp_demote(a);
1140 }
1141 
1142 inline auto get_mp_t(const integer_class &x) -> decltype(get_mpz_t(x))
1143 {
1144  return get_mpz_t(x);
1145 }
1146 
1147 inline auto get_mp_t(const rational_class &x) -> decltype(get_mpq_t(x))
1148 {
1149  return get_mpq_t(x);
1150 }
1151 
1152 inline int mp_cmpabs(const integer_class &a, const integer_class &b)
1153 {
1154  return mpz_cmpabs(get_mpz_t(a), get_mpz_t(b));
1155 }
1156 
1158 {
1159 public:
1160  // returns a uniformly distributed random integer between 0 and a-1,
1161  // inclusive
1162  void urandomint(integer_class &res, const integer_class &a)
1163  {
1164  mpz_urandomm(get_mpz_t(res), _state, get_mpz_t(a));
1165  mp_demote(res);
1166  }
1167 
1168  void seed(const integer_class &i)
1169  {
1170  gmp_randseed(_state, get_mpz_t(i));
1171  }
1172 
1173  mp_randstate()
1174  {
1175  gmp_randinit_default(_state);
1176  gmp_randseed_ui(_state, std::rand());
1177  }
1178 
1179  ~mp_randstate()
1180  {
1181  gmp_randclear(_state);
1182  }
1183 
1184 private:
1185  gmp_randstate_t _state;
1186 };
1187 
1188 #endif // SYMENGINE_INTEGER_CLASS == Piranha or Flint or GMP or GMPXX
1189 
1190 } // namespace SymEngine
1191 
1192 #if !defined(HAVE_SYMENGINE_GMP) && defined(HAVE_SYMENGINE_BOOST) \
1193  && BOOST_VERSION < 105900
1194 namespace boost
1195 {
1196 namespace detail
1197 {
1198 template <>
1199 struct make_unsigned_imp<SymEngine::integer_class> {
1200  typedef SymEngine::integer_class type;
1201 };
1202 } // namespace detail
1203 } // namespace boost
1204 #endif
1205 
1206 #endif // SYMENGINE_INTEGER_CLASS_H
Main namespace for SymEngine package.
Definition: add.cpp:19
std::enable_if< std::is_integral< T >::value, RCP< const Integer > >::type integer(T i)
Definition: integer.h:200
RCP< const Integer > gcd(const Integer &a, const Integer &b)
Greatest Common Divisor.
Definition: ntheory.cpp:29
RCP< const Number > rational(long n, long d)
convenience creator from two longs
Definition: rational.h:328
T rand(T... args)