carl  25.04
Computer ARithmetic Library
FloatExample.cpp
Go to the documentation of this file.
1 /**
2  * This file is for short testing - not intended for distribution.
3  * @file FloatExample.cpp
4  * @author Stefan Schupp <stefan.schupp@cs.rwth-aachen.de>
5  *
6  * @since 2013-12-16
7  * @version 2013-12-16
8  */
9 #include <carl-arith/interval/policies/rounding.h>
10 #include <carl-arith/numbers/numbers.h>
11 
12 #ifdef USE_MPFR_FLOAT
13 void to_int(mpz_t intRep, mpfr_t a) {
14  //std::cout << "Bits per limb " << mp_bits_per_limb << std::endl;
15  //std::cout << "Number limbs " << std::ceil(double(a->_mpfr_prec)/double(mp_bits_per_limb)) << std::endl;
16  //std::cout << "Precision is " << a->_mpfr_prec << std::endl;
17  //std::cout << "Sign is " << a->_mpfr_sign << std::endl;
18  //std::cout << "Exponent is " << carl::binary(a->_mpfr_exp) << std::endl;
19  //std::cout << "Exponent is " << a->_mpfr_exp << std::endl;
20  //std::cout << "Min Exponent is " << mpfr_get_emin() << std::endl;
21  //std::cout << "Min Exponent is " << carl::binary(mpfr_get_emin()) << std::endl;
22  //std::cout << "Scaled exponent: " << (a->_mpfr_exp + std::abs(mpfr_get_emin())) << std::endl;
23  //std::cout << "Scaled exponent: " << carl::binary((a->_mpfr_exp + std::abs(mpfr_get_emin()))) << std::endl;
24 
25  // mpfr mantissa is stored in limbs (usually 64-bit words) - the number of those depends on the precision.
26  int limbs = std::ceil(double(a->_mpfr_prec)/double(mp_bits_per_limb));
27  /*
28  std::cout << "Mantissa is ";
29  while( limbs > 0 ){
30  std::cout << carl::binary(h->_mpfr_d[limbs-1]) << " " << std::endl;
31  --limbs;
32  }
33  std::cout << std::endl;
34  limbs = std::ceil(double(h->_mpfr_prec)/double(mp_bits_per_limb));
35  */
36  mpz_t mant;
37  mpz_t tmp;
38  mpz_init(mant);
39  mpz_init(tmp);
40  mpz_set_ui(mant,0);
41  mpz_set_ui(tmp,0);
42  // as mpfr uses whole limbs (64-bit) we can cut away the additional zeroes (offset), if there are any
43  unsigned offset = mp_bits_per_limb - (a->_mpfr_prec % mp_bits_per_limb);
44  //std::cout << "Offset is " << offset << " bits" << std::endl;
45  //std::cout << "Mantissa is ";
46  //char outStr[1024];
47 
48  // assemble the integer representation of the mantissa. The limbs are stored in reversed order, least significant first.
49  while( limbs > 0 ){
50  mpz_set_ui(tmp, a->_mpfr_d[limbs-1]);
51  //std::cout << "Shift: " << (mp_bits_per_limb*(limbs-1)) << " bits" << std::endl;
52  mpz_mul_2exp(tmp, tmp, (mp_bits_per_limb*(limbs-1)));
53  mpz_add(mant, mant, tmp);
54  --limbs;
55  }
56  // cut away unnecessary zeroes at the end (divide by 2^offset -> left-shift).
57  mpz_cdiv_q_2exp(mant, mant, offset);
58 
59  //mpz_get_str(outStr, 2,mant);
60  //std::cout << "Mantissa: " << std::string(outStr) << std::endl;
61 
62  // set the exponent (8-bit), as it is in 2s complement, subtract the minimum to shift the exponent and get exponents ordered,
63  // right shift to put it before the mantissa later.
64  mpz_set_ui(tmp, (a->_mpfr_exp + std::abs(mpfr_get_emin())));
65  //std::cout << "Shift by " << (std::ceil(double(h->_mpfr_prec)/double(mp_bits_per_limb))*64+64-offset) << " bits" << std::endl;
66  mpz_mul_2exp(tmp,tmp,a->_mpfr_prec);
67 
68  // assemble the whole representation by addition and finally add sign.
69  mpz_add(mant,mant,tmp);
70  mpz_mul_si(mant,mant,a->_mpfr_sign);
71  mpz_set(intRep, mant);
72  /*
73  mpz_get_str(outStr, 2,mant);
74  std::cout << std::string(outStr) << std::endl;
75  mpz_get_str(outStr, 10,mant);
76  std::cout << std::string(outStr) << std::endl;
77  */
78 
79  // cleanup.
80  mpz_clear(mant);
81  mpz_clear(tmp);
82 }
83 
84 unsigned distance(mpfr_t a, mpfr_t b) {
85  // initialize variables
86  mpz_t intRepA;
87  mpz_t intRepB;
88  mpz_init(intRepA);
89  mpz_init(intRepB);
90  mpz_t distance;
91  mpz_init(distance);
92 
93  // the offset is used to cope with the exponent differences
94  long offset = a->_mpfr_exp - b->_mpfr_exp;
95  //std::cout << "Offset " << offset << std::endl;
96 
97  // get integer representations, we use absolute values for simplicity.
98  to_int(intRepA, a);
99  to_int(intRepB, b);
100  mpz_abs(intRepA, intRepA);
101  mpz_abs(intRepB, intRepB);
102 
103 
104 
105  // case distinction to cope with zero.
106  if(mpfr_zero_p(a) != 0) { // a is zero
107  if(mpfr_zero_p(b) != 0){ // b is also zero
108  mpz_clear(distance);
109  mpz_clear(intRepA);
110  mpz_clear(intRepB);
111  return 0;
112  }
113 
114  // b is not zero -> we compute the distance from close to zero to b and add 1.
115  mpfr_t zero;
116  mpz_t intRepZero;
117  mpfr_init2(zero,mpfr_get_prec(a));
118  mpz_init(intRepZero);
119 
120  mpfr_set_ui(zero,0, MPFR_RNDZ);
121  if(b->_mpfr_sign > 0) {
122  mpfr_nextabove(zero);
123  }else{
124  mpfr_nextbelow(zero);
125  }
126 
127  to_int(intRepZero, zero);
128  mpz_abs(intRepZero, intRepZero);
129  mpz_sub(distance, intRepB,intRepZero);
130  mpz_add_ui(distance,distance, 1);
131 
132  mpfr_clear(zero);
133  mpz_clear(intRepZero);
134  } else if(mpfr_zero_p(b) != 0) { // a is not zero, b is zero
135  mpfr_t zero;
136  mpz_t intRepZero;
137  mpfr_init2(zero,mpfr_get_prec(a));
138  mpz_init(intRepZero);
139 
140  mpfr_set_ui(zero,0, MPFR_RNDZ);
141  if(a->_mpfr_sign > 0) {
142  mpfr_nextabove(zero);
143  }else{
144  mpfr_nextbelow(zero);
145  }
146 
147 
148  to_int(intRepZero, zero);
149  mpz_abs(intRepZero, intRepZero);
150  mpz_sub(distance, intRepA,intRepZero);
151  mpz_add_ui(distance,distance, 1);
152 
153  mpfr_clear(zero);
154  mpz_clear(intRepZero);
155  } else if(a->_mpfr_sign == b->_mpfr_sign) { // both are not zero and at the same side
156  mpz_sub(distance, intRepA, intRepB);
157  mpz_abs(distance,distance);
158  } else { // both are not zero and one is larger, the other one is less zero, compute both distances to zero and add 2.
159  mpfr_t zeroA;
160  mpfr_init2(zeroA,mpfr_get_prec(a));
161  mpfr_t zeroB;
162  mpfr_init2(zeroB,mpfr_get_prec(a));
163 
164  mpfr_set_ui(zeroA,0, MPFR_RNDZ);
165  mpfr_set_ui(zeroB,0, MPFR_RNDZ);
166 
167  if(a->_mpfr_sign > 0) {
168  mpfr_nextabove(zeroA);
169  mpfr_nextbelow(zeroB);
170  }else{
171  mpfr_nextbelow(zeroA);
172  mpfr_nextabove(zeroB);
173  }
174  mpz_t intRepZeroA;
175  mpz_init(intRepZeroA);
176  mpz_t intRepZeroB;
177  mpz_init(intRepZeroB);
178  mpz_t d2;
179  mpz_init(d2);
180 
181  to_int(intRepZeroA, zeroA);
182  mpz_abs(intRepZeroA, intRepZeroA);
183  to_int(intRepZeroB, zeroB);
184  mpz_abs(intRepZeroB, intRepZeroB);
185 
186  mpz_sub(distance, intRepA,intRepZeroA);
187  mpz_sub(d2, intRepB,intRepZeroB);
188  mpz_add(distance, distance, d2);
189  mpz_add_ui(distance,distance, 2);
190 
191  mpfr_clear(zeroA);
192  mpfr_clear(zeroB);
193  mpz_clear(intRepZeroA);
194  mpz_clear(intRepZeroB);
195  mpz_clear(d2);
196  }
197  //std::cout << "Modify by " << 2*std::abs(offset)*a->_mpfr_prec << std::endl;
198 
199  // shift by offset (exponent differences).
200  unsigned result = mpz_get_ui(distance) - 2*std::abs(offset)*a->_mpfr_prec;
201 
202  // cleanup.
203  mpz_clear(distance);
204  mpz_clear(intRepA);
205  mpz_clear(intRepB);
206  return result;
207 }
208 #endif
209 
210 int main (int argc, char** argv)
211 {
212 #ifdef USE_MPFR_FLOAT
213  //carl::FLOAT_T<mpfr_t> eps = std::numeric_limits<carl::FLOAT_T<mpfr_t>>::epsilon();
214  carl::FLOAT_T<mpfr_t> eps = carl::FLOAT_T<mpfr_t>::maxval();
215  char out[1000];
216  std::string fString = "%." + std::to_string(eps.precision()+5) + "RNf";
217  //char format[fString.length()];
218  char* format = (char*)fString.c_str();
219  std::cout << "Precision is " << eps.precision() << std::endl;
220  mpfr_sprintf(out, "%.RDe", eps.value());
221  std::cout << std::string(out) << std::endl;
222 
223 
224  /*
225 #ifdef USE_MPFR_FLOAT
226  /*
227  mpfr_ptr numberA = new mpfr_t();
228  mpfr_init(numberA);
229  mpfr_set_ui(numberA, 5, MPFR_RNDN);
230 
231  mpfr_ptr numberB = new mpfr_t;
232  mpfr_init(numberB);
233  mpfr_set_ui(numberB, 5, MPFR_RNDN);
234 
235  mpfr_ptr result;
236  mpfr_init(result);
237 
238  //carl::convRnd<mpfr_ptr> r;
239  //mpfr_rnd_t tmp = r.operator() (carl::CARL_RND::CARL_RNDD);
240 
241  //result = carl::add_down(numberA, numberB);
242 
243  mpfr_clears(numberA, numberB, result);
244  */
245  /*
246  mpfr_t f;
247  mpfr_init2(f,100);
248  mpfr_set_d(f, 0, MPFR_RNDN);
249  //mpfr_nextabove(f);
250  //mpfr_nextbelow(f);
251  //mpfr_nextbelow(f);
252  //mpfr_nextbelow(f);
253  //mpfr_nextbelow(f);
254 
255 
256  char out1[120];
257  char out2[120];
258  for(unsigned target = 0; target < 400; ++target){
259  int targetDistance = target;
260  mpfr_t g;
261  mpfr_init2(g,100);
262  mpfr_set_d(g, 0, MPFR_RNDN);
263  while(targetDistance > 0){
264  mpfr_nextbelow(g);
265  --targetDistance;
266  }
267 
268  mpfr_sprintf(out1, "%.64RDe", f);
269  mpfr_sprintf(out2, "%.64RDe", g);
270  std::cout << "Distance (" << std::string(out1) << ", " << std::string(out2) << ") =" << distance(f,g) << std::endl;
271  mpfr_clear(g);
272  std::cout << "##############################" << std::endl;
273  }
274 
275  mpfr_clear(f);
276  */
277  /*
278  carl::FLOAT_T<mpfr_t> a(1);
279  carl::FLOAT_T<mpfr_t> b(2);
280 
281  carl::AlmostEqual2sComplement(a,b);
282  */
283  /*
284  std::cout << "F" << std::endl;
285  char out[120];
286  mpfr_sprintf(out, "%.100RDe", f);
287  std::cout << std::string(out) << std::endl;
288  mpfr_sprintf(out, "%.100RDb", f);
289  std::cout << std::string(out) << std::endl;
290 
291  std::cout << "G" << std::endl;
292  mpfr_sprintf(out, "%.100RDe", g);
293  std::cout << std::string(out) << std::endl;
294  mpfr_sprintf(out, "%.100RDb", g);
295  std::cout << std::string(out) << std::endl;
296 
297  mpz_t significantA;
298  mpz_init(significantA);
299  mpz_t significantB;
300  mpz_init(significantB);
301 
302  mpfr_exp_t expA = mpfr_get_z_exp(significantA,f);
303  std::cout << "ExpA: " << expA << std::endl;
304  mpfr_exp_t expB = mpfr_get_z_exp(significantB,g);
305  std::cout << "ExpB: " << expB << std::endl;
306  //std::cout << "Exp converted: " << carl::binary(expA) << std::endl;
307  std::cout << __func__ << ": 2^" << expA << " * " << mpz_get_si(significantA) << " and 2^" << expB << " * " << mpz_get_si(significantB) << std::endl;
308 
309  if(expA == expB) {
310  std::cout << "Equal exponents, compare mantissa" << std::endl;
311  mpz_t diff;
312  mpz_init(diff);
313  mpz_sub(diff, significantA, significantB);
314  mpz_abs(diff,diff);
315  std::cout << "difference: " << mpz_get_si(diff) << std::endl;
316  }else if (expA > expB) {
317  std::cout << "expA > expB" << std::endl;
318  mpz_mul_2exp(significantB,significantB,expA-expB);
319  std::cout << "Shifted: " << mpz_get_si(significantB) << std::endl;
320  mpz_t diff;
321  mpz_init(diff);
322  mpz_sub(diff, significantA, significantB);
323  mpz_abs(diff,diff);
324  std::cout << "difference: " << mpz_get_si(diff) << std::endl;
325  }else{
326  std::cout << "expB > expA" << std::endl;
327  mpz_mul_2exp(significantA,significantA,expB-expA);
328  std::cout << "Shifted: " << mpz_get_si(significantA) << std::endl;
329  mpz_t diff;
330  mpz_init(diff);
331  mpz_sub(diff, significantA, significantB);
332  mpz_abs(diff,diff);
333  std::cout << "difference: " << mpz_get_si(diff) << std::endl;
334  }
335  */
336  /*
337  std::cout << "DECONSTRUCT mpfr_t" << std::endl;
338  mpfr_t h;
339  mpfr_init2(h,65);
340  mpfr_set_d(h, 0, MPFR_RNDN);
341 
342  targetDistance = 2;
343  while(targetDistance > 0){
344  mpfr_nextbelow(h);
345  --targetDistance;
346  }
347 
348  std::cout << "Bits per limb " << mp_bits_per_limb << std::endl;
349  std::cout << "Number limbs " << std::ceil(double(h->_mpfr_prec)/double(mp_bits_per_limb)) << std::endl;
350  std::cout << "Precision is " << h->_mpfr_prec << std::endl;
351  std::cout << "Sign is " << h->_mpfr_sign << std::endl;
352  std::cout << "Exponent is " << carl::binary(h->_mpfr_exp) << std::endl;
353  std::cout << "Exponent is " << h->_mpfr_exp << std::endl;
354  std::cout << "Min Exponent is " << mpfr_get_emin() << std::endl;
355  std::cout << "Min Exponent is " << carl::binary(mpfr_get_emin()) << std::endl;
356 
357  std::cout << "Scaled exponent: " << (h->_mpfr_exp + std::abs(mpfr_get_emin())) << std::endl;
358  std::cout << "Scaled exponent: " << carl::binary((h->_mpfr_exp + std::abs(mpfr_get_emin()))) << std::endl;
359 
360 
361  int limbs = std::ceil(double(h->_mpfr_prec)/double(mp_bits_per_limb));
362  std::cout << "Mantissa is ";
363  while( limbs > 0 ){
364  std::cout << carl::binary(h->_mpfr_d[limbs-1]) << " " << std::endl;
365  --limbs;
366  }
367  std::cout << std::endl;
368  limbs = std::ceil(double(h->_mpfr_prec)/double(mp_bits_per_limb));
369 
370  mpz_t mant;
371  mpz_t tmp;
372  mpz_init(mant);
373  mpz_init(tmp);
374  unsigned offset = mp_bits_per_limb - (h->_mpfr_prec % mp_bits_per_limb);
375  std::cout << "Offset is " << offset << " bits" << std::endl;
376  std::cout << "Mantissa is ";
377  char outStr[1024];
378  while( limbs > 0 ){
379  mpz_set_ui(tmp, h->_mpfr_d[limbs-1]);
380  //std::cout << "Shift: " << (mp_bits_per_limb*(limbs-1)) << " bits" << std::endl;
381  mpz_mul_2exp(tmp, tmp, (mp_bits_per_limb*(limbs-1)));
382  mpz_add(mant, mant, tmp);
383  //mant = (h->_mpfr_d[limbs-1] << (mp_bits_per_limb*(limbs)));
384  --limbs;
385  }
386  mpz_cdiv_q_2exp(mant, mant, offset);
387  mpz_get_str(outStr, 2,mant);
388  std::cout << std::string(outStr) << std::endl;
389 
390  mpz_set_ui(tmp, (h->_mpfr_exp + std::abs(mpfr_get_emin())));
391  std::cout << "Shift by " << (std::ceil(double(h->_mpfr_prec)/double(mp_bits_per_limb))*64+64-offset) << " bits" << std::endl;
392  mpz_mul_2exp(tmp,tmp,h->_mpfr_prec);
393  mpz_add(mant,mant,tmp);
394 
395  mpz_get_str(outStr, 2,mant);
396  std::cout << std::string(outStr) << std::endl;
397  mpz_get_str(outStr, 10,mant);
398  std::cout << std::string(outStr) << std::endl;
399  */
400 
401  //mpfr_sprintf(out, "%.100RDb", h);
402  //std::cout << std::string(out) << std::endl;
403  //std::cout << carl::binary(h->_mpfr_exp) << std::endl;
404  /*
405  mpfr_t res;
406  mpfr_init2(res,100);
407 
408  mpfr_t mantissa;
409  mpfr_init2(mantissa,100);
410  mpfr_set_z(mantissa,significant, MPFR_RNDN);
411 
412  mpfr_t exponent;
413  mpfr_init2(exponent,100);
414  mpfr_set_si(exponent, expA, MPFR_RNDN);
415 
416  mpfr_mul_2si(res, mantissa, expA, MPFR_RNDN);
417  mpfr_sprintf(out, "%.100RDb", res);
418  std::cout << "RES: " << std::endl << std::string(out) << std::endl;
419 
420  mpfr_clear(f);
421 
422  #endif
423  */
424 #endif
425  return 0;
426 }
main
int main(int argc, char **argv)
Definition: FloatExample.cpp:210
carl::ceil
Interval< Number > ceil(const Interval< Number > &_in)
Method which returns the next larger integer of the passed number or the number itself,...
Definition: Interval.h:1535
carl::abs
Interval< Number > abs(const Interval< Number > &_in)
Method which returns the absolute value of the passed number.
Definition: Interval.h:1511
carl::to_int
Integer to_int(const Interval< Number > &_floatInterval)
Casts the Interval to an arbitrary integer type which has a constructor for a native int.
Definition: Interval.h:1500
carl::statistics::timing::zero
auto zero()
Return a zero duration.
Definition: Timing.h:25
carl::FLOAT_T
Templated wrapper class which allows universal usage of different IEEE 754 implementations.
Definition: FLOAT_T.h:114
carl::FLOAT_T::precision
precision_t precision() const
If precision is used, this getter returns the acutal precision (default: 53 bit).
Definition: FLOAT_T.h:222
carl::FLOAT_T::value
const FloatType & value() const
Getter for the raw value contained.
Definition: FLOAT_T.h:213