Daily bump.
[official-gcc.git] / gcc / sreal.c
blob5c429c5be3a0c9ff605a79ad331f0019cd2cfb1d
1 /* Simple data type for positive real numbers for the GNU compiler.
2 Copyright (C) 2002-2014 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 3, or (at your option) any later
9 version.
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
14 for more details.
16 You should have received a copy of the GNU General Public License
17 along with GCC; see the file COPYING3. If not see
18 <http://www.gnu.org/licenses/>. */
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.
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
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.
38 Invariant: The numbers are normalized before and after each call of sreal_*.
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
45 If the number would be too large, it is set to upper bounds of these
46 conditions.
48 If the number is zero or would be too small it meets following conditions:
49 sig == 0 && exp == -SREAL_MAX_EXP
52 #include "config.h"
53 #include "system.h"
54 #include "coretypes.h"
55 #include "sreal.h"
57 static inline void copy (sreal *, sreal *);
58 static inline void shift_right (sreal *, int);
59 static void normalize (sreal *);
61 /* Print the content of struct sreal. */
63 void
64 dump_sreal (FILE *file, sreal *x)
66 #if SREAL_PART_BITS < 32
67 fprintf (file, "((" HOST_WIDE_INT_PRINT_UNSIGNED " * 2^16 + "
68 HOST_WIDE_INT_PRINT_UNSIGNED ") * 2^%d)",
69 x->sig_hi, x->sig_lo, x->exp);
70 #else
71 fprintf (file, "(" HOST_WIDE_INT_PRINT_UNSIGNED " * 2^%d)", x->sig, x->exp);
72 #endif
75 DEBUG_FUNCTION void
76 debug (sreal &ref)
78 dump_sreal (stderr, &ref);
81 DEBUG_FUNCTION void
82 debug (sreal *ptr)
84 if (ptr)
85 debug (*ptr);
86 else
87 fprintf (stderr, "<nil>\n");
91 /* Copy the sreal number. */
93 static inline void
94 copy (sreal *r, sreal *a)
96 #if SREAL_PART_BITS < 32
97 r->sig_lo = a->sig_lo;
98 r->sig_hi = a->sig_hi;
99 #else
100 r->sig = a->sig;
101 #endif
102 r->exp = a->exp;
105 /* Shift X right by S bits. Needed: 0 < S <= SREAL_BITS.
106 When the most significant bit shifted out is 1, add 1 to X (rounding). */
108 static inline void
109 shift_right (sreal *x, int s)
111 gcc_assert (s > 0);
112 gcc_assert (s <= SREAL_BITS);
113 /* Exponent should never be so large because shift_right is used only by
114 sreal_add and sreal_sub ant thus the number cannot be shifted out from
115 exponent range. */
116 gcc_assert (x->exp + s <= SREAL_MAX_EXP);
118 x->exp += s;
120 #if SREAL_PART_BITS < 32
121 if (s > SREAL_PART_BITS)
123 s -= SREAL_PART_BITS;
124 x->sig_hi += (uhwi) 1 << (s - 1);
125 x->sig_lo = x->sig_hi >> s;
126 x->sig_hi = 0;
128 else
130 x->sig_lo += (uhwi) 1 << (s - 1);
131 if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
133 x->sig_hi++;
134 x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
136 x->sig_lo >>= s;
137 x->sig_lo |= (x->sig_hi & (((uhwi) 1 << s) - 1)) << (SREAL_PART_BITS - s);
138 x->sig_hi >>= s;
140 #else
141 x->sig += (uhwi) 1 << (s - 1);
142 x->sig >>= s;
143 #endif
146 /* Normalize *X. */
148 static void
149 normalize (sreal *x)
151 #if SREAL_PART_BITS < 32
152 int shift;
153 HOST_WIDE_INT mask;
155 if (x->sig_lo == 0 && x->sig_hi == 0)
157 x->exp = -SREAL_MAX_EXP;
159 else if (x->sig_hi < SREAL_MIN_SIG)
161 if (x->sig_hi == 0)
163 /* Move lower part of significant to higher part. */
164 x->sig_hi = x->sig_lo;
165 x->sig_lo = 0;
166 x->exp -= SREAL_PART_BITS;
168 shift = 0;
169 while (x->sig_hi < SREAL_MIN_SIG)
171 x->sig_hi <<= 1;
172 x->exp--;
173 shift++;
175 /* Check underflow. */
176 if (x->exp < -SREAL_MAX_EXP)
178 x->exp = -SREAL_MAX_EXP;
179 x->sig_hi = 0;
180 x->sig_lo = 0;
182 else if (shift)
184 mask = (1 << SREAL_PART_BITS) - (1 << (SREAL_PART_BITS - shift));
185 x->sig_hi |= (x->sig_lo & mask) >> (SREAL_PART_BITS - shift);
186 x->sig_lo = (x->sig_lo << shift) & (((uhwi) 1 << SREAL_PART_BITS) - 1);
189 else if (x->sig_hi > SREAL_MAX_SIG)
191 unsigned HOST_WIDE_INT tmp = x->sig_hi;
193 /* Find out how many bits will be shifted. */
194 shift = 0;
197 tmp >>= 1;
198 shift++;
200 while (tmp > SREAL_MAX_SIG);
202 /* Round the number. */
203 x->sig_lo += (uhwi) 1 << (shift - 1);
205 x->sig_lo >>= shift;
206 x->sig_lo += ((x->sig_hi & (((uhwi) 1 << shift) - 1))
207 << (SREAL_PART_BITS - shift));
208 x->sig_hi >>= shift;
209 x->exp += shift;
210 if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
212 x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
213 x->sig_hi++;
214 if (x->sig_hi > SREAL_MAX_SIG)
216 /* x->sig_hi was SREAL_MAX_SIG before increment
217 so now last bit is zero. */
218 x->sig_hi >>= 1;
219 x->sig_lo >>= 1;
220 x->exp++;
224 /* Check overflow. */
225 if (x->exp > SREAL_MAX_EXP)
227 x->exp = SREAL_MAX_EXP;
228 x->sig_hi = SREAL_MAX_SIG;
229 x->sig_lo = SREAL_MAX_SIG;
232 #else
233 if (x->sig == 0)
235 x->exp = -SREAL_MAX_EXP;
237 else if (x->sig < SREAL_MIN_SIG)
241 x->sig <<= 1;
242 x->exp--;
244 while (x->sig < SREAL_MIN_SIG);
246 /* Check underflow. */
247 if (x->exp < -SREAL_MAX_EXP)
249 x->exp = -SREAL_MAX_EXP;
250 x->sig = 0;
253 else if (x->sig > SREAL_MAX_SIG)
255 int last_bit;
258 last_bit = x->sig & 1;
259 x->sig >>= 1;
260 x->exp++;
262 while (x->sig > SREAL_MAX_SIG);
264 /* Round the number. */
265 x->sig += last_bit;
266 if (x->sig > SREAL_MAX_SIG)
268 x->sig >>= 1;
269 x->exp++;
272 /* Check overflow. */
273 if (x->exp > SREAL_MAX_EXP)
275 x->exp = SREAL_MAX_EXP;
276 x->sig = SREAL_MAX_SIG;
279 #endif
282 /* Set *R to SIG * 2 ^ EXP. Return R. */
284 sreal *
285 sreal_init (sreal *r, unsigned HOST_WIDE_INT sig, signed int exp)
287 #if SREAL_PART_BITS < 32
288 r->sig_lo = 0;
289 r->sig_hi = sig;
290 r->exp = exp - 16;
291 #else
292 r->sig = sig;
293 r->exp = exp;
294 #endif
295 normalize (r);
296 return r;
299 /* Return integer value of *R. */
301 HOST_WIDE_INT
302 sreal_to_int (sreal *r)
304 #if SREAL_PART_BITS < 32
305 if (r->exp <= -SREAL_BITS)
306 return 0;
307 if (r->exp >= 0)
308 return MAX_HOST_WIDE_INT;
309 return ((r->sig_hi << SREAL_PART_BITS) + r->sig_lo) >> -r->exp;
310 #else
311 if (r->exp <= -SREAL_BITS)
312 return 0;
313 if (r->exp >= SREAL_PART_BITS)
314 return MAX_HOST_WIDE_INT;
315 if (r->exp > 0)
316 return r->sig << r->exp;
317 if (r->exp < 0)
318 return r->sig >> -r->exp;
319 return r->sig;
320 #endif
323 /* Compare *A and *B. Return -1 if *A < *B, 1 if *A > *B and 0 if *A == *B. */
326 sreal_compare (sreal *a, sreal *b)
328 if (a->exp > b->exp)
329 return 1;
330 if (a->exp < b->exp)
331 return -1;
332 #if SREAL_PART_BITS < 32
333 if (a->sig_hi > b->sig_hi)
334 return 1;
335 if (a->sig_hi < b->sig_hi)
336 return -1;
337 if (a->sig_lo > b->sig_lo)
338 return 1;
339 if (a->sig_lo < b->sig_lo)
340 return -1;
341 #else
342 if (a->sig > b->sig)
343 return 1;
344 if (a->sig < b->sig)
345 return -1;
346 #endif
347 return 0;
350 /* *R = *A + *B. Return R. */
352 sreal *
353 sreal_add (sreal *r, sreal *a, sreal *b)
355 int dexp;
356 sreal tmp;
357 sreal *bb;
359 if (sreal_compare (a, b) < 0)
361 sreal *swap;
362 swap = a;
363 a = b;
364 b = swap;
367 dexp = a->exp - b->exp;
368 r->exp = a->exp;
369 if (dexp > SREAL_BITS)
371 #if SREAL_PART_BITS < 32
372 r->sig_hi = a->sig_hi;
373 r->sig_lo = a->sig_lo;
374 #else
375 r->sig = a->sig;
376 #endif
377 return r;
380 if (dexp == 0)
381 bb = b;
382 else
384 copy (&tmp, b);
385 shift_right (&tmp, dexp);
386 bb = &tmp;
389 #if SREAL_PART_BITS < 32
390 r->sig_hi = a->sig_hi + bb->sig_hi;
391 r->sig_lo = a->sig_lo + bb->sig_lo;
392 if (r->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
394 r->sig_hi++;
395 r->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
397 #else
398 r->sig = a->sig + bb->sig;
399 #endif
400 normalize (r);
401 return r;
404 /* *R = *A - *B. Return R. */
406 sreal *
407 sreal_sub (sreal *r, sreal *a, sreal *b)
409 int dexp;
410 sreal tmp;
411 sreal *bb;
413 gcc_assert (sreal_compare (a, b) >= 0);
415 dexp = a->exp - b->exp;
416 r->exp = a->exp;
417 if (dexp > SREAL_BITS)
419 #if SREAL_PART_BITS < 32
420 r->sig_hi = a->sig_hi;
421 r->sig_lo = a->sig_lo;
422 #else
423 r->sig = a->sig;
424 #endif
425 return r;
427 if (dexp == 0)
428 bb = b;
429 else
431 copy (&tmp, b);
432 shift_right (&tmp, dexp);
433 bb = &tmp;
436 #if SREAL_PART_BITS < 32
437 if (a->sig_lo < bb->sig_lo)
439 r->sig_hi = a->sig_hi - bb->sig_hi - 1;
440 r->sig_lo = a->sig_lo + ((uhwi) 1 << SREAL_PART_BITS) - bb->sig_lo;
442 else
444 r->sig_hi = a->sig_hi - bb->sig_hi;
445 r->sig_lo = a->sig_lo - bb->sig_lo;
447 #else
448 r->sig = a->sig - bb->sig;
449 #endif
450 normalize (r);
451 return r;
454 /* *R = *A * *B. Return R. */
456 sreal *
457 sreal_mul (sreal *r, sreal *a, sreal *b)
459 #if SREAL_PART_BITS < 32
460 if (a->sig_hi < SREAL_MIN_SIG || b->sig_hi < SREAL_MIN_SIG)
462 r->sig_lo = 0;
463 r->sig_hi = 0;
464 r->exp = -SREAL_MAX_EXP;
466 else
468 unsigned HOST_WIDE_INT tmp1, tmp2, tmp3;
469 if (sreal_compare (a, b) < 0)
471 sreal *swap;
472 swap = a;
473 a = b;
474 b = swap;
477 r->exp = a->exp + b->exp + SREAL_PART_BITS;
479 tmp1 = a->sig_lo * b->sig_lo;
480 tmp2 = a->sig_lo * b->sig_hi;
481 tmp3 = a->sig_hi * b->sig_lo + (tmp1 >> SREAL_PART_BITS);
483 r->sig_hi = a->sig_hi * b->sig_hi;
484 r->sig_hi += (tmp2 >> SREAL_PART_BITS) + (tmp3 >> SREAL_PART_BITS);
485 tmp2 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
486 tmp3 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
487 tmp1 = tmp2 + tmp3;
489 r->sig_lo = tmp1 & (((uhwi) 1 << SREAL_PART_BITS) - 1);
490 r->sig_hi += tmp1 >> SREAL_PART_BITS;
492 normalize (r);
494 #else
495 if (a->sig < SREAL_MIN_SIG || b->sig < SREAL_MIN_SIG)
497 r->sig = 0;
498 r->exp = -SREAL_MAX_EXP;
500 else
502 r->sig = a->sig * b->sig;
503 r->exp = a->exp + b->exp;
504 normalize (r);
506 #endif
507 return r;
510 /* *R = *A / *B. Return R. */
512 sreal *
513 sreal_div (sreal *r, sreal *a, sreal *b)
515 #if SREAL_PART_BITS < 32
516 unsigned HOST_WIDE_INT tmp, tmp1, tmp2;
518 gcc_assert (b->sig_hi >= SREAL_MIN_SIG);
519 if (a->sig_hi < SREAL_MIN_SIG)
521 r->sig_hi = 0;
522 r->sig_lo = 0;
523 r->exp = -SREAL_MAX_EXP;
525 else
527 /* Since division by the whole number is pretty ugly to write
528 we are dividing by first 3/4 of bits of number. */
530 tmp1 = (a->sig_hi << SREAL_PART_BITS) + a->sig_lo;
531 tmp2 = ((b->sig_hi << (SREAL_PART_BITS / 2))
532 + (b->sig_lo >> (SREAL_PART_BITS / 2)));
533 if (b->sig_lo & ((uhwi) 1 << ((SREAL_PART_BITS / 2) - 1)))
534 tmp2++;
536 r->sig_lo = 0;
537 tmp = tmp1 / tmp2;
538 tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
539 r->sig_hi = tmp << SREAL_PART_BITS;
541 tmp = tmp1 / tmp2;
542 tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
543 r->sig_hi += tmp << (SREAL_PART_BITS / 2);
545 tmp = tmp1 / tmp2;
546 r->sig_hi += tmp;
548 r->exp = a->exp - b->exp - SREAL_BITS - SREAL_PART_BITS / 2;
549 normalize (r);
551 #else
552 gcc_assert (b->sig != 0);
553 r->sig = (a->sig << SREAL_PART_BITS) / b->sig;
554 r->exp = a->exp - b->exp - SREAL_PART_BITS;
555 normalize (r);
556 #endif
557 return r;