1 /* IEC 559 floating point format related functions.
2    Copyright (C) 2001-2010 Roberto Bagnara <bagnara@cs.unipr.it>
3    Copyright (C) 2010-2016 BUGSENG srl (http://bugseng.com)
4 
5 This file is part of the Parma Polyhedra Library (PPL).
6 
7 The PPL is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 The PPL is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software Foundation,
19 Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA.
20 
21 For the most up-to-date information see the Parma Polyhedra Library
22 site: http://bugseng.com/products/ppl/ . */
23 
24 #ifndef PPL_Float_defs_hh
25 #define PPL_Float_defs_hh 1
26 
27 #include "globals_types.hh"
28 #include "meta_programming.hh"
29 #include "compiler.hh"
30 #include "Concrete_Expression_types.hh"
31 #include "Variable_types.hh"
32 #include "Linear_Form_types.hh"
33 #include <set>
34 #include <cmath>
35 #include <map>
36 #include <gmp.h>
37 
38 #ifdef NAN
39 #define PPL_NAN NAN
40 #else
41 #define PPL_NAN (HUGE_VAL - HUGE_VAL)
42 #endif
43 
44 namespace Parma_Polyhedra_Library {
45 
46 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
47 /*! \ingroup PPL_CXX_interface */
48 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
49 struct float_ieee754_half {
50   uint16_t word;
51   static const uint16_t SGN_MASK = 0x8000U;
52   static const uint16_t EXP_MASK = 0xfc00U;
53   static const uint16_t WRD_MAX = 0x7bffU;
54   static const uint16_t POS_INF = 0x7c00U;
55   static const uint16_t NEG_INF = 0xfc00U;
56   static const uint16_t POS_ZERO = 0x0000U;
57   static const uint16_t NEG_ZERO = 0x8000U;
58   static const unsigned int BASE = 2;
59   static const unsigned int EXPONENT_BITS = 5;
60   static const unsigned int MANTISSA_BITS = 10;
61   static const int EXPONENT_MAX = (1 << (EXPONENT_BITS - 1)) - 1;
62   static const int EXPONENT_BIAS = EXPONENT_MAX;
63   static const int EXPONENT_MIN = -EXPONENT_MAX + 1;
64   static const int EXPONENT_MIN_DENORM = EXPONENT_MIN
65                                         - static_cast<int>(MANTISSA_BITS);
66   static const Floating_Point_Format floating_point_format = IEEE754_HALF;
67   int inf_sign() const;
68   bool is_nan() const;
69   int zero_sign() const;
70   bool sign_bit() const;
71   void negate();
72   void dec();
73   void inc();
74   void set_max(bool negative);
75   void build(bool negative, mpz_t mantissa, int exponent);
76 
77 };
78 
79 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
80 /*! \ingroup PPL_CXX_interface */
81 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
82 struct float_ieee754_single {
83   uint32_t word;
84   static const uint32_t SGN_MASK = 0x80000000U;
85   static const uint32_t EXP_MASK = 0x7f800000U;
86   static const uint32_t WRD_MAX = 0x7f7fffffU;
87   static const uint32_t POS_INF = 0x7f800000U;
88   static const uint32_t NEG_INF = 0xff800000U;
89   static const uint32_t POS_ZERO = 0x00000000U;
90   static const uint32_t NEG_ZERO = 0x80000000U;
91   static const unsigned int BASE = 2;
92   static const unsigned int EXPONENT_BITS = 8;
93   static const unsigned int MANTISSA_BITS = 23;
94   static const int EXPONENT_MAX = (1 << (EXPONENT_BITS - 1)) - 1;
95   static const int EXPONENT_BIAS = EXPONENT_MAX;
96   static const int EXPONENT_MIN = -EXPONENT_MAX + 1;
97   static const int EXPONENT_MIN_DENORM = EXPONENT_MIN
98                                         - static_cast<int>(MANTISSA_BITS);
99   static const Floating_Point_Format floating_point_format = IEEE754_SINGLE;
100   int inf_sign() const;
101   bool is_nan() const;
102   int zero_sign() const;
103   bool sign_bit() const;
104   void negate();
105   void dec();
106   void inc();
107   void set_max(bool negative);
108   void build(bool negative, mpz_t mantissa, int exponent);
109 };
110 
111 #ifdef WORDS_BIGENDIAN
112 #ifndef PPL_WORDS_BIGENDIAN
113 #define PPL_WORDS_BIGENDIAN
114 #endif
115 #endif
116 
117 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
118 /*! \ingroup PPL_CXX_interface */
119 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
120 struct float_ieee754_double {
121 #ifdef PPL_WORDS_BIGENDIAN
122   uint32_t msp;
123   uint32_t lsp;
124 #else
125   uint32_t lsp;
126   uint32_t msp;
127 #endif
128   static const uint32_t MSP_SGN_MASK = 0x80000000U;
129   static const uint32_t MSP_POS_INF = 0x7ff00000U;
130   static const uint32_t MSP_NEG_INF = 0xfff00000U;
131   static const uint32_t MSP_POS_ZERO = 0x00000000U;
132   static const uint32_t MSP_NEG_ZERO = 0x80000000U;
133   static const uint32_t LSP_INF = 0;
134   static const uint32_t LSP_ZERO = 0;
135   static const uint32_t MSP_MAX = 0x7fefffffU;
136   static const uint32_t LSP_MAX = 0xffffffffU;
137   static const unsigned int BASE = 2;
138   static const unsigned int EXPONENT_BITS = 11;
139   static const unsigned int MANTISSA_BITS = 52;
140   static const int EXPONENT_MAX = (1 << (EXPONENT_BITS - 1)) - 1;
141   static const int EXPONENT_BIAS = EXPONENT_MAX;
142   static const int EXPONENT_MIN = -EXPONENT_MAX + 1;
143   static const int EXPONENT_MIN_DENORM = EXPONENT_MIN
144                                         - static_cast<int>(MANTISSA_BITS);
145   static const Floating_Point_Format floating_point_format = IEEE754_DOUBLE;
146   int inf_sign() const;
147   bool is_nan() const;
148   int zero_sign() const;
149   bool sign_bit() const;
150   void negate();
151   void dec();
152   void inc();
153   void set_max(bool negative);
154   void build(bool negative, mpz_t mantissa, int exponent);
155 };
156 
157 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
158 /*! \ingroup PPL_CXX_interface */
159 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
160 struct float_ibm_single {
161   uint32_t word;
162   static const uint32_t SGN_MASK = 0x80000000U;
163   static const uint32_t EXP_MASK = 0x7f000000U;
164   static const uint32_t WRD_MAX = 0x7fffffffU;
165   static const uint32_t POS_INF = 0x7f000000U;
166   static const uint32_t NEG_INF = 0xff000000U;
167   static const uint32_t POS_ZERO = 0x00000000U;
168   static const uint32_t NEG_ZERO = 0x80000000U;
169   static const unsigned int BASE = 16;
170   static const unsigned int EXPONENT_BITS = 7;
171   static const unsigned int MANTISSA_BITS = 24;
172   static const int EXPONENT_BIAS = 64;
173   static const int EXPONENT_MAX = (1 << (EXPONENT_BITS - 1)) - 1;
174   static const int EXPONENT_MIN = -EXPONENT_MAX + 1;
175   static const int EXPONENT_MIN_DENORM = EXPONENT_MIN
176                                         - static_cast<int>(MANTISSA_BITS);
177   static const Floating_Point_Format floating_point_format = IBM_SINGLE;
178   int inf_sign() const;
179   bool is_nan() const;
180   int zero_sign() const;
181   bool sign_bit() const;
182   void negate();
183   void dec();
184   void inc();
185   void set_max(bool negative);
186   void build(bool negative, mpz_t mantissa, int exponent);
187 };
188 
189 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
190 /*! \ingroup PPL_CXX_interface */
191 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
192 struct float_ibm_double {
193   static const unsigned int BASE = 16;
194   static const unsigned int EXPONENT_BITS = 7;
195   static const unsigned int MANTISSA_BITS = 56;
196   static const int EXPONENT_BIAS = 64;
197 };
198 
199 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
200 /*! \ingroup PPL_CXX_interface */
201 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
202 struct float_intel_double_extended {
203 #ifdef PPL_WORDS_BIGENDIAN
204   uint32_t msp;
205   uint64_t lsp;
206 #else
207   uint64_t lsp;
208   uint32_t msp;
209 #endif
210   static const uint32_t MSP_SGN_MASK = 0x00008000U;
211   static const uint32_t MSP_POS_INF = 0x00007fffU;
212   static const uint32_t MSP_NEG_INF = 0x0000ffffU;
213   static const uint32_t MSP_POS_ZERO = 0x00000000U;
214   static const uint32_t MSP_NEG_ZERO = 0x00008000U;
215   static const uint64_t LSP_INF = static_cast<uint64_t>(0x8000000000000000ULL);
216   static const uint64_t LSP_ZERO = 0;
217   static const uint32_t MSP_MAX = 0x00007ffeU;
218   static const uint64_t LSP_DMAX = static_cast<uint64_t>(0x7fffffffffffffffULL);
219   static const uint64_t LSP_NMAX = static_cast<uint64_t>(0xffffffffffffffffULL);
220   static const unsigned int BASE = 2;
221   static const unsigned int EXPONENT_BITS = 15;
222   static const unsigned int MANTISSA_BITS = 63;
223   static const int EXPONENT_MAX = (1 << (EXPONENT_BITS - 1)) - 1;
224   static const int EXPONENT_BIAS = EXPONENT_MAX;
225   static const int EXPONENT_MIN = -EXPONENT_MAX + 1;
226   static const int EXPONENT_MIN_DENORM = EXPONENT_MIN
227                                         - static_cast<int>(MANTISSA_BITS);
228   static const Floating_Point_Format floating_point_format =
229                                      INTEL_DOUBLE_EXTENDED;
230   int inf_sign() const;
231   bool is_nan() const;
232   int zero_sign() const;
233   bool sign_bit() const;
234   void negate();
235   void dec();
236   void inc();
237   void set_max(bool negative);
238   void build(bool negative, mpz_t mantissa, int exponent);
239 };
240 
241 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
242 /*! \ingroup PPL_CXX_interface */
243 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
244 struct float_ieee754_quad {
245 #ifdef PPL_WORDS_BIGENDIAN
246   uint64_t msp;
247   uint64_t lsp;
248 #else
249   uint64_t lsp;
250   uint64_t msp;
251 #endif
252   static const uint64_t MSP_SGN_MASK = static_cast<uint64_t>(0x8000000000000000ULL);
253   static const uint64_t MSP_POS_INF = static_cast<uint64_t>(0x7fff000000000000ULL);
254   static const uint64_t MSP_NEG_INF = static_cast<uint64_t>(0xffff000000000000ULL);
255   static const uint64_t MSP_POS_ZERO = static_cast<uint64_t>(0x0000000000000000ULL);
256   static const uint64_t MSP_NEG_ZERO = static_cast<uint64_t>(0x8000000000000000ULL);
257   static const uint64_t LSP_INF = 0;
258   static const uint64_t LSP_ZERO = 0;
259   static const uint64_t MSP_MAX = static_cast<uint64_t>(0x7ffeffffffffffffULL);
260   static const uint64_t LSP_MAX = static_cast<uint64_t>(0xffffffffffffffffULL);
261   static const unsigned int BASE = 2;
262   static const unsigned int EXPONENT_BITS = 15;
263   static const unsigned int MANTISSA_BITS = 112;
264   static const int EXPONENT_MAX = (1 << (EXPONENT_BITS - 1)) - 1;
265   static const int EXPONENT_BIAS = EXPONENT_MAX;
266   static const int EXPONENT_MIN = -EXPONENT_MAX + 1;
267   static const int EXPONENT_MIN_DENORM = EXPONENT_MIN
268                                         - static_cast<int>(MANTISSA_BITS);
269   static const Floating_Point_Format floating_point_format = IEEE754_QUAD;
270   int inf_sign() const;
271   bool is_nan() const;
272   int zero_sign() const;
273   bool sign_bit() const;
274   void negate();
275   void dec();
276   void inc();
277   void set_max(bool negative);
278   void build(bool negative, mpz_t mantissa, int exponent);
279 };
280 
281 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
282 /*! \ingroup PPL_CXX_interface */
283 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
284 template <typename T>
285 class Float : public False { };
286 
287 #if PPL_SUPPORTED_FLOAT
288 template <>
289 class Float<float> : public True {
290 public:
291 #if PPL_CXX_FLOAT_BINARY_FORMAT == PPL_FLOAT_IEEE754_HALF
292   typedef float_ieee754_half Binary;
293 #elif PPL_CXX_FLOAT_BINARY_FORMAT == PPL_FLOAT_IEEE754_SINGLE
294   typedef float_ieee754_single Binary;
295 #elif PPL_CXX_FLOAT_BINARY_FORMAT == PPL_FLOAT_IEEE754_DOUBLE
296   typedef float_ieee754_double Binary;
297 #elif PPL_CXX_FLOAT_BINARY_FORMAT == PPL_FLOAT_IBM_SINGLE
298   typedef float_ibm_single Binary;
299 #elif PPL_CXX_FLOAT_BINARY_FORMAT == PPL_FLOAT_IEEE754_QUAD
300   typedef float_ieee754_quad Binary;
301 #elif PPL_CXX_FLOAT_BINARY_FORMAT == PPL_FLOAT_INTEL_DOUBLE_EXTENDED
302   typedef float_intel_double_extended Binary;
303 #else
304 #error "Invalid value for PPL_CXX_FLOAT_BINARY_FORMAT"
305 #endif
306   union {
307     float number;
308     Binary binary;
309   } u;
310   Float();
311   Float(float v);
312   float value();
313 };
314 #endif
315 
316 #if PPL_SUPPORTED_DOUBLE
317 template <>
318 class Float<double> : public True {
319 public:
320 #if PPL_CXX_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IEEE754_HALF
321   typedef float_ieee754_half Binary;
322 #elif PPL_CXX_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IEEE754_SINGLE
323   typedef float_ieee754_single Binary;
324 #elif PPL_CXX_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IEEE754_DOUBLE
325   typedef float_ieee754_double Binary;
326 #elif PPL_CXX_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IBM_SINGLE
327   typedef float_ibm_single Binary;
328 #elif PPL_CXX_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IEEE754_QUAD
329   typedef float_ieee754_quad Binary;
330 #elif PPL_CXX_DOUBLE_BINARY_FORMAT == PPL_FLOAT_INTEL_DOUBLE_EXTENDED
331   typedef float_intel_double_extended Binary;
332 #else
333 #error "Invalid value for PPL_CXX_DOUBLE_BINARY_FORMAT"
334 #endif
335   union {
336     double number;
337     Binary binary;
338   } u;
339   Float();
340   Float(double v);
341   double value();
342 };
343 #endif
344 
345 #if PPL_SUPPORTED_LONG_DOUBLE
346 template <>
347 class Float<long double> : public True {
348 public:
349 #if PPL_CXX_LONG_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IEEE754_HALF
350   typedef float_ieee754_half Binary;
351 #elif PPL_CXX_LONG_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IEEE754_SINGLE
352   typedef float_ieee754_single Binary;
353 #elif PPL_CXX_LONG_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IEEE754_DOUBLE
354   typedef float_ieee754_double Binary;
355 #elif PPL_CXX_LONG_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IBM_SINGLE
356   typedef float_ibm_single Binary;
357 #elif PPL_CXX_LONG_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IEEE754_QUAD
358   typedef float_ieee754_quad Binary;
359 #elif PPL_CXX_LONG_DOUBLE_BINARY_FORMAT == PPL_FLOAT_INTEL_DOUBLE_EXTENDED
360   typedef float_intel_double_extended Binary;
361 #else
362 #error "Invalid value for PPL_CXX_LONG_DOUBLE_BINARY_FORMAT"
363 #endif
364   union {
365     long double number;
366     Binary binary;
367   } u;
368   Float();
369   Float(long double v);
370   long double value();
371 };
372 #endif
373 
374 // FIXME: is this the right place for this function?
375 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
376 /*! \brief
377   If \p v is nonzero, returns the position of the most significant bit
378   in \p a.
379 
380   The behavior is undefined if \p v is zero.
381 */
382 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
383 unsigned int msb_position(unsigned long long v);
384 
385 /*! \brief
386   An abstract class to be implemented by an external analyzer such
387   as ECLAIR in order to provide to the PPL the necessary information
388   for performing the analysis of floating point computations.
389 
390   \par Template type parameters
391 
392   - The class template parameter \p Target specifies the implementation
393   of Concrete_Expression to be used.
394   - The class template parameter \p FP_Interval_Type represents the type
395   of the intervals used in the abstract domain. The interval bounds
396   should have a floating point type.
397 */
398 template <typename Target, typename FP_Interval_Type>
399 class FP_Oracle {
400 public:
401   /*
402     FIXME: the const qualifiers on expressions may raise CLANG
403     compatibility issues. It may be necessary to omit them.
404   */
405 
406   /*! \brief
407     Asks the external analyzer for an interval that correctly
408     approximates the floating point entity referenced by \p dim.
409     Result is stored into \p result.
410 
411     \return <CODE>true</CODE> if the analyzer was able to find a correct
412     approximation, <CODE>false</CODE> otherwise.
413   */
414   virtual bool get_interval(dimension_type dim, FP_Interval_Type& result) const
415     = 0;
416 
417   /*! \brief
418     Asks the external analyzer for an interval that correctly
419     approximates the value of floating point constant \p expr.
420     Result is stored into \p result.
421 
422     \return <CODE>true</CODE> if the analyzer was able to find a correct
423     approximation, <CODE>false</CODE> otherwise.
424   */
425   virtual bool get_fp_constant_value(
426                const Floating_Point_Constant<Target>& expr,
427                      FP_Interval_Type& result) const = 0;
428 
429   /*! \brief
430     Asks the external analyzer for an interval that correctly approximates
431     the value of \p expr, which must be of integer type.
432     Result is stored into \p result.
433 
434     \return <CODE>true</CODE> if the analyzer was able to find a correct
435     approximation, <CODE>false</CODE> otherwise.
436   */
437   virtual bool get_integer_expr_value(const Concrete_Expression<Target>& expr,
438                                       FP_Interval_Type& result) const = 0;
439 
440   /*! \brief
441     Asks the external analyzer for the possible space dimensions that
442     are associated to the approximable reference \p expr.
443     Result is stored into \p result.
444 
445     \return <CODE>true</CODE> if the analyzer was able to return
446     the (possibly empty!) set, <CODE>false</CODE> otherwise.
447 
448     The resulting set MUST NOT contain <CODE>not_a_dimension()</CODE>.
449   */
450   virtual bool get_associated_dimensions(
451           const Approximable_Reference<Target>& expr,
452           std::set<dimension_type>& result) const = 0;
453 
454 };
455 
456 /* FIXME: some of the following  documentation should probably be
457    under PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS */
458 
459 /*! \brief \relates Float
460   Returns <CODE>true</CODE> if and only if there is some floating point
461   number that is representable by \p f2 but not by \p f1.
462 */
463 bool is_less_precise_than(Floating_Point_Format f1, Floating_Point_Format f2);
464 
465 /*! \brief \relates Float
466   Computes the absolute error of floating point computations.
467 
468   \par Template type parameters
469 
470   - The class template parameter \p FP_Interval_Type represents the type
471   of the intervals used in the abstract domain. The interval bounds
472   should have a floating point type.
473 
474   \param analyzed_format The floating point format used by the analyzed
475   program.
476 
477   \return The interval \f$[-\omega, \omega]\f$ where \f$\omega\f$ is the
478   smallest non-zero positive number in the less precise floating point
479   format between the analyzer format and the analyzed format.
480 */
481 template <typename FP_Interval_Type>
482 const FP_Interval_Type&
483 compute_absolute_error(Floating_Point_Format analyzed_format);
484 
485 /*! \brief \relates Linear_Form
486   Discards all linear forms containing variable \p var from the
487   linear form abstract store \p lf_store.
488 */
489 template <typename FP_Interval_Type>
490 void
491 discard_occurrences(std::map<dimension_type,
492                              Linear_Form<FP_Interval_Type> >& lf_store,
493                     Variable var);
494 
495 /*! \brief \relates Linear_Form
496   Assigns the linear form \p lf to \p var in the linear form abstract
497   store \p lf_store, then discards all occurrences of \p var from it.
498 */
499 template <typename FP_Interval_Type>
500 void
501 affine_form_image(std::map<dimension_type,
502                            Linear_Form<FP_Interval_Type> >& lf_store,
503                   Variable var,
504                   const Linear_Form<FP_Interval_Type>& lf);
505 
506 /*! \brief \relates Linear_Form
507   Discards from \p ls1 all linear forms but those that are associated
508   to the same variable in \p ls2.
509 */
510 template <typename FP_Interval_Type>
511 void
512 upper_bound_assign(std::map<dimension_type,
513                             Linear_Form<FP_Interval_Type> >& ls1,
514                    const std::map<dimension_type,
515                                   Linear_Form<FP_Interval_Type> >& ls2);
516 
517 } // namespace Parma_Polyhedra_Library
518 
519 #include "Float_inlines.hh"
520 #include "Float_templates.hh"
521 
522 #endif // !defined(PPL_Float_defs_hh)
523