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