1 /* IEC 559 floating point format related functions:
2    non-inline template functions.
3    Copyright (C) 2001-2010 Roberto Bagnara <bagnara@cs.unipr.it>
4    Copyright (C) 2010-2016 BUGSENG srl (http://bugseng.com)
5 
6 This file is part of the Parma Polyhedra Library (PPL).
7 
8 The PPL is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 The PPL is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with this program; if not, write to the Free Software Foundation,
20 Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA.
21 
22 For the most up-to-date information see the Parma Polyhedra Library
23 site: http://bugseng.com/products/ppl/ . */
24 
25 #ifndef PPL_Float_templates_hh
26 #define PPL_Float_templates_hh 1
27 
28 #include "Variable_defs.hh"
29 #include "Linear_Form_defs.hh"
30 #include <cmath>
31 
32 namespace Parma_Polyhedra_Library {
33 
34 template <typename FP_Interval_Type>
compute_absolute_error(const Floating_Point_Format analyzed_format)35 const FP_Interval_Type& compute_absolute_error(
36                         const Floating_Point_Format analyzed_format) {
37   typedef typename FP_Interval_Type::boundary_type analyzer_format;
38 
39   // FIXME: check if initializing caches with EMPTY is better.
40   static const FP_Interval_Type ZERO_INTERVAL = FP_Interval_Type(0);
41   // Cached results for each different analyzed format.
42   static FP_Interval_Type ieee754_half_result = ZERO_INTERVAL;
43   static FP_Interval_Type ieee754_single_result = ZERO_INTERVAL;
44   static FP_Interval_Type ieee754_double_result = ZERO_INTERVAL;
45   static FP_Interval_Type ibm_single_result = ZERO_INTERVAL;
46   static FP_Interval_Type ieee754_quad_result = ZERO_INTERVAL;
47   static FP_Interval_Type intel_double_extended_result = ZERO_INTERVAL;
48 
49   FP_Interval_Type* to_compute = NULL;
50   // Get the necessary information on the analyzed's format.
51   unsigned int f_base;
52   int f_exponent_bias;
53   unsigned int f_mantissa_bits;
54   switch (analyzed_format) {
55     case IEEE754_HALF:
56       if (ieee754_half_result != ZERO_INTERVAL) {
57         return ieee754_half_result;
58       }
59       to_compute = &ieee754_half_result;
60       f_base = float_ieee754_half::BASE;
61       f_exponent_bias = float_ieee754_half::EXPONENT_BIAS;
62       f_mantissa_bits = float_ieee754_half::MANTISSA_BITS;
63       break;
64     case IEEE754_SINGLE:
65       if (ieee754_single_result != ZERO_INTERVAL) {
66         return ieee754_single_result;
67       }
68 
69       to_compute = &ieee754_single_result;
70       f_base = float_ieee754_single::BASE;
71       f_exponent_bias = float_ieee754_single::EXPONENT_BIAS;
72       f_mantissa_bits = float_ieee754_single::MANTISSA_BITS;
73       break;
74     case IEEE754_DOUBLE:
75       if (ieee754_double_result != ZERO_INTERVAL) {
76         return ieee754_double_result;
77       }
78 
79       to_compute = &ieee754_double_result;
80       f_base = float_ieee754_double::BASE;
81       f_exponent_bias = float_ieee754_double::EXPONENT_BIAS;
82       f_mantissa_bits = float_ieee754_double::MANTISSA_BITS;
83       break;
84     case IBM_SINGLE:
85       if (ibm_single_result != ZERO_INTERVAL) {
86         return ibm_single_result;
87       }
88 
89       to_compute = &ibm_single_result;
90       f_base = float_ibm_single::BASE;
91       f_exponent_bias = float_ibm_single::EXPONENT_BIAS;
92       f_mantissa_bits = float_ibm_single::MANTISSA_BITS;
93       break;
94     case IEEE754_QUAD:
95       if (ieee754_quad_result != ZERO_INTERVAL) {
96         return ieee754_quad_result;
97       }
98 
99       to_compute = &ieee754_quad_result;
100       f_base = float_ieee754_quad::BASE;
101       f_exponent_bias = float_ieee754_quad::EXPONENT_BIAS;
102       f_mantissa_bits = float_ieee754_quad::MANTISSA_BITS;
103       break;
104     case INTEL_DOUBLE_EXTENDED:
105       if (intel_double_extended_result != ZERO_INTERVAL) {
106         return intel_double_extended_result;
107       }
108 
109       to_compute = &intel_double_extended_result;
110       f_base = float_intel_double_extended::BASE;
111       f_exponent_bias = float_intel_double_extended::EXPONENT_BIAS;
112       f_mantissa_bits = float_intel_double_extended::MANTISSA_BITS;
113       break;
114     default:
115       PPL_UNREACHABLE;
116       break;
117   }
118 
119   PPL_ASSERT(to_compute != NULL);
120 
121   // We assume that f_base is a power of 2.
122   analyzer_format omega;
123   int power = static_cast<int>(msb_position(f_base))
124     * ((1 - f_exponent_bias) - static_cast<int>(f_mantissa_bits));
125   omega = std::max(static_cast<analyzer_format>(ldexp(1.0, power)),
126                    std::numeric_limits<analyzer_format>::denorm_min());
127 
128   to_compute->build(i_constraint(GREATER_OR_EQUAL, -omega),
129                     i_constraint(LESS_OR_EQUAL, omega));
130   return *to_compute;
131 }
132 
133 template <typename FP_Interval_Type>
134 void
discard_occurrences(std::map<dimension_type,Linear_Form<FP_Interval_Type>> & lf_store,Variable var)135 discard_occurrences(std::map<dimension_type,
136                              Linear_Form<FP_Interval_Type> >& lf_store,
137                     Variable var) {
138   typedef Linear_Form<FP_Interval_Type> FP_Linear_Form;
139   typedef typename std::map<dimension_type, FP_Linear_Form>::iterator Iter;
140   for (Iter i = lf_store.begin(); i != lf_store.end(); ) {
141     if ((i->second).coefficient(var) != 0) {
142       i = lf_store.erase(i);
143     }
144     else {
145       ++i;
146     }
147   }
148 }
149 
150 /* FIXME: improve efficiency by adding the list of potentially conflicting
151    variables as an argument. */
152 template <typename FP_Interval_Type>
upper_bound_assign(std::map<dimension_type,Linear_Form<FP_Interval_Type>> & ls1,const std::map<dimension_type,Linear_Form<FP_Interval_Type>> & ls2)153 void upper_bound_assign(std::map<dimension_type,
154                                  Linear_Form<FP_Interval_Type> >& ls1,
155                         const std::map<dimension_type,
156                                        Linear_Form<FP_Interval_Type> >& ls2) {
157   typedef Linear_Form<FP_Interval_Type> FP_Linear_Form;
158   typedef typename std::map<dimension_type, FP_Linear_Form>::iterator Iter;
159   typedef typename std::map<dimension_type,
160                             FP_Linear_Form>::const_iterator Const_Iter;
161 
162   Const_Iter i2_end = ls2.end();
163   for (Iter i1 = ls1.begin(), i1_end = ls1.end(); i1 != i1_end; ) {
164     Const_Iter i2 = ls2.find(i1->first);
165     if ((i2 == i2_end) || (i1->second != i2->second)) {
166       i1 = ls1.erase(i1);
167     }
168     else {
169       ++i1;
170     }
171   }
172 }
173 
174 } // namespace Parma_Polyhedra_Library
175 
176 #endif // !defined(PPL_Float_templates_hh)
177