Loading...
Searching...
No Matches
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
23namespace 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
50typedef boost::multiprecision::number<boost::multiprecision::cpp_int_backend<>,
51 boost::multiprecision::et_off>
52 integer_class;
53typedef boost::multiprecision::number<
54 boost::multiprecision::cpp_rational_backend, boost::multiprecision::et_off>
55 rational_class;
56#elif SYMENGINE_INTEGER_CLASS == SYMENGINE_PIRANHA
57typedef piranha::integer integer_class;
58typedef piranha::rational rational_class;
59#elif SYMENGINE_INTEGER_CLASS == SYMENGINE_FLINT
60typedef fmpz_wrapper integer_class;
61typedef fmpq_wrapper rational_class;
62#elif SYMENGINE_INTEGER_CLASS == SYMENGINE_GMP
63typedef mpz_wrapper integer_class;
64typedef mpq_wrapper rational_class;
65#elif SYMENGINE_INTEGER_CLASS == SYMENGINE_GMPXX
66typedef mpz_class integer_class;
67typedef 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;
72inline namespace literals
73{
75inline integer_class operator"" _z(const char *str)
76{
77 return integer_class(str);
78}
79
80inline 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
89inline 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
96inline int mp_sign(const integer_class &i)
97{
98 return mpz_sgn(i.get_mpz_t());
99}
100
101inline 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
108inline double mp_get_d(const integer_class &i)
109{
110 return static_cast<double>(i.get_d());
111}
112
113inline void mp_set_d(integer_class &i, double a)
114{
115 mpz_set_d(i.get_mpz_t(), a);
116}
117
118inline void mp_demote(integer_class &i) {}
119
120inline bool mp_fits_ulong_p(const integer_class &i)
121{
122 return i.fits_ulong_p() != 0;
123}
124
125inline bool mp_fits_slong_p(const integer_class &i)
126{
127 return i.fits_slong_p() != 0;
128}
129
130inline unsigned long mp_get_ui(const integer_class &i)
131{
132 return i.get_ui();
133}
134
135inline long mp_get_si(const integer_class &i)
136{
137 return i.get_si();
138}
139
140inline mpz_srcptr get_mpz_t(const integer_class &i)
141{
142 return i.get_mpz_t();
143}
144
145inline mpz_ptr get_mpz_t(integer_class &i)
146{
147 return i.get_mpz_t();
148}
149
150inline 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
156inline 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
163inline 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
169inline 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
175inline 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
181inline 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
188inline 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
194inline 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
200inline 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
206inline 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
212inline 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
218inline 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
224inline 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
230inline 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
236inline 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
242inline 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
249inline const integer_class &get_den(const rational_class &i)
250{
251 return i.get_den();
252}
253
254inline const integer_class &get_num(const rational_class &i)
255{
256 return i.get_num();
257}
258
259inline integer_class &get_den(rational_class &i)
260{
261 return i.get_den();
262}
263
264inline integer_class &get_num(rational_class &i)
265{
266 return i.get_num();
267}
268
269inline mpq_srcptr get_mpq_t(const rational_class &i)
270{
271 return i.get_mpq_t();
272}
273
274inline void canonicalize(rational_class &i)
275{
276 i.canonicalize();
277}
278
279inline double mp_get_d(const rational_class &i)
280{
281 return i.get_d();
282}
283
284inline int mp_sign(const rational_class &i)
285{
286 return mpq_sgn(i.get_mpq_t());
287}
288
289inline 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
296inline 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
305inline piranha::integer mp_abs(const piranha::integer &i)
306{
307 return i.abs();
308}
309
310inline piranha::integer mp_sqrt(const piranha::integer &i)
311{
312 return i.sqrt();
313}
314
315inline void mp_demote(piranha::integer &i) {}
316
317inline mpz_ptr get_mpz_t(piranha::integer &i)
318{
319 return i._get_mpz_ptr();
320}
321
322inline auto get_mpz_t(const piranha::integer &i) -> decltype(i.get_mpz_view())
323{
324 return i.get_mpz_view();
325}
326
327inline void mp_pow_ui(piranha::integer &res, const piranha::integer &i,
328 unsigned long n)
329{
330 res = i.pow(n);
331}
332
333inline void mp_pow_ui(piranha::rational &res, const piranha::rational &i,
334 unsigned long n)
335{
336 res = i.pow(n);
337}
338
339inline 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
346inline 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
353inline 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
359inline 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
369inline 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
376inline 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
383inline 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
390inline 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
397inline 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
404inline 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
411inline 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
417inline 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
424inline 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
430inline int mp_sign(const piranha::integer &i)
431{
432 return i.sign();
433}
434
435inline long mp_get_si(const piranha::integer &i)
436{
437 return mpz_get_si(i.get_mpz_view());
438}
439
440inline unsigned long mp_get_ui(const piranha::integer &i)
441{
442 return mpz_get_ui(i.get_mpz_view());
443}
444
445inline double mp_get_d(const piranha::integer &i)
446{
447 return mpz_get_d(i.get_mpz_view());
448}
449
450inline void mp_set_d(piranha::integer &i, double a)
451{
452 i = a;
453}
454
455inline bool mp_fits_ulong_p(const piranha::integer &i)
456{
457 return mpz_fits_ulong_p(i.get_mpz_view()) != 0;
458}
459
460inline bool mp_fits_slong_p(const piranha::integer &i)
461{
462 return mpz_fits_slong_p(i.get_mpz_view()) != 0;
463}
464
465inline 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
471inline 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
480inline piranha::rational mp_abs(const piranha::rational &i)
481{
482 return i.abs();
483}
484
485inline const piranha::integer &get_den(const piranha::rational &i)
486{
487 return i.den();
488}
489
490inline const piranha::integer &get_num(const piranha::rational &i)
491{
492 return i.num();
493}
494
495inline piranha::integer &get_den(piranha::rational &i)
496{
497 return i._den();
498}
499
500inline piranha::integer &get_num(piranha::rational &i)
501{
502 return i._num();
503}
504
505inline void canonicalize(piranha::rational &i)
506{
507 i.canonicalise();
508}
509
510inline double mp_get_d(const piranha::rational &i)
511{
512 return mpq_get_d(i.get_mpq_view().get());
513}
514
515inline auto get_mpq_t(const piranha::rational &i) -> decltype(i.get_mpq_view())
516{
517 return i.get_mpq_view();
518}
519
520inline 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
528inline mpz_view_flint get_mpz_t(const fmpz_wrapper &i)
529{
530 return mpz_view_flint(i);
531}
532
533inline mpz_ptr get_mpz_t(fmpz_wrapper &i)
534{
535 return _fmpz_promote_val(i.get_fmpz_t());
536}
537
538inline void mp_demote(fmpz_wrapper &i)
539{
540 _fmpz_demote_val(i.get_fmpz_t());
541}
542
543inline int mp_sign(const fmpz_wrapper &i)
544{
545 return fmpz_sgn(i.get_fmpz_t());
546}
547
548inline long mp_get_si(const fmpz_wrapper &i)
549{
550 return fmpz_get_si(i.get_fmpz_t());
551}
552
553inline unsigned long mp_get_ui(const fmpz_wrapper &i)
554{
555 return fmpz_get_ui(i.get_fmpz_t());
556}
557
558inline bool mp_fits_slong_p(const fmpz_wrapper &i)
559{
560 return fmpz_fits_si(i.get_fmpz_t());
561}
562
563inline 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
568inline double mp_get_d(const fmpz_wrapper &i)
569{
570 return fmpz_get_d(i.get_fmpz_t());
571}
572
573inline void mp_set_d(fmpz_wrapper &i, double a)
574{
575 return fmpz_set_d(i.get_fmpz_t(), a);
576}
577
578inline 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
585inline 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
592inline 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
597inline 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
602inline 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
616inline 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
622inline 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
628inline 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
635inline 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
641inline 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
647inline 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
653inline 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
659inline 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
665inline 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
672inline 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
678inline 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
684inline 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
691inline 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
697inline 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
706inline const fmpz_wrapper &get_den(const fmpq_wrapper &i)
707{
708 return i.get_den();
709}
710
711inline const fmpz_wrapper &get_num(const fmpq_wrapper &i)
712{
713 return i.get_num();
714}
715
716inline fmpz_wrapper &get_den(fmpq_wrapper &i)
717{
718 return i.get_den();
719}
720
721inline fmpz_wrapper &get_num(fmpq_wrapper &i)
722{
723 return i.get_num();
724}
725
726inline mpq_view_flint get_mpq_t(const fmpq_wrapper &i)
727{
728 return mpq_view_flint(i);
729}
730
731inline void canonicalize(fmpq_wrapper &i)
732{
733 fmpq_canonicalise(i.get_fmpq_t());
734}
735
736inline 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
741inline int mp_sign(const fmpq_wrapper &i)
742{
743 return fmpq_sgn(i.get_fmpq_t());
744}
745
746inline 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
755inline 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
763inline 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
774inline int mp_sign(const integer_class &i)
775{
776 return boost::math::sign(i);
777}
778
779inline double mp_get_d(const integer_class &i)
780{
781 return i.convert_to<double>();
782}
783
784inline void mp_set_d(integer_class &i, double a)
785{
786 i.assign(a);
787}
788
789inline unsigned long mp_get_ui(const integer_class &i)
790{
791 return mp_abs(i).convert_to<unsigned long>();
792}
793
794inline long mp_get_si(const integer_class &i)
795{
796 return i.convert_to<long>();
797}
798
799inline bool mp_fits_ulong_p(const integer_class &i)
800{
801 return (i >= 0) && (i <= ULONG_MAX);
802}
803
804inline bool mp_fits_slong_p(const integer_class &i)
805{
806 return (i >= LONG_MIN) && (i <= LONG_MAX);
807}
808
809// bitwise and
810inline 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
816inline 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
822inline 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
828void mp_fdiv_qr(integer_class &q, integer_class &r, const integer_class &a,
829 const integer_class &b);
830
831void mp_cdiv_qr(integer_class &q, integer_class &r, const integer_class &a,
832 const integer_class &b);
833
834inline 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
842inline 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
850inline 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
858inline 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
864inline 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
871inline 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
878inline 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
884inline 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
893inline const integer_class get_den(const rational_class &i)
894{
895 return boost::multiprecision::denominator(i);
896}
897
898inline const integer_class get_num(const rational_class &i)
899{
900 return boost::multiprecision::numerator(i);
901}
902
903inline 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
911inline double mp_get_d(const rational_class &i)
912{
913 return i.convert_to<double>();
914}
915
916inline int mp_sign(const rational_class &i)
917{
918 return i.sign();
919}
920
921inline rational_class mp_abs(const rational_class &i)
922{
923 return boost::multiprecision::abs(i);
924}
925
926inline bool mp_divisible_p(const integer_class &a, const integer_class &b)
927{
928 return a % b == 0;
929}
930
931void mp_pow_ui(rational_class &res, const rational_class &i, unsigned long n);
932
933void 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 */
941void mp_gcdext(integer_class &gcd, integer_class &s, integer_class &t,
942 const integer_class &a, const integer_class &b);
943
944bool mp_invert(integer_class &res, const integer_class &a,
945 const integer_class &m);
946
947bool mp_root(integer_class &res, const integer_class &i, unsigned long n);
948
949integer_class mp_sqrt(const integer_class &i);
950
951void mp_rootrem(integer_class &a, integer_class &b, const integer_class &i,
952 unsigned long n);
953
954void mp_sqrtrem(integer_class &a, integer_class &b, const integer_class &i);
955
956int mp_probab_prime_p(const integer_class &i, unsigned retries);
957
958void mp_nextprime(integer_class &res, const integer_class &i);
959
960unsigned long mp_scan1(const integer_class &i);
961
962void mp_fib_ui(integer_class &res, unsigned long n);
963
964void mp_fib2_ui(integer_class &a, integer_class &b, unsigned long n);
965
966void mp_lucnum_ui(integer_class &res, unsigned long n);
967
968void mp_lucnum2_ui(integer_class &a, integer_class &b, unsigned long n);
969
970void mp_fac_ui(integer_class &res, unsigned long n);
971
972void mp_bin_ui(integer_class &res, const integer_class &n, unsigned long r);
973
974bool mp_perfect_power_p(const integer_class &i);
975
976bool mp_perfect_square_p(const integer_class &i);
977
978int mp_legendre(const integer_class &a, const integer_class &n);
979
980int mp_jacobi(const integer_class &a, const integer_class &n);
981
982int mp_kronecker(const integer_class &a, const integer_class &n);
983
984integer_class mp_primorial(unsigned long n);
985
986class mp_randstate
987{
988public:
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
1002private:
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
1019inline 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
1027inline 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
1034inline 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
1044inline 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
1054inline unsigned long mp_scan1(const integer_class &i)
1055{
1056 return mpz_scan1(get_mpz_t(i), 0);
1057}
1058
1059inline 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
1065inline 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
1072inline 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
1078inline 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
1085inline 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
1093inline 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
1099inline 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
1104inline 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
1109inline 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
1114inline bool mp_perfect_power_p(const integer_class &i)
1115{
1116 return mpz_perfect_power_p(get_mpz_t(i)) != 0;
1117}
1118
1119inline bool mp_perfect_square_p(const integer_class &i)
1120{
1121 return mpz_perfect_square_p(get_mpz_t(i)) != 0;
1122}
1123
1124inline 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
1129inline 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
1134inline 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
1142inline auto get_mp_t(const integer_class &x) -> decltype(get_mpz_t(x))
1143{
1144 return get_mpz_t(x);
1145}
1146
1147inline auto get_mp_t(const rational_class &x) -> decltype(get_mpq_t(x))
1148{
1149 return get_mpq_t(x);
1150}
1151
1152inline 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{
1159public:
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
1180 {
1181 gmp_randclear(_state);
1182 }
1183
1184private:
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
1194namespace boost
1195{
1196namespace detail
1197{
1198template <>
1199struct 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
RCP< const Integer > gcd(const Integer &a, const Integer &b)
Greatest Common Divisor.
Definition: ntheory.cpp:31
T rand(T... args)