Mercurial > hg > CbC > CbC_gcc

comparison gcc/sreal.c @ 0:a06113de4d67

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