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