gcc/sreal.c

   1 /* Simple data type for real numbers for the GNU compiler.
   2    Copyright (C) 2002-2014 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 "coretypes.h"
  53 #include "sreal.h"
  54
  55 /* Print the content of struct sreal.  */
  56
  57 void
  58 sreal::dump (FILE *file) const
  59 {
  60   fprintf (file, "(%" PRIu64 " * 2^%d)", m_sig, m_exp);
  61 }
  62
  63 DEBUG_FUNCTION void
  64 debug (sreal &ref)
  65 {
  66   ref.dump (stderr);
  67 }
  68
  69 DEBUG_FUNCTION void
  70 debug (sreal *ptr)
  71 {
  72   if (ptr)
  73     debug (*ptr);
  74   else
  75     fprintf (stderr, "<nil>\n");
  76 }
  77
  78 /* Shift this right by S bits.  Needed: 0 < S <= SREAL_BITS.
  79    When the most significant bit shifted out is 1, add 1 to this (rounding).
  80    */
  81
  82 void
  83 sreal::shift_right (int s)
  84 {
  85   gcc_checking_assert (s > 0);
  86   gcc_checking_assert (s <= SREAL_BITS);
  87   /* Exponent should never be so large because shift_right is used only by
  88      sreal_add and sreal_sub ant thus the number cannot be shifted out from
  89      exponent range.  */
  90   gcc_checking_assert (m_exp + s <= SREAL_MAX_EXP);
  91
  92   m_exp += s;
  93
  94   m_sig += (uint64_t) 1 << (s - 1);
  95   m_sig >>= s;
  96 }
  97
  98 /* Normalize *this.  */
  99
 100 void
 101 sreal::normalize ()
 102 {
 103   if (m_sig == 0)
 104     {
 105       m_negative = 0;
 106       m_exp = -SREAL_MAX_EXP;
 107     }
 108   else if (m_sig < SREAL_MIN_SIG)
 109     {
 110       do
 111         {
 112           m_sig <<= 1;
 113           m_exp--;
 114         }
 115       while (m_sig < SREAL_MIN_SIG);
 116
 117       /* Check underflow.  */
 118       if (m_exp < -SREAL_MAX_EXP)
 119         {
 120           m_exp = -SREAL_MAX_EXP;
 121           m_sig = 0;
 122         }
 123     }
 124   else if (m_sig > SREAL_MAX_SIG)
 125     {
 126       int last_bit;
 127       do
 128         {
 129           last_bit = m_sig & 1;
 130           m_sig >>= 1;
 131           m_exp++;
 132         }
 133       while (m_sig > SREAL_MAX_SIG);
 134
 135       /* Round the number.  */
 136       m_sig += last_bit;
 137       if (m_sig > SREAL_MAX_SIG)
 138         {
 139           m_sig >>= 1;
 140           m_exp++;
 141         }
 142
 143       /* Check overflow.  */
 144       if (m_exp > SREAL_MAX_EXP)
 145         {
 146           m_exp = SREAL_MAX_EXP;
 147           m_sig = SREAL_MAX_SIG;
 148         }
 149     }
 150 }
 151
 152 /* Return integer value of *this.  */
 153
 154 int64_t
 155 sreal::to_int () const
 156 {
 157   int64_t sign = m_negative ? -1 : 1;
 158
 159   if (m_exp <= -SREAL_BITS)
 160     return 0;
 161   if (m_exp >= SREAL_PART_BITS)
 162     return sign * INTTYPE_MAXIMUM (int64_t);
 163   if (m_exp > 0)
 164     return sign * (m_sig << m_exp);
 165   if (m_exp < 0)
 166     return sign * (m_sig >> -m_exp);
 167   return sign * m_sig;
 168 }
 169
 170 /* Return *this + other.  */
 171
 172 sreal
 173 sreal::operator+ (const sreal &other) const
 174 {
 175   const sreal *a_p = this, *b_p = &other;
 176
 177   if (a_p->m_negative && !b_p->m_negative)
 178     std::swap (a_p, b_p);
 179
 180   /* a + -b => a - b.  */
 181   if (!a_p->m_negative && b_p->m_negative)
 182     {
 183       sreal tmp = -(*b_p);
 184       if (*a_p < tmp)
 185         return signedless_minus (tmp, *a_p, true);
 186       else
 187         return signedless_minus (*a_p, tmp, false);
 188     }
 189
 190   gcc_checking_assert (a_p->m_negative == b_p->m_negative);
 191
 192   sreal r = signedless_plus (*a_p, *b_p, a_p->m_negative);
 193
 194   return r;
 195 }
 196
 197 sreal
 198 sreal::signedless_plus (const sreal &a, const sreal &b, bool negative)
 199 {
 200   const sreal *bb;
 201   sreal r, tmp;
 202   int dexp;
 203   const sreal *a_p = &a;
 204   const sreal *b_p = &b;
 205
 206   if (a_p->m_exp < b_p->m_exp)
 207     std::swap (a_p, b_p);
 208
 209   dexp = a_p->m_exp - b_p->m_exp;
 210   r.m_exp = a_p->m_exp;
 211   if (dexp > SREAL_BITS)
 212     {
 213       r.m_sig = a_p->m_sig;
 214       r.m_negative = negative;
 215       return r;
 216     }
 217
 218   if (dexp == 0)
 219     bb = b_p;
 220   else
 221     {
 222       tmp = *b_p;
 223       tmp.shift_right (dexp);
 224       bb = &tmp;
 225     }
 226
 227   r.m_sig = a_p->m_sig + bb->m_sig;
 228   r.normalize ();
 229
 230   r.m_negative = negative;
 231   return r;
 232 }
 233
 234 /* Return *this - other.  */
 235
 236 sreal
 237 sreal::operator- (const sreal &other) const
 238 {
 239   /* -a - b => -a + (-b).  */
 240   if (m_negative && !other.m_negative)
 241     return signedless_plus (*this, -other, true);
 242
 243   /* a - (-b) => a + b.  */
 244   if (!m_negative && other.m_negative)
 245     return signedless_plus (*this, -other, false);
 246
 247   gcc_checking_assert (m_negative == other.m_negative);
 248
 249   /* We want to substract a smaller number from bigger
 250     for nonegative numbers.  */
 251   if (!m_negative && *this < other)
 252     return signedless_minus (other, *this, true);
 253
 254   /* Example: -2 - (-3) => 3 - 2 */
 255   if (m_negative && *this > other)
 256     return signedless_minus (-other, -(*this), false);
 257
 258   sreal r = signedless_minus (*this, other, m_negative);
 259
 260   return r;
 261 }
 262
 263 sreal
 264 sreal::signedless_minus (const sreal &a, const sreal &b, bool negative)
 265 {
 266   int dexp;
 267   sreal tmp, r;
 268   const sreal *bb;
 269   const sreal *a_p = &a;
 270   const sreal *b_p = &b;
 271
 272   dexp = a_p->m_exp - b_p->m_exp;
 273
 274   r.m_exp = a_p->m_exp;
 275   if (dexp > SREAL_BITS)
 276     {
 277       r.m_sig = a_p->m_sig;
 278       r.m_negative = negative;
 279       return r;
 280     }
 281   if (dexp == 0)
 282     bb = b_p;
 283   else
 284     {
 285       tmp = *b_p;
 286       tmp.shift_right (dexp);
 287       bb = &tmp;
 288     }
 289
 290   r.m_sig = a_p->m_sig - bb->m_sig;
 291   r.normalize ();
 292
 293   r.m_negative = negative;
 294   return r;
 295 }
 296
 297 /* Return *this * other.  */
 298
 299 sreal
 300 sreal::operator* (const sreal &other) const
 301 {
 302   sreal r;
 303   if (m_sig < SREAL_MIN_SIG || other.m_sig < SREAL_MIN_SIG)
 304     {
 305       r.m_sig = 0;
 306       r.m_exp = -SREAL_MAX_EXP;
 307     }
 308   else
 309     {
 310       r.m_sig = m_sig * other.m_sig;
 311       r.m_exp = m_exp + other.m_exp;
 312       r.normalize ();
 313     }
 314
 315   r.m_negative = m_negative ^ other.m_negative;
 316   return r;
 317 }
 318
 319 /* Return *this / other.  */
 320
 321 sreal
 322 sreal::operator/ (const sreal &other) const
 323 {
 324   gcc_checking_assert (other.m_sig != 0);
 325   sreal r;
 326   r.m_sig = (m_sig << SREAL_PART_BITS) / other.m_sig;
 327   r.m_exp = m_exp - other.m_exp - SREAL_PART_BITS;
 328   r.m_negative = m_negative ^ other.m_negative;
 329   r.normalize ();
 330   return r;
 331 }