1/* Simple data type for positive real numbers for the GNU compiler.
2   Copyright (C) 2002, 2003, 2004, 2007 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 positive real numbers and 0;
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 HOST_WIDE_INT is used for the significant on 64-bit (and more than
32   64-bit) machines,
33   otherwise two HOST_WIDE_INTs are used for the significant.
34   Only a half of significant bits is used (in normalized sreals) so that we do
35   not have problems with overflow, for example when c->sig = a->sig * b->sig.
36   So the precision for 64-bit and 32-bit machines is 32-bit.
37
38   Invariant: The numbers are normalized before and after each call of sreal_*.
39
40   Normalized sreals:
41   All numbers (except zero) meet following conditions:
42	 SREAL_MIN_SIG <= sig && sig <= SREAL_MAX_SIG
43	-SREAL_MAX_EXP <= exp && exp <= SREAL_MAX_EXP
44
45   If the number would be too large, it is set to upper bounds of these
46   conditions.
47
48   If the number is zero or would be too small it meets following conditions:
49	sig == 0 && exp == -SREAL_MAX_EXP
50*/
51
52#include "config.h"
53#include "system.h"
54#include "coretypes.h"
55#include "tm.h"
56#include "sreal.h"
57
58static inline void copy (sreal *, sreal *);
59static inline void shift_right (sreal *, int);
60static void normalize (sreal *);
61
62/* Print the content of struct sreal.  */
63
64void
65dump_sreal (FILE *file, sreal *x)
66{
67#if SREAL_PART_BITS < 32
68  fprintf (file, "((" HOST_WIDE_INT_PRINT_UNSIGNED " * 2^16 + "
69	   HOST_WIDE_INT_PRINT_UNSIGNED ") * 2^%d)",
70	   x->sig_hi, x->sig_lo, x->exp);
71#else
72  fprintf (file, "(" HOST_WIDE_INT_PRINT_UNSIGNED " * 2^%d)", x->sig, x->exp);
73#endif
74}
75
76/* Copy the sreal number.  */
77
78static inline void
79copy (sreal *r, sreal *a)
80{
81#if SREAL_PART_BITS < 32
82  r->sig_lo = a->sig_lo;
83  r->sig_hi = a->sig_hi;
84#else
85  r->sig = a->sig;
86#endif
87  r->exp = a->exp;
88}
89
90/* Shift X right by S bits.  Needed: 0 < S <= SREAL_BITS.
91   When the most significant bit shifted out is 1, add 1 to X (rounding).  */
92
93static inline void
94shift_right (sreal *x, int s)
95{
96  gcc_assert (s > 0);
97  gcc_assert (s <= SREAL_BITS);
98  /* Exponent should never be so large because shift_right is used only by
99     sreal_add and sreal_sub ant thus the number cannot be shifted out from
100     exponent range.  */
101  gcc_assert (x->exp + s <= SREAL_MAX_EXP);
102
103  x->exp += s;
104
105#if SREAL_PART_BITS < 32
106  if (s > SREAL_PART_BITS)
107    {
108      s -= SREAL_PART_BITS;
109      x->sig_hi += (uhwi) 1 << (s - 1);
110      x->sig_lo = x->sig_hi >> s;
111      x->sig_hi = 0;
112    }
113  else
114    {
115      x->sig_lo += (uhwi) 1 << (s - 1);
116      if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
117	{
118	  x->sig_hi++;
119	  x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
120	}
121      x->sig_lo >>= s;
122      x->sig_lo |= (x->sig_hi & (((uhwi) 1 << s) - 1)) << (SREAL_PART_BITS - s);
123      x->sig_hi >>= s;
124    }
125#else
126  x->sig += (uhwi) 1 << (s - 1);
127  x->sig >>= s;
128#endif
129}
130
131/* Normalize *X.  */
132
133static void
134normalize (sreal *x)
135{
136#if SREAL_PART_BITS < 32
137  int shift;
138  HOST_WIDE_INT mask;
139
140  if (x->sig_lo == 0 && x->sig_hi == 0)
141    {
142      x->exp = -SREAL_MAX_EXP;
143    }
144  else if (x->sig_hi < SREAL_MIN_SIG)
145    {
146      if (x->sig_hi == 0)
147	{
148	  /* Move lower part of significant to higher part.  */
149	  x->sig_hi = x->sig_lo;
150	  x->sig_lo = 0;
151	  x->exp -= SREAL_PART_BITS;
152	}
153      shift = 0;
154      while (x->sig_hi < SREAL_MIN_SIG)
155	{
156	  x->sig_hi <<= 1;
157	  x->exp--;
158	  shift++;
159	}
160      /* Check underflow.  */
161      if (x->exp < -SREAL_MAX_EXP)
162	{
163	  x->exp = -SREAL_MAX_EXP;
164	  x->sig_hi = 0;
165	  x->sig_lo = 0;
166	}
167      else if (shift)
168	{
169	  mask = (1 << SREAL_PART_BITS) - (1 << (SREAL_PART_BITS - shift));
170	  x->sig_hi |= (x->sig_lo & mask) >> (SREAL_PART_BITS - shift);
171	  x->sig_lo = (x->sig_lo << shift) & (((uhwi) 1 << SREAL_PART_BITS) - 1);
172	}
173    }
174  else if (x->sig_hi > SREAL_MAX_SIG)
175    {
176      unsigned HOST_WIDE_INT tmp = x->sig_hi;
177
178      /* Find out how many bits will be shifted.  */
179      shift = 0;
180      do
181	{
182	  tmp >>= 1;
183	  shift++;
184	}
185      while (tmp > SREAL_MAX_SIG);
186
187      /* Round the number.  */
188      x->sig_lo += (uhwi) 1 << (shift - 1);
189
190      x->sig_lo >>= shift;
191      x->sig_lo += ((x->sig_hi & (((uhwi) 1 << shift) - 1))
192		    << (SREAL_PART_BITS - shift));
193      x->sig_hi >>= shift;
194      x->exp += shift;
195      if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
196	{
197	  x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
198	  x->sig_hi++;
199	  if (x->sig_hi > SREAL_MAX_SIG)
200	    {
201	      /* x->sig_hi was SREAL_MAX_SIG before increment
202		 so now last bit is zero.  */
203	      x->sig_hi >>= 1;
204	      x->sig_lo >>= 1;
205	      x->exp++;
206	    }
207	}
208
209      /* Check overflow.  */
210      if (x->exp > SREAL_MAX_EXP)
211	{
212	  x->exp = SREAL_MAX_EXP;
213	  x->sig_hi = SREAL_MAX_SIG;
214	  x->sig_lo = SREAL_MAX_SIG;
215	}
216    }
217#else
218  if (x->sig == 0)
219    {
220      x->exp = -SREAL_MAX_EXP;
221    }
222  else if (x->sig < SREAL_MIN_SIG)
223    {
224      do
225	{
226	  x->sig <<= 1;
227	  x->exp--;
228	}
229      while (x->sig < SREAL_MIN_SIG);
230
231      /* Check underflow.  */
232      if (x->exp < -SREAL_MAX_EXP)
233	{
234	  x->exp = -SREAL_MAX_EXP;
235	  x->sig = 0;
236	}
237    }
238  else if (x->sig > SREAL_MAX_SIG)
239    {
240      int last_bit;
241      do
242	{
243	  last_bit = x->sig & 1;
244	  x->sig >>= 1;
245	  x->exp++;
246	}
247      while (x->sig > SREAL_MAX_SIG);
248
249      /* Round the number.  */
250      x->sig += last_bit;
251      if (x->sig > SREAL_MAX_SIG)
252	{
253	  x->sig >>= 1;
254	  x->exp++;
255	}
256
257      /* Check overflow.  */
258      if (x->exp > SREAL_MAX_EXP)
259	{
260	  x->exp = SREAL_MAX_EXP;
261	  x->sig = SREAL_MAX_SIG;
262	}
263    }
264#endif
265}
266
267/* Set *R to SIG * 2 ^ EXP.  Return R.  */
268
269sreal *
270sreal_init (sreal *r, unsigned HOST_WIDE_INT sig, signed int exp)
271{
272#if SREAL_PART_BITS < 32
273  r->sig_lo = 0;
274  r->sig_hi = sig;
275  r->exp = exp - 16;
276#else
277  r->sig = sig;
278  r->exp = exp;
279#endif
280  normalize (r);
281  return r;
282}
283
284/* Return integer value of *R.  */
285
286HOST_WIDE_INT
287sreal_to_int (sreal *r)
288{
289#if SREAL_PART_BITS < 32
290  if (r->exp <= -SREAL_BITS)
291    return 0;
292  if (r->exp >= 0)
293    return MAX_HOST_WIDE_INT;
294  return ((r->sig_hi << SREAL_PART_BITS) + r->sig_lo) >> -r->exp;
295#else
296  if (r->exp <= -SREAL_BITS)
297    return 0;
298  if (r->exp >= SREAL_PART_BITS)
299    return MAX_HOST_WIDE_INT;
300  if (r->exp > 0)
301    return r->sig << r->exp;
302  if (r->exp < 0)
303    return r->sig >> -r->exp;
304  return r->sig;
305#endif
306}
307
308/* Compare *A and *B. Return -1 if *A < *B, 1 if *A > *B and 0 if *A == *B.  */
309
310int
311sreal_compare (sreal *a, sreal *b)
312{
313  if (a->exp > b->exp)
314    return 1;
315  if (a->exp < b->exp)
316    return -1;
317#if SREAL_PART_BITS < 32
318  if (a->sig_hi > b->sig_hi)
319    return 1;
320  if (a->sig_hi < b->sig_hi)
321    return -1;
322  if (a->sig_lo > b->sig_lo)
323    return 1;
324  if (a->sig_lo < b->sig_lo)
325    return -1;
326#else
327  if (a->sig > b->sig)
328    return 1;
329  if (a->sig < b->sig)
330    return -1;
331#endif
332  return 0;
333}
334
335/* *R = *A + *B.  Return R.  */
336
337sreal *
338sreal_add (sreal *r, sreal *a, sreal *b)
339{
340  int dexp;
341  sreal tmp;
342  sreal *bb;
343
344  if (sreal_compare (a, b) < 0)
345    {
346      sreal *swap;
347      swap = a;
348      a = b;
349      b = swap;
350    }
351
352  dexp = a->exp - b->exp;
353  r->exp = a->exp;
354  if (dexp > SREAL_BITS)
355    {
356#if SREAL_PART_BITS < 32
357      r->sig_hi = a->sig_hi;
358      r->sig_lo = a->sig_lo;
359#else
360      r->sig = a->sig;
361#endif
362      return r;
363    }
364
365  if (dexp == 0)
366    bb = b;
367  else
368    {
369      copy (&tmp, b);
370      shift_right (&tmp, dexp);
371      bb = &tmp;
372    }
373
374#if SREAL_PART_BITS < 32
375  r->sig_hi = a->sig_hi + bb->sig_hi;
376  r->sig_lo = a->sig_lo + bb->sig_lo;
377  if (r->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
378    {
379      r->sig_hi++;
380      r->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
381    }
382#else
383  r->sig = a->sig + bb->sig;
384#endif
385  normalize (r);
386  return r;
387}
388
389/* *R = *A - *B.  Return R.  */
390
391sreal *
392sreal_sub (sreal *r, sreal *a, sreal *b)
393{
394  int dexp;
395  sreal tmp;
396  sreal *bb;
397
398  gcc_assert (sreal_compare (a, b) >= 0);
399
400  dexp = a->exp - b->exp;
401  r->exp = a->exp;
402  if (dexp > SREAL_BITS)
403    {
404#if SREAL_PART_BITS < 32
405      r->sig_hi = a->sig_hi;
406      r->sig_lo = a->sig_lo;
407#else
408      r->sig = a->sig;
409#endif
410      return r;
411    }
412  if (dexp == 0)
413    bb = b;
414  else
415    {
416      copy (&tmp, b);
417      shift_right (&tmp, dexp);
418      bb = &tmp;
419    }
420
421#if SREAL_PART_BITS < 32
422  if (a->sig_lo < bb->sig_lo)
423    {
424      r->sig_hi = a->sig_hi - bb->sig_hi - 1;
425      r->sig_lo = a->sig_lo + ((uhwi) 1 << SREAL_PART_BITS) - bb->sig_lo;
426    }
427  else
428    {
429      r->sig_hi = a->sig_hi - bb->sig_hi;
430      r->sig_lo = a->sig_lo - bb->sig_lo;
431    }
432#else
433  r->sig = a->sig - bb->sig;
434#endif
435  normalize (r);
436  return r;
437}
438
439/* *R = *A * *B.  Return R.  */
440
441sreal *
442sreal_mul (sreal *r, sreal *a, sreal *b)
443{
444#if SREAL_PART_BITS < 32
445  if (a->sig_hi < SREAL_MIN_SIG || b->sig_hi < SREAL_MIN_SIG)
446    {
447      r->sig_lo = 0;
448      r->sig_hi = 0;
449      r->exp = -SREAL_MAX_EXP;
450    }
451  else
452    {
453      unsigned HOST_WIDE_INT tmp1, tmp2, tmp3;
454      if (sreal_compare (a, b) < 0)
455	{
456	  sreal *swap;
457	  swap = a;
458	  a = b;
459	  b = swap;
460	}
461
462      r->exp = a->exp + b->exp + SREAL_PART_BITS;
463
464      tmp1 = a->sig_lo * b->sig_lo;
465      tmp2 = a->sig_lo * b->sig_hi;
466      tmp3 = a->sig_hi * b->sig_lo + (tmp1 >> SREAL_PART_BITS);
467
468      r->sig_hi = a->sig_hi * b->sig_hi;
469      r->sig_hi += (tmp2 >> SREAL_PART_BITS) + (tmp3 >> SREAL_PART_BITS);
470      tmp2 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
471      tmp3 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
472      tmp1 = tmp2 + tmp3;
473
474      r->sig_lo = tmp1 & (((uhwi) 1 << SREAL_PART_BITS) - 1);
475      r->sig_hi += tmp1 >> SREAL_PART_BITS;
476
477      normalize (r);
478    }
479#else
480  if (a->sig < SREAL_MIN_SIG || b->sig < SREAL_MIN_SIG)
481    {
482      r->sig = 0;
483      r->exp = -SREAL_MAX_EXP;
484    }
485  else
486    {
487      r->sig = a->sig * b->sig;
488      r->exp = a->exp + b->exp;
489      normalize (r);
490    }
491#endif
492  return r;
493}
494
495/* *R = *A / *B.  Return R.  */
496
497sreal *
498sreal_div (sreal *r, sreal *a, sreal *b)
499{
500#if SREAL_PART_BITS < 32
501  unsigned HOST_WIDE_INT tmp, tmp1, tmp2;
502
503  gcc_assert (b->sig_hi >= SREAL_MIN_SIG);
504  if (a->sig_hi < SREAL_MIN_SIG)
505    {
506      r->sig_hi = 0;
507      r->sig_lo = 0;
508      r->exp = -SREAL_MAX_EXP;
509    }
510  else
511    {
512      /* Since division by the whole number is pretty ugly to write
513	 we are dividing by first 3/4 of bits of number.  */
514
515      tmp1 = (a->sig_hi << SREAL_PART_BITS) + a->sig_lo;
516      tmp2 = ((b->sig_hi << (SREAL_PART_BITS / 2))
517	      + (b->sig_lo >> (SREAL_PART_BITS / 2)));
518      if (b->sig_lo & ((uhwi) 1 << ((SREAL_PART_BITS / 2) - 1)))
519	tmp2++;
520
521      r->sig_lo = 0;
522      tmp = tmp1 / tmp2;
523      tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
524      r->sig_hi = tmp << SREAL_PART_BITS;
525
526      tmp = tmp1 / tmp2;
527      tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
528      r->sig_hi += tmp << (SREAL_PART_BITS / 2);
529
530      tmp = tmp1 / tmp2;
531      r->sig_hi += tmp;
532
533      r->exp = a->exp - b->exp - SREAL_BITS - SREAL_PART_BITS / 2;
534      normalize (r);
535    }
536#else
537  gcc_assert (b->sig != 0);
538  r->sig = (a->sig << SREAL_PART_BITS) / b->sig;
539  r->exp = a->exp - b->exp - SREAL_PART_BITS;
540  normalize (r);
541#endif
542  return r;
543}
544