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