cprover
float_utils.cpp
Go to the documentation of this file.
1 /*******************************************************************\
2 
3 Module:
4 
5 Author: Daniel Kroening, kroening@kroening.com
6 
7 \*******************************************************************/
8 
9 #include "float_utils.h"
10 
11 #include <cassert>
12 #include <algorithm>
13 
14 #include <util/arith_tools.h>
15 
17 {
18  bvt round_to_even=
20  bvt round_to_plus_inf=
22  bvt round_to_minus_inf=
24  bvt round_to_zero=
26 
27  rounding_mode_bits.round_to_even=bv_utils.equal(src, round_to_even);
28  rounding_mode_bits.round_to_plus_inf=bv_utils.equal(src, round_to_plus_inf);
29  rounding_mode_bits.round_to_minus_inf=bv_utils.equal(src, round_to_minus_inf);
30  rounding_mode_bits.round_to_zero=bv_utils.equal(src, round_to_zero);
31 }
32 
34 {
35  unbiased_floatt result;
36 
37  // we need to convert negative integers
38  result.sign=sign_bit(src);
39 
40  result.fraction=bv_utils.absolute_value(src);
41 
42  // build an exponent (unbiased) -- this is signed!
43  result.exponent=
45  src.size()-1,
46  address_bits(src.size()-1).to_long()+1);
47 
48  return rounder(result);
49 }
50 
52 {
53  unbiased_floatt result;
54 
55  result.fraction=src;
56 
57  // build an exponent (unbiased) -- this is signed!
58  result.exponent=
60  src.size()-1,
61  address_bits(src.size()-1).to_long()+1);
62 
63  result.sign=const_literal(false);
64 
65  return rounder(result);
66 }
67 
69  const bvt &src,
70  std::size_t dest_width)
71 {
72  return to_integer(src, dest_width, true);
73 }
74 
76  const bvt &src,
77  std::size_t dest_width)
78 {
79  return to_integer(src, dest_width, false);
80 }
81 
83  const bvt &src,
84  std::size_t dest_width,
85  bool is_signed)
86 {
87  assert(src.size()==spec.width());
88 
89  const unbiased_floatt unpacked=unpack(src);
90 
91  // The following is the usual case in ANSI-C, and we optimize for that.
93  {
94  bvt fraction=unpacked.fraction;
95 
96  if(dest_width>fraction.size())
97  {
98  bvt lsb_extension=bv_utils.build_constant(0U, dest_width-fraction.size());
99  fraction.insert(fraction.begin(),
100  lsb_extension.begin(),
101  lsb_extension.end());
102  }
103 
104  // if the exponent is positive, shift right
105  bvt offset=bv_utils.build_constant(fraction.size()-1,
106  unpacked.exponent.size());
107  bvt distance=bv_utils.sub(offset, unpacked.exponent);
108  bvt shift_result=bv_utils.shift(
109  fraction, bv_utilst::shiftt::LRIGHT, distance);
110 
111  // if the exponent is negative, we have zero anyways
112  bvt result=shift_result;
113  literalt exponent_sign=unpacked.exponent[unpacked.exponent.size()-1];
114 
115  for(std::size_t i=0; i<result.size(); i++)
116  result[i]=prop.land(result[i], !exponent_sign);
117 
118  // chop out the right number of bits from the result
119  if(result.size()>dest_width)
120  {
121  result.resize(dest_width);
122  }
123 
124  assert(result.size()==dest_width);
125 
126  // if signed, apply sign.
127  if(is_signed)
128  result=bv_utils.cond_negate(result, unpacked.sign);
129  else
130  {
131  // It's unclear what the behaviour for negative floats
132  // to integer shall be.
133  }
134 
135  return result;
136  }
137  else
138  throw "unsupported rounding mode";
139 }
140 
142 {
143  unbiased_floatt result;
144 
145  result.sign=const_literal(src.get_sign());
146  result.NaN=const_literal(src.is_NaN());
147  result.infinity=const_literal(src.is_infinity());
150 
151  return pack(bias(result));
152 }
153 
155  const bvt &src,
156  const ieee_float_spect &dest_spec)
157 {
158  assert(src.size()==spec.width());
159 
160  #if 1
161  // Catch the special case in which we extend,
162  // e.g. single to double.
163  // In this case, rounding can be avoided,
164  // but a denormal number may be come normal.
165  // Be careful to exclude the difficult case
166  // when denormalised numbers in the old format
167  // can be converted to denormalised numbers in the
168  // new format. Note that this is rare and will only
169  // happen with very non-standard formats.
170 
171  int sourceSmallestNormalExponent=-((1 << (spec.e - 1)) - 1);
172  int sourceSmallestDenormalExponent =
173  sourceSmallestNormalExponent - spec.f;
174 
175  // Using the fact that f doesn't include the hidden bit
176 
177  int destSmallestNormalExponent=-((1 << (dest_spec.e - 1)) - 1);
178 
179  if(dest_spec.e>=spec.e &&
180  dest_spec.f>=spec.f &&
181  !(sourceSmallestDenormalExponent < destSmallestNormalExponent))
182  {
183  unbiased_floatt unpacked_src=unpack(src);
184  unbiased_floatt result;
185 
186  // the fraction gets zero-padded
187  std::size_t padding=dest_spec.f-spec.f;
188  result.fraction=
189  bv_utils.concatenate(bv_utils.zeros(padding), unpacked_src.fraction);
190 
191  // the exponent gets sign-extended
192  result.exponent=
193  bv_utils.sign_extension(unpacked_src.exponent, dest_spec.e);
194 
195  // if the number was denormal and is normal in the new format,
196  // normalise it!
197  if(dest_spec.e > spec.e)
198  {
199  normalization_shift(result.fraction, result.exponent);
200  }
201 
202  // the flags get copied
203  result.sign=unpacked_src.sign;
204  result.NaN=unpacked_src.NaN;
205  result.infinity=unpacked_src.infinity;
206 
207  // no rounding needed!
208  spec=dest_spec;
209  return pack(bias(result));
210  }
211  else // NOLINT(readability/braces)
212  #endif
213  {
214  // we actually need to round
215  unbiased_floatt result=unpack(src);
216  spec=dest_spec;
217  return rounder(result);
218  }
219 }
220 
222 {
223  return prop.land(
224  !exponent_all_zeros(src),
225  !exponent_all_ones(src));
226 }
227 
230  const unbiased_floatt &src1,
231  const unbiased_floatt &src2)
232 {
233  // extend both
234  bvt extended_exponent1=
235  bv_utils.sign_extension(src1.exponent, src1.exponent.size()+1);
236  bvt extended_exponent2=
237  bv_utils.sign_extension(src2.exponent, src2.exponent.size()+1);
238 
239  assert(extended_exponent1.size()==extended_exponent2.size());
240 
241  // compute shift distance (here is the subtraction)
242  return bv_utils.sub(extended_exponent1, extended_exponent2);
243 }
244 
246  const bvt &src1,
247  const bvt &src2,
248  bool subtract)
249 {
250  unbiased_floatt unpacked1=unpack(src1);
251  unbiased_floatt unpacked2=unpack(src2);
252 
253  // subtract?
254  if(subtract)
255  unpacked2.sign=!unpacked2.sign;
256 
257  // figure out which operand has the bigger exponent
258  const bvt exponent_difference=subtract_exponents(unpacked1, unpacked2);
259  literalt src2_bigger=exponent_difference.back();
260 
261  const bvt bigger_exponent=
262  bv_utils.select(src2_bigger, unpacked2.exponent, unpacked1.exponent);
263 
264  // swap fractions as needed
265  const bvt new_fraction1=
266  bv_utils.select(src2_bigger, unpacked2.fraction, unpacked1.fraction);
267 
268  const bvt new_fraction2=
269  bv_utils.select(src2_bigger, unpacked1.fraction, unpacked2.fraction);
270 
271  // compute distance
272  const bvt distance=bv_utils.absolute_value(exponent_difference);
273 
274  // limit the distance: shifting more than f+3 bits is unnecessary
275  const bvt limited_dist=limit_distance(distance, spec.f+3);
276 
277  // pad fractions with 2 zeros from below
278  const bvt fraction1_padded=
279  bv_utils.concatenate(bv_utils.zeros(3), new_fraction1);
280  const bvt fraction2_padded=
281  bv_utils.concatenate(bv_utils.zeros(3), new_fraction2);
282 
283  // shift new_fraction2
284  literalt sticky_bit;
285  const bvt fraction1_shifted=fraction1_padded;
286  const bvt fraction2_shifted=sticky_right_shift(
287  fraction2_padded, limited_dist, sticky_bit);
288 
289  // sticky bit: or of the bits lost by the right-shift
290  bvt fraction2_stickied=fraction2_shifted;
291  fraction2_stickied[0]=prop.lor(fraction2_shifted[0], sticky_bit);
292 
293  // need to have two extra fraction bits for addition and rounding
294  const bvt fraction1_ext=
295  bv_utils.zero_extension(fraction1_shifted, fraction1_shifted.size()+2);
296  const bvt fraction2_ext=
297  bv_utils.zero_extension(fraction2_stickied, fraction2_stickied.size()+2);
298 
299  unbiased_floatt result;
300 
301  // now add/sub them
302  literalt subtract_lit=prop.lxor(unpacked1.sign, unpacked2.sign);
303  result.fraction=
304  bv_utils.add_sub(fraction1_ext, fraction2_ext, subtract_lit);
305 
306  // sign of result
307  literalt fraction_sign=result.fraction.back();
308  result.fraction=bv_utils.absolute_value(result.fraction);
309 
310  result.exponent=bigger_exponent;
311 
312  // adjust the exponent for the fact that we added two bits to the fraction
313  result.exponent=
314  bv_utils.add(
315  bv_utils.sign_extension(result.exponent, result.exponent.size()+1),
316  bv_utils.build_constant(2, result.exponent.size()+1));
317 
318  // NaN?
319  result.NaN=prop.lor(
320  prop.land(prop.land(unpacked1.infinity, unpacked2.infinity),
321  prop.lxor(unpacked1.sign, unpacked2.sign)),
322  prop.lor(unpacked1.NaN, unpacked2.NaN));
323 
324  // infinity?
325  result.infinity=prop.land(
326  !result.NaN,
327  prop.lor(unpacked1.infinity, unpacked2.infinity));
328 
329  // zero?
330  // Note that:
331  // 1. The zero flag isn't used apart from in divide and
332  // is only set on unpack
333  // 2. Subnormals mean that addition or subtraction can't round to 0,
334  // thus we can perform this test now
335  // 3. The rules for sign are different for zero
336  result.zero=prop.land(
337  !prop.lor(result.infinity, result.NaN),
338  !prop.lor(result.fraction));
339 
340 
341  // sign
342  literalt add_sub_sign=
343  prop.lxor(prop.lselect(src2_bigger, unpacked2.sign, unpacked1.sign),
344  fraction_sign);
345 
346  literalt infinity_sign=
347  prop.lselect(unpacked1.infinity, unpacked1.sign, unpacked2.sign);
348 
349  #if 1
350  literalt zero_sign=
352  prop.lor(unpacked1.sign, unpacked2.sign),
353  prop.land(unpacked1.sign, unpacked2.sign));
354 
355  result.sign=prop.lselect(
356  result.infinity,
357  infinity_sign,
358  prop.lselect(result.zero,
359  zero_sign,
360  add_sub_sign));
361  #else
362  result.sign=prop.lselect(
363  result.infinity,
364  infinity_sign,
365  add_sub_sign);
366  #endif
367 
368  #if 0
369  result.sign=const_literal(false);
370  result.fraction.resize(spec.f+1, const_literal(true));
371  result.exponent.resize(spec.e, const_literal(false));
372  result.NaN=const_literal(false);
373  result.infinity=const_literal(false);
374  // for(std::size_t i=0; i<result.fraction.size(); i++)
375  // result.fraction[i]=const_literal(true);
376 
377  for(std::size_t i=0; i<result.fraction.size(); i++)
378  result.fraction[i]=new_fraction2[i];
379 
380  return pack(bias(result));
381  #endif
382 
383  return rounder(result);
384 }
385 
388  const bvt &dist,
389  mp_integer limit)
390 {
391  std::size_t nb_bits=integer2unsigned(address_bits(limit));
392 
393  bvt upper_bits=dist;
394  upper_bits.erase(upper_bits.begin(), upper_bits.begin()+nb_bits);
395  literalt or_upper_bits=prop.lor(upper_bits);
396 
397  bvt lower_bits=dist;
398  lower_bits.resize(nb_bits);
399 
400  bvt result;
401  result.resize(lower_bits.size());
402 
403  // bitwise or with or_upper_bits
404  for(std::size_t i=0; i<result.size(); i++)
405  result[i]=prop.lor(lower_bits[i], or_upper_bits);
406 
407  return result;
408 }
409 
410 bvt float_utilst::mul(const bvt &src1, const bvt &src2)
411 {
412  // unpack
413  const unbiased_floatt unpacked1=unpack(src1);
414  const unbiased_floatt unpacked2=unpack(src2);
415 
416  // zero-extend the fractions
417  const bvt fraction1=
418  bv_utils.zero_extension(unpacked1.fraction, unpacked1.fraction.size()*2);
419  const bvt fraction2=
420  bv_utils.zero_extension(unpacked2.fraction, unpacked2.fraction.size()*2);
421 
422  // multiply fractions
423  unbiased_floatt result;
424  result.fraction=bv_utils.unsigned_multiplier(fraction1, fraction2);
425 
426  // extend exponents to account for overflow
427  // add two bits, as we do extra arithmetic on it later
428  const bvt exponent1=
429  bv_utils.sign_extension(unpacked1.exponent, unpacked1.exponent.size()+2);
430  const bvt exponent2=
431  bv_utils.sign_extension(unpacked2.exponent, unpacked2.exponent.size()+2);
432 
433  bvt added_exponent=bv_utils.add(exponent1, exponent2);
434 
435  // adjust, we are thowing in an extra fraction bit
436  // it has been extended above
437  result.exponent=bv_utils.inc(added_exponent);
438 
439  // new sign
440  result.sign=prop.lxor(unpacked1.sign, unpacked2.sign);
441 
442  // infinity?
443  result.infinity=prop.lor(unpacked1.infinity, unpacked2.infinity);
444 
445  // NaN?
446  {
447  bvt NaN_cond;
448 
449  NaN_cond.push_back(is_NaN(src1));
450  NaN_cond.push_back(is_NaN(src2));
451 
452  // infinity * 0 is NaN!
453  NaN_cond.push_back(prop.land(unpacked1.zero, unpacked2.infinity));
454  NaN_cond.push_back(prop.land(unpacked2.zero, unpacked1.infinity));
455 
456  result.NaN=prop.lor(NaN_cond);
457  }
458 
459  return rounder(result);
460 }
461 
462 bvt float_utilst::div(const bvt &src1, const bvt &src2)
463 {
464  // unpack
465  const unbiased_floatt unpacked1=unpack(src1);
466  const unbiased_floatt unpacked2=unpack(src2);
467 
468  std::size_t div_width=unpacked1.fraction.size()*2+1;
469 
470  // pad fraction1 with zeros
471  bvt fraction1=unpacked1.fraction;
472  fraction1.reserve(div_width);
473  while(fraction1.size()<div_width)
474  fraction1.insert(fraction1.begin(), const_literal(false));
475 
476  // zero-extend fraction2
477  const bvt fraction2=
478  bv_utils.zero_extension(unpacked2.fraction, div_width);
479 
480  // divide fractions
481  unbiased_floatt result;
482  bvt rem;
483  bv_utils.unsigned_divider(fraction1, fraction2, result.fraction, rem);
484 
485  // is there a remainder?
486  literalt have_remainder=bv_utils.is_not_zero(rem);
487 
488  // we throw this into the result, as one additional bit,
489  // to get the right rounding decision
490  result.fraction.insert(
491  result.fraction.begin(), have_remainder);
492 
493  // We will subtract the exponents;
494  // to account for overflow, we add a bit.
495  // we add a second bit for the adjust by extra fraction bits
496  const bvt exponent1=
497  bv_utils.sign_extension(unpacked1.exponent, unpacked1.exponent.size()+2);
498  const bvt exponent2=
499  bv_utils.sign_extension(unpacked2.exponent, unpacked2.exponent.size()+2);
500 
501  // subtract exponents
502  bvt added_exponent=bv_utils.sub(exponent1, exponent2);
503 
504  // adjust, as we have thown in extra fraction bits
505  result.exponent=bv_utils.add(
506  added_exponent,
507  bv_utils.build_constant(spec.f, added_exponent.size()));
508 
509  // new sign
510  result.sign=prop.lxor(unpacked1.sign, unpacked2.sign);
511 
512  // Infinity? This happens when
513  // 1) dividing a non-nan/non-zero by zero, or
514  // 2) first operand is inf and second is non-nan and non-zero
515  // In particular, inf/0=inf.
516  result.infinity=
517  prop.lor(
518  prop.land(!unpacked1.zero,
519  prop.land(!unpacked1.NaN,
520  unpacked2.zero)),
521  prop.land(unpacked1.infinity,
522  prop.land(!unpacked2.NaN,
523  !unpacked2.zero)));
524 
525  // NaN?
526  result.NaN=prop.lor(unpacked1.NaN,
527  prop.lor(unpacked2.NaN,
528  prop.lor(prop.land(unpacked1.zero, unpacked2.zero),
529  prop.land(unpacked1.infinity, unpacked2.infinity))));
530 
531  // Division by infinity produces zero, unless we have NaN
532  literalt force_zero=
533  prop.land(!unpacked1.NaN, unpacked2.infinity);
534 
535  result.fraction=bv_utils.select(force_zero,
536  bv_utils.zeros(result.fraction.size()), result.fraction);
537 
538  return rounder(result);
539 }
540 
541 bvt float_utilst::rem(const bvt &src1, const bvt &src2)
542 {
543  /* The semantics of floating-point remainder implemented as below
544  is the sensible one. Unfortunately this is not the one required
545  by IEEE-754 or fmod / remainder. Martin has discussed the
546  'correct' semantics with Christoph and Alberto at length as
547  well as talking to various hardware designers and we still
548  hasn't found a good way to implement them in a solver.
549  We have some approaches that are correct but they really
550  don't scale. */
551 
552  // stub: do src1-(src1/src2)*src2
553  return sub(src1, mul(div(src1, src2), src2));
554 }
555 
557 {
558  bvt result=src;
559  assert(!src.empty());
560  literalt &sign_bit=result[result.size()-1];
562  return result;
563 }
564 
566 {
567  bvt result=src;
568  assert(!src.empty());
569  result[result.size()-1]=const_literal(false);
570  return result;
571 }
572 
574  const bvt &src1,
575  relt rel,
576  const bvt &src2)
577 {
578  if(rel==relt::GT)
579  return relation(src2, relt::LT, src1); // swapped
580  else if(rel==relt::GE)
581  return relation(src2, relt::LE, src1); // swapped
582 
583  assert(rel==relt::EQ || rel==relt::LT || rel==relt::LE);
584 
585  // special cases: -0 and 0 are equal
586  literalt is_zero1=is_zero(src1);
587  literalt is_zero2=is_zero(src2);
588  literalt both_zero=prop.land(is_zero1, is_zero2);
589 
590  // NaN compares to nothing
591  literalt is_NaN1=is_NaN(src1);
592  literalt is_NaN2=is_NaN(src2);
593  literalt NaN=prop.lor(is_NaN1, is_NaN2);
594 
595  if(rel==relt::LT || rel==relt::LE)
596  {
597  literalt bitwise_equal=bv_utils.equal(src1, src2);
598 
599  // signs different? trivial! Unless Zero.
600 
601  literalt signs_different=
602  prop.lxor(sign_bit(src1), sign_bit(src2));
603 
604  // as long as the signs match: compare like unsigned numbers
605 
606  // this works due to the BIAS
607  literalt less_than1=bv_utils.unsigned_less_than(src1, src2);
608 
609  // if both are negative (and not the same), need to turn around!
610  literalt less_than2=
611  prop.lxor(less_than1, prop.land(sign_bit(src1), sign_bit(src2)));
612 
613  literalt less_than3=
614  prop.lselect(signs_different,
615  sign_bit(src1),
616  less_than2);
617 
618  if(rel==relt::LT)
619  {
620  bvt and_bv;
621  and_bv.push_back(less_than3);
622  and_bv.push_back(!bitwise_equal); // for the case of two negative numbers
623  and_bv.push_back(!both_zero);
624  and_bv.push_back(!NaN);
625 
626  return prop.land(and_bv);
627  }
628  else if(rel==relt::LE)
629  {
630  bvt or_bv;
631  or_bv.push_back(less_than3);
632  or_bv.push_back(both_zero);
633  or_bv.push_back(bitwise_equal);
634 
635  return prop.land(prop.lor(or_bv), !NaN);
636  }
637  else
638  assert(false);
639  }
640  else if(rel==relt::EQ)
641  {
642  literalt bitwise_equal=bv_utils.equal(src1, src2);
643 
644  return prop.land(
645  prop.lor(bitwise_equal, both_zero),
646  !NaN);
647  }
648  else
649  assert(0);
650 
651  // not reached
652  return const_literal(false);
653 }
654 
656 {
657  assert(!src.empty());
658  bvt all_but_sign;
659  all_but_sign=src;
660  all_but_sign.resize(all_but_sign.size()-1);
661  return bv_utils.is_zero(all_but_sign);
662 }
663 
665 {
666  bvt and_bv;
667  and_bv.push_back(!sign_bit(src));
668  and_bv.push_back(exponent_all_ones(src));
669  and_bv.push_back(fraction_all_zeros(src));
670  return prop.land(and_bv);
671 }
672 
674 {
675  return prop.land(
676  exponent_all_ones(src),
677  fraction_all_zeros(src));
678 }
679 
682 {
683  return bv_utils.extract(src, spec.f, spec.f+spec.e-1);
684 }
685 
688 {
689  return bv_utils.extract(src, 0, spec.f-1);
690 }
691 
693 {
694  bvt and_bv;
695  and_bv.push_back(sign_bit(src));
696  and_bv.push_back(exponent_all_ones(src));
697  and_bv.push_back(fraction_all_zeros(src));
698  return prop.land(and_bv);
699 }
700 
702 {
703  return prop.land(exponent_all_ones(src),
704  !fraction_all_zeros(src));
705 }
706 
708 {
709  bvt exponent=src;
710 
711  // removes the fractional part
712  exponent.erase(exponent.begin(), exponent.begin()+spec.f);
713 
714  // removes the sign
715  exponent.resize(spec.e);
716 
717  return bv_utils.is_all_ones(exponent);
718 }
719 
721 {
722  bvt exponent=src;
723 
724  // removes the fractional part
725  exponent.erase(exponent.begin(), exponent.begin()+spec.f);
726 
727  // removes the sign
728  exponent.resize(spec.e);
729 
730  return bv_utils.is_zero(exponent);
731 }
732 
734 {
735  // does not include hidden bit
736  bvt tmp=src;
737  assert(src.size()==spec.width());
738  tmp.resize(spec.f);
739  return bv_utils.is_zero(tmp);
740 }
741 
743 void float_utilst::normalization_shift(bvt &fraction, bvt &exponent)
744 {
745  #if 0
746  // this thing is quadratic!
747 
748  bvt new_fraction=prop.new_variables(fraction.size());
749  bvt new_exponent=prop.new_variables(exponent.size());
750 
751  // i is the shift distance
752  for(std::size_t i=0; i<fraction.size(); i++)
753  {
754  bvt equal;
755 
756  // the bits above need to be zero
757  for(std::size_t j=0; j<i; j++)
758  equal.push_back(
759  !fraction[fraction.size()-1-j]);
760 
761  // this one needs to be one
762  equal.push_back(fraction[fraction.size()-1-i]);
763 
764  // iff all of that holds, we shift here!
765  literalt shift=prop.land(equal);
766 
767  // build shifted value
768  bvt shifted_fraction=bv_utils.shift(fraction, bv_utilst::LEFT, i);
769  bv_utils.cond_implies_equal(shift, shifted_fraction, new_fraction);
770 
771  // build new exponent
772  bvt adjustment=bv_utils.build_constant(-i, exponent.size());
773  bvt added_exponent=bv_utils.add(exponent, adjustment);
774  bv_utils.cond_implies_equal(shift, added_exponent, new_exponent);
775  }
776 
777  // Fraction all zero? It stays zero.
778  // The exponent is undefined in that case.
779  literalt fraction_all_zero=bv_utils.is_zero(fraction);
780  bvt zero_fraction;
781  zero_fraction.resize(fraction.size(), const_literal(false));
782  bv_utils.cond_implies_equal(fraction_all_zero, zero_fraction, new_fraction);
783 
784  fraction=new_fraction;
785  exponent=new_exponent;
786 
787  #else
788 
789  // n-log-n alignment shifter.
790  // The worst-case shift is the number of fraction
791  // bits minus one, in case the faction is one exactly.
792  assert(!fraction.empty());
793  unsigned depth=integer2unsigned(address_bits(fraction.size()-1));
794 
795  if(exponent.size()<depth)
796  exponent=bv_utils.sign_extension(exponent, depth);
797 
798  bvt exponent_delta=bv_utils.zeros(exponent.size());
799 
800  for(int d=depth-1; d>=0; d--)
801  {
802  std::size_t distance=(1<<d);
803  assert(fraction.size()>distance);
804 
805  // check if first 'distance'-many bits are zeros
806  const bvt prefix=bv_utils.extract_msb(fraction, distance);
807  literalt prefix_is_zero=bv_utils.is_zero(prefix);
808 
809  // If so, shift the zeros out left by 'distance'.
810  // Otherwise, leave as is.
811  const bvt shifted=
812  bv_utils.shift(fraction, bv_utilst::shiftt::LEFT, distance);
813 
814  fraction=
815  bv_utils.select(prefix_is_zero, shifted, fraction);
816 
817  // add corresponding weight to exponent
818  assert(d<(signed)exponent_delta.size());
819  exponent_delta[d]=prefix_is_zero;
820  }
821 
822  exponent=bv_utils.sub(exponent, exponent_delta);
823 
824  #endif
825 }
826 
828 void float_utilst::denormalization_shift(bvt &fraction, bvt &exponent)
829 {
831 
832  // Is the exponent strictly less than -bias+1, i.e., exponent<-bias+1?
833  // This is transformed to distance=(-bias+1)-exponent
834  // i.e., distance>0
835  // Note that 1-bias is the exponent represented by 0...01,
836  // i.e. the exponent of the smallest normal number and thus the 'base'
837  // exponent for subnormal numbers.
838 
839  assert(exponent.size()>=spec.e);
840 
841 #if 1
842  // Need to sign extend to avoid overflow. Note that this is a
843  // relatively rare problem as the value needs to be close to the top
844  // of the exponent range and then range must not have been
845  // previously extended as add, multiply, etc. do. This is primarily
846  // to handle casting down from larger ranges.
847  exponent=bv_utils.sign_extension(exponent, exponent.size() + 1);
848 #endif
849 
850  bvt distance=bv_utils.sub(
851  bv_utils.build_constant(-bias+1, exponent.size()), exponent);
852 
853  // use sign bit
854  literalt denormal=prop.land(
855  !distance.back(),
856  !bv_utils.is_zero(distance));
857 
858 #if 1
859  // Care must be taken to not loose information required for the
860  // guard and sticky bits. +3 is for the hidden, guard and sticky bits.
861  if(fraction.size() < (spec.f + 3))
862  {
863  // Add zeros at the LSB end for the guard bit to shift into
864  fraction=
865  bv_utils.concatenate(bv_utils.zeros((spec.f + 3) - fraction.size()),
866  fraction);
867  }
868 
869  bvt denormalisedFraction=fraction;
870 
871  literalt sticky_bit=const_literal(false);
872  denormalisedFraction =
873  sticky_right_shift(fraction, distance, sticky_bit);
874  denormalisedFraction[0]=prop.lor(denormalisedFraction[0], sticky_bit);
875 
876  fraction=
878  denormal,
879  denormalisedFraction,
880  fraction);
881 
882 #else
883  fraction=
885  denormal,
886  bv_utils.shift(fraction, bv_utilst::LRIGHT, distance),
887  fraction);
888 #endif
889 
890  exponent=
891  bv_utils.select(denormal,
892  bv_utils.build_constant(-bias, exponent.size()),
893  exponent);
894 }
895 
897 {
898  // incoming: some fraction (with explicit 1),
899  // some exponent without bias
900  // outgoing: rounded, with right size, with hidden bit, bias
901 
902  bvt aligned_fraction=src.fraction,
903  aligned_exponent=src.exponent;
904 
905  {
906  std::size_t exponent_bits=
907  std::max((std::size_t)integer2size_t(address_bits(spec.f)),
908  (std::size_t)spec.e)+1;
909 
910  // before normalization, make sure exponent is large enough
911  if(aligned_exponent.size()<exponent_bits)
912  {
913  // sign extend
914  aligned_exponent=
915  bv_utils.sign_extension(aligned_exponent, exponent_bits);
916  }
917  }
918 
919  // align it!
920  normalization_shift(aligned_fraction, aligned_exponent);
921  denormalization_shift(aligned_fraction, aligned_exponent);
922 
923  unbiased_floatt result;
924  result.fraction=aligned_fraction;
925  result.exponent=aligned_exponent;
926  result.sign=src.sign;
927  result.NaN=src.NaN;
928  result.infinity=src.infinity;
929 
930  round_fraction(result);
931  round_exponent(result);
932 
933  return pack(bias(result));
934 }
935 
938  const std::size_t dest_bits,
939  const literalt sign,
940  const bvt &fraction)
941 {
942  assert(dest_bits<fraction.size());
943 
944  // we have too many fraction bits
945  std::size_t extra_bits=fraction.size()-dest_bits;
946 
947  // more than two extra bits are superflus, and are
948  // turned into a sticky bit
949 
950  literalt sticky_bit=const_literal(false);
951 
952  if(extra_bits>=2)
953  {
954  // We keep most-significant bits, and thus the tail is made
955  // of least-significant bits.
956  bvt tail=bv_utils.extract(fraction, 0, extra_bits-2);
957  sticky_bit=prop.lor(tail);
958  }
959 
960  // the rounding bit is the last extra bit
961  assert(extra_bits>=1);
962  literalt rounding_bit=fraction[extra_bits-1];
963 
964  // we get one bit of the fraction for some rounding decisions
965  literalt rounding_least=fraction[extra_bits];
966 
967  // round-to-nearest (ties to even)
968  literalt round_to_even=
969  prop.land(rounding_bit,
970  prop.lor(rounding_least, sticky_bit));
971 
972  // round up
973  literalt round_to_plus_inf=
974  prop.land(!sign,
975  prop.lor(rounding_bit, sticky_bit));
976 
977  // round down
978  literalt round_to_minus_inf=
979  prop.land(sign,
980  prop.lor(rounding_bit, sticky_bit));
981 
982  // round to zero
983  literalt round_to_zero=
984  const_literal(false);
985 
986  // now select appropriate one
987  return prop.lselect(rounding_mode_bits.round_to_even, round_to_even,
991  prop.new_variable())))); // otherwise non-det
992 }
993 
995 {
996  std::size_t fraction_size=spec.f+1;
997 
998  // do we need to enlarge the fraction?
999  if(result.fraction.size()<fraction_size)
1000  {
1001  // pad with zeros at bottom
1002  std::size_t padding=fraction_size-result.fraction.size();
1003 
1004  result.fraction=bv_utils.concatenate(
1005  bv_utils.zeros(padding),
1006  result.fraction);
1007 
1008  assert(result.fraction.size()==fraction_size);
1009  }
1010  else if(result.fraction.size()==fraction_size) // it stays
1011  {
1012  // do nothing
1013  }
1014  else // fraction gets smaller -- rounding
1015  {
1016  std::size_t extra_bits=result.fraction.size()-fraction_size;
1017  assert(extra_bits>=1);
1018 
1019  // this computes the rounding decision
1021  fraction_size, result.sign, result.fraction);
1022 
1023  // chop off all the extra bits
1024  result.fraction=bv_utils.extract(
1025  result.fraction, extra_bits, result.fraction.size()-1);
1026 
1027  assert(result.fraction.size()==fraction_size);
1028 
1029 #if 0
1030  // *** does not catch when the overflow goes subnormal -> normal ***
1031  // incrementing the fraction might result in an overflow
1032  result.fraction=
1033  bv_utils.zero_extension(result.fraction, result.fraction.size()+1);
1034 
1035  result.fraction=bv_utils.incrementer(result.fraction, increment);
1036 
1037  literalt overflow=result.fraction.back();
1038 
1039  // In case of an overflow, the exponent has to be incremented.
1040  // "Post normalization" is then required.
1041  result.exponent=
1042  bv_utils.incrementer(result.exponent, overflow);
1043 
1044  // post normalization of the fraction
1045  literalt integer_part1=result.fraction.back();
1046  literalt integer_part0=result.fraction[result.fraction.size()-2];
1047  literalt new_integer_part=prop.lor(integer_part1, integer_part0);
1048 
1049  result.fraction.resize(result.fraction.size()-1);
1050  result.fraction.back()=new_integer_part;
1051 
1052 #else
1053  // When incrementing due to rounding there are two edge
1054  // cases we need to be aware of:
1055  // 1. If the number is normal, the increment can overflow.
1056  // In this case we need to increment the exponent and
1057  // set the MSB of the fraction to 1.
1058  // 2. If the number is the largest subnormal, the increment
1059  // can change the MSB making it normal. Thus the exponent
1060  // must be incremented but the fraction will be OK.
1061  literalt oldMSB=result.fraction.back();
1062 
1063  result.fraction=bv_utils.incrementer(result.fraction, increment);
1064 
1065  // Normal overflow when old MSB == 1 and new MSB == 0
1066  literalt overflow=prop.land(oldMSB, neg(result.fraction.back()));
1067 
1068  // Subnormal to normal transition when old MSB == 0 and new MSB == 1
1069  literalt subnormal_to_normal=
1070  prop.land(neg(oldMSB), result.fraction.back());
1071 
1072  // In case of an overflow or subnormal to normal conversion,
1073  // the exponent has to be incremented.
1074  result.exponent=
1075  bv_utils.incrementer(result.exponent,
1076  prop.lor(overflow, subnormal_to_normal));
1077 
1078  // post normalization of the fraction
1079  // In the case of overflow, set the MSB to 1
1080  // The subnormal case will have (only) the MSB set to 1
1081  result.fraction.back()=prop.lor(result.fraction.back(), overflow);
1082 #endif
1083  }
1084 }
1085 
1087 {
1088  // do we need to enlarge the exponent?
1089  if(result.exponent.size()<spec.e)
1090  {
1091  // should have been done before
1092  assert(false);
1093  }
1094  else if(result.exponent.size()==spec.e) // it stays
1095  {
1096  // do nothing
1097  }
1098  else // exponent gets smaller -- chop off top bits
1099  {
1100  bvt old_exponent=result.exponent;
1101  result.exponent.resize(spec.e);
1102 
1103  // max_exponent is the maximum representable
1104  // i.e. 1 higher than the maximum possible for a normal number
1105  bvt max_exponent=
1107  spec.max_exponent()-spec.bias(), old_exponent.size());
1108 
1109  // the exponent is garbage if the fractional is zero
1110 
1111  literalt exponent_too_large=
1112  prop.land(
1113  !bv_utils.signed_less_than(old_exponent, max_exponent),
1114  !bv_utils.is_zero(result.fraction));
1115 
1116 #if 1
1117  // Directed rounding modes round overflow to the maximum normal
1118  // depending on the particular mode and the sign
1119  literalt overflow_to_inf=
1122  !result.sign),
1124  result.sign)));
1125 
1126  literalt set_to_max=
1127  prop.land(exponent_too_large, !overflow_to_inf);
1128 
1129 
1130  bvt largest_normal_exponent=
1132  spec.max_exponent()-(spec.bias() + 1), result.exponent.size());
1133 
1134  result.exponent=
1135  bv_utils.select(set_to_max, largest_normal_exponent, result.exponent);
1136 
1137  result.fraction=
1138  bv_utils.select(set_to_max,
1139  bv_utils.inverted(bv_utils.zeros(result.fraction.size())),
1140  result.fraction);
1141 
1142  result.infinity=prop.lor(result.infinity,
1143  prop.land(exponent_too_large,
1144  overflow_to_inf));
1145 #else
1146  result.infinity=prop.lor(result.infinity, exponent_too_large);
1147 #endif
1148  }
1149 }
1150 
1153 {
1154  biased_floatt result;
1155 
1156  result.sign=src.sign;
1157  result.NaN=src.NaN;
1158  result.infinity=src.infinity;
1159 
1160  // we need to bias the new exponent
1161  result.exponent=add_bias(src.exponent);
1162 
1163  // strip off hidden bit
1164  assert(src.fraction.size()==spec.f+1);
1165 
1166  literalt hidden_bit=src.fraction[src.fraction.size()-1];
1167  literalt denormal=!hidden_bit;
1168 
1169  result.fraction=src.fraction;
1170  result.fraction.resize(spec.f);
1171 
1172  // make exponent zero if its denormal
1173  // (includes zero)
1174  for(std::size_t i=0; i<result.exponent.size(); i++)
1175  result.exponent[i]=
1176  prop.land(result.exponent[i], !denormal);
1177 
1178  return result;
1179 }
1180 
1182 {
1183  assert(src.size()==spec.e);
1184 
1185  return bv_utils.add(
1186  src,
1188 }
1189 
1191 {
1192  assert(src.size()==spec.e);
1193 
1194  return bv_utils.sub(
1195  src,
1197 }
1198 
1200 {
1201  assert(src.size()==spec.width());
1202 
1203  unbiased_floatt result;
1204 
1205  result.sign=sign_bit(src);
1206 
1207  result.fraction=get_fraction(src);
1208  result.fraction.push_back(is_normal(src)); // add hidden bit
1209 
1210  result.exponent=get_exponent(src);
1211  assert(result.exponent.size()==spec.e);
1212 
1213  // unbias the exponent
1214  literalt denormal=bv_utils.is_zero(result.exponent);
1215 
1216  result.exponent=
1217  bv_utils.select(denormal,
1219  sub_bias(result.exponent));
1220 
1221  result.infinity=is_infinity(src);
1222  result.zero=is_zero(src);
1223  result.NaN=is_NaN(src);
1224 
1225  return result;
1226 }
1227 
1229 {
1230  assert(src.fraction.size()==spec.f);
1231  assert(src.exponent.size()==spec.e);
1232 
1233  bvt result;
1234  result.resize(spec.width());
1235 
1236  // do sign
1237  // we make this 'false' for NaN
1238  result[result.size()-1]=
1239  prop.lselect(src.NaN, const_literal(false), src.sign);
1240 
1241  literalt infinity_or_NaN=
1242  prop.lor(src.NaN, src.infinity);
1243 
1244  // just copy fraction
1245  for(std::size_t i=0; i<spec.f; i++)
1246  result[i]=prop.land(src.fraction[i], !infinity_or_NaN);
1247 
1248  result[0]=prop.lor(result[0], src.NaN);
1249 
1250  // do exponent
1251  for(std::size_t i=0; i<spec.e; i++)
1252  result[i+spec.f]=prop.lor(
1253  src.exponent[i],
1254  infinity_or_NaN);
1255 
1256  return result;
1257 }
1258 
1260 {
1261  mp_integer int_value=0;
1262 
1263  for(std::size_t i=0; i<src.size(); i++)
1264  int_value+=power(2, i)*prop.l_get(src[i]).is_true();
1265 
1266  ieee_floatt result;
1267  result.spec=spec;
1268  result.unpack(int_value);
1269 
1270  return result;
1271 }
1272 
1274  const bvt &op,
1275  const bvt &dist,
1276  literalt &sticky)
1277 {
1278  std::size_t d=1;
1279  bvt result=op;
1280  sticky=const_literal(false);
1281 
1282  for(std::size_t stage=0; stage<dist.size(); stage++)
1283  {
1284  if(dist[stage]!=const_literal(false))
1285  {
1286  bvt tmp=bv_utils.shift(result, bv_utilst::shiftt::LRIGHT, d);
1287 
1288  bvt lost_bits;
1289 
1290  if(d<=result.size())
1291  lost_bits=bv_utils.extract(result, 0, d-1);
1292  else
1293  lost_bits=result;
1294 
1295  sticky=prop.lor(
1296  prop.land(dist[stage], prop.lor(lost_bits)),
1297  sticky);
1298 
1299  result=bv_utils.select(dist[stage], tmp, result);
1300  }
1301 
1302  d=d<<1;
1303  }
1304 
1305  return result;
1306 }
1307 
1309  const bvt &src1,
1310  const bvt &src2)
1311 {
1312  return src1;
1313 }
1314 
1316  const bvt &op0,
1317  const bvt &op1)
1318 {
1319  return op0;
1320 }
bool get_sign() const
Definition: ieee_float.h:225
bool is_signed(const typet &t)
Convenience function – is the type signed?
Definition: util.cpp:45
ieee_floatt get(const bvt &) const
bvt inverted(const bvt &op)
Definition: bv_utils.cpp:561
BigInt mp_integer
Definition: mp_arith.h:19
static bvt extract(const bvt &a, std::size_t first, std::size_t last)
Definition: bv_utils.cpp:42
literalt is_minus_inf(const bvt &)
void round_exponent(unbiased_floatt &result)
bvt cond_negate(const bvt &bv, const literalt cond)
Definition: bv_utils.cpp:741
literalt equal(const bvt &op0, const bvt &op1)
Bit-blasting ID_equal and use in other encodings.
Definition: bv_utils.cpp:1089
bvt new_variables(std::size_t width)
generates a bitvector of given width with new variables
Definition: prop.cpp:39
bvt sub_bias(const bvt &exponent)
virtual void normalization_shift(bvt &fraction, bvt &exponent)
normalize fraction/exponent pair returns &#39;zero&#39; if fraction is zero
bvt to_signed_integer(const bvt &src, std::size_t int_width)
Definition: float_utils.cpp:68
virtual literalt lor(literalt a, literalt b)=0
mp_integer address_bits(const mp_integer &size)
ceil(log2(size))
bvt sub(const bvt &op0, const bvt &op1)
Definition: bv_utils.h:65
mp_integer max_exponent() const
Definition: ieee_float.cpp:36
bool is_NaN() const
Definition: ieee_float.h:226
literalt is_zero(const bvt &op)
Definition: bv_utils.h:137
const mp_integer & get_exponent() const
Definition: ieee_float.h:230
void denormalization_shift(bvt &fraction, bvt &exponent)
make sure exponent is not too small; the exponent is unbiased
virtual bvt rem(const bvt &src1, const bvt &src2)
literalt is_normal(const bvt &)
bool is_infinity() const
Definition: ieee_float.h:227
literalt exponent_all_ones(const bvt &)
virtual literalt land(literalt a, literalt b)=0
literalt is_all_ones(const bvt &op)
Definition: bv_utils.h:152
virtual literalt new_variable()=0
void unpack(const mp_integer &i)
Definition: ieee_float.cpp:312
bvt concatenate(const bvt &a, const bvt &b) const
Definition: bv_utils.cpp:80
bvt to_unsigned_integer(const bvt &src, std::size_t int_width)
Definition: float_utils.cpp:75
bvt add_bias(const bvt &exponent)
bvt negate(const bvt &)
propt & prop
Definition: float_utils.h:152
virtual literalt lxor(literalt a, literalt b)=0
virtual bvt mul(const bvt &src1, const bvt &src2)
bvt debug1(const bvt &op0, const bvt &op1)
bool is_true() const
Definition: literal.h:155
bvt conversion(const bvt &src, const ieee_float_spect &dest_spec)
literalt is_not_zero(const bvt &op)
Definition: bv_utils.h:140
static bvt extract_msb(const bvt &a, std::size_t n)
Definition: bv_utils.cpp:58
bvt get_exponent(const bvt &)
Gets the unbiased exponent in a floating-point bit-vector.
literalt is_zero(const bvt &)
static literalt sign_bit(const bvt &src)
Definition: float_utils.h:90
bvt absolute_value(const bvt &op)
Definition: bv_utils.cpp:754
ieee_float_spect spec
Definition: ieee_float.h:123
literalt is_plus_inf(const bvt &)
bvt unsigned_multiplier(const bvt &op0, const bvt &op1)
Definition: bv_utils.cpp:614
literalt signed_less_than(const bvt &bv0, const bvt &bv1)
Definition: bv_utils.cpp:1255
bvt select(literalt s, const bvt &a, const bvt &b)
If s is true, selects a otherwise selects b.
Definition: bv_utils.cpp:96
unbiased_floatt unpack(const bvt &)
virtual bvt div(const bvt &src1, const bvt &src2)
literalt fraction_rounding_decision(const std::size_t dest_bits, const literalt sign, const bvt &fraction)
rounding decision for fraction using sticky bit
mp_integer bias() const
Definition: ieee_float.cpp:21
bvt get_fraction(const bvt &)
Gets the fraction without hidden bit in a floating-point bit-vector src.
bvt zero_extension(const bvt &bv, std::size_t new_size)
Definition: bv_utils.h:180
virtual bvt add_sub(const bvt &src1, const bvt &src2, bool subtract)
bvt debug2(const bvt &op0, const bvt &op1)
void unsigned_divider(const bvt &op0, const bvt &op1, bvt &res, bvt &rem)
Definition: bv_utils.cpp:873
bool is_true() const
Definition: threeval.h:25
bvt subtract_exponents(const unbiased_floatt &src1, const unbiased_floatt &src2)
Subtracts the exponents.
bvt incrementer(const bvt &op, literalt carry_in)
Definition: bv_utils.cpp:553
literalt exponent_all_zeros(const bvt &)
void set_rounding_mode(const bvt &)
Definition: float_utils.cpp:16
bvt shift(const bvt &op, const shiftt shift, std::size_t distance)
Definition: bv_utils.cpp:480
literalt unsigned_less_than(const bvt &bv0, const bvt &bv1)
Definition: bv_utils.cpp:1243
literalt is_infinity(const bvt &)
unsigned integer2unsigned(const mp_integer &n)
Definition: mp_arith.cpp:203
literalt const_literal(bool value)
Definition: literal.h:187
literalt is_NaN(const bvt &)
bvt add(const bvt &op0, const bvt &op1)
Definition: bv_utils.h:64
ieee_float_spect spec
Definition: float_utils.h:86
bvt from_unsigned_integer(const bvt &)
Definition: float_utils.cpp:51
rounding_mode_bitst rounding_mode_bits
Definition: float_utils.h:65
bvt to_integer(const bvt &src, std::size_t int_width, bool is_signed)
Definition: float_utils.cpp:82
bvt abs(const bvt &)
bv_utilst bv_utils
Definition: float_utils.h:153
biased_floatt bias(const unbiased_floatt &)
takes an unbiased float, and applies the bias
bvt inc(const bvt &op)
Definition: bv_utils.h:36
virtual tvt l_get(literalt a) const =0
literalt neg(literalt a)
Definition: literal.h:192
bvt limit_distance(const bvt &dist, mp_integer limit)
Limits the shift distance.
void cond_implies_equal(literalt cond, const bvt &a, const bvt &b)
Definition: bv_utils.cpp:1293
bvt from_signed_integer(const bvt &)
Definition: float_utils.cpp:33
bvt sign_extension(const bvt &bv, std::size_t new_size)
Definition: bv_utils.h:175
std::size_t width() const
Definition: ieee_float.h:50
std::size_t f
Definition: ieee_float.h:26
bvt sub(const bvt &src1, const bvt &src2)
Definition: float_utils.h:114
virtual literalt lselect(literalt a, literalt b, literalt c)=0
std::size_t integer2size_t(const mp_integer &n)
Definition: mp_arith.cpp:195
bvt sticky_right_shift(const bvt &op, const bvt &dist, literalt &sticky)
bvt add_sub(const bvt &op0, const bvt &op1, bool subtract)
Definition: bv_utils.cpp:339
std::size_t e
Definition: ieee_float.h:26
std::vector< literalt > bvt
Definition: literal.h:197
bvt build_constant(const mp_integer &i, std::size_t width)
Definition: bv_utils.cpp:15
bvt zeros(std::size_t new_size) const
Definition: bv_utils.h:185
mp_integer power(const mp_integer &base, const mp_integer &exponent)
A multi-precision implementation of the power operator.
const mp_integer & get_fraction() const
Definition: ieee_float.h:231
virtual bvt rounder(const unbiased_floatt &)
bvt pack(const biased_floatt &)
void round_fraction(unbiased_floatt &result)
bvt build_constant(const ieee_floatt &)
literalt fraction_all_zeros(const bvt &)
literalt relation(const bvt &src1, relt rel, const bvt &src2)