PR other/22202
[official-gcc.git] / gcc / sreal.c
blob5982d265b522d89b4bf41c5b729019ab50af9075
1 /* Simple data type for positive real numbers for the GNU compiler.
2 Copyright (C) 2002, 2003, 2004 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
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 COPYING. If not, write to the Free
18 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
19 02110-1301, USA. */
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.
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
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.
39 Invariant: The numbers are normalized before and after each call of sreal_*.
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
46 If the number would be too large, it is set to upper bounds of these
47 conditions.
49 If the number is zero or would be too small it meets following conditions:
50 sig == 0 && exp == -SREAL_MAX_EXP
53 #include "config.h"
54 #include "system.h"
55 #include "coretypes.h"
56 #include "tm.h"
57 #include "sreal.h"
59 static inline void copy (sreal *, sreal *);
60 static inline void shift_right (sreal *, int);
61 static void normalize (sreal *);
63 /* Print the content of struct sreal. */
65 void
66 dump_sreal (FILE *file, sreal *x)
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
77 /* Copy the sreal number. */
79 static inline void
80 copy (sreal *r, sreal *a)
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;
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). */
94 static inline void
95 shift_right (sreal *x, int s)
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);
104 x->exp += s;
106 #if SREAL_PART_BITS < 32
107 if (s > SREAL_PART_BITS)
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;
114 else
116 x->sig_lo += (uhwi) 1 << (s - 1);
117 if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
119 x->sig_hi++;
120 x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
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;
126 #else
127 x->sig += (uhwi) 1 << (s - 1);
128 x->sig >>= s;
129 #endif
132 /* Normalize *X. */
134 static void
135 normalize (sreal *x)
137 #if SREAL_PART_BITS < 32
138 int shift;
139 HOST_WIDE_INT mask;
141 if (x->sig_lo == 0 && x->sig_hi == 0)
143 x->exp = -SREAL_MAX_EXP;
145 else if (x->sig_hi < SREAL_MIN_SIG)
147 if (x->sig_hi == 0)
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;
154 shift = 0;
155 while (x->sig_hi < SREAL_MIN_SIG)
157 x->sig_hi <<= 1;
158 x->exp--;
159 shift++;
161 /* Check underflow. */
162 if (x->exp < -SREAL_MAX_EXP)
164 x->exp = -SREAL_MAX_EXP;
165 x->sig_hi = 0;
166 x->sig_lo = 0;
168 else if (shift)
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);
175 else if (x->sig_hi > SREAL_MAX_SIG)
177 unsigned HOST_WIDE_INT tmp = x->sig_hi;
179 /* Find out how many bits will be shifted. */
180 shift = 0;
183 tmp >>= 1;
184 shift++;
186 while (tmp > SREAL_MAX_SIG);
188 /* Round the number. */
189 x->sig_lo += (uhwi) 1 << (shift - 1);
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))
198 x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
199 x->sig_hi++;
200 if (x->sig_hi > SREAL_MAX_SIG)
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++;
210 /* Check overflow. */
211 if (x->exp > SREAL_MAX_EXP)
213 x->exp = SREAL_MAX_EXP;
214 x->sig_hi = SREAL_MAX_SIG;
215 x->sig_lo = SREAL_MAX_SIG;
218 #else
219 if (x->sig == 0)
221 x->exp = -SREAL_MAX_EXP;
223 else if (x->sig < SREAL_MIN_SIG)
227 x->sig <<= 1;
228 x->exp--;
230 while (x->sig < SREAL_MIN_SIG);
232 /* Check underflow. */
233 if (x->exp < -SREAL_MAX_EXP)
235 x->exp = -SREAL_MAX_EXP;
236 x->sig = 0;
239 else if (x->sig > SREAL_MAX_SIG)
241 int last_bit;
244 last_bit = x->sig & 1;
245 x->sig >>= 1;
246 x->exp++;
248 while (x->sig > SREAL_MAX_SIG);
250 /* Round the number. */
251 x->sig += last_bit;
252 if (x->sig > SREAL_MAX_SIG)
254 x->sig >>= 1;
255 x->exp++;
258 /* Check overflow. */
259 if (x->exp > SREAL_MAX_EXP)
261 x->exp = SREAL_MAX_EXP;
262 x->sig = SREAL_MAX_SIG;
265 #endif
268 /* Set *R to SIG * 2 ^ EXP. Return R. */
270 sreal *
271 sreal_init (sreal *r, unsigned HOST_WIDE_INT sig, signed int exp)
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;
285 /* Return integer value of *R. */
287 HOST_WIDE_INT
288 sreal_to_int (sreal *r)
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
309 /* Compare *A and *B. Return -1 if *A < *B, 1 if *A > *B and 0 if *A == *B. */
312 sreal_compare (sreal *a, sreal *b)
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;
336 /* *R = *A + *B. Return R. */
338 sreal *
339 sreal_add (sreal *r, sreal *a, sreal *b)
341 int dexp;
342 sreal tmp;
343 sreal *bb;
345 if (sreal_compare (a, b) < 0)
347 sreal *swap;
348 swap = a;
349 a = b;
350 b = swap;
353 dexp = a->exp - b->exp;
354 r->exp = a->exp;
355 if (dexp > SREAL_BITS)
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;
366 if (dexp == 0)
367 bb = b;
368 else
370 copy (&tmp, b);
371 shift_right (&tmp, dexp);
372 bb = &tmp;
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))
380 r->sig_hi++;
381 r->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
383 #else
384 r->sig = a->sig + bb->sig;
385 #endif
386 normalize (r);
387 return r;
390 /* *R = *A - *B. Return R. */
392 sreal *
393 sreal_sub (sreal *r, sreal *a, sreal *b)
395 int dexp;
396 sreal tmp;
397 sreal *bb;
399 gcc_assert (sreal_compare (a, b) >= 0);
401 dexp = a->exp - b->exp;
402 r->exp = a->exp;
403 if (dexp > SREAL_BITS)
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;
413 if (dexp == 0)
414 bb = b;
415 else
417 copy (&tmp, b);
418 shift_right (&tmp, dexp);
419 bb = &tmp;
422 #if SREAL_PART_BITS < 32
423 if (a->sig_lo < bb->sig_lo)
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;
428 else
430 r->sig_hi = a->sig_hi - bb->sig_hi;
431 r->sig_lo = a->sig_lo - bb->sig_lo;
433 #else
434 r->sig = a->sig - bb->sig;
435 #endif
436 normalize (r);
437 return r;
440 /* *R = *A * *B. Return R. */
442 sreal *
443 sreal_mul (sreal *r, sreal *a, sreal *b)
445 #if SREAL_PART_BITS < 32
446 if (a->sig_hi < SREAL_MIN_SIG || b->sig_hi < SREAL_MIN_SIG)
448 r->sig_lo = 0;
449 r->sig_hi = 0;
450 r->exp = -SREAL_MAX_EXP;
452 else
454 unsigned HOST_WIDE_INT tmp1, tmp2, tmp3;
455 if (sreal_compare (a, b) < 0)
457 sreal *swap;
458 swap = a;
459 a = b;
460 b = swap;
463 r->exp = a->exp + b->exp + SREAL_PART_BITS;
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);
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;
475 r->sig_lo = tmp1 & (((uhwi) 1 << SREAL_PART_BITS) - 1);
476 r->sig_hi += tmp1 >> SREAL_PART_BITS;
478 normalize (r);
480 #else
481 if (a->sig < SREAL_MIN_SIG || b->sig < SREAL_MIN_SIG)
483 r->sig = 0;
484 r->exp = -SREAL_MAX_EXP;
486 else
488 r->sig = a->sig * b->sig;
489 r->exp = a->exp + b->exp;
490 normalize (r);
492 #endif
493 return r;
496 /* *R = *A / *B. Return R. */
498 sreal *
499 sreal_div (sreal *r, sreal *a, sreal *b)
501 #if SREAL_PART_BITS < 32
502 unsigned HOST_WIDE_INT tmp, tmp1, tmp2;
504 gcc_assert (b->sig_hi >= SREAL_MIN_SIG);
505 if (a->sig_hi < SREAL_MIN_SIG)
507 r->sig_hi = 0;
508 r->sig_lo = 0;
509 r->exp = -SREAL_MAX_EXP;
511 else
513 /* Since division by the whole number is pretty ugly to write
514 we are dividing by first 3/4 of bits of number. */
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++;
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;
527 tmp = tmp1 / tmp2;
528 tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
529 r->sig_hi += tmp << (SREAL_PART_BITS / 2);
531 tmp = tmp1 / tmp2;
532 r->sig_hi += tmp;
534 r->exp = a->exp - b->exp - SREAL_BITS - SREAL_PART_BITS / 2;
535 normalize (r);
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;