1 /* Simple data type for positive real numbers for the GNU compiler.
2 Copyright (C) 2002, 2003 Free Software Foundation, Inc.
4 This file is part of GCC.
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 2, or (at your option) any later
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
16 You should have received a copy of the GNU General Public License
17 along with GCC; see the file COPYING. If not, write to the Free
18 Software Foundation, 59 Temple Place - Suite 330, Boston, MA
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.
29 (for < 64-bit machines sig = sig_lo + sig_hi * 2 ^ SREAL_PART_BITS)
32 One HOST_WIDE_INT is used for the significant on 64-bit (and more than
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 precission for 64-bit and 32-bit machines is 32-bit.
39 Invariant: The numbers are normalized before and after each call of sreal_*.
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
46 If the number would be too large, it is set to upper bounds of these
49 If the number is zero or would be too small it meets following conditions:
50 sig == 0 && exp == -SREAL_MAX_EXP
55 #include "coretypes.h"
59 void dump_sreal
PARAMS ((FILE *, sreal
*));
60 static inline void copy
PARAMS ((sreal
*, sreal
*));
61 static inline void shift_right
PARAMS ((sreal
*, int));
62 static void normalize
PARAMS ((sreal
*));
64 /* Print the content of struct sreal. */
71 #if SREAL_PART_BITS < 32
72 fprintf (file
, "((" HOST_WIDE_INT_PRINT_UNSIGNED
" * 2^16 + "
73 HOST_WIDE_INT_PRINT_UNSIGNED
") * 2^%d)",
74 x
->sig_hi
, x
->sig_lo
, x
->exp
);
76 fprintf (file
, "(" HOST_WIDE_INT_PRINT_UNSIGNED
" * 2^%d)", x
->sig
, x
->exp
);
80 /* Copy the sreal number. */
87 #if SREAL_PART_BITS < 32
88 r
->sig_lo
= a
->sig_lo
;
89 r
->sig_hi
= a
->sig_hi
;
96 /* Shift X right by S bits. Needed: 0 < S <= SREAL_BITS.
97 When the most significant bit shifted out is 1, add 1 to X (rounding). */
104 #ifdef ENABLE_CHECKING
105 if (s
<= 0 || s
> SREAL_BITS
)
107 if (x
->exp
+ s
> SREAL_MAX_EXP
)
109 /* Exponent should never be so large because shift_right is used only by
110 sreal_add and sreal_sub ant thus the number cannot be shifted out from
118 #if SREAL_PART_BITS < 32
119 if (s
> SREAL_PART_BITS
)
121 s
-= SREAL_PART_BITS
;
122 x
->sig_hi
+= (uhwi
) 1 << (s
- 1);
123 x
->sig_lo
= x
->sig_hi
>> s
;
128 x
->sig_lo
+= (uhwi
) 1 << (s
- 1);
129 if (x
->sig_lo
& ((uhwi
) 1 << SREAL_PART_BITS
))
132 x
->sig_lo
-= (uhwi
) 1 << SREAL_PART_BITS
;
135 x
->sig_lo
|= (x
->sig_hi
& (((uhwi
) 1 << s
) - 1)) << (SREAL_PART_BITS
- s
);
139 x
->sig
+= (uhwi
) 1 << (s
- 1);
150 #if SREAL_PART_BITS < 32
154 if (x
->sig_lo
== 0 && x
->sig_hi
== 0)
156 x
->exp
= -SREAL_MAX_EXP
;
158 else if (x
->sig_hi
< SREAL_MIN_SIG
)
162 /* Move lower part of significant to higher part. */
163 x
->sig_hi
= x
->sig_lo
;
165 x
->exp
-= SREAL_PART_BITS
;
168 while (x
->sig_hi
< SREAL_MIN_SIG
)
174 /* Check underflow. */
175 if (x
->exp
< -SREAL_MAX_EXP
)
177 x
->exp
= -SREAL_MAX_EXP
;
183 mask
= (1 << SREAL_PART_BITS
) - (1 << (SREAL_PART_BITS
- shift
));
184 x
->sig_hi
|= (x
->sig_lo
& mask
) >> (SREAL_PART_BITS
- shift
);
185 x
->sig_lo
= (x
->sig_lo
<< shift
) & (((uhwi
) 1 << SREAL_PART_BITS
) - 1);
188 else if (x
->sig_hi
> SREAL_MAX_SIG
)
190 unsigned HOST_WIDE_INT tmp
= x
->sig_hi
;
192 /* Find out how many bits will be shifted. */
199 while (tmp
> SREAL_MAX_SIG
);
201 /* Round the number. */
202 x
->sig_lo
+= (uhwi
) 1 << (shift
- 1);
205 x
->sig_lo
+= ((x
->sig_hi
& (((uhwi
) 1 << shift
) - 1))
206 << (SREAL_PART_BITS
- shift
));
209 if (x
->sig_lo
& ((uhwi
) 1 << SREAL_PART_BITS
))
211 x
->sig_lo
-= (uhwi
) 1 << SREAL_PART_BITS
;
213 if (x
->sig_hi
> SREAL_MAX_SIG
)
215 /* x->sig_hi was SREAL_MAX_SIG before increment
216 so now last bit is zero. */
223 /* Check overflow. */
224 if (x
->exp
> SREAL_MAX_EXP
)
226 x
->exp
= SREAL_MAX_EXP
;
227 x
->sig_hi
= SREAL_MAX_SIG
;
228 x
->sig_lo
= SREAL_MAX_SIG
;
234 x
->exp
= -SREAL_MAX_EXP
;
236 else if (x
->sig
< SREAL_MIN_SIG
)
243 while (x
->sig
< SREAL_MIN_SIG
);
245 /* Check underflow. */
246 if (x
->exp
< -SREAL_MAX_EXP
)
248 x
->exp
= -SREAL_MAX_EXP
;
252 else if (x
->sig
> SREAL_MAX_SIG
)
257 last_bit
= x
->sig
& 1;
261 while (x
->sig
> SREAL_MAX_SIG
);
263 /* Round the number. */
265 if (x
->sig
> SREAL_MAX_SIG
)
271 /* Check overflow. */
272 if (x
->exp
> SREAL_MAX_EXP
)
274 x
->exp
= SREAL_MAX_EXP
;
275 x
->sig
= SREAL_MAX_SIG
;
281 /* Set *R to SIG * 2 ^ EXP. Return R. */
284 sreal_init (r
, sig
, exp
)
286 unsigned HOST_WIDE_INT sig
;
289 #if SREAL_PART_BITS < 32
301 /* Return integer value of *R. */
307 #if SREAL_PART_BITS < 32
308 if (r
->exp
<= -SREAL_BITS
)
311 return MAX_HOST_WIDE_INT
;
312 return ((r
->sig_hi
<< SREAL_PART_BITS
) + r
->sig_lo
) >> -r
->exp
;
314 if (r
->exp
<= -SREAL_BITS
)
316 if (r
->exp
>= SREAL_PART_BITS
)
317 return MAX_HOST_WIDE_INT
;
319 return r
->sig
<< r
->exp
;
321 return r
->sig
>> -r
->exp
;
326 /* Compare *A and *B. Return -1 if *A < *B, 1 if *A > *B and 0 if *A == *B. */
337 #if SREAL_PART_BITS < 32
338 if (a
->sig_hi
> b
->sig_hi
)
340 if (a
->sig_hi
< b
->sig_hi
)
342 if (a
->sig_lo
> b
->sig_lo
)
344 if (a
->sig_lo
< b
->sig_lo
)
355 /* *R = *A + *B. Return R. */
367 if (sreal_compare (a
, b
) < 0)
375 dexp
= a
->exp
- b
->exp
;
377 if (dexp
> SREAL_BITS
)
379 #if SREAL_PART_BITS < 32
380 r
->sig_hi
= a
->sig_hi
;
381 r
->sig_lo
= a
->sig_lo
;
393 shift_right (&tmp
, dexp
);
397 #if SREAL_PART_BITS < 32
398 r
->sig_hi
= a
->sig_hi
+ bb
->sig_hi
;
399 r
->sig_lo
= a
->sig_lo
+ bb
->sig_lo
;
400 if (r
->sig_lo
& ((uhwi
) 1 << SREAL_PART_BITS
))
403 r
->sig_lo
-= (uhwi
) 1 << SREAL_PART_BITS
;
406 r
->sig
= a
->sig
+ bb
->sig
;
412 /* *R = *A - *B. Return R. */
424 if (sreal_compare (a
, b
) < 0)
429 dexp
= a
->exp
- b
->exp
;
431 if (dexp
> SREAL_BITS
)
433 #if SREAL_PART_BITS < 32
434 r
->sig_hi
= a
->sig_hi
;
435 r
->sig_lo
= a
->sig_lo
;
446 shift_right (&tmp
, dexp
);
450 #if SREAL_PART_BITS < 32
451 if (a
->sig_lo
< bb
->sig_lo
)
453 r
->sig_hi
= a
->sig_hi
- bb
->sig_hi
- 1;
454 r
->sig_lo
= a
->sig_lo
+ ((uhwi
) 1 << SREAL_PART_BITS
) - bb
->sig_lo
;
458 r
->sig_hi
= a
->sig_hi
- bb
->sig_hi
;
459 r
->sig_lo
= a
->sig_lo
- bb
->sig_lo
;
462 r
->sig
= a
->sig
- bb
->sig
;
468 /* *R = *A * *B. Return R. */
476 #if SREAL_PART_BITS < 32
477 if (a
->sig_hi
< SREAL_MIN_SIG
|| b
->sig_hi
< SREAL_MIN_SIG
)
481 r
->exp
= -SREAL_MAX_EXP
;
485 unsigned HOST_WIDE_INT tmp1
, tmp2
, tmp3
;
486 if (sreal_compare (a
, b
) < 0)
494 r
->exp
= a
->exp
+ b
->exp
+ SREAL_PART_BITS
;
496 tmp1
= a
->sig_lo
* b
->sig_lo
;
497 tmp2
= a
->sig_lo
* b
->sig_hi
;
498 tmp3
= a
->sig_hi
* b
->sig_lo
+ (tmp1
>> SREAL_PART_BITS
);
500 r
->sig_hi
= a
->sig_hi
* b
->sig_hi
;
501 r
->sig_hi
+= (tmp2
>> SREAL_PART_BITS
) + (tmp3
>> SREAL_PART_BITS
);
502 tmp2
&= ((uhwi
) 1 << SREAL_PART_BITS
) - 1;
503 tmp3
&= ((uhwi
) 1 << SREAL_PART_BITS
) - 1;
506 r
->sig_lo
= tmp1
& (((uhwi
) 1 << SREAL_PART_BITS
) - 1);
507 r
->sig_hi
+= tmp1
>> SREAL_PART_BITS
;
512 if (a
->sig
< SREAL_MIN_SIG
|| b
->sig
< SREAL_MIN_SIG
)
515 r
->exp
= -SREAL_MAX_EXP
;
519 r
->sig
= a
->sig
* b
->sig
;
520 r
->exp
= a
->exp
+ b
->exp
;
527 /* *R = *A / *B. Return R. */
535 #if SREAL_PART_BITS < 32
536 unsigned HOST_WIDE_INT tmp
, tmp1
, tmp2
;
538 if (b
->sig_hi
< SREAL_MIN_SIG
)
542 else if (a
->sig_hi
< SREAL_MIN_SIG
)
546 r
->exp
= -SREAL_MAX_EXP
;
550 /* Since division by the whole number is pretty ugly to write
551 we are dividing by first 3/4 of bits of number. */
553 tmp1
= (a
->sig_hi
<< SREAL_PART_BITS
) + a
->sig_lo
;
554 tmp2
= ((b
->sig_hi
<< (SREAL_PART_BITS
/ 2))
555 + (b
->sig_lo
>> (SREAL_PART_BITS
/ 2)));
556 if (b
->sig_lo
& ((uhwi
) 1 << ((SREAL_PART_BITS
/ 2) - 1)))
561 tmp1
= (tmp1
% tmp2
) << (SREAL_PART_BITS
/ 2);
562 r
->sig_hi
= tmp
<< SREAL_PART_BITS
;
565 tmp1
= (tmp1
% tmp2
) << (SREAL_PART_BITS
/ 2);
566 r
->sig_hi
+= tmp
<< (SREAL_PART_BITS
/ 2);
571 r
->exp
= a
->exp
- b
->exp
- SREAL_BITS
- SREAL_PART_BITS
/ 2;
581 r
->sig
= (a
->sig
<< SREAL_PART_BITS
) / b
->sig
;
582 r
->exp
= a
->exp
- b
->exp
- SREAL_PART_BITS
;