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, false);
 186       else
 187         return signedless_minus (*a_p, tmp, true);
 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 < *b_p)
 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       return r;
 215     }
 216
 217   if (dexp == 0)
 218     bb = b_p;
 219   else
 220     {
 221       tmp = *b_p;
 222       tmp.shift_right (dexp);
 223       bb = &tmp;
 224     }
 225
 226   r.m_sig = a_p->m_sig + bb->m_sig;
 227   r.normalize ();
 228
 229   r.m_negative = negative;
 230   return r;
 231 }
 232
 233 /* Return *this - other.  */
 234
 235 sreal
 236 sreal::operator- (const sreal &other) const
 237 {
 238   /* -a - b => -a + (-b).  */
 239   if (m_negative && !other.m_negative)
 240     return signedless_plus (*this, -other, true);
 241
 242   /* a - (-b) => a + b.  */
 243   if (!m_negative && other.m_negative)
 244     return signedless_plus (*this, -other, false);
 245
 246   gcc_checking_assert (m_negative == other.m_negative);
 247
 248   /* We want to substract a smaller number from bigger
 249     for nonegative numbers.  */
 250   if (!m_negative && *this < other)
 251     return -signedless_minus (other, *this, true);
 252
 253   /* Example: -2 - (-3) => 3 - 2 */
 254   if (m_negative && *this > other)
 255     return signedless_minus (-other, -(*this), true);
 256
 257   sreal r = signedless_minus (*this, other, m_negative);
 258
 259   return r;
 260 }
 261
 262 sreal
 263 sreal::signedless_minus (const sreal &a, const sreal &b, bool negative)
 264 {
 265   int dexp;
 266   sreal tmp, r;
 267   const sreal *bb;
 268   const sreal *a_p = &a;
 269   const sreal *b_p = &b;
 270
 271   dexp = a_p->m_exp - b_p->m_exp;
 272
 273   r.m_exp = a_p->m_exp;
 274   if (dexp > SREAL_BITS)
 275     {
 276       r.m_sig = a_p->m_sig;
 277       return r;
 278     }
 279   if (dexp == 0)
 280     bb = b_p;
 281   else
 282     {
 283       tmp = *b_p;
 284       tmp.shift_right (dexp);
 285       bb = &tmp;
 286     }
 287
 288   r.m_sig = a_p->m_sig - bb->m_sig;
 289   r.normalize ();
 290
 291   r.m_negative = negative;
 292   return r;
 293 }
 294
 295 /* Return *this * other.  */
 296
 297 sreal
 298 sreal::operator* (const sreal &other) const
 299 {
 300   sreal r;
 301   if (m_sig < SREAL_MIN_SIG || other.m_sig < SREAL_MIN_SIG)
 302     {
 303       r.m_sig = 0;
 304       r.m_exp = -SREAL_MAX_EXP;
 305     }
 306   else
 307     {
 308       r.m_sig = m_sig * other.m_sig;
 309       r.m_exp = m_exp + other.m_exp;
 310       r.normalize ();
 311     }
 312
 313   r.m_negative = m_negative ^ other.m_negative;
 314   return r;
 315 }
 316
 317 /* Return *this / other.  */
 318
 319 sreal
 320 sreal::operator/ (const sreal &other) const
 321 {
 322   gcc_checking_assert (other.m_sig != 0);
 323   sreal r;
 324   r.m_sig = (m_sig << SREAL_PART_BITS) / other.m_sig;
 325   r.m_exp = m_exp - other.m_exp - SREAL_PART_BITS;
 326   r.m_negative = m_negative ^ other.m_negative;
 327   r.normalize ();
 328   return r;
 329 }