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 static inline void copy
PARAMS ((sreal
*, sreal
*));
60 static inline void shift_right
PARAMS ((sreal
*, int));
61 static void normalize
PARAMS ((sreal
*));
63 /* Print the content of struct sreal. */
70 #if SREAL_PART_BITS < 32
71 fprintf (file
, "((" HOST_WIDE_INT_PRINT_UNSIGNED
" * 2^16 + "
72 HOST_WIDE_INT_PRINT_UNSIGNED
") * 2^%d)",
73 x
->sig_hi
, x
->sig_lo
, x
->exp
);
75 fprintf (file
, "(" HOST_WIDE_INT_PRINT_UNSIGNED
" * 2^%d)", x
->sig
, x
->exp
);
79 /* Copy the sreal number. */
86 #if SREAL_PART_BITS < 32
87 r
->sig_lo
= a
->sig_lo
;
88 r
->sig_hi
= a
->sig_hi
;
95 /* Shift X right by S bits. Needed: 0 < S <= SREAL_BITS.
96 When the most significant bit shifted out is 1, add 1 to X (rounding). */
103 #ifdef ENABLE_CHECKING
104 if (s
<= 0 || s
> SREAL_BITS
)
106 if (x
->exp
+ s
> SREAL_MAX_EXP
)
108 /* Exponent should never be so large because shift_right is used only by
109 sreal_add and sreal_sub ant thus the number cannot be shifted out from
117 #if SREAL_PART_BITS < 32
118 if (s
> SREAL_PART_BITS
)
120 s
-= SREAL_PART_BITS
;
121 x
->sig_hi
+= (uhwi
) 1 << (s
- 1);
122 x
->sig_lo
= x
->sig_hi
>> s
;
127 x
->sig_lo
+= (uhwi
) 1 << (s
- 1);
128 if (x
->sig_lo
& ((uhwi
) 1 << SREAL_PART_BITS
))
131 x
->sig_lo
-= (uhwi
) 1 << SREAL_PART_BITS
;
134 x
->sig_lo
|= (x
->sig_hi
& (((uhwi
) 1 << s
) - 1)) << (SREAL_PART_BITS
- s
);
138 x
->sig
+= (uhwi
) 1 << (s
- 1);
149 #if SREAL_PART_BITS < 32
153 if (x
->sig_lo
== 0 && x
->sig_hi
== 0)
155 x
->exp
= -SREAL_MAX_EXP
;
157 else if (x
->sig_hi
< SREAL_MIN_SIG
)
161 /* Move lower part of significant to higher part. */
162 x
->sig_hi
= x
->sig_lo
;
164 x
->exp
-= SREAL_PART_BITS
;
167 while (x
->sig_hi
< SREAL_MIN_SIG
)
173 /* Check underflow. */
174 if (x
->exp
< -SREAL_MAX_EXP
)
176 x
->exp
= -SREAL_MAX_EXP
;
182 mask
= (1 << SREAL_PART_BITS
) - (1 << (SREAL_PART_BITS
- shift
));
183 x
->sig_hi
|= (x
->sig_lo
& mask
) >> (SREAL_PART_BITS
- shift
);
184 x
->sig_lo
= (x
->sig_lo
<< shift
) & (((uhwi
) 1 << SREAL_PART_BITS
) - 1);
187 else if (x
->sig_hi
> SREAL_MAX_SIG
)
189 unsigned HOST_WIDE_INT tmp
= x
->sig_hi
;
191 /* Find out how many bits will be shifted. */
198 while (tmp
> SREAL_MAX_SIG
);
200 /* Round the number. */
201 x
->sig_lo
+= (uhwi
) 1 << (shift
- 1);
204 x
->sig_lo
+= ((x
->sig_hi
& (((uhwi
) 1 << shift
) - 1))
205 << (SREAL_PART_BITS
- shift
));
208 if (x
->sig_lo
& ((uhwi
) 1 << SREAL_PART_BITS
))
210 x
->sig_lo
-= (uhwi
) 1 << SREAL_PART_BITS
;
212 if (x
->sig_hi
> SREAL_MAX_SIG
)
214 /* x->sig_hi was SREAL_MAX_SIG before increment
215 so now last bit is zero. */
222 /* Check overflow. */
223 if (x
->exp
> SREAL_MAX_EXP
)
225 x
->exp
= SREAL_MAX_EXP
;
226 x
->sig_hi
= SREAL_MAX_SIG
;
227 x
->sig_lo
= SREAL_MAX_SIG
;
233 x
->exp
= -SREAL_MAX_EXP
;
235 else if (x
->sig
< SREAL_MIN_SIG
)
242 while (x
->sig
< SREAL_MIN_SIG
);
244 /* Check underflow. */
245 if (x
->exp
< -SREAL_MAX_EXP
)
247 x
->exp
= -SREAL_MAX_EXP
;
251 else if (x
->sig
> SREAL_MAX_SIG
)
256 last_bit
= x
->sig
& 1;
260 while (x
->sig
> SREAL_MAX_SIG
);
262 /* Round the number. */
264 if (x
->sig
> SREAL_MAX_SIG
)
270 /* Check overflow. */
271 if (x
->exp
> SREAL_MAX_EXP
)
273 x
->exp
= SREAL_MAX_EXP
;
274 x
->sig
= SREAL_MAX_SIG
;
280 /* Set *R to SIG * 2 ^ EXP. Return R. */
283 sreal_init (r
, sig
, exp
)
285 unsigned HOST_WIDE_INT sig
;
288 #if SREAL_PART_BITS < 32
300 /* Return integer value of *R. */
306 #if SREAL_PART_BITS < 32
307 if (r
->exp
<= -SREAL_BITS
)
310 return MAX_HOST_WIDE_INT
;
311 return ((r
->sig_hi
<< SREAL_PART_BITS
) + r
->sig_lo
) >> -r
->exp
;
313 if (r
->exp
<= -SREAL_BITS
)
315 if (r
->exp
>= SREAL_PART_BITS
)
316 return MAX_HOST_WIDE_INT
;
318 return r
->sig
<< r
->exp
;
320 return r
->sig
>> -r
->exp
;
325 /* Compare *A and *B. Return -1 if *A < *B, 1 if *A > *B and 0 if *A == *B. */
336 #if SREAL_PART_BITS < 32
337 if (a
->sig_hi
> b
->sig_hi
)
339 if (a
->sig_hi
< b
->sig_hi
)
341 if (a
->sig_lo
> b
->sig_lo
)
343 if (a
->sig_lo
< b
->sig_lo
)
354 /* *R = *A + *B. Return R. */
366 if (sreal_compare (a
, b
) < 0)
374 dexp
= a
->exp
- b
->exp
;
376 if (dexp
> SREAL_BITS
)
378 #if SREAL_PART_BITS < 32
379 r
->sig_hi
= a
->sig_hi
;
380 r
->sig_lo
= a
->sig_lo
;
392 shift_right (&tmp
, dexp
);
396 #if SREAL_PART_BITS < 32
397 r
->sig_hi
= a
->sig_hi
+ bb
->sig_hi
;
398 r
->sig_lo
= a
->sig_lo
+ bb
->sig_lo
;
399 if (r
->sig_lo
& ((uhwi
) 1 << SREAL_PART_BITS
))
402 r
->sig_lo
-= (uhwi
) 1 << SREAL_PART_BITS
;
405 r
->sig
= a
->sig
+ bb
->sig
;
411 /* *R = *A - *B. Return R. */
423 if (sreal_compare (a
, b
) < 0)
428 dexp
= a
->exp
- b
->exp
;
430 if (dexp
> SREAL_BITS
)
432 #if SREAL_PART_BITS < 32
433 r
->sig_hi
= a
->sig_hi
;
434 r
->sig_lo
= a
->sig_lo
;
445 shift_right (&tmp
, dexp
);
449 #if SREAL_PART_BITS < 32
450 if (a
->sig_lo
< bb
->sig_lo
)
452 r
->sig_hi
= a
->sig_hi
- bb
->sig_hi
- 1;
453 r
->sig_lo
= a
->sig_lo
+ ((uhwi
) 1 << SREAL_PART_BITS
) - bb
->sig_lo
;
457 r
->sig_hi
= a
->sig_hi
- bb
->sig_hi
;
458 r
->sig_lo
= a
->sig_lo
- bb
->sig_lo
;
461 r
->sig
= a
->sig
- bb
->sig
;
467 /* *R = *A * *B. Return R. */
475 #if SREAL_PART_BITS < 32
476 if (a
->sig_hi
< SREAL_MIN_SIG
|| b
->sig_hi
< SREAL_MIN_SIG
)
480 r
->exp
= -SREAL_MAX_EXP
;
484 unsigned HOST_WIDE_INT tmp1
, tmp2
, tmp3
;
485 if (sreal_compare (a
, b
) < 0)
493 r
->exp
= a
->exp
+ b
->exp
+ SREAL_PART_BITS
;
495 tmp1
= a
->sig_lo
* b
->sig_lo
;
496 tmp2
= a
->sig_lo
* b
->sig_hi
;
497 tmp3
= a
->sig_hi
* b
->sig_lo
+ (tmp1
>> SREAL_PART_BITS
);
499 r
->sig_hi
= a
->sig_hi
* b
->sig_hi
;
500 r
->sig_hi
+= (tmp2
>> SREAL_PART_BITS
) + (tmp3
>> SREAL_PART_BITS
);
501 tmp2
&= ((uhwi
) 1 << SREAL_PART_BITS
) - 1;
502 tmp3
&= ((uhwi
) 1 << SREAL_PART_BITS
) - 1;
505 r
->sig_lo
= tmp1
& (((uhwi
) 1 << SREAL_PART_BITS
) - 1);
506 r
->sig_hi
+= tmp1
>> SREAL_PART_BITS
;
511 if (a
->sig
< SREAL_MIN_SIG
|| b
->sig
< SREAL_MIN_SIG
)
514 r
->exp
= -SREAL_MAX_EXP
;
518 r
->sig
= a
->sig
* b
->sig
;
519 r
->exp
= a
->exp
+ b
->exp
;
526 /* *R = *A / *B. Return R. */
534 #if SREAL_PART_BITS < 32
535 unsigned HOST_WIDE_INT tmp
, tmp1
, tmp2
;
537 if (b
->sig_hi
< SREAL_MIN_SIG
)
541 else if (a
->sig_hi
< SREAL_MIN_SIG
)
545 r
->exp
= -SREAL_MAX_EXP
;
549 /* Since division by the whole number is pretty ugly to write
550 we are dividing by first 3/4 of bits of number. */
552 tmp1
= (a
->sig_hi
<< SREAL_PART_BITS
) + a
->sig_lo
;
553 tmp2
= ((b
->sig_hi
<< (SREAL_PART_BITS
/ 2))
554 + (b
->sig_lo
>> (SREAL_PART_BITS
/ 2)));
555 if (b
->sig_lo
& ((uhwi
) 1 << ((SREAL_PART_BITS
/ 2) - 1)))
560 tmp1
= (tmp1
% tmp2
) << (SREAL_PART_BITS
/ 2);
561 r
->sig_hi
= tmp
<< SREAL_PART_BITS
;
564 tmp1
= (tmp1
% tmp2
) << (SREAL_PART_BITS
/ 2);
565 r
->sig_hi
+= tmp
<< (SREAL_PART_BITS
/ 2);
570 r
->exp
= a
->exp
- b
->exp
- SREAL_BITS
- SREAL_PART_BITS
/ 2;
580 r
->sig
= (a
->sig
<< SREAL_PART_BITS
) / b
->sig
;
581 r
->exp
= a
->exp
- b
->exp
- SREAL_PART_BITS
;