-/* Simple data type for positive real numbers for the GNU compiler.
- Copyright (C) 2002, 2003 Free Software Foundation, Inc.
+/* Simple data type for real numbers for the GNU compiler.
+ Copyright (C) 2002-2015 Free Software Foundation, Inc.
This file is part of GCC.
GCC is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free
-Software Foundation; either version 2, or (at your option) any later
+Software Foundation; either version 3, or (at your option) any later
version.
GCC is distributed in the hope that it will be useful, but WITHOUT ANY
for more details.
You should have received a copy of the GNU General Public License
-along with GCC; see the file COPYING. If not, write to the Free
-Software Foundation, 59 Temple Place - Suite 330, Boston, MA
-02111-1307, USA. */
+along with GCC; see the file COPYING3. If not see
+<http://www.gnu.org/licenses/>. */
-/* This library supports positive real numbers and 0;
+/* This library supports real numbers;
inf and nan are NOT supported.
It is written to be simple and fast.
(for < 64-bit machines sig = sig_lo + sig_hi * 2 ^ SREAL_PART_BITS)
exp = exponent
- One HOST_WIDE_INT is used for the significant on 64-bit (and more than
- 64-bit) machines,
- otherwise two HOST_WIDE_INTs are used for the significant.
+ One uint64_t is used for the significant.
Only a half of significant bits is used (in normalized sreals) so that we do
not have problems with overflow, for example when c->sig = a->sig * b->sig.
- So the precision for 64-bit and 32-bit machines is 32-bit.
+ So the precision is 32-bit.
Invariant: The numbers are normalized before and after each call of sreal_*.
#include "config.h"
#include "system.h"
+#include <math.h>
#include "coretypes.h"
-#include "tm.h"
#include "sreal.h"
-static inline void copy (sreal *, sreal *);
-static inline void shift_right (sreal *, int);
-static void normalize (sreal *);
-
/* Print the content of struct sreal. */
void
-dump_sreal (FILE *file, sreal *x)
+sreal::dump (FILE *file) const
{
-#if SREAL_PART_BITS < 32
- fprintf (file, "((" HOST_WIDE_INT_PRINT_UNSIGNED " * 2^16 + "
- HOST_WIDE_INT_PRINT_UNSIGNED ") * 2^%d)",
- x->sig_hi, x->sig_lo, x->exp);
-#else
- fprintf (file, "(" HOST_WIDE_INT_PRINT_UNSIGNED " * 2^%d)", x->sig, x->exp);
-#endif
+ fprintf (file, "(%" PRIi64 " * 2^%d)", m_sig, m_exp);
}
-/* Copy the sreal number. */
-
-static inline void
-copy (sreal *r, sreal *a)
+DEBUG_FUNCTION void
+debug (const sreal &ref)
{
-#if SREAL_PART_BITS < 32
- r->sig_lo = a->sig_lo;
- r->sig_hi = a->sig_hi;
-#else
- r->sig = a->sig;
-#endif
- r->exp = a->exp;
+ ref.dump (stderr);
}
-/* Shift X right by S bits. Needed: 0 < S <= SREAL_BITS.
- When the most significant bit shifted out is 1, add 1 to X (rounding). */
-
-static inline void
-shift_right (sreal *x, int s)
+DEBUG_FUNCTION void
+debug (const sreal *ptr)
{
-#ifdef ENABLE_CHECKING
- if (s <= 0 || s > SREAL_BITS)
- abort ();
- if (x->exp + s > SREAL_MAX_EXP)
- {
- /* Exponent should never be so large because shift_right is used only by
- sreal_add and sreal_sub ant thus the number cannot be shifted out from
- exponent range. */
- abort ();
- }
-#endif
-
- x->exp += s;
-
-#if SREAL_PART_BITS < 32
- if (s > SREAL_PART_BITS)
- {
- s -= SREAL_PART_BITS;
- x->sig_hi += (uhwi) 1 << (s - 1);
- x->sig_lo = x->sig_hi >> s;
- x->sig_hi = 0;
- }
+ if (ptr)
+ debug (*ptr);
else
- {
- x->sig_lo += (uhwi) 1 << (s - 1);
- if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
- {
- x->sig_hi++;
- x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
- }
- x->sig_lo >>= s;
- x->sig_lo |= (x->sig_hi & (((uhwi) 1 << s) - 1)) << (SREAL_PART_BITS - s);
- x->sig_hi >>= s;
- }
-#else
- x->sig += (uhwi) 1 << (s - 1);
- x->sig >>= s;
-#endif
+ fprintf (stderr, "<nil>\n");
}
-/* Normalize *X. */
+/* Shift this right by S bits. Needed: 0 < S <= SREAL_BITS.
+ When the most significant bit shifted out is 1, add 1 to this (rounding).
+ */
-static void
-normalize (sreal *x)
+void
+sreal::shift_right (int s)
{
-#if SREAL_PART_BITS < 32
- int shift;
- HOST_WIDE_INT mask;
+ gcc_checking_assert (s > 0);
+ gcc_checking_assert (s <= SREAL_BITS);
+ /* Exponent should never be so large because shift_right is used only by
+ sreal_add and sreal_sub ant thus the number cannot be shifted out from
+ exponent range. */
+ gcc_checking_assert (m_exp + s <= SREAL_MAX_EXP);
- if (x->sig_lo == 0 && x->sig_hi == 0)
- {
- x->exp = -SREAL_MAX_EXP;
- }
- else if (x->sig_hi < SREAL_MIN_SIG)
- {
- if (x->sig_hi == 0)
- {
- /* Move lower part of significant to higher part. */
- x->sig_hi = x->sig_lo;
- x->sig_lo = 0;
- x->exp -= SREAL_PART_BITS;
- }
- shift = 0;
- while (x->sig_hi < SREAL_MIN_SIG)
- {
- x->sig_hi <<= 1;
- x->exp--;
- shift++;
- }
- /* Check underflow. */
- if (x->exp < -SREAL_MAX_EXP)
- {
- x->exp = -SREAL_MAX_EXP;
- x->sig_hi = 0;
- x->sig_lo = 0;
- }
- else if (shift)
- {
- mask = (1 << SREAL_PART_BITS) - (1 << (SREAL_PART_BITS - shift));
- x->sig_hi |= (x->sig_lo & mask) >> (SREAL_PART_BITS - shift);
- x->sig_lo = (x->sig_lo << shift) & (((uhwi) 1 << SREAL_PART_BITS) - 1);
- }
- }
- else if (x->sig_hi > SREAL_MAX_SIG)
- {
- unsigned HOST_WIDE_INT tmp = x->sig_hi;
-
- /* Find out how many bits will be shifted. */
- shift = 0;
- do
- {
- tmp >>= 1;
- shift++;
- }
- while (tmp > SREAL_MAX_SIG);
-
- /* Round the number. */
- x->sig_lo += (uhwi) 1 << (shift - 1);
-
- x->sig_lo >>= shift;
- x->sig_lo += ((x->sig_hi & (((uhwi) 1 << shift) - 1))
- << (SREAL_PART_BITS - shift));
- x->sig_hi >>= shift;
- x->exp += shift;
- if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
- {
- x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
- x->sig_hi++;
- if (x->sig_hi > SREAL_MAX_SIG)
- {
- /* x->sig_hi was SREAL_MAX_SIG before increment
- so now last bit is zero. */
- x->sig_hi >>= 1;
- x->sig_lo >>= 1;
- x->exp++;
- }
- }
-
- /* Check overflow. */
- if (x->exp > SREAL_MAX_EXP)
- {
- x->exp = SREAL_MAX_EXP;
- x->sig_hi = SREAL_MAX_SIG;
- x->sig_lo = SREAL_MAX_SIG;
- }
- }
-#else
- if (x->sig == 0)
- {
- x->exp = -SREAL_MAX_EXP;
- }
- else if (x->sig < SREAL_MIN_SIG)
- {
- do
- {
- x->sig <<= 1;
- x->exp--;
- }
- while (x->sig < SREAL_MIN_SIG);
-
- /* Check underflow. */
- if (x->exp < -SREAL_MAX_EXP)
- {
- x->exp = -SREAL_MAX_EXP;
- x->sig = 0;
- }
- }
- else if (x->sig > SREAL_MAX_SIG)
- {
- int last_bit;
- do
- {
- last_bit = x->sig & 1;
- x->sig >>= 1;
- x->exp++;
- }
- while (x->sig > SREAL_MAX_SIG);
-
- /* Round the number. */
- x->sig += last_bit;
- if (x->sig > SREAL_MAX_SIG)
- {
- x->sig >>= 1;
- x->exp++;
- }
-
- /* Check overflow. */
- if (x->exp > SREAL_MAX_EXP)
- {
- x->exp = SREAL_MAX_EXP;
- x->sig = SREAL_MAX_SIG;
- }
- }
-#endif
-}
+ m_exp += s;
-/* Set *R to SIG * 2 ^ EXP. Return R. */
-
-sreal *
-sreal_init (sreal *r, unsigned HOST_WIDE_INT sig, signed int exp)
-{
-#if SREAL_PART_BITS < 32
- r->sig_lo = 0;
- r->sig_hi = sig;
- r->exp = exp - 16;
-#else
- r->sig = sig;
- r->exp = exp;
-#endif
- normalize (r);
- return r;
+ m_sig += (int64_t) 1 << (s - 1);
+ m_sig >>= s;
}
-/* Return integer value of *R. */
+/* Return integer value of *this. */
-HOST_WIDE_INT
-sreal_to_int (sreal *r)
+int64_t
+sreal::to_int () const
{
-#if SREAL_PART_BITS < 32
- if (r->exp <= -SREAL_BITS)
- return 0;
- if (r->exp >= 0)
- return MAX_HOST_WIDE_INT;
- return ((r->sig_hi << SREAL_PART_BITS) + r->sig_lo) >> -r->exp;
-#else
- if (r->exp <= -SREAL_BITS)
+ int64_t sign = m_sig < 0 ? -1 : 1;
+
+ if (m_exp <= -SREAL_BITS)
return 0;
- if (r->exp >= SREAL_PART_BITS)
- return MAX_HOST_WIDE_INT;
- if (r->exp > 0)
- return r->sig << r->exp;
- if (r->exp < 0)
- return r->sig >> -r->exp;
- return r->sig;
-#endif
+ if (m_exp >= SREAL_PART_BITS)
+ return sign * INTTYPE_MAXIMUM (int64_t);
+ if (m_exp > 0)
+ return m_sig << m_exp;
+ if (m_exp < 0)
+ return m_sig >> -m_exp;
+ return m_sig;
}
-/* Compare *A and *B. Return -1 if *A < *B, 1 if *A > *B and 0 if *A == *B. */
+/* Return value of *this as double.
+ This should be used for debug output only. */
-int
-sreal_compare (sreal *a, sreal *b)
+double
+sreal::to_double () const
{
- if (a->exp > b->exp)
- return 1;
- if (a->exp < b->exp)
- return -1;
-#if SREAL_PART_BITS < 32
- if (a->sig_hi > b->sig_hi)
- return 1;
- if (a->sig_hi < b->sig_hi)
- return -1;
- if (a->sig_lo > b->sig_lo)
- return 1;
- if (a->sig_lo < b->sig_lo)
- return -1;
-#else
- if (a->sig > b->sig)
- return 1;
- if (a->sig < b->sig)
- return -1;
-#endif
- return 0;
+ double val = m_sig;
+ if (m_exp)
+ val = ldexp (val, m_exp);
+ return val;
}
-/* *R = *A + *B. Return R. */
+/* Return *this + other. */
-sreal *
-sreal_add (sreal *r, sreal *a, sreal *b)
+sreal
+sreal::operator+ (const sreal &other) const
{
int dexp;
- sreal tmp;
- sreal *bb;
+ sreal tmp, r;
- if (sreal_compare (a, b) < 0)
- {
- sreal *swap;
- swap = a;
- a = b;
- b = swap;
- }
+ const sreal *a_p = this, *b_p = &other, *bb;
+
+ if (a_p->m_exp < b_p->m_exp)
+ std::swap (a_p, b_p);
- dexp = a->exp - b->exp;
- r->exp = a->exp;
+ dexp = a_p->m_exp - b_p->m_exp;
+ r.m_exp = a_p->m_exp;
if (dexp > SREAL_BITS)
{
-#if SREAL_PART_BITS < 32
- r->sig_hi = a->sig_hi;
- r->sig_lo = a->sig_lo;
-#else
- r->sig = a->sig;
-#endif
+ r.m_sig = a_p->m_sig;
return r;
}
if (dexp == 0)
- bb = b;
+ bb = b_p;
else
{
- copy (&tmp, b);
- shift_right (&tmp, dexp);
+ tmp = *b_p;
+ tmp.shift_right (dexp);
bb = &tmp;
}
-#if SREAL_PART_BITS < 32
- r->sig_hi = a->sig_hi + bb->sig_hi;
- r->sig_lo = a->sig_lo + bb->sig_lo;
- if (r->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
- {
- r->sig_hi++;
- r->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
- }
-#else
- r->sig = a->sig + bb->sig;
-#endif
- normalize (r);
+ r.m_sig = a_p->m_sig + bb->m_sig;
+ r.normalize ();
return r;
}
-/* *R = *A - *B. Return R. */
-sreal *
-sreal_sub (sreal *r, sreal *a, sreal *b)
+/* Return *this - other. */
+
+sreal
+sreal::operator- (const sreal &other) const
{
int dexp;
- sreal tmp;
- sreal *bb;
+ sreal tmp, r;
+ const sreal *bb;
+ const sreal *a_p = this, *b_p = &other;
- if (sreal_compare (a, b) < 0)
+ int64_t sign = 1;
+ if (a_p->m_exp < b_p->m_exp)
{
- abort ();
+ sign = -1;
+ std::swap (a_p, b_p);
}
- dexp = a->exp - b->exp;
- r->exp = a->exp;
+ dexp = a_p->m_exp - b_p->m_exp;
+ r.m_exp = a_p->m_exp;
if (dexp > SREAL_BITS)
{
-#if SREAL_PART_BITS < 32
- r->sig_hi = a->sig_hi;
- r->sig_lo = a->sig_lo;
-#else
- r->sig = a->sig;
-#endif
+ r.m_sig = sign * a_p->m_sig;
return r;
}
if (dexp == 0)
- bb = b;
+ bb = b_p;
else
{
- copy (&tmp, b);
- shift_right (&tmp, dexp);
+ tmp = *b_p;
+ tmp.shift_right (dexp);
bb = &tmp;
}
-#if SREAL_PART_BITS < 32
- if (a->sig_lo < bb->sig_lo)
- {
- r->sig_hi = a->sig_hi - bb->sig_hi - 1;
- r->sig_lo = a->sig_lo + ((uhwi) 1 << SREAL_PART_BITS) - bb->sig_lo;
- }
- else
- {
- r->sig_hi = a->sig_hi - bb->sig_hi;
- r->sig_lo = a->sig_lo - bb->sig_lo;
- }
-#else
- r->sig = a->sig - bb->sig;
-#endif
- normalize (r);
+ r.m_sig = sign * (a_p->m_sig - bb->m_sig);
+ r.normalize ();
return r;
}
-/* *R = *A * *B. Return R. */
+/* Return *this * other. */
-sreal *
-sreal_mul (sreal *r, sreal *a, sreal *b)
+sreal
+sreal::operator* (const sreal &other) const
{
-#if SREAL_PART_BITS < 32
- if (a->sig_hi < SREAL_MIN_SIG || b->sig_hi < SREAL_MIN_SIG)
+ sreal r;
+ if (absu_hwi (m_sig) < SREAL_MIN_SIG || absu_hwi (other.m_sig) < SREAL_MIN_SIG)
{
- r->sig_lo = 0;
- r->sig_hi = 0;
- r->exp = -SREAL_MAX_EXP;
+ r.m_sig = 0;
+ r.m_exp = -SREAL_MAX_EXP;
}
else
{
- unsigned HOST_WIDE_INT tmp1, tmp2, tmp3;
- if (sreal_compare (a, b) < 0)
- {
- sreal *swap;
- swap = a;
- a = b;
- b = swap;
- }
-
- r->exp = a->exp + b->exp + SREAL_PART_BITS;
-
- tmp1 = a->sig_lo * b->sig_lo;
- tmp2 = a->sig_lo * b->sig_hi;
- tmp3 = a->sig_hi * b->sig_lo + (tmp1 >> SREAL_PART_BITS);
-
- r->sig_hi = a->sig_hi * b->sig_hi;
- r->sig_hi += (tmp2 >> SREAL_PART_BITS) + (tmp3 >> SREAL_PART_BITS);
- tmp2 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
- tmp3 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
- tmp1 = tmp2 + tmp3;
-
- r->sig_lo = tmp1 & (((uhwi) 1 << SREAL_PART_BITS) - 1);
- r->sig_hi += tmp1 >> SREAL_PART_BITS;
-
- normalize (r);
+ r.m_sig = m_sig * other.m_sig;
+ r.m_exp = m_exp + other.m_exp;
+ r.normalize ();
}
-#else
- if (a->sig < SREAL_MIN_SIG || b->sig < SREAL_MIN_SIG)
- {
- r->sig = 0;
- r->exp = -SREAL_MAX_EXP;
- }
- else
- {
- r->sig = a->sig * b->sig;
- r->exp = a->exp + b->exp;
- normalize (r);
- }
-#endif
+
return r;
}
-/* *R = *A / *B. Return R. */
+/* Return *this / other. */
-sreal *
-sreal_div (sreal *r, sreal *a, sreal *b)
+sreal
+sreal::operator/ (const sreal &other) const
{
-#if SREAL_PART_BITS < 32
- unsigned HOST_WIDE_INT tmp, tmp1, tmp2;
-
- if (b->sig_hi < SREAL_MIN_SIG)
- {
- abort ();
- }
- else if (a->sig_hi < SREAL_MIN_SIG)
- {
- r->sig_hi = 0;
- r->sig_lo = 0;
- r->exp = -SREAL_MAX_EXP;
- }
- else
- {
- /* Since division by the whole number is pretty ugly to write
- we are dividing by first 3/4 of bits of number. */
-
- tmp1 = (a->sig_hi << SREAL_PART_BITS) + a->sig_lo;
- tmp2 = ((b->sig_hi << (SREAL_PART_BITS / 2))
- + (b->sig_lo >> (SREAL_PART_BITS / 2)));
- if (b->sig_lo & ((uhwi) 1 << ((SREAL_PART_BITS / 2) - 1)))
- tmp2++;
-
- r->sig_lo = 0;
- tmp = tmp1 / tmp2;
- tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
- r->sig_hi = tmp << SREAL_PART_BITS;
-
- tmp = tmp1 / tmp2;
- tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
- r->sig_hi += tmp << (SREAL_PART_BITS / 2);
-
- tmp = tmp1 / tmp2;
- r->sig_hi += tmp;
-
- r->exp = a->exp - b->exp - SREAL_BITS - SREAL_PART_BITS / 2;
- normalize (r);
- }
-#else
- if (b->sig == 0)
- {
- abort ();
- }
- else
- {
- r->sig = (a->sig << SREAL_PART_BITS) / b->sig;
- r->exp = a->exp - b->exp - SREAL_PART_BITS;
- normalize (r);
- }
-#endif
+ gcc_checking_assert (other.m_sig != 0);
+ sreal r;
+ r.m_sig = (m_sig << SREAL_PART_BITS) / other.m_sig;
+ r.m_exp = m_exp - other.m_exp - SREAL_PART_BITS;
+ r.normalize ();
return r;
}
projects / gcc.git / blobdiff