carl  24.04
Computer ARithmetic Library
operations.cpp
Go to the documentation of this file.
1 #include "../numbers.h"
2 #include <limits>
3 
4 #include "parser.h"
5 
6 namespace carl
7 {
8 
9  bool sqrt_exact(const mpq_class& a, mpq_class& b)
10  {
11  if( mpq_sgn(a.__get_mp()) < 0 ) return false;
12  mpz_class den = a.get_den();
13  mpz_class num = a.get_num();
14  mpz_class root_den;
15  mpz_class root_den_rem;
16  mpz_sqrtrem(root_den.__get_mp(), root_den_rem.__get_mp(), den.__get_mp());
17  if( !carl::is_zero( root_den_rem ) )
18  return false;
19 
20  mpz_class root_num;
21  mpz_class root_num_rem;
22  mpz_sqrtrem(root_num.__get_mp(), root_num_rem.__get_mp(), num.__get_mp());
23  if( !carl::is_zero( root_num_rem ) )
24  return false;
25 
26  mpq_class resNum;
27  mpq_set_z(resNum.get_mpq_t(), root_num.get_mpz_t());
28  mpq_class resDen;
29  mpq_set_z(resDen.get_mpq_t(), root_den.get_mpz_t());
30  mpq_div(b.get_mpq_t(), resNum.get_mpq_t(), resDen.get_mpq_t());
31  return true;
32  }
33 
34  mpq_class sqrt(const mpq_class& a) {
35  auto r = sqrt_safe(a);
36  return (r.first + r.second) / 2;
37  }
38 
39  std::pair<mpq_class,mpq_class> sqrt_safe(const mpq_class& a)
40  {
41  if (mpq_sgn(a.__get_mp()) == 0) return std::make_pair(mpq_class(0),mpq_class(0));
42  assert( mpq_sgn(a.__get_mp()) > 0 );
43  mpz_class den = a.get_den();
44  mpz_class num = a.get_num();
45  mpz_class root_den;
46  mpz_class root_den_rem;
47  mpz_sqrtrem(root_den.__get_mp(), root_den_rem.__get_mp(), den.__get_mp());
48 
49  mpz_class root_num;
50  mpz_class root_num_rem;
51  mpz_sqrtrem(root_num.__get_mp(), root_num_rem.__get_mp(), num.__get_mp());
52 
53  mpq_class lower;
54  mpq_class upper;
55 
56  lower = root_num;
57  if(root_den_rem == 0)
58  lower /= root_den;
59  else
60  lower /= root_den+1;
61 
62  if(root_num_rem == 0)
63  upper = root_num;
64  else
65  upper = root_num+1;
66 
67  upper /= root_den;
68 
69  return std::make_pair(lower,upper);
70  }
71 
72  std::pair<mpq_class,mpq_class> root_safe(const mpq_class& a, uint n)
73  {
74  assert(mpq_sgn(a.__get_mp()) > 0);
75  mpz_class den = a.get_den();
76  mpz_class num = a.get_num();
77  mpz_class root_den;
78  int den_exact = mpz_root(root_den.__get_mp(), den.__get_mp(), n);
79 
80  mpz_class root_num;
81  int num_exact = mpz_root(root_num.__get_mp(), num.__get_mp(), n);
82 
83  mpq_class lower;
84  mpq_class upper;
85 
86  lower = root_num;
87  if(den_exact)
88  lower /= root_den;
89  else
90  lower /= root_den+1;
91 
92  if(num_exact)
93  upper = root_num;
94  else
95  upper = root_num+1;
96 
97  upper /= root_den;
98  return std::make_pair(lower,upper);
99  }
100 
101  std::pair<mpq_class, mpq_class> sqrt_fast(const mpq_class& a)
102  {
103  assert(a >= 0);
104 #if 1
105  mpz_class num;
106  mpz_class num_rem;
107  mpz_sqrtrem(num.__get_mp(), num_rem.__get_mp(), a.get_num().__get_mp());
108  mpz_class den;
109  mpz_class den_rem;
110  mpz_sqrtrem(den.__get_mp(), den_rem.__get_mp(), a.get_den().__get_mp());
111 
112  if (carl::is_zero(num_rem)) {
113  if (carl::is_zero(den_rem)) {
114  mpq_class exact_root = mpq_class(num) / den;
115  return std::make_pair(exact_root, exact_root);
116  } else {
117  mpq_class lower = mpq_class(num) / (den + 1);
118  mpq_class upper = mpq_class(num) / den;
119  return std::make_pair(lower, upper);
120  }
121  } else {
122  if (carl::is_zero(den_rem)) {
123  mpq_class lower = mpq_class(num) / den;
124  mpq_class upper = (mpq_class(num) + 1) / den;
125  return std::make_pair(lower, upper);
126  } else {
127  mpq_class lower = mpq_class(num) / (den + 1);
128  mpq_class upper = (mpq_class(num) + 1) / den;
129  return std::make_pair(lower, upper);
130  }
131  }
132 #else
133  mpq_class exact_root;
134  if (carl::sqrtp(a, exact_root)) {
135  // root can be computed exactly.
136  return std::make_pair(exact_root, exact_root);
137  } else {
138  // compute an approximation with sqrt(). we can assume that the surrounding integers contain the actual root.
139  mpf_class af = sqrt(mpf_class(a));
140  mpq_class lower(af - carl::constant_one<mpf_class>::get());
141  mpq_class upper(af + carl::constant_one<mpf_class>::get());
142  assert(lower * lower < a);
143  assert(upper * upper > a);
144  return std::make_pair(lower, upper);
145  }
146 #endif
147  }
148 
149  template<>
150  mpz_class parse<mpz_class>(const std::string& n) {
151  mpz_class res;
152  bool success = parser::parseDecimal(n, res);
153  assert(success);
154  return res;
155  }
156 
157  template<>
158  bool try_parse<mpz_class>(const std::string& n, mpz_class& res) {
159  return parser::parseDecimal(n, res);
160  }
161 
162  template<>
163  mpq_class parse<mpq_class>(const std::string& n) {
164  mpq_class res;
165  bool success = parser::parseRational(n, res);
166  assert(success);
167  return res;
168  }
169 
170  template<>
171  bool try_parse<mpq_class>(const std::string& n, mpq_class& res) {
172  return parser::parseRational(n, res);
173  }
174 
175  std::string toString(const mpq_class& _number, bool _infix)
176  {
177  std::stringstream s;
178  bool negative = (_number < mpq_class(0));
179  if(negative) s << "(-" << (_infix ? "" : " ");
180  if(_infix) s << carl::abs(_number);
181  else
182  {
183  mpz_class d = carl::get_denom(_number);
184  if(d != mpz_class(1)) s << "(/ " << carl::abs(carl::get_num(_number)) << " " << carl::abs(d) << ")";
185  else s << carl::abs(_number);
186  }
187  if(negative)
188  s << ")";
189  return s.str();
190  }
191 
192  std::string toString(const mpz_class& _number, bool _infix)
193  {
194  std::stringstream s;
195  bool negative = (_number < mpz_class(0));
196  if(negative) s << "(-" << (_infix ? "" : " ");
197  s << carl::abs(_number);
198  if(negative)
199  s << ")";
200  return s.str();
201  }
202 }
carl
carl is the main namespace for the library.
carl::uint
std::uint64_t uint
Definition: numbers.h:16
carl::sqrt_fast
std::pair< cln::cl_RA, cln::cl_RA > sqrt_fast(const cln::cl_RA &a)
Compute square root in a fast but less precise way.
carl::toString
std::string toString(Relation r)
Definition: Relation.h:73
carl::abs
Interval< Number > abs(const Interval< Number > &_in)
Method which returns the absolute value of the passed number.
Definition: Interval.h:1511
carl::sqrt_exact
bool sqrt_exact(const cln::cl_RA &a, cln::cl_RA &b)
Calculate the square root of a fraction if possible.
carl::root_safe
std::pair< cln::cl_RA, cln::cl_RA > root_safe(const cln::cl_RA &a, uint n)
carl::is_zero
bool is_zero(const Interval< Number > &i)
Check if this interval is a point-interval containing 0.
Definition: Interval.h:1453
carl::get_num
cln::cl_I get_num(const cln::cl_RA &n)
Extract the numerator from a fraction.
Definition: operations.h:60
carl::sqrt
Interval< Number > sqrt(const Interval< Number > &i)
Definition: Power.h:51
carl::try_parse< mpq_class >
bool try_parse< mpq_class >(const std::string &n, mpq_class &res)
Definition: operations.cpp:171
carl::get_denom
cln::cl_I get_denom(const cln::cl_RA &n)
Extract the denominator from a fraction.
Definition: operations.h:69
carl::parse< mpz_class >
mpz_class parse< mpz_class >(const std::string &n)
Definition: operations.cpp:150
carl::sqrt_safe
std::pair< cln::cl_RA, cln::cl_RA > sqrt_safe(const cln::cl_RA &a)
Calculate the square root of a fraction.
carl::try_parse< mpz_class >
bool try_parse< mpz_class >(const std::string &n, mpz_class &res)
Definition: operations.cpp:158
carl::parse< mpq_class >
mpq_class parse< mpq_class >(const std::string &n)
Definition: operations.cpp:163
carl::parser::parseDecimal
bool parseDecimal(const std::string &input, T &output)
Definition: parser.h:146
carl::parser::parseRational
bool parseRational(const std::string &input, T &output)
Definition: parser.h:150