1 /* Simple data type for positive real numbers for the GNU compiler.
2 Copyright (C) 2002, 2003, 2004, 2007, 2010 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 3, 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 COPYING3. If not see
18 <http://www.gnu.org/licenses/>. */
19
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.
23
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
30
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.
37
38 Invariant: The numbers are normalized before and after each call of sreal_*.
39
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
44
45 If the number would be too large, it is set to upper bounds of these
46 conditions.
47
48 If the number is zero or would be too small it meets following conditions:
49 sig == 0 && exp == -SREAL_MAX_EXP
50 */
51
52 #include "config.h"
53 #include "system.h"
54 #include "coretypes.h"
55 #include "sreal.h"
56
57 static inline void copy (sreal *, sreal *);
58 static inline void shift_right (sreal *, int);
59 static void normalize (sreal *);
60
61 /* Print the content of struct sreal. */
62
63 void
dump_sreal(FILE * file,sreal * x)64 dump_sreal (FILE *file, sreal *x)
65 {
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
73 }
74
75 /* Copy the sreal number. */
76
77 static inline void
copy(sreal * r,sreal * a)78 copy (sreal *r, sreal *a)
79 {
80 #if SREAL_PART_BITS < 32
81 r->sig_lo = a->sig_lo;
82 r->sig_hi = a->sig_hi;
83 #else
84 r->sig = a->sig;
85 #endif
86 r->exp = a->exp;
87 }
88
89 /* Shift X right by S bits. Needed: 0 < S <= SREAL_BITS.
90 When the most significant bit shifted out is 1, add 1 to X (rounding). */
91
92 static inline void
shift_right(sreal * x,int s)93 shift_right (sreal *x, int s)
94 {
95 gcc_assert (s > 0);
96 gcc_assert (s <= SREAL_BITS);
97 /* Exponent should never be so large because shift_right is used only by
98 sreal_add and sreal_sub ant thus the number cannot be shifted out from
99 exponent range. */
100 gcc_assert (x->exp + s <= SREAL_MAX_EXP);
101
102 x->exp += s;
103
104 #if SREAL_PART_BITS < 32
105 if (s > SREAL_PART_BITS)
106 {
107 s -= SREAL_PART_BITS;
108 x->sig_hi += (uhwi) 1 << (s - 1);
109 x->sig_lo = x->sig_hi >> s;
110 x->sig_hi = 0;
111 }
112 else
113 {
114 x->sig_lo += (uhwi) 1 << (s - 1);
115 if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
116 {
117 x->sig_hi++;
118 x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
119 }
120 x->sig_lo >>= s;
121 x->sig_lo |= (x->sig_hi & (((uhwi) 1 << s) - 1)) << (SREAL_PART_BITS - s);
122 x->sig_hi >>= s;
123 }
124 #else
125 x->sig += (uhwi) 1 << (s - 1);
126 x->sig >>= s;
127 #endif
128 }
129
130 /* Normalize *X. */
131
132 static void
normalize(sreal * x)133 normalize (sreal *x)
134 {
135 #if SREAL_PART_BITS < 32
136 int shift;
137 HOST_WIDE_INT mask;
138
139 if (x->sig_lo == 0 && x->sig_hi == 0)
140 {
141 x->exp = -SREAL_MAX_EXP;
142 }
143 else if (x->sig_hi < SREAL_MIN_SIG)
144 {
145 if (x->sig_hi == 0)
146 {
147 /* Move lower part of significant to higher part. */
148 x->sig_hi = x->sig_lo;
149 x->sig_lo = 0;
150 x->exp -= SREAL_PART_BITS;
151 }
152 shift = 0;
153 while (x->sig_hi < SREAL_MIN_SIG)
154 {
155 x->sig_hi <<= 1;
156 x->exp--;
157 shift++;
158 }
159 /* Check underflow. */
160 if (x->exp < -SREAL_MAX_EXP)
161 {
162 x->exp = -SREAL_MAX_EXP;
163 x->sig_hi = 0;
164 x->sig_lo = 0;
165 }
166 else if (shift)
167 {
168 mask = (1 << SREAL_PART_BITS) - (1 << (SREAL_PART_BITS - shift));
169 x->sig_hi |= (x->sig_lo & mask) >> (SREAL_PART_BITS - shift);
170 x->sig_lo = (x->sig_lo << shift) & (((uhwi) 1 << SREAL_PART_BITS) - 1);
171 }
172 }
173 else if (x->sig_hi > SREAL_MAX_SIG)
174 {
175 unsigned HOST_WIDE_INT tmp = x->sig_hi;
176
177 /* Find out how many bits will be shifted. */
178 shift = 0;
179 do
180 {
181 tmp >>= 1;
182 shift++;
183 }
184 while (tmp > SREAL_MAX_SIG);
185
186 /* Round the number. */
187 x->sig_lo += (uhwi) 1 << (shift - 1);
188
189 x->sig_lo >>= shift;
190 x->sig_lo += ((x->sig_hi & (((uhwi) 1 << shift) - 1))
191 << (SREAL_PART_BITS - shift));
192 x->sig_hi >>= shift;
193 x->exp += shift;
194 if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
195 {
196 x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
197 x->sig_hi++;
198 if (x->sig_hi > SREAL_MAX_SIG)
199 {
200 /* x->sig_hi was SREAL_MAX_SIG before increment
201 so now last bit is zero. */
202 x->sig_hi >>= 1;
203 x->sig_lo >>= 1;
204 x->exp++;
205 }
206 }
207
208 /* Check overflow. */
209 if (x->exp > SREAL_MAX_EXP)
210 {
211 x->exp = SREAL_MAX_EXP;
212 x->sig_hi = SREAL_MAX_SIG;
213 x->sig_lo = SREAL_MAX_SIG;
214 }
215 }
216 #else
217 if (x->sig == 0)
218 {
219 x->exp = -SREAL_MAX_EXP;
220 }
221 else if (x->sig < SREAL_MIN_SIG)
222 {
223 do
224 {
225 x->sig <<= 1;
226 x->exp--;
227 }
228 while (x->sig < SREAL_MIN_SIG);
229
230 /* Check underflow. */
231 if (x->exp < -SREAL_MAX_EXP)
232 {
233 x->exp = -SREAL_MAX_EXP;
234 x->sig = 0;
235 }
236 }
237 else if (x->sig > SREAL_MAX_SIG)
238 {
239 int last_bit;
240 do
241 {
242 last_bit = x->sig & 1;
243 x->sig >>= 1;
244 x->exp++;
245 }
246 while (x->sig > SREAL_MAX_SIG);
247
248 /* Round the number. */
249 x->sig += last_bit;
250 if (x->sig > SREAL_MAX_SIG)
251 {
252 x->sig >>= 1;
253 x->exp++;
254 }
255
256 /* Check overflow. */
257 if (x->exp > SREAL_MAX_EXP)
258 {
259 x->exp = SREAL_MAX_EXP;
260 x->sig = SREAL_MAX_SIG;
261 }
262 }
263 #endif
264 }
265
266 /* Set *R to SIG * 2 ^ EXP. Return R. */
267
268 sreal *
sreal_init(sreal * r,unsigned HOST_WIDE_INT sig,signed int exp)269 sreal_init (sreal *r, unsigned HOST_WIDE_INT sig, signed int exp)
270 {
271 #if SREAL_PART_BITS < 32
272 r->sig_lo = 0;
273 r->sig_hi = sig;
274 r->exp = exp - 16;
275 #else
276 r->sig = sig;
277 r->exp = exp;
278 #endif
279 normalize (r);
280 return r;
281 }
282
283 /* Return integer value of *R. */
284
285 HOST_WIDE_INT
sreal_to_int(sreal * r)286 sreal_to_int (sreal *r)
287 {
288 #if SREAL_PART_BITS < 32
289 if (r->exp <= -SREAL_BITS)
290 return 0;
291 if (r->exp >= 0)
292 return MAX_HOST_WIDE_INT;
293 return ((r->sig_hi << SREAL_PART_BITS) + r->sig_lo) >> -r->exp;
294 #else
295 if (r->exp <= -SREAL_BITS)
296 return 0;
297 if (r->exp >= SREAL_PART_BITS)
298 return MAX_HOST_WIDE_INT;
299 if (r->exp > 0)
300 return r->sig << r->exp;
301 if (r->exp < 0)
302 return r->sig >> -r->exp;
303 return r->sig;
304 #endif
305 }
306
307 /* Compare *A and *B. Return -1 if *A < *B, 1 if *A > *B and 0 if *A == *B. */
308
309 int
sreal_compare(sreal * a,sreal * b)310 sreal_compare (sreal *a, sreal *b)
311 {
312 if (a->exp > b->exp)
313 return 1;
314 if (a->exp < b->exp)
315 return -1;
316 #if SREAL_PART_BITS < 32
317 if (a->sig_hi > b->sig_hi)
318 return 1;
319 if (a->sig_hi < b->sig_hi)
320 return -1;
321 if (a->sig_lo > b->sig_lo)
322 return 1;
323 if (a->sig_lo < b->sig_lo)
324 return -1;
325 #else
326 if (a->sig > b->sig)
327 return 1;
328 if (a->sig < b->sig)
329 return -1;
330 #endif
331 return 0;
332 }
333
334 /* *R = *A + *B. Return R. */
335
336 sreal *
sreal_add(sreal * r,sreal * a,sreal * b)337 sreal_add (sreal *r, sreal *a, sreal *b)
338 {
339 int dexp;
340 sreal tmp;
341 sreal *bb;
342
343 if (sreal_compare (a, b) < 0)
344 {
345 sreal *swap;
346 swap = a;
347 a = b;
348 b = swap;
349 }
350
351 dexp = a->exp - b->exp;
352 r->exp = a->exp;
353 if (dexp > SREAL_BITS)
354 {
355 #if SREAL_PART_BITS < 32
356 r->sig_hi = a->sig_hi;
357 r->sig_lo = a->sig_lo;
358 #else
359 r->sig = a->sig;
360 #endif
361 return r;
362 }
363
364 if (dexp == 0)
365 bb = b;
366 else
367 {
368 copy (&tmp, b);
369 shift_right (&tmp, dexp);
370 bb = &tmp;
371 }
372
373 #if SREAL_PART_BITS < 32
374 r->sig_hi = a->sig_hi + bb->sig_hi;
375 r->sig_lo = a->sig_lo + bb->sig_lo;
376 if (r->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
377 {
378 r->sig_hi++;
379 r->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
380 }
381 #else
382 r->sig = a->sig + bb->sig;
383 #endif
384 normalize (r);
385 return r;
386 }
387
388 /* *R = *A - *B. Return R. */
389
390 sreal *
sreal_sub(sreal * r,sreal * a,sreal * b)391 sreal_sub (sreal *r, sreal *a, sreal *b)
392 {
393 int dexp;
394 sreal tmp;
395 sreal *bb;
396
397 gcc_assert (sreal_compare (a, b) >= 0);
398
399 dexp = a->exp - b->exp;
400 r->exp = a->exp;
401 if (dexp > SREAL_BITS)
402 {
403 #if SREAL_PART_BITS < 32
404 r->sig_hi = a->sig_hi;
405 r->sig_lo = a->sig_lo;
406 #else
407 r->sig = a->sig;
408 #endif
409 return r;
410 }
411 if (dexp == 0)
412 bb = b;
413 else
414 {
415 copy (&tmp, b);
416 shift_right (&tmp, dexp);
417 bb = &tmp;
418 }
419
420 #if SREAL_PART_BITS < 32
421 if (a->sig_lo < bb->sig_lo)
422 {
423 r->sig_hi = a->sig_hi - bb->sig_hi - 1;
424 r->sig_lo = a->sig_lo + ((uhwi) 1 << SREAL_PART_BITS) - bb->sig_lo;
425 }
426 else
427 {
428 r->sig_hi = a->sig_hi - bb->sig_hi;
429 r->sig_lo = a->sig_lo - bb->sig_lo;
430 }
431 #else
432 r->sig = a->sig - bb->sig;
433 #endif
434 normalize (r);
435 return r;
436 }
437
438 /* *R = *A * *B. Return R. */
439
440 sreal *
sreal_mul(sreal * r,sreal * a,sreal * b)441 sreal_mul (sreal *r, sreal *a, sreal *b)
442 {
443 #if SREAL_PART_BITS < 32
444 if (a->sig_hi < SREAL_MIN_SIG || b->sig_hi < SREAL_MIN_SIG)
445 {
446 r->sig_lo = 0;
447 r->sig_hi = 0;
448 r->exp = -SREAL_MAX_EXP;
449 }
450 else
451 {
452 unsigned HOST_WIDE_INT tmp1, tmp2, tmp3;
453 if (sreal_compare (a, b) < 0)
454 {
455 sreal *swap;
456 swap = a;
457 a = b;
458 b = swap;
459 }
460
461 r->exp = a->exp + b->exp + SREAL_PART_BITS;
462
463 tmp1 = a->sig_lo * b->sig_lo;
464 tmp2 = a->sig_lo * b->sig_hi;
465 tmp3 = a->sig_hi * b->sig_lo + (tmp1 >> SREAL_PART_BITS);
466
467 r->sig_hi = a->sig_hi * b->sig_hi;
468 r->sig_hi += (tmp2 >> SREAL_PART_BITS) + (tmp3 >> SREAL_PART_BITS);
469 tmp2 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
470 tmp3 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
471 tmp1 = tmp2 + tmp3;
472
473 r->sig_lo = tmp1 & (((uhwi) 1 << SREAL_PART_BITS) - 1);
474 r->sig_hi += tmp1 >> SREAL_PART_BITS;
475
476 normalize (r);
477 }
478 #else
479 if (a->sig < SREAL_MIN_SIG || b->sig < SREAL_MIN_SIG)
480 {
481 r->sig = 0;
482 r->exp = -SREAL_MAX_EXP;
483 }
484 else
485 {
486 r->sig = a->sig * b->sig;
487 r->exp = a->exp + b->exp;
488 normalize (r);
489 }
490 #endif
491 return r;
492 }
493
494 /* *R = *A / *B. Return R. */
495
496 sreal *
sreal_div(sreal * r,sreal * a,sreal * b)497 sreal_div (sreal *r, sreal *a, sreal *b)
498 {
499 #if SREAL_PART_BITS < 32
500 unsigned HOST_WIDE_INT tmp, tmp1, tmp2;
501
502 gcc_assert (b->sig_hi >= SREAL_MIN_SIG);
503 if (a->sig_hi < SREAL_MIN_SIG)
504 {
505 r->sig_hi = 0;
506 r->sig_lo = 0;
507 r->exp = -SREAL_MAX_EXP;
508 }
509 else
510 {
511 /* Since division by the whole number is pretty ugly to write
512 we are dividing by first 3/4 of bits of number. */
513
514 tmp1 = (a->sig_hi << SREAL_PART_BITS) + a->sig_lo;
515 tmp2 = ((b->sig_hi << (SREAL_PART_BITS / 2))
516 + (b->sig_lo >> (SREAL_PART_BITS / 2)));
517 if (b->sig_lo & ((uhwi) 1 << ((SREAL_PART_BITS / 2) - 1)))
518 tmp2++;
519
520 r->sig_lo = 0;
521 tmp = tmp1 / tmp2;
522 tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
523 r->sig_hi = tmp << SREAL_PART_BITS;
524
525 tmp = tmp1 / tmp2;
526 tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
527 r->sig_hi += tmp << (SREAL_PART_BITS / 2);
528
529 tmp = tmp1 / tmp2;
530 r->sig_hi += tmp;
531
532 r->exp = a->exp - b->exp - SREAL_BITS - SREAL_PART_BITS / 2;
533 normalize (r);
534 }
535 #else
536 gcc_assert (b->sig != 0);
537 r->sig = (a->sig << SREAL_PART_BITS) / b->sig;
538 r->exp = a->exp - b->exp - SREAL_PART_BITS;
539 normalize (r);
540 #endif
541 return r;
542 }
543