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 (const sreal &ref)
  65 {
  66   ref.dump (stderr);
  67 }
  68
  69 DEBUG_FUNCTION void
  70 debug (const 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 += (int64_t) 1 << (s - 1);
  95   m_sig >>= s;
  96 }
  97
  98 /* Normalize *this.  */
  99
 100 void
 101 sreal::normalize ()
 102 {
 103   int64_t s = m_sig < 0 ? -1 : 1;
 104   unsigned HOST_WIDE_INT sig = absu_hwi (m_sig);
 105
 106   if (sig == 0)
 107     {
 108       m_exp = -SREAL_MAX_EXP;
 109     }
 110   else if (sig < SREAL_MIN_SIG)
 111     {
 112       do
 113         {
 114           sig <<= 1;
 115           m_exp--;
 116         }
 117       while (sig < SREAL_MIN_SIG);
 118
 119       /* Check underflow.  */
 120       if (m_exp < -SREAL_MAX_EXP)
 121         {
 122           m_exp = -SREAL_MAX_EXP;
 123           sig = 0;
 124         }
 125     }
 126   else if (sig > SREAL_MAX_SIG)
 127     {
 128       int last_bit;
 129       do
 130         {
 131           last_bit = sig & 1;
 132           sig >>= 1;
 133           m_exp++;
 134         }
 135       while (sig > SREAL_MAX_SIG);
 136
 137       /* Round the number.  */
 138       sig += last_bit;
 139       if (sig > SREAL_MAX_SIG)
 140         {
 141           sig >>= 1;
 142           m_exp++;
 143         }
 144
 145       /* Check overflow.  */
 146       if (m_exp > SREAL_MAX_EXP)
 147         {
 148           m_exp = SREAL_MAX_EXP;
 149           sig = SREAL_MAX_SIG;
 150         }
 151     }
 152
 153   m_sig = s * sig;
 154 }
 155
 156 /* Return integer value of *this.  */
 157
 158 int64_t
 159 sreal::to_int () const
 160 {
 161   int64_t sign = m_sig < 0 ? -1 : 1;
 162
 163   if (m_exp <= -SREAL_BITS)
 164     return 0;
 165   if (m_exp >= SREAL_PART_BITS)
 166     return sign * INTTYPE_MAXIMUM (int64_t);
 167   if (m_exp > 0)
 168     return m_sig << m_exp;
 169   if (m_exp < 0)
 170     return m_sig >> -m_exp;
 171   return m_sig;
 172 }
 173
 174 /* Return *this + other.  */
 175
 176 sreal
 177 sreal::operator+ (const sreal &other) const
 178 {
 179   int dexp;
 180   sreal tmp, r;
 181
 182   const sreal *a_p = this, *b_p = &other, *bb;
 183
 184   if (a_p->m_exp < b_p->m_exp)
 185     std::swap (a_p, b_p);
 186
 187   dexp = a_p->m_exp - b_p->m_exp;
 188   r.m_exp = a_p->m_exp;
 189   if (dexp > SREAL_BITS)
 190     {
 191       r.m_sig = a_p->m_sig;
 192       return r;
 193     }
 194
 195   if (dexp == 0)
 196     bb = b_p;
 197   else
 198     {
 199       tmp = *b_p;
 200       tmp.shift_right (dexp);
 201       bb = &tmp;
 202     }
 203
 204   r.m_sig = a_p->m_sig + bb->m_sig;
 205   r.normalize ();
 206   return r;
 207 }
 208
 209
 210 /* Return *this - other.  */
 211
 212 sreal
 213 sreal::operator- (const sreal &other) const
 214 {
 215   int dexp;
 216   sreal tmp, r;
 217   const sreal *bb;
 218   const sreal *a_p = this, *b_p = &other;
 219
 220   int64_t sign = 1;
 221   if (a_p->m_exp < b_p->m_exp)
 222     {
 223       sign = -1;
 224       std::swap (a_p, b_p);
 225     }
 226
 227   dexp = a_p->m_exp - b_p->m_exp;
 228   r.m_exp = a_p->m_exp;
 229   if (dexp > SREAL_BITS)
 230     {
 231       r.m_sig = sign * a_p->m_sig;
 232       return r;
 233     }
 234   if (dexp == 0)
 235     bb = b_p;
 236   else
 237     {
 238       tmp = *b_p;
 239       tmp.shift_right (dexp);
 240       bb = &tmp;
 241     }
 242
 243   r.m_sig = sign * (a_p->m_sig - bb->m_sig);
 244   r.normalize ();
 245   return r;
 246 }
 247
 248 /* Return *this * other.  */
 249
 250 sreal
 251 sreal::operator* (const sreal &other) const
 252 {
 253   sreal r;
 254   if (std::abs (m_sig) < SREAL_MIN_SIG || std::abs (other.m_sig) < SREAL_MIN_SIG)
 255     {
 256       r.m_sig = 0;
 257       r.m_exp = -SREAL_MAX_EXP;
 258     }
 259   else
 260     {
 261       r.m_sig = m_sig * other.m_sig;
 262       r.m_exp = m_exp + other.m_exp;
 263       r.normalize ();
 264     }
 265
 266   return r;
 267 }
 268
 269 /* Return *this / other.  */
 270
 271 sreal
 272 sreal::operator/ (const sreal &other) const
 273 {
 274   gcc_checking_assert (other.m_sig != 0);
 275   sreal r;
 276   r.m_sig = (m_sig << SREAL_PART_BITS) / other.m_sig;
 277   r.m_exp = m_exp - other.m_exp - SREAL_PART_BITS;
 278   r.normalize ();
 279   return r;
 280 }