1 /* Simple data type for real numbers for the GNU compiler. 2 Copyright (C) 2002-2015 Free Software Foundation, Inc. 3 4 This file is part of GCC. 5 6 GCC is free software; you can redistribute it and/or modify it under 7 the terms of the GNU General Public License as published by the Free 8 Software Foundation; either version 3, or (at your option) any later 9 version. 10 11 GCC is distributed in the hope that it will be useful, but WITHOUT ANY 12 WARRANTY; without even the implied warranty of MERCHANTABILITY or 13 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 14 for more details. 15 16 You should have received a copy of the GNU General Public License 17 along with GCC; see the file COPYING3. If not see 18 <http://www.gnu.org/licenses/>. */ 19 20 /* This library supports real numbers; 21 inf and nan are NOT supported. 22 It is written to be simple and fast. 23 24 Value of sreal is 25 x = sig * 2 ^ exp 26 where 27 sig = significant 28 (for < 64-bit machines sig = sig_lo + sig_hi * 2 ^ SREAL_PART_BITS) 29 exp = exponent 30 31 One uint64_t is used for the significant. 32 Only a half of significant bits is used (in normalized sreals) so that we do 33 not have problems with overflow, for example when c->sig = a->sig * b->sig. 34 So the precision is 32-bit. 35 36 Invariant: The numbers are normalized before and after each call of sreal_*. 37 38 Normalized sreals: 39 All numbers (except zero) meet following conditions: 40 SREAL_MIN_SIG <= sig && sig <= SREAL_MAX_SIG 41 -SREAL_MAX_EXP <= exp && exp <= SREAL_MAX_EXP 42 43 If the number would be too large, it is set to upper bounds of these 44 conditions. 45 46 If the number is zero or would be too small it meets following conditions: 47 sig == 0 && exp == -SREAL_MAX_EXP 48 */ 49 50 #include "config.h" 51 #include "system.h" 52 #include <math.h> 53 #include "coretypes.h" 54 #include "sreal.h" 55 56 /* Print the content of struct sreal. */ 57 58 void 59 sreal::dump (FILE *file) const 60 { 61 fprintf (file, "(%" PRIi64 " * 2^%d)", m_sig, m_exp); 62 } 63 64 DEBUG_FUNCTION void 65 debug (const sreal &ref) 66 { 67 ref.dump (stderr); 68 } 69 70 DEBUG_FUNCTION void 71 debug (const sreal *ptr) 72 { 73 if (ptr) 74 debug (*ptr); 75 else 76 fprintf (stderr, "<nil>\n"); 77 } 78 79 /* Shift this right by S bits. Needed: 0 < S <= SREAL_BITS. 80 When the most significant bit shifted out is 1, add 1 to this (rounding). 81 */ 82 83 void 84 sreal::shift_right (int s) 85 { 86 gcc_checking_assert (s > 0); 87 gcc_checking_assert (s <= SREAL_BITS); 88 /* Exponent should never be so large because shift_right is used only by 89 sreal_add and sreal_sub ant thus the number cannot be shifted out from 90 exponent range. */ 91 gcc_checking_assert (m_exp + s <= SREAL_MAX_EXP); 92 93 m_exp += s; 94 95 m_sig += (int64_t) 1 << (s - 1); 96 m_sig >>= s; 97 } 98 99 /* Return integer value of *this. */ 100 101 int64_t 102 sreal::to_int () const 103 { 104 int64_t sign = m_sig < 0 ? -1 : 1; 105 106 if (m_exp <= -SREAL_BITS) 107 return 0; 108 if (m_exp >= SREAL_PART_BITS) 109 return sign * INTTYPE_MAXIMUM (int64_t); 110 if (m_exp > 0) 111 return m_sig << m_exp; 112 if (m_exp < 0) 113 return m_sig >> -m_exp; 114 return m_sig; 115 } 116 117 /* Return value of *this as double. 118 This should be used for debug output only. */ 119 120 double 121 sreal::to_double () const 122 { 123 double val = m_sig; 124 if (m_exp) 125 val = ldexp (val, m_exp); 126 return val; 127 } 128 129 /* Return *this + other. */ 130 131 sreal 132 sreal::operator+ (const sreal &other) const 133 { 134 int dexp; 135 sreal tmp, r; 136 137 const sreal *a_p = this, *b_p = &other, *bb; 138 139 if (a_p->m_exp < b_p->m_exp) 140 std::swap (a_p, b_p); 141 142 dexp = a_p->m_exp - b_p->m_exp; 143 r.m_exp = a_p->m_exp; 144 if (dexp > SREAL_BITS) 145 { 146 r.m_sig = a_p->m_sig; 147 return r; 148 } 149 150 if (dexp == 0) 151 bb = b_p; 152 else 153 { 154 tmp = *b_p; 155 tmp.shift_right (dexp); 156 bb = &tmp; 157 } 158 159 r.m_sig = a_p->m_sig + bb->m_sig; 160 r.normalize (); 161 return r; 162 } 163 164 165 /* Return *this - other. */ 166 167 sreal 168 sreal::operator- (const sreal &other) const 169 { 170 int dexp; 171 sreal tmp, r; 172 const sreal *bb; 173 const sreal *a_p = this, *b_p = &other; 174 175 int64_t sign = 1; 176 if (a_p->m_exp < b_p->m_exp) 177 { 178 sign = -1; 179 std::swap (a_p, b_p); 180 } 181 182 dexp = a_p->m_exp - b_p->m_exp; 183 r.m_exp = a_p->m_exp; 184 if (dexp > SREAL_BITS) 185 { 186 r.m_sig = sign * a_p->m_sig; 187 return r; 188 } 189 if (dexp == 0) 190 bb = b_p; 191 else 192 { 193 tmp = *b_p; 194 tmp.shift_right (dexp); 195 bb = &tmp; 196 } 197 198 r.m_sig = sign * (a_p->m_sig - bb->m_sig); 199 r.normalize (); 200 return r; 201 } 202 203 /* Return *this * other. */ 204 205 sreal 206 sreal::operator* (const sreal &other) const 207 { 208 sreal r; 209 if (absu_hwi (m_sig) < SREAL_MIN_SIG || absu_hwi (other.m_sig) < SREAL_MIN_SIG) 210 { 211 r.m_sig = 0; 212 r.m_exp = -SREAL_MAX_EXP; 213 } 214 else 215 { 216 r.m_sig = m_sig * other.m_sig; 217 r.m_exp = m_exp + other.m_exp; 218 r.normalize (); 219 } 220 221 return r; 222 } 223 224 /* Return *this / other. */ 225 226 sreal 227 sreal::operator/ (const sreal &other) const 228 { 229 gcc_checking_assert (other.m_sig != 0); 230 sreal r; 231 r.m_sig = (m_sig << SREAL_PART_BITS) / other.m_sig; 232 r.m_exp = m_exp - other.m_exp - SREAL_PART_BITS; 233 r.normalize (); 234 return r; 235 } 236