1/* Simple data type for real numbers for the GNU compiler. 2 Copyright (C) 2002-2015 Free Software Foundation, Inc. 3 4This file is part of GCC. 5 6GCC is free software; you can redistribute it and/or modify it under 7the terms of the GNU General Public License as published by the Free 8Software Foundation; either version 3, or (at your option) any later 9version. 10 11GCC is distributed in the hope that it will be useful, but WITHOUT ANY 12WARRANTY; without even the implied warranty of MERCHANTABILITY or 13FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 14for more details. 15 16You should have received a copy of the GNU General Public License 17along 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 58void 59sreal::dump (FILE *file) const 60{ 61 fprintf (file, "(%" PRIi64 " * 2^%d)", m_sig, m_exp); 62} 63 64DEBUG_FUNCTION void 65debug (const sreal &ref) 66{ 67 ref.dump (stderr); 68} 69 70DEBUG_FUNCTION void 71debug (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 83void 84sreal::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 101int64_t 102sreal::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 120double 121sreal::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 131sreal 132sreal::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 167sreal 168sreal::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 205sreal 206sreal::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 226sreal 227sreal::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