1 /* Simple data type for positive real numbers for the GNU compiler.
2 Copyright (C) 2002, 2003, 2004 Free Software Foundation, Inc.
3
4 This file is part of GCC.
5
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.
10
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.
15
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. */
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
59 static inline void copy (sreal *, sreal *);
60 static inline void shift_right (sreal *, int);
61 static void normalize (sreal *);
62
63 /* Print the content of struct sreal. */
64
65 void
dump_sreal(FILE * file,sreal * x)66 dump_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
79 static inline void
copy(sreal * r,sreal * a)80 copy (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
94 static inline void
shift_right(sreal * x,int s)95 shift_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
134 static void
normalize(sreal * x)135 normalize (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
270 sreal *
sreal_init(sreal * r,unsigned HOST_WIDE_INT sig,signed int exp)271 sreal_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
287 HOST_WIDE_INT
sreal_to_int(sreal * r)288 sreal_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
311 int
sreal_compare(sreal * a,sreal * b)312 sreal_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
338 sreal *
sreal_add(sreal * r,sreal * a,sreal * b)339 sreal_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
392 sreal *
sreal_sub(sreal * r,sreal * a,sreal * b)393 sreal_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
442 sreal *
sreal_mul(sreal * r,sreal * a,sreal * b)443 sreal_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
498 sreal *
sreal_div(sreal * r,sreal * a,sreal * b)499 sreal_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