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) + 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) + 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::SHIFT_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 = 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  std::size_t depth = 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::SHIFT_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 = std::max(address_bits(spec.f), spec.e) + 1;
907 
908  // before normalization, make sure exponent is large enough
909  if(aligned_exponent.size()<exponent_bits)
910  {
911  // sign extend
912  aligned_exponent=
913  bv_utils.sign_extension(aligned_exponent, exponent_bits);
914  }
915  }
916 
917  // align it!
918  normalization_shift(aligned_fraction, aligned_exponent);
919  denormalization_shift(aligned_fraction, aligned_exponent);
920 
921  unbiased_floatt result;
922  result.fraction=aligned_fraction;
923  result.exponent=aligned_exponent;
924  result.sign=src.sign;
925  result.NaN=src.NaN;
926  result.infinity=src.infinity;
927 
928  round_fraction(result);
929  round_exponent(result);
930 
931  return pack(bias(result));
932 }
933 
936  const std::size_t dest_bits,
937  const literalt sign,
938  const bvt &fraction)
939 {
940  assert(dest_bits<fraction.size());
941 
942  // we have too many fraction bits
943  std::size_t extra_bits=fraction.size()-dest_bits;
944 
945  // more than two extra bits are superflus, and are
946  // turned into a sticky bit
947 
948  literalt sticky_bit=const_literal(false);
949 
950  if(extra_bits>=2)
951  {
952  // We keep most-significant bits, and thus the tail is made
953  // of least-significant bits.
954  bvt tail=bv_utils.extract(fraction, 0, extra_bits-2);
955  sticky_bit=prop.lor(tail);
956  }
957 
958  // the rounding bit is the last extra bit
959  assert(extra_bits>=1);
960  literalt rounding_bit=fraction[extra_bits-1];
961 
962  // we get one bit of the fraction for some rounding decisions
963  literalt rounding_least=fraction[extra_bits];
964 
965  // round-to-nearest (ties to even)
966  literalt round_to_even=
967  prop.land(rounding_bit,
968  prop.lor(rounding_least, sticky_bit));
969 
970  // round up
971  literalt round_to_plus_inf=
972  prop.land(!sign,
973  prop.lor(rounding_bit, sticky_bit));
974 
975  // round down
976  literalt round_to_minus_inf=
977  prop.land(sign,
978  prop.lor(rounding_bit, sticky_bit));
979 
980  // round to zero
981  literalt round_to_zero=
982  const_literal(false);
983 
984  // now select appropriate one
985  return prop.lselect(rounding_mode_bits.round_to_even, round_to_even,
989  prop.new_variable())))); // otherwise non-det
990 }
991 
993 {
994  std::size_t fraction_size=spec.f+1;
995 
996  // do we need to enlarge the fraction?
997  if(result.fraction.size()<fraction_size)
998  {
999  // pad with zeros at bottom
1000  std::size_t padding=fraction_size-result.fraction.size();
1001 
1002  result.fraction=bv_utils.concatenate(
1003  bv_utils.zeros(padding),
1004  result.fraction);
1005 
1006  assert(result.fraction.size()==fraction_size);
1007  }
1008  else if(result.fraction.size()==fraction_size) // it stays
1009  {
1010  // do nothing
1011  }
1012  else // fraction gets smaller -- rounding
1013  {
1014  std::size_t extra_bits=result.fraction.size()-fraction_size;
1015  assert(extra_bits>=1);
1016 
1017  // this computes the rounding decision
1019  fraction_size, result.sign, result.fraction);
1020 
1021  // chop off all the extra bits
1022  result.fraction=bv_utils.extract(
1023  result.fraction, extra_bits, result.fraction.size()-1);
1024 
1025  assert(result.fraction.size()==fraction_size);
1026 
1027 #if 0
1028  // *** does not catch when the overflow goes subnormal -> normal ***
1029  // incrementing the fraction might result in an overflow
1030  result.fraction=
1031  bv_utils.zero_extension(result.fraction, result.fraction.size()+1);
1032 
1033  result.fraction=bv_utils.incrementer(result.fraction, increment);
1034 
1035  literalt overflow=result.fraction.back();
1036 
1037  // In case of an overflow, the exponent has to be incremented.
1038  // "Post normalization" is then required.
1039  result.exponent=
1040  bv_utils.incrementer(result.exponent, overflow);
1041 
1042  // post normalization of the fraction
1043  literalt integer_part1=result.fraction.back();
1044  literalt integer_part0=result.fraction[result.fraction.size()-2];
1045  literalt new_integer_part=prop.lor(integer_part1, integer_part0);
1046 
1047  result.fraction.resize(result.fraction.size()-1);
1048  result.fraction.back()=new_integer_part;
1049 
1050 #else
1051  // When incrementing due to rounding there are two edge
1052  // cases we need to be aware of:
1053  // 1. If the number is normal, the increment can overflow.
1054  // In this case we need to increment the exponent and
1055  // set the MSB of the fraction to 1.
1056  // 2. If the number is the largest subnormal, the increment
1057  // can change the MSB making it normal. Thus the exponent
1058  // must be incremented but the fraction will be OK.
1059  literalt oldMSB=result.fraction.back();
1060 
1061  result.fraction=bv_utils.incrementer(result.fraction, increment);
1062 
1063  // Normal overflow when old MSB == 1 and new MSB == 0
1064  literalt overflow=prop.land(oldMSB, neg(result.fraction.back()));
1065 
1066  // Subnormal to normal transition when old MSB == 0 and new MSB == 1
1067  literalt subnormal_to_normal=
1068  prop.land(neg(oldMSB), result.fraction.back());
1069 
1070  // In case of an overflow or subnormal to normal conversion,
1071  // the exponent has to be incremented.
1072  result.exponent=
1073  bv_utils.incrementer(result.exponent,
1074  prop.lor(overflow, subnormal_to_normal));
1075 
1076  // post normalization of the fraction
1077  // In the case of overflow, set the MSB to 1
1078  // The subnormal case will have (only) the MSB set to 1
1079  result.fraction.back()=prop.lor(result.fraction.back(), overflow);
1080 #endif
1081  }
1082 }
1083 
1085 {
1086  // do we need to enlarge the exponent?
1087  if(result.exponent.size()<spec.e)
1088  {
1089  // should have been done before
1090  assert(false);
1091  }
1092  else if(result.exponent.size()==spec.e) // it stays
1093  {
1094  // do nothing
1095  }
1096  else // exponent gets smaller -- chop off top bits
1097  {
1098  bvt old_exponent=result.exponent;
1099  result.exponent.resize(spec.e);
1100 
1101  // max_exponent is the maximum representable
1102  // i.e. 1 higher than the maximum possible for a normal number
1103  bvt max_exponent=
1105  spec.max_exponent()-spec.bias(), old_exponent.size());
1106 
1107  // the exponent is garbage if the fractional is zero
1108 
1109  literalt exponent_too_large=
1110  prop.land(
1111  !bv_utils.signed_less_than(old_exponent, max_exponent),
1112  !bv_utils.is_zero(result.fraction));
1113 
1114 #if 1
1115  // Directed rounding modes round overflow to the maximum normal
1116  // depending on the particular mode and the sign
1117  literalt overflow_to_inf=
1120  !result.sign),
1122  result.sign)));
1123 
1124  literalt set_to_max=
1125  prop.land(exponent_too_large, !overflow_to_inf);
1126 
1127 
1128  bvt largest_normal_exponent=
1130  spec.max_exponent()-(spec.bias() + 1), result.exponent.size());
1131 
1132  result.exponent=
1133  bv_utils.select(set_to_max, largest_normal_exponent, result.exponent);
1134 
1135  result.fraction=
1136  bv_utils.select(set_to_max,
1137  bv_utils.inverted(bv_utils.zeros(result.fraction.size())),
1138  result.fraction);
1139 
1140  result.infinity=prop.lor(result.infinity,
1141  prop.land(exponent_too_large,
1142  overflow_to_inf));
1143 #else
1144  result.infinity=prop.lor(result.infinity, exponent_too_large);
1145 #endif
1146  }
1147 }
1148 
1151 {
1152  biased_floatt result;
1153 
1154  result.sign=src.sign;
1155  result.NaN=src.NaN;
1156  result.infinity=src.infinity;
1157 
1158  // we need to bias the new exponent
1159  result.exponent=add_bias(src.exponent);
1160 
1161  // strip off hidden bit
1162  assert(src.fraction.size()==spec.f+1);
1163 
1164  literalt hidden_bit=src.fraction[src.fraction.size()-1];
1165  literalt denormal=!hidden_bit;
1166 
1167  result.fraction=src.fraction;
1168  result.fraction.resize(spec.f);
1169 
1170  // make exponent zero if its denormal
1171  // (includes zero)
1172  for(std::size_t i=0; i<result.exponent.size(); i++)
1173  result.exponent[i]=
1174  prop.land(result.exponent[i], !denormal);
1175 
1176  return result;
1177 }
1178 
1180 {
1181  assert(src.size()==spec.e);
1182 
1183  return bv_utils.add(
1184  src,
1186 }
1187 
1189 {
1190  assert(src.size()==spec.e);
1191 
1192  return bv_utils.sub(
1193  src,
1195 }
1196 
1198 {
1199  assert(src.size()==spec.width());
1200 
1201  unbiased_floatt result;
1202 
1203  result.sign=sign_bit(src);
1204 
1205  result.fraction=get_fraction(src);
1206  result.fraction.push_back(is_normal(src)); // add hidden bit
1207 
1208  result.exponent=get_exponent(src);
1209  assert(result.exponent.size()==spec.e);
1210 
1211  // unbias the exponent
1212  literalt denormal=bv_utils.is_zero(result.exponent);
1213 
1214  result.exponent=
1215  bv_utils.select(denormal,
1217  sub_bias(result.exponent));
1218 
1219  result.infinity=is_infinity(src);
1220  result.zero=is_zero(src);
1221  result.NaN=is_NaN(src);
1222 
1223  return result;
1224 }
1225 
1227 {
1228  assert(src.fraction.size()==spec.f);
1229  assert(src.exponent.size()==spec.e);
1230 
1231  bvt result;
1232  result.resize(spec.width());
1233 
1234  // do sign
1235  // we make this 'false' for NaN
1236  result[result.size()-1]=
1237  prop.lselect(src.NaN, const_literal(false), src.sign);
1238 
1239  literalt infinity_or_NaN=
1240  prop.lor(src.NaN, src.infinity);
1241 
1242  // just copy fraction
1243  for(std::size_t i=0; i<spec.f; i++)
1244  result[i]=prop.land(src.fraction[i], !infinity_or_NaN);
1245 
1246  result[0]=prop.lor(result[0], src.NaN);
1247 
1248  // do exponent
1249  for(std::size_t i=0; i<spec.e; i++)
1250  result[i+spec.f]=prop.lor(
1251  src.exponent[i],
1252  infinity_or_NaN);
1253 
1254  return result;
1255 }
1256 
1258 {
1259  mp_integer int_value=0;
1260 
1261  for(std::size_t i=0; i<src.size(); i++)
1262  int_value+=power(2, i)*prop.l_get(src[i]).is_true();
1263 
1264  ieee_floatt result;
1265  result.spec=spec;
1266  result.unpack(int_value);
1267 
1268  return result;
1269 }
1270 
1272  const bvt &op,
1273  const bvt &dist,
1274  literalt &sticky)
1275 {
1276  std::size_t d=1;
1277  bvt result=op;
1278  sticky=const_literal(false);
1279 
1280  for(std::size_t stage=0; stage<dist.size(); stage++)
1281  {
1282  if(dist[stage]!=const_literal(false))
1283  {
1285 
1286  bvt lost_bits;
1287 
1288  if(d<=result.size())
1289  lost_bits=bv_utils.extract(result, 0, d-1);
1290  else
1291  lost_bits=result;
1292 
1293  sticky=prop.lor(
1294  prop.land(dist[stage], prop.lor(lost_bits)),
1295  sticky);
1296 
1297  result=bv_utils.select(dist[stage], tmp, result);
1298  }
1299 
1300  d=d<<1;
1301  }
1302 
1303  return result;
1304 }
1305 
1307  const bvt &src1,
1308  const bvt &src2)
1309 {
1310  return src1;
1311 }
1312 
1314  const bvt &op0,
1315  const bvt &op1)
1316 {
1317  return op0;
1318 }
bool get_sign() const
Definition: ieee_float.h:236
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:580
BigInt mp_integer
Definition: mp_arith.h:22
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:760
literalt equal(const bvt &op0, const bvt &op1)
Bit-blasting ID_equal and use in other encodings.
Definition: bv_utils.cpp:1108
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
std::size_t 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:237
literalt is_zero(const bvt &op)
Definition: bv_utils.h:141
const mp_integer & get_exponent() const
Definition: ieee_float.h:241
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:238
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:156
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:144
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:773
ieee_float_spect spec
Definition: ieee_float.h:134
literalt is_plus_inf(const bvt &)
bvt unsigned_multiplier(const bvt &op0, const bvt &op1)
Definition: bv_utils.cpp:633
literalt signed_less_than(const bvt &bv0, const bvt &bv1)
Definition: bv_utils.cpp:1274
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:184
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:892
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:572
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:1262
literalt is_infinity(const bvt &)
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:1312
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:179
std::size_t width() const
Definition: ieee_float.h:54
std::size_t f
Definition: ieee_float.h:30
bvt sub(const bvt &src1, const bvt &src2)
Definition: float_utils.h:114
virtual literalt lselect(literalt a, literalt b, literalt c)=0
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:30
std::vector< literalt > bvt
Definition: literal.h:200
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:189
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:242
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)