1/* Simple data type for real numbers for the GNU compiler.
2   Copyright (C) 2002-2015 Free Software Foundation, Inc.
3
4This file is part of GCC.
5
6GCC is free software; you can redistribute it and/or modify it under
7the terms of the GNU General Public License as published by the Free
8Software Foundation; either version 3, or (at your option) any later
9version.
10
11GCC is distributed in the hope that it will be useful, but WITHOUT ANY
12WARRANTY; without even the implied warranty of MERCHANTABILITY or
13FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
14for more details.
15
16You should have received a copy of the GNU General Public License
17along 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
58void
59sreal::dump (FILE *file) const
60{
61  fprintf (file, "(%" PRIi64 " * 2^%d)", m_sig, m_exp);
62}
63
64DEBUG_FUNCTION void
65debug (const sreal &ref)
66{
67  ref.dump (stderr);
68}
69
70DEBUG_FUNCTION void
71debug (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
83void
84sreal::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
101int64_t
102sreal::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
120double
121sreal::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
131sreal
132sreal::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
167sreal
168sreal::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
205sreal
206sreal::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
226sreal
227sreal::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}
236