carl  24.04
Computer ARithmetic Library
mpfr_float.tpp
Go to the documentation of this file.
1 #include "FLOAT_T.h"
2 
3 
4 #ifdef USE_MPFR_FLOAT
5 namespace carl {
6 
7 
8 template<>
9 class FLOAT_T<mpfr_t>
10 {
11  private:
12  mpfr_t mValue;
13  static precision_t mDefaultPrecision;
14 
15  public:
16 
17  /**
18  * Constructors & Destructors
19  */
20 
21  FLOAT_T()
22  {
23  mpfr_init2(mValue, mDefaultPrecision);
24  mpfr_set_zero(mValue, 1);
25  }
26 
27  // Default precision is initially set to 53 bits in mpfr implementation
28  FLOAT_T(const double _double, const CARL_RND _rnd=CARL_RND::N, precision_t _prec=mDefaultPrecision)
29  {
30  mpfr_init2(mValue,_prec);
31  mpfr_set_d(mValue,_double,mpfr_rnd_t(_rnd));
32  }
33 
34  // Default precision is initially set to 53 bits in mpfr implementation
35  FLOAT_T(const float _float, const CARL_RND _rnd=CARL_RND::N, precision_t _prec=mDefaultPrecision)
36  {
37  mpfr_init2(mValue,_prec);
38  mpfr_set_flt(mValue, _float, mpfr_rnd_t(_rnd));
39  }
40 
41  // Default precision is initially set to 53 bits in mpfr implementation
42  FLOAT_T(const int _int, const CARL_RND _rnd=CARL_RND::N, precision_t _prec=mDefaultPrecision)
43  {
44  mpfr_init2(mValue,_prec);
45  mpfr_set_si(mValue,_int,mpfr_rnd_t(_rnd));
46  }
47 
48  // Default precision is initially set to 53 bits in mpfr implementation
49  FLOAT_T(const unsigned _int, const CARL_RND _rnd=CARL_RND::N, precision_t _prec=mDefaultPrecision)
50  {
51  mpfr_init2(mValue,_prec);
52  mpfr_set_ui(mValue,_int,mpfr_rnd_t(_rnd));
53  }
54 
55  // Default precision is initially set to 53 bits in mpfr implementation
56 
57  FLOAT_T(const long _long, const CARL_RND _rnd=CARL_RND::N, precision_t _prec=mDefaultPrecision)
58  {
59  mpfr_init2(mValue,_prec);
60  mpfr_set_si(mValue,_long,mpfr_rnd_t(_rnd));
61  }
62 
63  // Default precision is initially set to 53 bits in mpfr implementation
64  FLOAT_T(const unsigned long _long, const CARL_RND _rnd=CARL_RND::N, precision_t _prec=mDefaultPrecision)
65  {
66  mpfr_init2(mValue,_prec);
67  mpfr_set_ui(mValue,_long,mpfr_rnd_t(_rnd));
68  }
69 
70  FLOAT_T(const mpfr_t& _mpfrNumber)
71  {
72  mpfr_init2(mValue,mpfr_get_prec(_mpfrNumber));
73  mpfr_set(mValue, _mpfrNumber, MPFR_RNDN);
74  }
75 
76  FLOAT_T(const FLOAT_T<mpfr_t>& _float)
77  {
78  mpfr_init2(mValue, mpfr_get_prec(_float.mValue));
79  mpfr_set(mValue, _float.mValue, MPFR_RNDN);
80  }
81 
82  FLOAT_T(FLOAT_T<mpfr_t>&& _float)
83  {
84  if(this->mValue == _float.value()){
85  mpfr_swap(mValue,_float.mValue);
86  } else {
87  mpfr_init2(mValue, mpfr_get_prec(_float.value()));
88  mpfr_swap(mValue,_float.mValue);
89  }
90  }
91 
92  FLOAT_T(const std::string& _string)
93  {
94  mpfr_init_set_str(mValue, _string.c_str(), 10, MPFR_RNDN);
95  }
96 
97  template<typename F, DisableIf< std::is_same<F, mpfr_t> > = dummy>
98  FLOAT_T(const FLOAT_T<F>& _float, const CARL_RND _rnd=CARL_RND::N, precision_t _prec=mDefaultPrecision)
99  {
100  mpfr_init2(mValue,_prec);
101  mpfr_set_d(mValue,_float.to_double(),mpfr_rnd_t(_rnd));
102  }
103 
104  ~FLOAT_T()
105  {
106  mpfr_clear(mValue);
107  }
108 
109  static inline const FLOAT_T<mpfr_t> const_infinity(){
110  mpfr_t tmp;
111  mpfr_init2(tmp, FLOAT_T<mpfr_t>::defaultPrecision());
112  mpfr_set_inf(tmp,1);
113  return FLOAT_T<mpfr_t>(tmp);
114  }
115 
116  inline const FLOAT_T<mpfr_t> setNan(){
117  mpfr_set_nan(mValue);
118  return *this;
119  }
120 
121  inline FLOAT_T<mpfr_t> increment() const {
122  FLOAT_T<mpfr_t> tmp(*this);
123  mpfr_nextabove(tmp.mValue);
124  return tmp;
125  }
126 
127  inline FLOAT_T<mpfr_t> decrement() const {
128  FLOAT_T<mpfr_t> tmp(*this);
129  mpfr_nextbelow(tmp.mValue);
130  return tmp;
131  }
132 
133  static inline FLOAT_T<mpfr_t> machine_epsilon(precision_t prec)
134  {
135  /* the smallest eps such that 1 + eps != 1 */
136  return machine_epsilon(FLOAT_T<mpfr_t>(1,CARL_RND::N, prec));
137  }
138 
139  static inline FLOAT_T<mpfr_t> machine_epsilon(const FLOAT_T<mpfr_t>& x)
140  {
141  /* the smallest eps such that x + eps != x */
142  FLOAT_T<mpfr_t> tmp(x);
143  if( x < 0)
144  {
145  return -tmp.increment() + tmp;
146  }else{
147  return tmp.increment() - tmp;
148  }
149  }
150 
151  static inline carl::FLOAT_T<mpfr_t> minval(const CARL_RND _rnd=CARL_RND::N, precision_t prec=FLOAT_T<mpfr_t>::defaultPrecision()) {
152  /* min = 1/2 * 2^emin = 2^(emin - 1) */
153  mpfr_t tmp;
154  mpfr_init2(tmp,prec);
155  mpfr_set_ui(tmp,1,mpfr_rnd_t(_rnd));
156  mpfr_mul_2ui(tmp,tmp,(unsigned)mpfr_get_emin()-1,mpfr_rnd_t(_rnd));
157  return FLOAT_T<mpfr_t>(tmp);
158  }
159 
160  static inline FLOAT_T<mpfr_t> maxval(const CARL_RND _rnd=CARL_RND::N, precision_t prec=FLOAT_T<mpfr_t>::defaultPrecision())
161  {
162  /* max = (1 - eps) * 2^emax, eps is machine epsilon */
163  mpfr_t tmp;
164  mpfr_init2(tmp,prec);
165  mpfr_set_ui(tmp,1, mpfr_rnd_t(_rnd));
166  mpfr_sub(tmp,tmp,machine_epsilon(prec).value(), mpfr_rnd_t(_rnd));
167  mpfr_mul_2ui(tmp,tmp,(unsigned)mpfr_get_emax(),mpfr_rnd_t(_rnd));
168  return (FLOAT_T<mpfr_t>(tmp));
169  }
170 
171  /*******************
172  * Getter & Setter *
173  *******************/
174 
175  const mpfr_t& value() const
176  {
177  return mValue;
178  }
179 
180  static void setDefaultPrecision(precision_t _prec) {
181  mDefaultPrecision = _prec;
182  }
183 
184  static precision_t defaultPrecision () {
185  return mDefaultPrecision;
186  }
187 
188  precision_t precision() const
189  {
190  return mpfr_get_prec(mValue);
191  }
192 
193  FLOAT_T<mpfr_t>& setPrecision( const precision_t& _prec, const CARL_RND _rnd=CARL_RND::N )
194  {
195  mpfr_prec_round(mValue, convPrec(_prec), mpfr_rnd_t(_rnd));
196  return *this;
197  }
198 
199  /*************
200  * Operators *
201  *************/
202 
203  FLOAT_T<mpfr_t>& operator = (const FLOAT_T<mpfr_t>& _rhs)
204  {
205  if(this->mValue == _rhs.value())
206  return *this;
207 
208  // Note: This is a workaround to get the limb size correct. Instead use free and reallocate.
209  mpfr_set_prec(mValue, mpfr_get_prec(_rhs.mValue));
210  mpfr_set(mValue, _rhs.mValue, MPFR_RNDN);
211 
212  /*
213  mValue->_mpfr_prec = _rhs.mValue->_mpfr_prec;
214  mValue->_mpfr_sign = _rhs.mValue->_mpfr_sign;
215  mValue->_mpfr_exp = _rhs.mValue->_mpfr_exp;
216 
217  // deep-copy limbs
218  std::size_t limbs = (std::size_t)std::ceil(double(_rhs.mValue->_mpfr_prec)/double(mp_bits_per_limb));
219 
220  while( limbs > 0 ){
221  mValue->_mpfr_d[limbs-1] = _rhs.mValue->_mpfr_d[limbs-1];
222  --limbs;
223  }
224  */
225  return *this;
226  }
227 
228  FLOAT_T<mpfr_t>& safeSet (const FLOAT_T<mpfr_t>& _rhs, const CARL_RND _rnd=CARL_RND::N)
229  {
230  mpfr_set_prec(mValue, mpfr_get_prec(_rhs.mValue));
231  mpfr_set(mValue, _rhs.mValue, mpfr_rnd_t(_rnd));
232  return *this;
233  }
234 
235  /**
236  * Boolean operators
237  */
238 
239  bool operator == ( const FLOAT_T<mpfr_t>& _rhs) const
240  {
241  return mpfr_cmp(mValue,_rhs.mValue) == 0;
242  }
243 
244  bool operator != ( const FLOAT_T<mpfr_t> & _rhs) const
245  {
246  return mpfr_equal_p(mValue, _rhs.mValue) == 0;
247  }
248 
249  bool operator > ( const FLOAT_T<mpfr_t> & _rhs) const
250  {
251  return mpfr_greater_p(mValue, _rhs.mValue) != 0;
252  }
253 
254  bool operator > ( int _rhs) const
255  {
256  return mpfr_cmp_si(mValue, _rhs) > 0;
257  }
258 
259  bool operator > ( unsigned _rhs) const
260  {
261  return mpfr_cmp_ui(mValue, _rhs) > 0;
262  }
263 
264  bool operator < ( const FLOAT_T<mpfr_t> & _rhs) const
265  {
266  return mpfr_less_p(mValue, _rhs.mValue) != 0;
267  }
268 
269  bool operator < ( int _rhs) const
270  {
271  return mpfr_cmp_si(mValue, _rhs) < 0;
272  }
273 
274  bool operator < ( unsigned _rhs) const
275  {
276  return mpfr_cmp_ui(mValue, _rhs) < 0;
277  }
278 
279  bool operator <= ( const FLOAT_T<mpfr_t> & _rhs) const
280  {
281  return mpfr_lessequal_p(mValue, _rhs.mValue) != 0;
282  }
283 
284  bool operator >= ( const FLOAT_T<mpfr_t> & _rhs) const
285  {
286  return mpfr_greaterequal_p(mValue, _rhs.mValue) != 0;
287  }
288 
289  /**
290  * arithmetic operations
291  */
292 
293  FLOAT_T<mpfr_t>& add_assign( const FLOAT_T<mpfr_t>& _op2, CARL_RND _rnd=CARL_RND::N )
294  {
295  mpfr_add(mValue, mValue, _op2.mValue, mpfr_rnd_t(_rnd));
296  return *this;
297  }
298 
299  void add( FLOAT_T<mpfr_t>& _result, const FLOAT_T<mpfr_t>& _op2, CARL_RND _rnd=CARL_RND::N ) const
300  {
301  mpfr_add(_result.mValue, this->mValue, _op2.mValue, mpfr_rnd_t(_rnd));
302  }
303 
304  FLOAT_T<mpfr_t>& sub_assign( const FLOAT_T<mpfr_t>& _op2, CARL_RND _rnd=CARL_RND::N )
305  {
306  mpfr_sub(mValue, mValue, _op2.mValue, mpfr_rnd_t(_rnd));
307  return *this;
308  }
309 
310  void sub( FLOAT_T<mpfr_t>& _result, const FLOAT_T<mpfr_t>& _op2, CARL_RND _rnd=CARL_RND::N ) const
311  {
312  mpfr_sub(_result.mValue, this->mValue, _op2.mValue, mpfr_rnd_t(_rnd));
313  }
314 
315  FLOAT_T<mpfr_t>& mul_assign(const FLOAT_T<mpfr_t>& _op2, CARL_RND _rnd=CARL_RND::N)
316  {
317  mpfr_mul(mValue, mValue, _op2.mValue, mpfr_rnd_t(_rnd));
318  return *this;
319  }
320 
321  void mul( FLOAT_T<mpfr_t>& _result, const FLOAT_T<mpfr_t>& _op2, CARL_RND _rnd=CARL_RND::N) const
322  {
323  mpfr_mul(_result.mValue, this->mValue, _op2.mValue, mpfr_rnd_t(_rnd));
324  }
325 
326  FLOAT_T<mpfr_t>& div_assign(const FLOAT_T<mpfr_t>& _op2, CARL_RND _rnd=CARL_RND::N)
327  {
328  assert( _op2 != 0 );
329  mpfr_div(mValue, mValue, _op2.mValue, mpfr_rnd_t(_rnd));
330  return *this;
331  }
332 
333  void div( FLOAT_T<mpfr_t>& _result, const FLOAT_T<mpfr_t>& _op2, CARL_RND _rnd=CARL_RND::N) const
334  {
335  assert( _op2 != 0 );
336  mpfr_div(_result.mValue, this->mValue, _op2.mValue, mpfr_rnd_t(_rnd));
337  }
338 
339  /**
340  * special operators
341  */
342 
343  FLOAT_T<mpfr_t>& sqrt_assign(CARL_RND _rnd = CARL_RND::N)
344  {
345  assert( *this >= 0);
346  mpfr_sqrt(mValue, mValue, mpfr_rnd_t(_rnd));
347  return *this;
348  }
349 
350  void sqrt(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
351  {
352  assert( *this >= 0);
353  mpfr_sqrt(_result.mValue, mValue, mpfr_rnd_t(_rnd));
354  }
355 
356  FLOAT_T<mpfr_t>& cbrt_assign(CARL_RND _rnd = CARL_RND::N)
357  {
358  assert( *this >= 0);
359  mpfr_cbrt(mValue, mValue, mpfr_rnd_t(_rnd));
360  return *this;
361  }
362 
363  void cbrt(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
364  {
365  assert( *this >= 0);
366  mpfr_cbrt(_result.mValue, mValue, mpfr_rnd_t(_rnd));
367  }
368 
369  FLOAT_T<mpfr_t>& root_assign(unsigned long int _k, CARL_RND _rnd = CARL_RND::N)
370  {
371  assert( *this >= 0);
372  mpfr_root(mValue, mValue, _k, mpfr_rnd_t(_rnd));
373  return *this;
374  }
375 
376  void root(FLOAT_T<mpfr_t>& _result, unsigned long int _k, CARL_RND _rnd = CARL_RND::N) const
377  {
378  assert( *this >= 0);
379  mpfr_root(_result.mValue, mValue, _k, mpfr_rnd_t(_rnd));
380  }
381 
382  FLOAT_T<mpfr_t>& pow_assign(unsigned long int _exp, CARL_RND _rnd = CARL_RND::N)
383  {
384  mpfr_pow_ui(mValue, mValue, _exp, mpfr_rnd_t(_rnd));
385  return *this;
386  }
387 
388  void pow(FLOAT_T<mpfr_t>& _result, unsigned long int _exp, CARL_RND _rnd = CARL_RND::N) const
389  {
390  mpfr_pow_ui(_result.mValue, mValue, _exp, mpfr_rnd_t(_rnd));
391  }
392 
393  FLOAT_T<mpfr_t>& abs_assign(CARL_RND _rnd = CARL_RND::N)
394  {
395  mpfr_abs(mValue, mValue, mpfr_rnd_t(_rnd));
396  return *this;
397  }
398 
399  void abs(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
400  {
401  mpfr_abs(_result.mValue, mValue, mpfr_rnd_t(_rnd));
402  }
403 
404  FLOAT_T<mpfr_t>& exp_assign(CARL_RND _rnd = CARL_RND::N)
405  {
406  mpfr_exp(mValue, mValue, mpfr_rnd_t(_rnd));
407  return *this;
408  }
409 
410  void exp(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
411  {
412  mpfr_exp(_result.mValue, mValue, mpfr_rnd_t(_rnd));
413  }
414 
415  FLOAT_T<mpfr_t>& sin_assign(CARL_RND _rnd = CARL_RND::N)
416  {
417  mpfr_sin(mValue, mValue, mpfr_rnd_t(_rnd));
418  return *this;
419  }
420 
421  void sin(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
422  {
423  mpfr_sin(_result.mValue, mValue, mpfr_rnd_t(_rnd));
424  }
425 
426  FLOAT_T<mpfr_t>& cos_assign(CARL_RND _rnd = CARL_RND::N)
427  {
428  mpfr_cos(mValue, mValue, mpfr_rnd_t(_rnd));
429  return *this;
430  }
431 
432  void cos(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
433  {
434  mpfr_cos(_result.mValue, mValue, mpfr_rnd_t(_rnd));
435  }
436 
437  FLOAT_T<mpfr_t>& log_assign(CARL_RND _rnd = CARL_RND::N)
438  {
439  mpfr_log(mValue, mValue, mpfr_rnd_t(_rnd));
440  return *this;
441  }
442 
443  void log(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
444  {
445  mpfr_log(_result.mValue, mValue, mpfr_rnd_t(_rnd));
446  }
447 
448  FLOAT_T<mpfr_t>& tan_assign(CARL_RND _rnd = CARL_RND::N)
449  {
450  mpfr_tan(mValue, mValue, mpfr_rnd_t(_rnd));
451  return *this;
452  }
453 
454  void tan(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
455  {
456  mpfr_tan(_result.mValue, mValue, mpfr_rnd_t(_rnd));
457  }
458 
459  FLOAT_T<mpfr_t>& asin_assign(CARL_RND _rnd = CARL_RND::N)
460  {
461  mpfr_asin(mValue, mValue, mpfr_rnd_t(_rnd));
462  return *this;
463  }
464 
465  void asin(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
466  {
467  mpfr_asin(_result.mValue, mValue, mpfr_rnd_t(_rnd));
468  }
469 
470  FLOAT_T<mpfr_t>& acos_assign(CARL_RND _rnd = CARL_RND::N)
471  {
472  mpfr_acos(mValue, mValue, mpfr_rnd_t(_rnd));
473  return *this;
474  }
475 
476  void acos(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
477  {
478  mpfr_acos(_result.mValue, mValue, mpfr_rnd_t(_rnd));
479  }
480 
481  FLOAT_T<mpfr_t>& atan_assign(CARL_RND _rnd = CARL_RND::N)
482  {
483  mpfr_atan(mValue, mValue, mpfr_rnd_t(_rnd));
484  return *this;
485  }
486 
487  void atan(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
488  {
489  mpfr_atan(_result.mValue, mValue, mpfr_rnd_t(_rnd));
490  }
491 
492  FLOAT_T<mpfr_t>& sinh_assign(CARL_RND _rnd = CARL_RND::N)
493  {
494  mpfr_sinh(mValue, mValue, mpfr_rnd_t(_rnd));
495  return *this;
496  }
497 
498  void sinh(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
499  {
500  mpfr_sinh(_result.mValue, mValue, mpfr_rnd_t(_rnd));
501  }
502 
503  FLOAT_T<mpfr_t>& cosh_assign(CARL_RND _rnd = CARL_RND::N)
504  {
505  mpfr_cosh(mValue, mValue, mpfr_rnd_t(_rnd));
506  return *this;
507  }
508 
509  void cosh(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
510  {
511  mpfr_cosh(_result.mValue, mValue, mpfr_rnd_t(_rnd));
512  }
513 
514  FLOAT_T<mpfr_t>& tanh_assign(CARL_RND _rnd = CARL_RND::N)
515  {
516  mpfr_tanh(mValue, mValue, mpfr_rnd_t(_rnd));
517  return *this;
518  }
519 
520  void tanh(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
521  {
522  mpfr_tanh(_result.mValue, mValue, mpfr_rnd_t(_rnd));
523  }
524 
525  FLOAT_T<mpfr_t>& asinh_assign(CARL_RND _rnd = CARL_RND::N)
526  {
527  mpfr_asinh(mValue, mValue, mpfr_rnd_t(_rnd));
528  return *this;
529  }
530 
531  void asinh(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
532  {
533  mpfr_asinh(_result.mValue, mValue, mpfr_rnd_t(_rnd));
534  }
535 
536  FLOAT_T<mpfr_t>& acosh_assign(CARL_RND _rnd = CARL_RND::N)
537  {
538  mpfr_acosh(mValue, mValue, mpfr_rnd_t(_rnd));
539  return *this;
540  }
541 
542  void acosh(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
543  {
544  mpfr_acosh(_result.mValue, mValue, mpfr_rnd_t(_rnd));
545  }
546 
547  FLOAT_T<mpfr_t>& atanh_assign(CARL_RND _rnd = CARL_RND::N)
548  {
549  mpfr_atanh(mValue, mValue, mpfr_rnd_t(_rnd));
550  return *this;
551  }
552 
553  void atanh(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
554  {
555  mpfr_atanh(_result.mValue, mValue, mpfr_rnd_t(_rnd));
556  }
557 
558  void floor(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
559  {
560  mpfr_init(_result.mValue);
561  mpfr_set_zero(_result.mValue, 1);
562  mpfr_rint_floor(_result.mValue, mValue, mpfr_rnd_t(_rnd));
563  //_result = mpfr_integer_p(result);
564  //mpfr_clear(result);
565  }
566 
567  FLOAT_T& floor_assign(CARL_RND _rnd = CARL_RND::N)
568  {
569  mpfr_rint_floor(mValue, mValue, mpfr_rnd_t(_rnd));
570  return *this;
571  }
572 
573  void ceil(FLOAT_T<mpfr_t>& _result, CARL_RND _rnd = CARL_RND::N) const
574  {
575  mpfr_init(_result.mValue);
576  mpfr_set_zero(_result.mValue, 1);
577  mpfr_rint_ceil(_result.mValue, mValue, mpfr_rnd_t(_rnd));
578  //_result = mpfr_integer_p(result);
579  //mpfr_clear(result);
580  }
581 
582  FLOAT_T& ceil_assign(CARL_RND _rnd = CARL_RND::N)
583  {
584  mpfr_rint_ceil(mValue, mValue, mpfr_rnd_t(_rnd));
585  return *this;
586  }
587 
588  /**
589  * conversion operators
590  */
591 
592  double to_double(CARL_RND _rnd=CARL_RND::N) const
593  {
594  return mpfr_get_d(mValue, mpfr_rnd_t(_rnd));
595  }
596 
597  /**
598  * Explicit typecast operator to integer.
599  * @return Integer representation of this.
600  */
601  explicit operator int() const
602  {
603  return (int)this->to_double();
604  }
605 
606  /**
607  * Explicit typecast operator to long.
608  * @return Long representation of this.
609  */
610  explicit operator long() const
611  {
612  return (long)(this->to_double());
613  }
614 
615  /**
616  * Explicit typecast operator to double.
617  * @return Double representation of this.
618  */
619  explicit operator double() const
620  {
621  return this->to_double();
622  }
623 
624 
625  friend std::ostream & operator<< (std::ostream& ostr, const FLOAT_T<mpfr_t> & p) {
626  ostr << p.toString();
627  return ostr;
628  }
629 
630 
631 
632  friend bool operator== (const FLOAT_T<mpfr_t>& _lhs, const int _rhs)
633  {
634  return mpfr_cmp_si(_lhs.mValue, _rhs) == 0;
635  }
636 
637  friend bool operator== (const int _lhs, const FLOAT_T<mpfr_t>& _rhs)
638  {
639  return _rhs == _lhs;
640  }
641 
642  friend bool operator== (const FLOAT_T<mpfr_t>& _lhs, const double _rhs)
643  {
644  return mpfr_cmp_d(_lhs.mValue,_rhs) == 0;
645  }
646 
647  friend bool operator== (const double _lhs, const FLOAT_T<mpfr_t>& _rhs)
648  {
649  return _rhs == _lhs;
650  }
651 
652  friend bool operator== (const FLOAT_T<mpfr_t>& _lhs, const float _rhs)
653  {
654  return mpfr_cmp_d(_lhs.mValue, _rhs);
655  }
656 
657  friend bool operator== (const float _lhs, const FLOAT_T<mpfr_t>& _rhs)
658  {
659  return _rhs == _lhs;
660  }
661 
662  /**
663  * Operators
664  */
665 
666  friend FLOAT_T<mpfr_t> operator +(const FLOAT_T<mpfr_t>& _lhs, const FLOAT_T<mpfr_t>& _rhs)
667  {
668  FLOAT_T<mpfr_t> res;
669  mpfr_add(res.mValue, _lhs.mValue, _rhs.mValue, MPFR_RNDN);
670  return res;
671  }
672 
673  friend FLOAT_T<mpfr_t> operator +(const mpfr_t& _lhs, const FLOAT_T<mpfr_t>& _rhs)
674  {
675  FLOAT_T<mpfr_t> res;
676  mpfr_add(res.mValue, _lhs, _rhs.mValue, MPFR_RNDN);
677  return res;
678  }
679 
680  friend FLOAT_T<mpfr_t> operator +(const FLOAT_T<mpfr_t>& _lhs, const mpfr_t& _rhs)
681  {
682  FLOAT_T<mpfr_t> res;
683  mpfr_add(res.mValue, _lhs.mValue, _rhs, MPFR_RNDN);
684  return res;
685  }
686 
687  template<typename Other, EnableIf< is_number_type<Other>, Not<is_interval_type<Other>> > = dummy>
688  friend FLOAT_T<mpfr_t> operator +(const Other& _lhs, const FLOAT_T<mpfr_t>& _rhs)
689  {
690  FLOAT_T<mpfr_t> res;
691  mpfr_add(res.mValue, FLOAT_T<mpfr_t>(_lhs).mValue, _rhs.mValue, MPFR_RNDN);
692  return res;
693  }
694 
695  template<typename Other, EnableIf< is_number_type<Other>, Not<is_interval_type<Other>> > = dummy>
696  friend FLOAT_T<mpfr_t> operator +(const FLOAT_T<mpfr_t>& _lhs, const Other& _rhs)
697  {
698  return _rhs+_lhs;
699  }
700 
701  friend FLOAT_T<mpfr_t> operator -(const FLOAT_T<mpfr_t>& _lhs, const FLOAT_T<mpfr_t>& _rhs)
702  {
703  FLOAT_T<mpfr_t> res;
704  mpfr_sub(res.mValue, _lhs.mValue, _rhs.mValue, MPFR_RNDN);
705  return res;
706  }
707 
708  friend FLOAT_T<mpfr_t> operator -(const mpfr_t& _lhs, const FLOAT_T<mpfr_t>& _rhs)
709  {
710  FLOAT_T<mpfr_t> res;
711  mpfr_sub(res.mValue, _lhs, _rhs.mValue, MPFR_RNDN);
712  return res;
713  }
714 
715  friend FLOAT_T<mpfr_t> operator -(const FLOAT_T<mpfr_t>& _lhs, const mpfr_t& _rhs)
716  {
717  FLOAT_T<mpfr_t> res;
718  mpfr_sub(res.mValue, _lhs.mValue, _rhs, MPFR_RNDN);
719  return res;
720  }
721 
722  template<typename Other, EnableIf< is_number_type<Other>, Not<is_interval_type<Other>> > = dummy>
723  friend FLOAT_T<mpfr_t> operator -(const Other& _lhs, const FLOAT_T<mpfr_t>& _rhs)
724  {
725  FLOAT_T<mpfr_t> res;
726  mpfr_sub(res.mValue, FLOAT_T<mpfr_t>(_lhs).mValue, _rhs.mValue, MPFR_RNDN);
727  return res;
728  }
729 
730  template<typename Other, EnableIf< is_number_type<Other>, Not<is_interval_type<Other>> > = dummy>
731  friend FLOAT_T<mpfr_t> operator -(const FLOAT_T<mpfr_t>& _lhs, const Other& _rhs)
732  {
733  FLOAT_T<mpfr_t> res;
734  mpfr_sub(res.mValue, _lhs.mValue, FLOAT_T<mpfr_t>(_rhs).mValue, MPFR_RNDN);
735  return res;
736  }
737 
738  friend FLOAT_T<mpfr_t> operator -(const FLOAT_T<mpfr_t>& _lhs)
739  {
740  FLOAT_T<mpfr_t> res;
741  mpfr_mul_si(res.mValue, _lhs.mValue, -1, MPFR_RNDN);
742  return res;
743  }
744 
745  friend FLOAT_T<mpfr_t> operator *(const FLOAT_T<mpfr_t>& _lhs, const FLOAT_T<mpfr_t>& _rhs)
746  {
747  FLOAT_T<mpfr_t> res;
748  mpfr_mul(res.mValue, _lhs.mValue, _rhs.mValue, MPFR_RNDN);
749  return res;
750  }
751 
752  template<typename Other, EnableIf< is_number_type<Other>, Not<is_interval_type<Other>> > = dummy>
753  friend FLOAT_T<mpfr_t> operator *(const Other& _lhs, const FLOAT_T<mpfr_t>& _rhs)
754  {
755  FLOAT_T<mpfr_t> res;
756  FLOAT_T<mpfr_t> lhs = FLOAT_T<mpfr_t>(_lhs);
757  mpfr_mul(res.mValue, lhs.mValue, _rhs.mValue, MPFR_RNDN);
758  return res;
759  }
760 
761  template<typename Other, EnableIf< is_number_type<Other>, Not<is_interval_type<Other>> > = dummy>
762  friend FLOAT_T<mpfr_t> operator *(const FLOAT_T<mpfr_t>& _lhs, const Other& _rhs)
763  {
764  return _rhs*_lhs;
765  }
766 
767  friend FLOAT_T<mpfr_t> operator *(const mpfr_t& _lhs, const FLOAT_T<mpfr_t>& _rhs)
768  {
769  FLOAT_T<mpfr_t> res;
770  mpfr_mul(res.mValue, _lhs, _rhs.mValue, MPFR_RNDN);
771  return res;
772  }
773 
774  friend FLOAT_T<mpfr_t> operator *(const FLOAT_T<mpfr_t>& _lhs, const mpfr_t& _rhs)
775  {
776  FLOAT_T<mpfr_t> res;
777  mpfr_mul(res.mValue, _lhs.mValue, _rhs, MPFR_RNDN);
778  return res;
779  }
780 
781  friend FLOAT_T<mpfr_t> operator /(const FLOAT_T<mpfr_t>& _lhs, const FLOAT_T<mpfr_t>& _rhs)
782  {
783  // TODO: mpfr_div results in infty when dividing by zero, although this should not be defined.
784  FLOAT_T<mpfr_t> res;
785  mpfr_div(res.mValue, _lhs.mValue, _rhs.mValue, MPFR_RNDN);
786  return res;
787  }
788 
789  friend FLOAT_T<mpfr_t> operator /(const mpfr_t& _lhs, const FLOAT_T<mpfr_t>& _rhs)
790  {
791  // TODO: mpfr_div results in infty when dividing by zero, although this should not be defined.
792  FLOAT_T<mpfr_t> res;
793  mpfr_div(res.mValue, _lhs, _rhs.mValue, MPFR_RNDN);
794  return res;
795  }
796 
797  friend FLOAT_T<mpfr_t> operator /(const FLOAT_T<mpfr_t>& _lhs, const mpfr_t& _rhs)
798  {
799  // TODO: mpfr_div results in infty when dividing by zero, although this should not be defined.
800  FLOAT_T<mpfr_t> res;
801  mpfr_div(res.mValue, _lhs.mValue, _rhs, MPFR_RNDN);
802  return res;
803  }
804 
805  template<typename Other, EnableIf< is_number_type<Other>, Not<is_interval_type<Other>> > = dummy>
806  friend FLOAT_T<mpfr_t> operator /(const Other& _lhs, const FLOAT_T<mpfr_t>& _rhs)
807  {
808  // TODO: mpfr_div results in infty when dividing by zero, although this should not be defined.
809  FLOAT_T<mpfr_t> res;
810  mpfr_div(res.mValue, FLOAT_T<mpfr_t>(_lhs).mValue, _rhs.mValue, MPFR_RNDN);
811  return res;
812  }
813 
814  template<typename Other, EnableIf< is_number_type<Other>, Not<is_interval_type<Other>> > = dummy>
815  friend FLOAT_T<mpfr_t> operator /(const FLOAT_T<mpfr_t>& _lhs, const Other& _rhs)
816  {
817  // TODO: mpfr_div results in infty when dividing by zero, although this should not be defined.
818  FLOAT_T<mpfr_t> res;
819  mpfr_div(res.mValue, _lhs.mValue, FLOAT_T<mpfr_t>(_rhs).mValue, MPFR_RNDN);
820  return res;
821  }
822 
823  friend FLOAT_T<mpfr_t>& operator++(FLOAT_T<mpfr_t>& _num)
824  {
825  mpfr_add_ui(_num.mValue, _num.mValue, 1, MPFR_RNDN);
826  return _num;
827  }
828 
829  friend FLOAT_T<mpfr_t>& operator--(FLOAT_T<mpfr_t>& _num)
830  {
831  mpfr_sub_ui(_num.mValue, _num.mValue, 1, MPFR_RNDN);
832  return _num;
833  }
834 
835  FLOAT_T<mpfr_t>& operator +=(const FLOAT_T<mpfr_t>& _rhs)
836  {
837  mpfr_add(mValue, mValue, _rhs.mValue, MPFR_RNDN);
838  return *this;
839  }
840 
841  FLOAT_T<mpfr_t>& operator +=(const mpfr_t& _rhs)
842  {
843  mpfr_add(mValue, mValue, _rhs, MPFR_RNDN);
844  return *this;
845  }
846 
847  FLOAT_T<mpfr_t>& operator -=(const FLOAT_T<mpfr_t>& _rhs)
848  {
849  mpfr_sub(mValue,mValue, _rhs.mValue, MPFR_RNDN);
850  return *this;
851  }
852 
853  FLOAT_T<mpfr_t>& operator -=(const mpfr_t& _rhs)
854  {
855  mpfr_sub(mValue,mValue, _rhs, MPFR_RNDN);
856  return *this;
857  }
858 
859  FLOAT_T<mpfr_t> operator-()
860  {
861  FLOAT_T<mpfr_t> res;
862  mpfr_mul_si(res.mValue, this->mValue, -1, MPFR_RNDN);
863  return res;
864  }
865 
866  FLOAT_T<mpfr_t>& operator *=(const FLOAT_T<mpfr_t>& _rhs)
867  {
868  mpfr_mul(mValue, mValue, _rhs.mValue, MPFR_RNDN);
869  return *this;
870  }
871 
872  FLOAT_T<mpfr_t>& operator *=(const mpfr_t& _rhs)
873  {
874  mpfr_mul(mValue, mValue, _rhs, MPFR_RNDN);
875  return *this;
876  }
877 
878  FLOAT_T<mpfr_t>& operator /=(const FLOAT_T<mpfr_t>& _rhs)
879  {
880  // TODO: mpfr_div results in infty when dividing by zero, although this should not be defined.
881  mpfr_div(mValue, mValue, _rhs.mValue, MPFR_RNDN);
882  return *this;
883  }
884 
885  FLOAT_T<mpfr_t>& operator /=(const mpfr_t& _rhs)
886  {
887  // TODO: mpfr_div results in infty when dividing by zero, although this should not be defined.
888  mpfr_div(mValue, mValue, _rhs, MPFR_RNDN);
889  return *this;
890  }
891 
892  /**
893  * Auxiliary Functions
894  */
895 
896  std::string toString() const
897  {
898  // TODO: Better rounding mode?
899 // std::string out;
900  char out[80];
901 // str << mpfr_get_d(mValue, MPFR_RNDN);
902  mpfr_sprintf(out, "%.50RDe", mValue);
903  return std::string(out);
904  }
905 
906  static void integerDistance(const FLOAT_T<mpfr_t>& a, const FLOAT_T<mpfr_t>& b, mpz_t& dist) {
907  distance(a.value(), b.value(), dist);
908  }
909 
910  inline static std::size_t bits2digits(precision_t b) {
911  const double LOG10_2 = 0.30102999566398119;
912  return std::size_t(std::floor( double(b) * LOG10_2 ));
913  }
914 
915  private:
916 
917  void clear() {
918  mpfr_clear(mValue);
919  }
920 
921  mpfr_prec_t convPrec(precision_t _prec) const
922  {
923  return _prec;
924  }
925 
926  /**
927  * @brief Converts a given float to its integer representation
928  * @details This function converts a mpfr float to a corresponding integer representation. This is done as follows:
929  * We can use a transformation similar to the one used for c++ doubles (with modifications). We use the mantissa and extend
930  * it by the exponent (we left-shift the exponent in front of the mantissa) and add the sign. This can be interpreted as an
931  * integer. Note that zero is a special case, which should be handled separately.
932  *
933  * @param intRep Parameter where we write the integer to.
934  * @param a input mpfr float.
935  */
936  static void to_int(mpz_t& intRep, const mpfr_t a) {
937 
938  //std::cout << "Bits per limb " << mp_bits_per_limb << std::endl;
939  //std::cout << "Number limbs " << std::ceil(double(a->_mpfr_prec)/double(mp_bits_per_limb)) << std::endl;
940  //std::cout << "Precision is " << a->_mpfr_prec << std::endl;
941  //std::cout << "Sign is " << a->_mpfr_sign << std::endl;
942  //std::cout << "Exponent is " << carl::binary(a->_mpfr_exp) << std::endl;
943  //std::cout << "Exponent is " << a->_mpfr_exp << std::endl;
944  //std::cout << "Min Exponent is " << mpfr_get_emin() << std::endl;
945  //std::cout << "Min Exponent is " << carl::binary(mpfr_get_emin()) << std::endl;
946  //std::cout << "Scaled exponent: " << (a->_mpfr_exp + std::abs(mpfr_get_emin())) << std::endl;
947  //std::cout << "Scaled exponent: " << carl::binary((a->_mpfr_exp + std::abs(mpfr_get_emin()))) << std::endl;
948 
949  // mpfr mantissa is stored in limbs (usually 64-bit words) - the number of those depends on the precision.
950  std::size_t limbs = (std::size_t)std::ceil(double(a->_mpfr_prec)/double(mp_bits_per_limb));
951  /*
952  std::cout << "Mantissa is ";
953  while( limbs > 0 ){
954  std::cout << carl::binary(h->_mpfr_d[limbs-1]) << " " << std::endl;
955  --limbs;
956  }
957  std::cout << std::endl;
958  limbs = std::ceil(double(h->_mpfr_prec)/double(mp_bits_per_limb));
959  */
960  mpz_t mant;
961  mpz_t tmp;
962  mpz_init(mant);
963  mpz_init(tmp);
964  mpz_set_ui(mant,0);
965  mpz_set_ui(tmp,0);
966  // as mpfr uses whole limbs (64-bit) we can cut away the additional zeroes (offset), if there are any
967  long offset = mp_bits_per_limb - (a->_mpfr_prec % mp_bits_per_limb);
968  //std::cout << "Offset is " << offset << " bits" << std::endl;
969  //std::cout << "Mantissa is ";
970  //char outStr[1024];
971 
972  // assemble the integer representation of the mantissa. The limbs are stored in reversed order, least significant first.
973  while( limbs > 0 ){
974  mpz_set_ui(tmp, a->_mpfr_d[limbs-1]);
975  //std::cout << "Shift: " << (mp_bits_per_limb*(limbs-1)) << " bits" << std::endl;
976  mpz_mul_2exp(tmp, tmp, ((std::size_t)mp_bits_per_limb*(limbs-1)));
977  mpz_add(mant, mant, tmp);
978  --limbs;
979  }
980  // cut away unnecessary zeroes at the end (divide by 2^offset -> left-shift).
981  mpz_cdiv_q_2exp(mant, mant, offset);
982 
983  //mpz_get_str(outStr, 2,mant);
984  //std::cout << "Mantissa: " << std::string(outStr) << std::endl;
985 
986  // set the exponent (8-bit), as it is in 2s complement, subtract the minimum to shift the exponent and get exponents ordered,
987  // right shift to put it before the mantissa later.
988  mpz_set_ui(tmp, (a->_mpfr_exp + std::abs(mpfr_get_emin())));
989  //std::cout << "Shift by " << (std::ceil(double(h->_mpfr_prec)/double(mp_bits_per_limb))*64+64-offset) << " bits" << std::endl;
990  mpz_mul_2exp(tmp,tmp,a->_mpfr_prec);
991 
992  // assemble the whole representation by addition and finally add sign.
993  mpz_add(mant,mant,tmp);
994  mpz_mul_si(mant,mant,a->_mpfr_sign);
995  mpz_set(intRep, mant);
996  /*
997  mpz_get_str(outStr, 2,mant);
998  std::cout << std::string(outStr) << std::endl;
999  mpz_get_str(outStr, 10,mant);
1000  std::cout << std::string(outStr) << std::endl;
1001  */
1002 
1003  // cleanup.
1004  mpz_clear(mant);
1005  mpz_clear(tmp);
1006  }
1007 
1008  /**
1009  * @brief Calculates the distance in ULPs between two mpfr floats a and b.
1010  * @details We use an integer representation for calculating the distance. Special cases are whenever we reach or cross
1011  * zero. TODO: Include cases for NaN and INFTY.
1012  *
1013  * @param a
1014  * @param b
1015  *
1016  * @return [description]
1017  */
1018  static void distance(const mpfr_t& a, const mpfr_t& b, mpz_t& dist) {
1019  // initialize variables
1020  mpz_t intRepA;
1021  mpz_t intRepB;
1022  mpz_init(intRepA);
1023  mpz_init(intRepB);
1024 
1025  // the offset is used to cope with the exponent differences
1026  long offset = a->_mpfr_exp - b->_mpfr_exp;
1027  //std::cout << "Offset " << offset << std::endl;
1028 
1029  // get integer representations, we use absolute values for simplicity.
1030  to_int(intRepA, a);
1031  to_int(intRepB, b);
1032  mpz_abs(intRepA, intRepA);
1033  mpz_abs(intRepB, intRepB);
1034 
1035  // case distinction to cope with zero.
1036  if(mpfr_zero_p(a) != 0) { // a is zero
1037  if(mpfr_zero_p(b) != 0){ // b is also zero
1038  mpz_set_ui(dist,0);
1039  }else{
1040 
1041  // b is not zero -> we compute the distance from close to zero to b and add 1.
1042  mpfr_t zero;
1043  mpz_t intRepZero;
1044  mpfr_init2(zero,mpfr_get_prec(a));
1045  mpz_init(intRepZero);
1046 
1047  mpfr_set_ui(zero,0, MPFR_RNDZ);
1048  if(b->_mpfr_sign > 0) {
1049  mpfr_nextabove(zero);
1050  }else{
1051  mpfr_nextbelow(zero);
1052  }
1053 
1054  to_int(intRepZero, zero);
1055  mpz_abs(intRepZero, intRepZero);
1056  mpz_sub(dist, intRepB,intRepZero);
1057  mpz_add_ui(dist,dist, 1);
1058 
1059  mpfr_clear(zero);
1060  mpz_clear(intRepZero);
1061  }
1062  } else if(mpfr_zero_p(b) != 0) { // a is not zero, b is zero
1063  mpfr_t zero;
1064  mpz_t intRepZero;
1065  mpfr_init2(zero,mpfr_get_prec(a));
1066  mpz_init(intRepZero);
1067 
1068  mpfr_set_ui(zero,0, MPFR_RNDZ);
1069  if(a->_mpfr_sign > 0) {
1070  mpfr_nextabove(zero);
1071  }else{
1072  mpfr_nextbelow(zero);
1073  }
1074 
1075  to_int(intRepZero, zero);
1076  mpz_abs(intRepZero, intRepZero);
1077  mpz_sub(dist, intRepA,intRepZero);
1078  mpz_add_ui(dist,dist, 1);
1079 
1080  mpfr_clear(zero);
1081  mpz_clear(intRepZero);
1082  } else if(a->_mpfr_sign == b->_mpfr_sign) { // both are not zero and at the same side
1083  mpz_sub(dist, intRepA, intRepB);
1084  mpz_abs(dist,dist);
1085  } else { // both are not zero and one is larger, the other one is less zero, compute both dists to zero and add 2.
1086  mpfr_t zeroA;
1087  mpfr_init2(zeroA,mpfr_get_prec(a));
1088  mpfr_t zeroB;
1089  mpfr_init2(zeroB,mpfr_get_prec(a));
1090 
1091  mpfr_set_ui(zeroA,0, MPFR_RNDZ);
1092  mpfr_set_ui(zeroB,0, MPFR_RNDZ);
1093 
1094  if(a->_mpfr_sign > 0) {
1095  mpfr_nextabove(zeroA);
1096  mpfr_nextbelow(zeroB);
1097  }else{
1098  mpfr_nextbelow(zeroA);
1099  mpfr_nextabove(zeroB);
1100  }
1101  mpz_t intRepZeroA;
1102  mpz_init(intRepZeroA);
1103  mpz_t intRepZeroB;
1104  mpz_init(intRepZeroB);
1105  mpz_t d2;
1106  mpz_init(d2);
1107 
1108  to_int(intRepZeroA, zeroA);
1109  mpz_abs(intRepZeroA, intRepZeroA);
1110  to_int(intRepZeroB, zeroB);
1111  mpz_abs(intRepZeroB, intRepZeroB);
1112 
1113  mpz_sub(dist, intRepA,intRepZeroA);
1114  mpz_sub(d2, intRepB,intRepZeroB);
1115  mpz_add(dist, dist, d2);
1116  mpz_add_ui(dist,dist, 2);
1117 
1118  mpfr_clear(zeroA);
1119  mpfr_clear(zeroB);
1120  mpz_clear(intRepZeroA);
1121  mpz_clear(intRepZeroB);
1122  mpz_clear(d2);
1123  }
1124  //std::cout << "Modify by " << 2*std::abs(offset)*a->_mpfr_prec << std::endl;
1125 
1126  // shift by offset (exponent differences).
1127  long shift = 2*std::abs(offset)*unsigned(a->_mpfr_prec);
1128  mpz_sub_ui(dist, dist, shift);
1129  // cleanup.
1130  mpz_clear(intRepA);
1131  mpz_clear(intRepB);
1132  }
1133 };
1134 
1135 template<>
1136 inline bool isInfinity(const FLOAT_T<mpfr_t>& _in) {
1137  return (mpfr_inf_p(_in.value()) != 0);
1138 }
1139 
1140 template<>
1141 inline bool isNan(const FLOAT_T<mpfr_t>& _in) {
1142  return (mpfr_nan_p(_in.value()) != 0);
1143 }
1144 
1145 template<>
1146 inline bool AlmostEqual2sComplement<FLOAT_T<mpfr_t>>(const FLOAT_T<mpfr_t>& A, const FLOAT_T<mpfr_t>& B, unsigned maxUlps)
1147 {
1148  //std::cout << "Distance: " << FLOAT_T<mpfr_t>::integerDistance(A,B) << std::endl;
1149  mpz_t distance;
1150  mpz_init(distance);
1151  FLOAT_T<mpfr_t>::integerDistance(A,B,distance);
1152  mpz_sub_ui(distance, distance, maxUlps);
1153  bool less = (mpz_sgn(distance) <= 0);
1154  mpz_clear(distance);
1155  return less;
1156 }
1157 
1158 }// namespace
1159 
1160 namespace std{
1161 
1162  template<>
1163  struct hash<carl::FLOAT_T<mpfr_t>> {
1164  size_t operator()(const carl::FLOAT_T<mpfr_t>& _in) const {
1165 
1166  __mpfr_struct numStruct = *_in.value();
1167  std::size_t limbs = (std::size_t)std::ceil(double(numStruct._mpfr_prec)/double(mp_bits_per_limb));
1168 
1169  size_t seed = 0;
1170  //std::cout << "start hash" << std::endl;
1171  //std::cout << "seed: " << seed << std::endl;
1172  if(mpfr_number_p(_in.value()) && mpfr_zero_p(_in.value()) == 0){
1173  while(limbs > 0) {
1174  carl::hash_add(seed,numStruct._mpfr_d[limbs-1]);
1175  //std::cout << "seed: " << seed << std::endl;
1176  --limbs;
1177  }
1178  }
1179  //std::cout << "seed: " << seed << std::endl;
1180  if(mpfr_nan_p(_in.value()) == 0){
1181  carl::hash_add(seed, size_t(numStruct._mpfr_sign));
1182  }
1183  //std::cout << "seed: " << seed << std::endl;
1184  carl::hash_add(seed, size_t(numStruct._mpfr_prec));
1185  //std::cout << "seed: " << seed << std::endl;
1186  carl::hash_add(seed, size_t(numStruct._mpfr_exp));
1187  //std::cout << "seed: " << seed << std::endl;
1188  //std::cout << "end hash" << std::endl;
1189  return seed;
1190  }
1191  };
1192 
1193  template<>
1194  class numeric_limits<carl::FLOAT_T<mpfr_t>>
1195  {
1196  public:
1197  static const bool is_specialized = true;
1198  static const bool is_signed = true;
1199  static const bool is_integer = false;
1200  static const bool is_exact = false;
1201  static const int radix = 2;
1202 
1203  static const bool has_infinity = true;
1204  static const bool has_quiet_NaN = true;
1205  static const bool has_signaling_NaN = true;
1206 
1207  static const bool is_iec559 = true; // = IEEE 754
1208  static const bool is_bounded = true;
1209  static const bool is_modulo = false;
1210  static const bool traps = true;
1211  static const bool tinyness_before = true;
1212 
1213  static const float_denorm_style has_denorm = denorm_absent;
1214 
1215 
1216  inline static carl::FLOAT_T<mpfr_t> (min) (carl::CARL_RND _rnd = carl::CARL_RND::N, carl::precision_t precision = carl::FLOAT_T<mpfr_t>::defaultPrecision()) { return carl::FLOAT_T<mpfr_t>::minval(_rnd, precision); }
1217  inline static carl::FLOAT_T<mpfr_t> (max) (carl::CARL_RND _rnd = carl::CARL_RND::N, carl::precision_t precision = carl::FLOAT_T<mpfr_t>::defaultPrecision()) { return carl::FLOAT_T<mpfr_t>::maxval(_rnd, precision); }
1218  inline static carl::FLOAT_T<mpfr_t> lowest (carl::CARL_RND _rnd = carl::CARL_RND::N, carl::precision_t precision = carl::FLOAT_T<mpfr_t>::defaultPrecision()) { return -carl::FLOAT_T<mpfr_t>::maxval(_rnd, precision); }
1219 
1220  // Returns smallest eps such that 1 + eps != 1 (classic machine epsilon)
1221  inline static carl::FLOAT_T<mpfr_t> epsilon(carl::precision_t precision = carl::FLOAT_T<mpfr_t>::defaultPrecision()) { return carl::FLOAT_T<mpfr_t>::machine_epsilon(precision); }
1222 
1223  // Returns smallest eps such that x + eps != x (relative machine epsilon)
1224  inline static carl::FLOAT_T<mpfr_t> epsilon(const carl::FLOAT_T<mpfr_t>& x) { return carl::FLOAT_T<mpfr_t>::machine_epsilon(x); }
1225 
1226  inline static carl::FLOAT_T<mpfr_t> round_error(carl::CARL_RND _rnd = carl::CARL_RND::N, carl::precision_t precision = carl::FLOAT_T<mpfr_t>::defaultPrecision()) { return carl::FLOAT_T<mpfr_t>(0.5, _rnd, precision); }
1227 
1228  inline static const carl::FLOAT_T<mpfr_t> infinity() { return carl::FLOAT_T<mpfr_t>::const_infinity(); }
1229 
1230  inline static const carl::FLOAT_T<mpfr_t> quiet_NaN() { return carl::FLOAT_T<mpfr_t>().setNan(); }
1231  inline static const carl::FLOAT_T<mpfr_t> signaling_NaN() { return carl::FLOAT_T<mpfr_t>().setNan(); }
1232  inline static const carl::FLOAT_T<mpfr_t> denorm_min() { return (min)(); }
1233 
1234  // Please note, exponent range is not fixed in MPFR
1235  static const int min_exponent = MPFR_EMIN_DEFAULT;
1236  static const int max_exponent = MPFR_EMAX_DEFAULT;
1237  static const int min_exponent10 = (int) (MPFR_EMIN_DEFAULT * 0.3010299956639811);
1238  static const int max_exponent10 = (int) (MPFR_EMAX_DEFAULT * 0.3010299956639811);
1239 
1240  // Following members should be constant according to standard, but they can be variable in MPFR
1241  // So we define them as functions here.
1242  //
1243  // This is preferable way for std::numeric_limits<carl::FLOAT_T<mpfr_t>> specialization.
1244  // But it is incompatible with standard std::numeric_limits and might not work with other libraries, e.g. boost.
1245  // See below for compatible implementation.
1246  inline static float_round_style round_style() { return round_to_nearest; }
1247 
1248  inline static std::size_t digits() { return std::size_t(carl::FLOAT_T<mpfr_t>::defaultPrecision()); }
1249  inline static std::size_t digits( const carl::FLOAT_T<mpfr_t>& x ) { return std::size_t(x.precision()); }
1250 
1251  inline static std::size_t digits10( carl::precision_t precision = carl::FLOAT_T<mpfr_t>::defaultPrecision() ) { return carl::FLOAT_T<mpfr_t>::bits2digits(precision); }
1252  inline static std::size_t digits10( const carl::FLOAT_T<mpfr_t>& x ) { return carl::FLOAT_T<mpfr_t>::bits2digits(x.precision()); }
1253  inline static std::size_t max_digits10( carl::precision_t precision = carl::FLOAT_T<mpfr_t>::defaultPrecision() ) { return digits10(precision); }
1254  };
1255 }
1256 
1257 #endif