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