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