CBMC
ieee_float.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 "ieee_float.h"
10 
11 #include "arith_tools.h"
12 #include "bitvector_types.h"
13 #include "floatbv_expr.h"
14 #include "invariant.h"
15 #include "std_expr.h"
16 
17 #include <cstdint>
18 #include <limits>
19 
21 {
22  return power(2, e-1)-1;
23 }
24 
26 {
27  floatbv_typet result;
28  result.set_f(f);
29  result.set_width(width());
30  if(x86_extended)
31  result.set(ID_x86_extended, true);
32  return result;
33 }
34 
36 {
37  return power(2, e)-1;
38 }
39 
41 {
42  return power(2, f)-1;
43 }
44 
46 {
47  std::size_t width=type.get_width();
48  f=type.get_f();
49  DATA_INVARIANT(f != 0, "mantissa must be at least 1 bit");
51  f < width,
52  "mantissa bits must be less than "
53  "originating type width");
54  e=width-f-1;
55  x86_extended=type.get_bool(ID_x86_extended);
56  if(x86_extended)
57  e=e-1; // no hidden bit
58 }
59 
61 {
62  return floatbv_rounding_mode(static_cast<unsigned>(rm));
63 }
64 
65 void ieee_floatt::print(std::ostream &out) const
66 {
67  out << to_ansi_c_string();
68 }
69 
70 std::string ieee_floatt::format(const format_spect &format_spec) const
71 {
72  std::string result;
73 
74  switch(format_spec.style)
75  {
77  result+=to_string_decimal(format_spec.precision);
78  break;
79 
81  result+=to_string_scientific(format_spec.precision);
82  break;
83 
85  {
86  // On Linux, the man page says:
87  // "Style e is used if the exponent from its conversion
88  // is less than -4 or greater than or equal to the precision."
89  //
90  // On BSD, it's
91  // "The argument is printed in style f (F) or in style e (E)
92  // whichever gives full precision in minimum space."
93 
94  mp_integer _exponent, _fraction;
95  extract_base10(_fraction, _exponent);
96 
97  mp_integer adjusted_exponent=base10_digits(_fraction)+_exponent;
98 
99  if(adjusted_exponent>=format_spec.precision ||
100  adjusted_exponent<-4)
101  {
102  result+=to_string_scientific(format_spec.precision);
103  }
104  else
105  {
106  result+=to_string_decimal(format_spec.precision);
107 
108  // Implementations tested also appear to suppress trailing zeros
109  // and trailing dots.
110 
111  {
112  std::size_t trunc_pos=result.find_last_not_of('0');
113  if(trunc_pos!=std::string::npos)
114  result.resize(trunc_pos+1);
115  }
116 
117  if(!result.empty() && result.back()=='.')
118  result.resize(result.size()-1);
119  }
120  }
121  break;
122  }
123 
124  while(result.size()<format_spec.min_width)
125  result=" "+result;
126 
127  return result;
128 }
129 
131 {
132  PRECONDITION(src >= 0);
133  mp_integer tmp=src;
134  mp_integer result=0;
135  while(tmp!=0) { ++result; tmp/=10; }
136  return result;
137 }
138 
139 std::string ieee_floatt::to_string_decimal(std::size_t precision) const
140 {
141  std::string result;
142 
143  if(sign_flag)
144  result+='-';
145 
146  if((NaN_flag || infinity_flag) && !sign_flag)
147  result+='+';
148 
149  // special cases
150  if(NaN_flag)
151  result+="NaN";
152  else if(infinity_flag)
153  result+="inf";
154  else if(is_zero())
155  {
156  result+='0';
157 
158  // add zeros, if needed
159  if(precision>0)
160  {
161  result+='.';
162  for(std::size_t i=0; i<precision; i++)
163  result+='0';
164  }
165  }
166  else
167  {
168  mp_integer _exponent, _fraction;
169  extract_base2(_fraction, _exponent);
170 
171  // convert to base 10
172  if(_exponent>=0)
173  {
174  result+=integer2string(_fraction*power(2, _exponent));
175 
176  // add dot and zeros, if needed
177  if(precision>0)
178  {
179  result+='.';
180  for(std::size_t i=0; i<precision; i++)
181  result+='0';
182  }
183  }
184  else
185  {
186  #if 1
187  mp_integer position=-_exponent;
188 
189  // 10/2=5 -- this makes it base 10
190  _fraction*=power(5, position);
191 
192  // apply rounding
193  if(position>precision)
194  {
195  mp_integer r=power(10, position-precision);
196  mp_integer remainder=_fraction%r;
197  _fraction/=r;
198  // not sure if this is the right kind of rounding here
199  if(remainder>=r/2)
200  ++_fraction;
201  position=precision;
202  }
203 
204  std::string tmp=integer2string(_fraction);
205 
206  // pad with zeros from the front, if needed
207  while(mp_integer(tmp.size())<=position) tmp="0"+tmp;
208 
209  const std::size_t dot =
210  tmp.size() - numeric_cast_v<std::size_t>(position);
211  result+=std::string(tmp, 0, dot)+'.';
212  result+=std::string(tmp, dot, std::string::npos);
213 
214  // append zeros if needed
215  for(mp_integer i=position; i<precision; ++i)
216  result+='0';
217  #else
218 
219  result+=integer2string(_fraction);
220 
221  if(_exponent!=0)
222  result+="*2^"+integer2string(_exponent);
223 
224  #endif
225  }
226  }
227 
228  return result;
229 }
230 
233 std::string ieee_floatt::to_string_scientific(std::size_t precision) const
234 {
235  std::string result;
236 
237  if(sign_flag)
238  result+='-';
239 
240  if((NaN_flag || infinity_flag) && !sign_flag)
241  result+='+';
242 
243  // special cases
244  if(NaN_flag)
245  result+="NaN";
246  else if(infinity_flag)
247  result+="inf";
248  else if(is_zero())
249  {
250  result+='0';
251 
252  // add zeros, if needed
253  if(precision>0)
254  {
255  result+='.';
256  for(std::size_t i=0; i<precision; i++)
257  result+='0';
258  }
259 
260  result+="e0";
261  }
262  else
263  {
264  mp_integer _exponent, _fraction;
265  extract_base10(_fraction, _exponent);
266 
267  // C99 appears to say that conversion to decimal should
268  // use the currently selected IEEE rounding mode.
269  if(base10_digits(_fraction)>precision+1)
270  {
271  // re-align
272  mp_integer distance=base10_digits(_fraction)-(precision+1);
273  mp_integer p=power(10, distance);
274  mp_integer remainder=_fraction%p;
275  _fraction/=p;
276  _exponent+=distance;
277 
278  if(remainder==p/2)
279  {
280  // need to do rounding mode here
281  ++_fraction;
282  }
283  else if(remainder>p/2)
284  ++_fraction;
285  }
286 
287  std::string decimals=integer2string(_fraction);
288 
289  CHECK_RETURN(!decimals.empty());
290 
291  // First add top digit to result.
292  result+=decimals[0];
293 
294  // Now add dot and further zeros, if needed.
295  if(precision>0)
296  {
297  result+='.';
298 
299  while(decimals.size()<precision+1)
300  decimals+='0';
301 
302  result+=decimals.substr(1, precision);
303  }
304 
305  // add exponent
306  result+='e';
307 
308  std::string exponent_str=
309  integer2string(base10_digits(_fraction)+_exponent-1);
310 
311  if(!exponent_str.empty() && exponent_str[0]!='-')
312  result+='+';
313 
314  result+=exponent_str;
315  }
316 
317  return result;
318 }
319 
321 {
322  PRECONDITION(spec.f != 0);
323  PRECONDITION(spec.e != 0);
324 
325  {
326  mp_integer tmp=i;
327 
328  // split this apart
329  mp_integer pf=power(2, spec.f);
330  fraction=tmp%pf;
331  tmp/=pf;
332 
333  mp_integer pe=power(2, spec.e);
334  exponent=tmp%pe;
335  tmp/=pe;
336 
337  sign_flag=(tmp!=0);
338  }
339 
340  // NaN?
341  if(exponent==spec.max_exponent() && fraction!=0)
342  {
343  make_NaN();
344  }
345  else if(exponent==spec.max_exponent() && fraction==0) // Infinity
346  {
347  NaN_flag=false;
348  infinity_flag=true;
349  }
350  else if(exponent==0 && fraction==0) // zero
351  {
352  NaN_flag=false;
353  infinity_flag=false;
354  }
355  else if(exponent==0) // denormal?
356  {
357  NaN_flag=false;
358  infinity_flag=false;
359  exponent=-spec.bias()+1; // NOT -spec.bias()!
360  }
361  else // normal
362  {
363  NaN_flag=false;
364  infinity_flag=false;
365  fraction+=power(2, spec.f); // hidden bit!
366  exponent-=spec.bias(); // un-bias
367  }
368 }
369 
371 {
372  return fraction>=power(2, spec.f);
373 }
374 
376 {
377  mp_integer result=0;
378 
379  // sign bit
380  if(sign_flag)
381  result+=power(2, spec.e+spec.f);
382 
383  if(NaN_flag)
384  {
385  result+=power(2, spec.f)*spec.max_exponent();
386  result+=1;
387  }
388  else if(infinity_flag)
389  {
390  result+=power(2, spec.f)*spec.max_exponent();
391  }
392  else if(fraction==0 && exponent==0)
393  {
394  // zero
395  }
396  else if(is_normal()) // normal?
397  {
398  // fraction -- need to hide hidden bit
399  result+=fraction-power(2, spec.f); // hidden bit
400 
401  // exponent -- bias!
402  result+=power(2, spec.f)*(exponent+spec.bias());
403  }
404  else // denormal
405  {
406  result+=fraction; // denormal -- no hidden bit
407  // the exponent is zero
408  }
409 
410  return result;
411 }
412 
414  mp_integer &_fraction,
415  mp_integer &_exponent) const
416 {
417  if(is_zero() || is_NaN() || is_infinity())
418  {
419  _fraction=_exponent=0;
420  return;
421  }
422 
423  _exponent=exponent;
424  _fraction=fraction;
425 
426  // adjust exponent
427  _exponent-=spec.f;
428 
429  // try to integer-ize
430  while((_fraction%2)==0)
431  {
432  _fraction/=2;
433  ++_exponent;
434  }
435 }
436 
438  mp_integer &_fraction,
439  mp_integer &_exponent) const
440 {
441  if(is_zero() || is_NaN() || is_infinity())
442  {
443  _fraction=_exponent=0;
444  return;
445  }
446 
447  _exponent=exponent;
448  _fraction=fraction;
449 
450  // adjust exponent
451  _exponent-=spec.f;
452 
453  // now make it base 10
454  if(_exponent>=0)
455  {
456  _fraction*=power(2, _exponent);
457  _exponent=0;
458  }
459  else // _exponent<0
460  {
461  // 10/2=5 -- this makes it base 10
462  _fraction*=power(5, -_exponent);
463  }
464 
465  // try to re-normalize
466  while((_fraction%10)==0)
467  {
468  _fraction/=10;
469  ++_exponent;
470  }
471 }
472 
474  const mp_integer &_fraction,
475  const mp_integer &_exponent)
476 {
477  sign_flag=_fraction<0;
478  fraction=_fraction;
479  if(sign_flag)
481  exponent=_exponent;
482  exponent+=spec.f;
483  align();
484 }
485 
488  const mp_integer &_fraction,
489  const mp_integer &_exponent)
490 {
491  NaN_flag=infinity_flag=false;
492  sign_flag=_fraction<0;
493  fraction=_fraction;
494  if(sign_flag)
496  exponent=spec.f;
497  exponent+=_exponent;
498 
499  if(_exponent<0)
500  {
501  // bring to max. precision
502  mp_integer e_power=power(2, spec.e);
503  fraction*=power(2, e_power);
504  exponent-=e_power;
505  fraction/=power(5, -_exponent);
506  }
507  else if(_exponent>0)
508  {
509  // fix base
510  fraction*=power(5, _exponent);
511  }
512 
513  align();
514 }
515 
517 {
519  exponent=spec.f;
520  fraction=i;
521  align();
522 }
523 
525 {
526  // NaN?
527  if(NaN_flag)
528  {
529  fraction=0;
530  exponent=0;
531  sign_flag=false;
532  return;
533  }
534 
535  // do sign
536  if(fraction<0)
537  {
540  }
541 
542  // zero?
543  if(fraction==0)
544  {
545  exponent=0;
546  return;
547  }
548 
549  // 'usual case'
550  mp_integer f_power=power(2, spec.f);
551  mp_integer f_power_next=power(2, spec.f+1);
552 
553  std::size_t lowPower2=fraction.floorPow2();
554  mp_integer exponent_offset=0;
555 
556  if(lowPower2<spec.f) // too small
557  {
558  exponent_offset-=(spec.f-lowPower2);
559 
560  INVARIANT(
561  fraction * power(2, (spec.f - lowPower2)) >= f_power,
562  "normalisation result must be >= lower bound");
563  INVARIANT(
564  fraction * power(2, (spec.f - lowPower2)) < f_power_next,
565  "normalisation result must be < upper bound");
566  }
567  else if(lowPower2>spec.f) // too large
568  {
569  exponent_offset+=(lowPower2-spec.f);
570 
571  INVARIANT(
572  fraction / power(2, (lowPower2 - spec.f)) >= f_power,
573  "normalisation result must be >= lower bound");
574  INVARIANT(
575  fraction / power(2, (lowPower2 - spec.f)) < f_power_next,
576  "normalisation result must be < upper bound");
577  }
578 
579  mp_integer biased_exponent=exponent+exponent_offset+spec.bias();
580 
581  // exponent too large (infinity)?
582  if(biased_exponent>=spec.max_exponent())
583  {
584  // we need to consider the rounding mode here
585  switch(rounding_mode)
586  {
587  case UNKNOWN:
588  case NONDETERMINISTIC:
589  case ROUND_TO_EVEN:
590  infinity_flag=true;
591  break;
592 
593  case ROUND_TO_MINUS_INF:
594  // the result of the rounding is never larger than the argument
595  if(sign_flag)
596  infinity_flag=true;
597  else
598  make_fltmax();
599  break;
600 
601  case ROUND_TO_PLUS_INF:
602  // the result of the rounding is never smaller than the argument
603  if(sign_flag)
604  {
605  make_fltmax();
606  sign_flag=true; // restore sign
607  }
608  else
609  infinity_flag=true;
610  break;
611 
612  case ROUND_TO_ZERO:
613  if(sign_flag)
614  {
615  make_fltmax();
616  sign_flag=true; // restore sign
617  }
618  else
619  make_fltmax(); // positive
620  break;
621  }
622 
623  return; // done
624  }
625  else if(biased_exponent<=0) // exponent too small?
626  {
627  // produce a denormal (or zero)
628  mp_integer new_exponent=mp_integer(1)-spec.bias();
629  exponent_offset=new_exponent-exponent;
630  }
631 
632  exponent+=exponent_offset;
633 
634  if(exponent_offset>0)
635  {
636  divide_and_round(fraction, power(2, exponent_offset));
637 
638  // rounding might make the fraction too big!
639  if(fraction==f_power_next)
640  {
641  fraction=f_power;
642  ++exponent;
643  }
644  }
645  else if(exponent_offset<0)
646  fraction*=power(2, -exponent_offset);
647 
648  if(fraction==0)
649  exponent=0;
650 }
651 
653  mp_integer &dividend,
654  const mp_integer &divisor)
655 {
656  const mp_integer remainder = dividend % divisor;
657  dividend /= divisor;
658 
659  if(remainder!=0)
660  {
661  switch(rounding_mode)
662  {
663  case ROUND_TO_EVEN:
664  {
665  mp_integer divisor_middle = divisor / 2;
666  if(remainder < divisor_middle)
667  {
668  // crop
669  }
670  else if(remainder > divisor_middle)
671  {
672  ++dividend;
673  }
674  else // exactly in the middle -- go to even
675  {
676  if((dividend % 2) != 0)
677  ++dividend;
678  }
679  }
680  break;
681 
682  case ROUND_TO_ZERO:
683  // this means just crop
684  break;
685 
686  case ROUND_TO_MINUS_INF:
687  if(sign_flag)
688  ++dividend;
689  break;
690 
691  case ROUND_TO_PLUS_INF:
692  if(!sign_flag)
693  ++dividend;
694  break;
695 
696  case NONDETERMINISTIC:
697  case UNKNOWN:
698  UNREACHABLE;
699  }
700  }
701 }
702 
704 {
706 }
707 
709 {
710  PRECONDITION(other.spec.f == spec.f);
711 
712  // NaN/x = NaN
713  if(NaN_flag)
714  return *this;
715 
716  // x/NaN = NaN
717  if(other.NaN_flag)
718  {
719  make_NaN();
720  return *this;
721  }
722 
723  // 0/0 = NaN
724  if(is_zero() && other.is_zero())
725  {
726  make_NaN();
727  return *this;
728  }
729 
730  // x/0 = +-inf
731  if(other.is_zero())
732  {
733  infinity_flag=true;
734  if(other.sign_flag)
735  negate();
736  return *this;
737  }
738 
739  // x/inf = NaN
740  if(other.infinity_flag)
741  {
742  if(infinity_flag)
743  {
744  make_NaN();
745  return *this;
746  }
747 
748  bool old_sign=sign_flag;
749  make_zero();
750  sign_flag=old_sign;
751 
752  if(other.sign_flag)
753  negate();
754 
755  return *this;
756  } // inf/x = inf
757  else if(infinity_flag)
758  {
759  if(other.sign_flag)
760  negate();
761 
762  return *this;
763  }
764 
765  exponent-=other.exponent;
766  fraction*=power(2, spec.f);
767 
768  // to account for error
769  fraction*=power(2, spec.f);
770  exponent-=spec.f;
771 
772  fraction/=other.fraction;
773 
774  if(other.sign_flag)
775  negate();
776 
777  align();
778 
779  return *this;
780 }
781 
783 {
784  PRECONDITION(other.spec.f == spec.f);
785 
786  if(other.NaN_flag)
787  make_NaN();
788  if(NaN_flag)
789  return *this;
790 
791  if(infinity_flag || other.infinity_flag)
792  {
793  if(is_zero() || other.is_zero())
794  {
795  // special case Inf * 0 is NaN
796  make_NaN();
797  return *this;
798  }
799 
800  if(other.sign_flag)
801  negate();
802  infinity_flag=true;
803  return *this;
804  }
805 
806  exponent+=other.exponent;
807  exponent-=spec.f;
808  fraction*=other.fraction;
809 
810  if(other.sign_flag)
811  negate();
812 
813  align();
814 
815  return *this;
816 }
817 
819 {
820  PRECONDITION(other.spec == spec);
821  ieee_floatt _other=other;
822 
823  if(other.NaN_flag)
824  make_NaN();
825  if(NaN_flag)
826  return *this;
827 
828  if(infinity_flag && other.infinity_flag)
829  {
830  if(sign_flag==other.sign_flag)
831  return *this;
832  make_NaN();
833  return *this;
834  }
835  else if(infinity_flag)
836  return *this;
837  else if(other.infinity_flag)
838  {
839  infinity_flag=true;
840  sign_flag=other.sign_flag;
841  return *this;
842  }
843 
844  // 0 + 0 needs special treatment for the signs
845  if(is_zero() && other.is_zero())
846  {
847  if(get_sign()==other.get_sign())
848  return *this;
849  else
850  {
852  {
853  set_sign(true);
854  return *this;
855  }
856  else
857  {
858  set_sign(false);
859  return *this;
860  }
861  }
862  }
863 
864  // get smaller exponent
865  if(_other.exponent<exponent)
866  {
867  fraction*=power(2, exponent-_other.exponent);
868  exponent=_other.exponent;
869  }
870  else if(exponent<_other.exponent)
871  {
872  _other.fraction*=power(2, _other.exponent-exponent);
873  _other.exponent=exponent;
874  }
875 
876  INVARIANT(
877  exponent == _other.exponent,
878  "prior block equalises the exponents by setting both to the "
879  "minimum of their previous values while adjusting the mantissa");
880 
881  if(sign_flag)
882  fraction.negate();
883  if(_other.sign_flag)
884  _other.fraction.negate();
885 
886  fraction+=_other.fraction;
887 
888  // if the result is zero,
889  // there is some set of rules to get the sign
890  if(fraction==0)
891  {
893  sign_flag=true;
894  else
895  sign_flag=false;
896  }
897  else // fraction!=0
898  {
899  sign_flag=(fraction<0);
900  if(sign_flag)
901  fraction.negate();
902  }
903 
904  align();
905 
906  return *this;
907 }
908 
910 {
911  ieee_floatt _other=other;
912  _other.sign_flag=!_other.sign_flag;
913  return (*this)+=_other;
914 }
915 
916 bool ieee_floatt::operator<(const ieee_floatt &other) const
917 {
918  if(NaN_flag || other.NaN_flag)
919  return false;
920 
921  // check both zero?
922  if(is_zero() && other.is_zero())
923  return false;
924 
925  // one of them zero?
926  if(is_zero())
927  return !other.sign_flag;
928  else if(other.is_zero())
929  return sign_flag;
930 
931  // check sign
932  if(sign_flag!=other.sign_flag)
933  return sign_flag;
934 
935  // handle infinity
936  if(infinity_flag)
937  {
938  if(other.infinity_flag)
939  return false;
940  else
941  return sign_flag;
942  }
943  else if(other.infinity_flag)
944  return !sign_flag;
945 
946  // check exponent
947  if(exponent!=other.exponent)
948  {
949  if(sign_flag) // both negative
950  return exponent>other.exponent;
951  else
952  return exponent<other.exponent;
953  }
954 
955  // check significand
956  if(sign_flag) // both negative
957  return fraction>other.fraction;
958  else
959  return fraction<other.fraction;
960 }
961 
962 bool ieee_floatt::operator<=(const ieee_floatt &other) const
963 {
964  if(NaN_flag || other.NaN_flag)
965  return false;
966 
967  // check zero
968  if(is_zero() && other.is_zero())
969  return true;
970 
971  // handle infinity
972  if(infinity_flag && other.infinity_flag &&
973  sign_flag==other.sign_flag)
974  return true;
975 
976  if(!infinity_flag && !other.infinity_flag &&
977  sign_flag==other.sign_flag &&
978  exponent==other.exponent &&
979  fraction==other.fraction)
980  return true;
981 
982  return *this<other;
983 }
984 
985 bool ieee_floatt::operator>(const ieee_floatt &other) const
986 {
987  return other<*this;
988 }
989 
990 bool ieee_floatt::operator>=(const ieee_floatt &other) const
991 {
992  return other<=*this;
993 }
994 
995 bool ieee_floatt::operator==(const ieee_floatt &other) const
996 {
997  // packed equality!
998  if(NaN_flag && other.NaN_flag)
999  return true;
1000  else if(NaN_flag || other.NaN_flag)
1001  return false;
1002 
1003  if(infinity_flag && other.infinity_flag &&
1004  sign_flag==other.sign_flag)
1005  return true;
1006  else if(infinity_flag || other.infinity_flag)
1007  return false;
1008 
1009  // if(a.is_zero() && b.is_zero()) return true;
1010 
1011  return
1012  exponent==other.exponent &&
1013  fraction==other.fraction &&
1014  sign_flag==other.sign_flag;
1015 }
1016 
1017 bool ieee_floatt::ieee_equal(const ieee_floatt &other) const
1018 {
1019  if(NaN_flag || other.NaN_flag)
1020  return false;
1021  if(is_zero() && other.is_zero())
1022  return true;
1023  PRECONDITION(spec == other.spec);
1024  return *this==other;
1025 }
1026 
1027 bool ieee_floatt::operator==(int i) const
1028 {
1029  ieee_floatt other(spec);
1030  other.from_integer(i);
1031  return *this==other;
1032 }
1033 
1034 bool ieee_floatt::operator!=(const ieee_floatt &other) const
1035 {
1036  return !(*this==other);
1037 }
1038 
1040 {
1041  if(NaN_flag || other.NaN_flag)
1042  return true; // !!!
1043  if(is_zero() && other.is_zero())
1044  return false;
1045  PRECONDITION(spec == other.spec);
1046  return *this!=other;
1047 }
1048 
1050 {
1051  mp_integer _exponent=exponent-spec.f;
1052  mp_integer _fraction=fraction;
1053 
1054  if(sign_flag)
1055  _fraction.negate();
1056 
1057  spec=dest_spec;
1058 
1059  if(_fraction==0)
1060  {
1061  // We have a zero. It stays a zero.
1062  // Don't call build to preserve sign.
1063  }
1064  else
1065  build(_fraction, _exponent);
1066 }
1067 
1069 {
1071  unpack(bvrep2integer(expr.get_value(), spec.width(), false));
1072 }
1073 
1075 {
1076  if(NaN_flag || infinity_flag || is_zero())
1077  return 0;
1078 
1079  mp_integer result=fraction;
1080 
1081  mp_integer new_exponent=exponent-spec.f;
1082 
1083  // if the exponent is negative, divide
1084  if(new_exponent<0)
1085  result/=power(2, -new_exponent);
1086  else
1087  result*=power(2, new_exponent); // otherwise, multiply
1088 
1089  if(sign_flag)
1090  result.negate();
1091 
1092  return result;
1093 }
1094 
1095 void ieee_floatt::from_double(const double d)
1096 {
1097  spec.f=52;
1098  spec.e=11;
1099  DATA_INVARIANT(spec.width() == 64, "width must be 64 bits");
1100 
1101  static_assert(
1102  std::numeric_limits<double>::is_iec559,
1103  "this requires the double layout is according to the ieee754 "
1104  "standard");
1105  static_assert(
1106  sizeof(double) == sizeof(std::uint64_t), "ensure double has 64 bit width");
1107 
1108  union
1109  {
1110  double d;
1111  std::uint64_t i;
1112  } u;
1113 
1114  u.d=d;
1115 
1116  unpack(u.i);
1117 }
1118 
1119 void ieee_floatt::from_float(const float f)
1120 {
1121  spec.f=23;
1122  spec.e=8;
1123  DATA_INVARIANT(spec.width() == 32, "width must be 32 bits");
1124 
1125  static_assert(
1126  std::numeric_limits<float>::is_iec559,
1127  "this requires the float layout is according to the ieee754 "
1128  "standard");
1129  static_assert(
1130  sizeof(float) == sizeof(std::uint32_t), "ensure float has 32 bit width");
1131 
1132  union
1133  {
1134  float f;
1135  std::uint32_t i;
1136  } u;
1137 
1138  u.f=f;
1139 
1140  unpack(u.i);
1141 }
1142 
1144 {
1145  NaN_flag=true;
1146  // sign=false;
1147  exponent=0;
1148  fraction=0;
1149  infinity_flag=false;
1150 }
1151 
1153 {
1154  mp_integer bit_pattern=
1155  power(2, spec.e+spec.f)-1-power(2, spec.f);
1156  unpack(bit_pattern);
1157 }
1158 
1160 {
1161  unpack(power(2, spec.f));
1162 }
1163 
1165 {
1166  NaN_flag=false;
1167  sign_flag=false;
1168  exponent=0;
1169  fraction=0;
1170  infinity_flag=true;
1171 }
1172 
1174 {
1176  sign_flag=true;
1177 }
1178 
1180 {
1181  return spec.f==52 && spec.e==11;
1182 }
1183 
1185 {
1186  return spec.f==23 && spec.e==8;
1187 }
1188 
1192 {
1194  union { double f; uint64_t i; } a;
1195 
1196  if(infinity_flag)
1197  {
1198  if(sign_flag)
1199  return -std::numeric_limits<double>::infinity();
1200  else
1201  return std::numeric_limits<double>::infinity();
1202  }
1203 
1204  if(NaN_flag)
1205  {
1206  if(sign_flag)
1207  return -std::numeric_limits<double>::quiet_NaN();
1208  else
1209  return std::numeric_limits<double>::quiet_NaN();
1210  }
1211 
1212  mp_integer i=pack();
1213  CHECK_RETURN(i.is_ulong());
1214  CHECK_RETURN(i <= std::numeric_limits<std::uint64_t>::max());
1215 
1216  a.i=i.to_ulong();
1217  return a.f;
1218 }
1219 
1223 {
1224  // if false - ieee_floatt::to_float not supported on this architecture
1225  static_assert(
1226  std::numeric_limits<float>::is_iec559,
1227  "this requires the float layout is according to the IEC 559/IEEE 754 "
1228  "standard");
1229  static_assert(
1230  sizeof(float) == sizeof(uint32_t), "a 32 bit float type is required");
1231 
1232  union { float f; uint32_t i; } a;
1233 
1234  if(infinity_flag)
1235  {
1236  if(sign_flag)
1237  return -std::numeric_limits<float>::infinity();
1238  else
1239  return std::numeric_limits<float>::infinity();
1240  }
1241 
1242  if(NaN_flag)
1243  {
1244  if(sign_flag)
1245  return -std::numeric_limits<float>::quiet_NaN();
1246  else
1247  return std::numeric_limits<float>::quiet_NaN();
1248  }
1249 
1250  a.i = numeric_cast_v<uint32_t>(pack());
1251  return a.f;
1252 }
1253 
1257 {
1258  if(is_NaN())
1259  return;
1260 
1261  bool old_sign=get_sign();
1262 
1263  if(is_zero())
1264  {
1265  unpack(1);
1266  set_sign(!greater);
1267 
1268  return;
1269  }
1270 
1271  if(is_infinity())
1272  {
1273  if(get_sign()==greater)
1274  {
1275  make_fltmax();
1276  set_sign(old_sign);
1277  }
1278  return;
1279  }
1280 
1281  bool dir;
1282  if(greater)
1283  dir=!get_sign();
1284  else
1285  dir=get_sign();
1286 
1287  set_sign(false);
1288 
1289  mp_integer old=pack();
1290  if(dir)
1291  ++old;
1292  else
1293  --old;
1294 
1295  unpack(old);
1296 
1297  // sign change impossible (zero case caught earlier)
1298  set_sign(old_sign);
1299 }
mp_integer bvrep2integer(const irep_idt &src, std::size_t width, bool is_signed)
convert a bit-vector representation (possibly signed) to integer
irep_idt integer2bvrep(const mp_integer &src, std::size_t width)
convert an integer to bit-vector representation with given width This uses two's complement for negat...
mp_integer power(const mp_integer &base, const mp_integer &exponent)
A multi-precision implementation of the power operator.
Pre-defined bitvector types.
const floatbv_typet & to_floatbv_type(const typet &type)
Cast a typet to a floatbv_typet.
void set_width(std::size_t width)
Definition: std_types.h:932
std::size_t get_width() const
Definition: std_types.h:925
A constant literal expression.
Definition: std_expr.h:2990
const irep_idt & get_value() const
Definition: std_expr.h:2998
typet & type()
Return the type of the expression.
Definition: expr.h:84
Fixed-width bit-vector with IEEE floating-point interpretation.
void set_f(std::size_t b)
std::size_t get_f() const
stylet style
Definition: format_spec.h:28
unsigned precision
Definition: format_spec.h:19
unsigned min_width
Definition: format_spec.h:18
mp_integer max_fraction() const
Definition: ieee_float.cpp:40
class floatbv_typet to_type() const
Definition: ieee_float.cpp:25
mp_integer bias() const
Definition: ieee_float.cpp:20
mp_integer max_exponent() const
Definition: ieee_float.cpp:35
void from_type(const floatbv_typet &type)
Definition: ieee_float.cpp:45
static ieee_float_spect double_precision()
Definition: ieee_float.h:76
std::size_t f
Definition: ieee_float.h:26
std::size_t width() const
Definition: ieee_float.h:50
std::size_t e
Definition: ieee_float.h:26
bool is_double() const
ieee_floatt & operator*=(const ieee_floatt &other)
Definition: ieee_float.cpp:782
void extract_base10(mp_integer &_exponent, mp_integer &_fraction) const
Definition: ieee_float.cpp:437
void make_minus_infinity()
void make_fltmax()
void extract_base2(mp_integer &_exponent, mp_integer &_fraction) const
Definition: ieee_float.cpp:413
float to_float() const
Note that calling from_float -> to_float can return different bit patterns for NaN.
void unpack(const mp_integer &i)
Definition: ieee_float.cpp:320
ieee_float_spect spec
Definition: ieee_float.h:134
mp_integer exponent
Definition: ieee_float.h:321
mp_integer to_integer() const
std::string to_ansi_c_string() const
Definition: ieee_float.h:280
bool is_NaN() const
Definition: ieee_float.h:248
void make_plus_infinity()
bool operator!=(const ieee_floatt &other) const
bool ieee_equal(const ieee_floatt &other) const
ieee_floatt & operator+=(const ieee_floatt &other)
Definition: ieee_float.cpp:818
constant_exprt to_expr() const
Definition: ieee_float.cpp:703
void divide_and_round(mp_integer &dividend, const mp_integer &divisor)
Definition: ieee_float.cpp:652
bool is_float() const
bool get_sign() const
Definition: ieee_float.h:247
void make_zero()
Definition: ieee_float.h:186
void set_sign(bool _sign)
Definition: ieee_float.h:183
void from_float(const float f)
bool NaN_flag
Definition: ieee_float.h:323
void from_expr(const constant_exprt &expr)
bool operator==(const ieee_floatt &other) const
Definition: ieee_float.cpp:995
static constant_exprt rounding_mode_expr(rounding_modet)
Definition: ieee_float.cpp:60
bool operator<=(const ieee_floatt &other) const
Definition: ieee_float.cpp:962
std::string to_string_scientific(std::size_t precision) const
format as [-]d.ddde+-d Note that printf always produces at least two digits for the exponent.
Definition: ieee_float.cpp:233
std::string to_string_decimal(std::size_t precision) const
Definition: ieee_float.cpp:139
void negate()
Definition: ieee_float.h:178
void make_NaN()
bool sign_flag
Definition: ieee_float.h:320
ieee_floatt & operator-=(const ieee_floatt &other)
Definition: ieee_float.cpp:909
bool operator>(const ieee_floatt &other) const
Definition: ieee_float.cpp:985
void next_representable(bool greater)
Sets *this to the next representable number closer to plus infinity (greater = true) or minus infinit...
bool ieee_not_equal(const ieee_floatt &other) const
bool infinity_flag
Definition: ieee_float.h:323
bool is_zero() const
Definition: ieee_float.h:243
double to_double() const
Note that calling from_double -> to_double can return different bit patterns for NaN.
ieee_floatt & operator/=(const ieee_floatt &other)
Definition: ieee_float.cpp:708
bool operator>=(const ieee_floatt &other) const
Definition: ieee_float.cpp:990
std::string format(const format_spect &format_spec) const
Definition: ieee_float.cpp:70
bool is_infinity() const
Definition: ieee_float.h:249
void make_fltmin()
void from_integer(const mp_integer &i)
Definition: ieee_float.cpp:516
rounding_modet rounding_mode
Definition: ieee_float.h:132
static mp_integer base10_digits(const mp_integer &src)
Definition: ieee_float.cpp:130
bool is_normal() const
Definition: ieee_float.cpp:370
mp_integer pack() const
Definition: ieee_float.cpp:375
mp_integer fraction
Definition: ieee_float.h:322
@ NONDETERMINISTIC
Definition: ieee_float.h:126
@ ROUND_TO_PLUS_INF
Definition: ieee_float.h:125
@ ROUND_TO_MINUS_INF
Definition: ieee_float.h:124
void build(const mp_integer &exp, const mp_integer &frac)
Definition: ieee_float.cpp:473
void from_base10(const mp_integer &exp, const mp_integer &frac)
compute f * (10^e)
Definition: ieee_float.cpp:487
void from_double(const double d)
bool operator<(const ieee_floatt &other) const
Definition: ieee_float.cpp:916
void change_spec(const ieee_float_spect &dest_spec)
void align()
Definition: ieee_float.cpp:524
void print(std::ostream &out) const
Definition: ieee_float.cpp:65
bool get_bool(const irep_idt &name) const
Definition: irep.cpp:57
void set(const irep_idt &name, const irep_idt &value)
Definition: irep.h:412
void dot(const goto_modelt &src, std::ostream &out)
Definition: dot.cpp:359
constant_exprt floatbv_rounding_mode(unsigned rm)
returns the a rounding mode expression for a given IEEE rounding mode, encoded using the recommendati...
API to expression classes for floating-point arithmetic.
static int8_t r
Definition: irep_hash.h:60
double remainder(double x, double y)
Definition: math.c:2396
const std::string integer2string(const mp_integer &n, unsigned base)
Definition: mp_arith.cpp:103
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 DATA_INVARIANT(CONDITION, REASON)
This condition should be used to document that assumptions that are made on goto_functions,...
Definition: invariant.h:534
#define PRECONDITION(CONDITION)
Definition: invariant.h:463
API to expression classes.