Searching...
No Matches
rational.cpp
2#include <symengine/pow.h>
3#include <symengine/symengine_exception.h>
4
5namespace SymEngine
6{
7
8bool Rational::is_canonical(const rational_class &i) const
9{
10 rational_class x = i;
11 canonicalize(x);
12 // If 'x' is an integer, it should not be Rational:
13 if (SymEngine::get_den(x) == 1)
14 return false;
15 // if 'i' is not in canonical form:
16 if (SymEngine::get_num(x) != SymEngine::get_num(i))
17 return false;
18 if (SymEngine::get_den(x) != SymEngine::get_den(i))
19 return false;
20 return true;
21}
22
23RCP<const Number> Rational::from_mpq(const rational_class &i)
24{
25 // If the result is an Integer, return an Integer:
26 if (SymEngine::get_den(i) == 1) {
27 return integer(SymEngine::get_num(i));
28 } else {
29 rational_class j(i);
30 return make_rcp<const Rational>(std::move(j));
31 }
32}
33
34RCP<const Number> Rational::from_mpq(rational_class &&i)
35{
36 // If the result is an Integer, return an Integer:
37 if (SymEngine::get_den(i) == 1) {
38 return integer(SymEngine::get_num(i));
39 } else {
40 return make_rcp<const Rational>(std::move(i));
41 }
42}
43
44RCP<const Number> Rational::from_two_ints(const Integer &n, const Integer &d)
45{
46 if (d.as_integer_class() == 0) {
47 if (n.as_integer_class() == 0) {
48 return Nan;
49 } else {
50 return ComplexInf;
51 }
52 }
53#if SYMENGINE_INTEGER_CLASS == SYMENGINE_BOOSTMP
54 // workaround https://github.com/boostorg/rational/issues/27
55 rational_class q;
56 if (d.as_integer_class() < 0) {
57 q = rational_class(-n.as_integer_class(), -d.as_integer_class());
58 } else {
59 q = rational_class(n.as_integer_class(), d.as_integer_class());
60 }
61#else
62 rational_class q(n.as_integer_class(), d.as_integer_class());
63#endif
64 // This is potentially slow, but has to be done, since 'n/d' might not be
65 // in canonical form.
66 canonicalize(q);
67
69}
70
71RCP<const Number> Rational::from_two_ints(long n, long d)
72{
73 if (d == 0) {
74 if (n == 0) {
75 return Nan;
76 } else {
77 return ComplexInf;
78 }
79 }
80#if SYMENGINE_INTEGER_CLASS == SYMENGINE_BOOSTMP
81 // workaround https://github.com/boostorg/rational/issues/27
82 if (d < 0) {
83 d *= -1;
84 n *= -1;
85 }
86#endif
87 rational_class q(n, d);
88
89 // This is potentially slow, but has to be done, since 'n/d' might not be
90 // in canonical form.
91 canonicalize(q);
92
93 return Rational::from_mpq(q);
94}
95
96hash_t Rational::__hash__() const
97{
98 // only the least significant bits that fit into "signed long int" are
99 // hashed:
100 hash_t seed = SYMENGINE_RATIONAL;
101 hash_combine<long long int>(seed, mp_get_si(SymEngine::get_num(this->i)));
102 hash_combine<long long int>(seed, mp_get_si(SymEngine::get_den(this->i)));
103 return seed;
104}
105
106bool Rational::__eq__(const Basic &o) const
107{
108 if (is_a<Rational>(o)) {
109 const Rational &s = down_cast<const Rational &>(o);
110 return this->i == s.i;
111 }
112 return false;
113}
114
115int Rational::compare(const Basic &o) const
116{
117 if (is_a<Rational>(o)) {
118 const Rational &s = down_cast<const Rational &>(o);
119 if (i == s.i)
120 return 0;
121 return i < s.i ? -1 : 1;
122 }
123 if (is_a<Integer>(o)) {
124 const Integer &s = down_cast<const Integer &>(o);
125 return i < s.as_integer_class() ? -1 : 1;
126 }
127 throw NotImplementedError("unhandled comparison of Rational");
128}
129
130void get_num_den(const Rational &rat, const Ptr<RCP<const Integer>> &num,
131 const Ptr<RCP<const Integer>> &den)
132{
133 *num = integer(SymEngine::get_num(rat.as_rational_class()));
134 *den = integer(SymEngine::get_den(rat.as_rational_class()));
135}
136
137bool Rational::is_perfect_power(bool is_expected) const
138{
139 const integer_class &num = SymEngine::get_num(i);
140 if (num == 1)
141 return mp_perfect_power_p(SymEngine::get_den(i));
142
143 const integer_class &den = SymEngine::get_den(i);
144 // TODO: fix this
145 if (not is_expected) {
146 if (mp_cmpabs(num, den) > 0) {
147 if (!mp_perfect_power_p(den))
148 return false;
149 } else {
150 if (!mp_perfect_power_p(num))
151 return false;
152 }
153 }
154 integer_class prod = num * den;
155 return mp_perfect_power_p(prod);
156}
157
158bool Rational::nth_root(const Ptr<RCP<const Number>> &the_rat,
159 unsigned long n) const
160{
161 if (n == 0)
162 throw SymEngineException("i_nth_root: Can not find Zeroth root");
163
164#if SYMENGINE_INTEGER_CLASS != SYMENGINE_BOOSTMP
165 rational_class r;
166 int ret = mp_root(SymEngine::get_num(r), SymEngine::get_num(i), n);
167 if (ret == 0)
168 return false;
169 ret = mp_root(SymEngine::get_den(r), SymEngine::get_den(i), n);
170 if (ret == 0)
171 return false;
172#else
173 // boost::multiprecision::cpp_rational doesn't provide
174 // non-const get_num and get_den
175 integer_class num, den;
176 int ret = mp_root(num, SymEngine::get_num(i), n);
177 if (ret == 0)
178 return false;
179 ret = mp_root(den, SymEngine::get_den(i), n);
180 if (ret == 0)
181 return false;
182 rational_class r(num, den);
183#endif
184 // No need to canonicalize since this is in canonical form
185 *the_rat = make_rcp<const Rational>(std::move(r));
186 return true;
187}
188
189RCP<const Basic> Rational::powrat(const Rational &other) const
190{
191 return SymEngine::mul(other.rpowrat(*this->get_num()),
192 other.neg()->rpowrat(*this->get_den()));
193}
194
195RCP<const Basic> Rational::rpowrat(const Integer &other) const
196{
197 if (not(mp_fits_ulong_p(SymEngine::get_den(i))))
198 throw SymEngineException("powrat: den of 'exp' does not fit ulong.");
199 unsigned long exp = mp_get_ui(SymEngine::get_den(i));
200 RCP<const Integer> res;
201 if (other.is_negative()) {
202 if (i_nth_root(outArg(res), *other.neg(), exp)) {
203 if (exp % 2 == 0) {
204 return I->pow(*get_num())->mul(*res->powint(*get_num()));
205 } else {
206 return SymEngine::neg(res->powint(*get_num()));
207 }
208 }
209 } else {
210 if (i_nth_root(outArg(res), other, exp)) {
211 return res->powint(*get_num());
212 }
213 }
214 integer_class q, r;
215 auto num = SymEngine::get_num(i);
216 auto den = SymEngine::get_den(i);
217
218 mp_fdiv_qr(q, r, num, den);
219 // Here we make the exponent postive and a fraction between
220 // 0 and 1. We multiply numerator and denominator appropriately
221 // to achieve this
222 RCP<const Number> coef = other.powint(*integer(q));
223 map_basic_basic surd;
224
225 if ((other.is_negative()) and den == 2) {
226 imulnum(outArg(coef), I);
227 // if other.neg() is one, no need to add it to dict
228 if (other.as_integer_class() != -1)
229 insert(surd, other.neg(),
230 Rational::from_mpq(rational_class(r, den)));
231 } else {
232 insert(surd, other.rcp_from_this(),
233 Rational::from_mpq(rational_class(r, den)));
234 }
235 return Mul::from_dict(coef, std::move(surd));
236}
237
238} // namespace SymEngine
The lowest unit of symbolic representation.
Definition: basic.h:97
RCP< T > rcp_from_this()
Get RCP<T> pointer to self (it will cast the pointer to T)
Integer Class.
Definition: integer.h:19
const integer_class & as_integer_class() const
Convert to integer_class.
Definition: integer.h:45
RCP< const Integer > neg() const
Definition: integer.h:116
RCP< const Number > powint(const Integer &other) const
Fast Power Evaluation.
Definition: integer.h:102
bool is_negative() const override
Definition: integer.h:70
static RCP< const Basic > from_dict(const RCP< const Number > &coef, map_basic_basic &&d)
Create a Mul from a dict.
Definition: mul.cpp:115
Rational Class.
Definition: rational.h:16
RCP< const Basic > rpowrat(const Integer &other) const
Definition: rational.cpp:195
static RCP< const Number > from_mpq(const rational_class &i)
Definition: rational.cpp:23
hash_t __hash__() const override
Definition: rational.cpp:96
int compare(const Basic &o) const override
Definition: rational.cpp:115
rational_class i
i : object of rational_class
Definition: rational.h:19
RCP< const Number > powrat(const Integer &other) const
Definition: rational.h:193
const rational_class & as_rational_class() const
Convert to rational_class.
Definition: rational.h:50
static RCP< const Number > from_two_ints(const Integer &n, const Integer &d)
Definition: rational.cpp:44
bool is_canonical(const rational_class &i) const
Definition: rational.cpp:8
RCP< const Rational > neg() const
Definition: rational.h:92
bool __eq__(const Basic &o) const override
Definition: rational.cpp:106
T move(T... args)
Main namespace for SymEngine package.
int i_nth_root(const Ptr< RCP< const Integer > > &r, const Integer &a, unsigned long int n)
Integer nth root.
Definition: integer.cpp:128
RCP< const Basic > exp(const RCP< const Basic > &x)
Returns the natural exponential function E**x = pow(E, x)
Definition: pow.cpp:271
void get_num_den(const Rational &rat, const Ptr< RCP< const Integer > > &num, const Ptr< RCP< const Integer > > &den)
returns the num and den of rational rat as RCP<const Integer>
Definition: rational.cpp:130
RCP< const Basic > mul(const RCP< const Basic > &a, const RCP< const Basic > &b)
Multiplication.
Definition: mul.cpp:352
void insert(T1 &m, const T2 &first, const T3 &second)
Definition: dict.h:83
std::enable_if< std::is_integral< T >::value, RCP< constInteger > >::type integer(T i)
Definition: integer.h:197
RCP< const Basic > neg(const RCP< const Basic > &a)
Negation.
Definition: mul.cpp:443