gcc/sreal.c

   1 /* Simple data type for real numbers for the GNU compiler.
   2    Copyright (C) 2002-2015 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, "(%" PRIi64 " * 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 /* Return integer value of *this.  */
 100
 101 int64_t
 102 sreal::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
 120 double
 121 sreal::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
 131 sreal
 132 sreal::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
 167 sreal
 168 sreal::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
 205 sreal
 206 sreal::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
 226 sreal
 227 sreal::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 }