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