rational.cpp
1 #include <symengine/rational.h>
2 #include <symengine/pow.h>
3 #include <symengine/symengine_exception.h>
4 
5 namespace SymEngine
6 {
7 
8 bool 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 
23 RCP<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 
34 RCP<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 
44 RCP<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 
68  return Rational::from_mpq(std::move(q));
69 }
70 
71 RCP<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 
96 hash_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 
106 bool 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 
115 int 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 
130 void 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 
137 bool 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 
158 bool 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 
189 RCP<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 
195 RCP<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:95
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:48
RCP< const Integer > neg() const
Definition: integer.h:119
RCP< const Number > powint(const Integer &other) const
Fast Power Evaluation.
Definition: integer.h:105
virtual bool is_negative() const
Definition: integer.h:73
static RCP< const Basic > from_dict(const RCP< const Number > &coef, map_basic_basic &&d)
Create a Mul from a dict.
Definition: mul.cpp:116
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
const rational_class & as_rational_class() const
Convert to rational_class.
Definition: rational.h:50
virtual int compare(const Basic &o) const
Definition: rational.cpp:115
rational_class i
i : object of rational_class
Definition: rational.h:19
virtual bool __eq__(const Basic &o) const
Definition: rational.cpp:106
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
RCP< const Number > powrat(const Integer &other) const
Definition: rational.h:193
virtual hash_t __hash__() const
Definition: rational.cpp:96
T move(T... args)
Main namespace for SymEngine package.
Definition: add.cpp:19
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
std::enable_if< std::is_integral< T >::value, RCP< const Integer > >::type integer(T i)
Definition: integer.h:200
RCP< const Basic > exp(const RCP< const Basic > &x)
Returns the natural exponential function E**x = pow(E, x)
Definition: pow.cpp:270
RCP< const Basic > mul(const RCP< const Basic > &a, const RCP< const Basic > &b)
Multiplication.
Definition: mul.cpp:347
void insert(T1 &m, const T2 &first, const T3 &second)
Definition: dict.h:83
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 > neg(const RCP< const Basic > &a)
Negation.
Definition: mul.cpp:438