diff options
author | Josef Zlomek <zlomekj@suse.cz> | 2003-02-10 13:34:24 +0100 |
---|---|---|
committer | Josef Zlomek <zlomek@gcc.gnu.org> | 2003-02-10 12:34:24 +0000 |
commit | ac5e69da5659335b8f8d942b29bd01dc7384b6ad (patch) | |
tree | 8de9fc53d7c5e025943cb8ddb54de48abc65a501 /gcc/sreal.c | |
parent | 637df4f4ba4066794c61f17133d820400fa34c8a (diff) | |
download | gcc-ac5e69da5659335b8f8d942b29bd01dc7384b6ad.tar.gz |
Makefile.in (sreal.o): Added.
* Makefile.in (sreal.o): Added.
(predict.o): Depends on sreal.h instead of real.h.
* sreal.c: New file.
* sreal.h: New file.
* predict.c: Use sreal.c instead of real.c.
From-SVN: r62630
Diffstat (limited to 'gcc/sreal.c')
-rw-r--r-- | gcc/sreal.c | 591 |
1 files changed, 591 insertions, 0 deletions
diff --git a/gcc/sreal.c b/gcc/sreal.c new file mode 100644 index 00000000000..54648eafeec --- /dev/null +++ b/gcc/sreal.c @@ -0,0 +1,591 @@ +/* Simple data type for positive real numbers for the GNU compiler. + Copyright (C) 2002, 2003 Free Software Foundation, Inc. + +This file is part of GCC. + +GCC is free software; you can redistribute it and/or modify it under +the terms of the GNU General Public License as published by the Free +Software Foundation; either version 2, or (at your option) any later +version. + +GCC is distributed in the hope that it will be useful, but WITHOUT ANY +WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with GCC; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA +02111-1307, USA. */ + +/* This library supports positive real numbers and 0; + inf and nan are NOT supported. + It is written to be simple and fast. + + Value of sreal is + x = sig * 2 ^ exp + where + sig = significant + (for < 64-bit machines sig = sig_lo + sig_hi * 2 ^ SREAL_PART_BITS) + exp = exponent + + One HOST_WIDE_INT is used for the significant on 64-bit (and more than + 64-bit) machines, + otherwise two HOST_WIDE_INTs are used for the significant. + Only a half of significant bits is used (in normalized sreals) so that we do + not have problems with overflow, for example when c->sig = a->sig * b->sig. + So the precission for 64-bit and 32-bit machines is 32-bit. + + Invariant: The numbers are normalized before and after each call of sreal_*. + + Normalized sreals: + All numbers (except zero) meet following conditions: + SREAL_MIN_SIG <= sig && sig <= SREAL_MAX_SIG + -SREAL_MAX_EXP <= exp && exp <= SREAL_MAX_EXP + + If the number would be too large, it is set to upper bounds of these + conditions. + + If the number is zero or would be too small it meets following conditions: + sig == 0 && exp == -SREAL_MAX_EXP +*/ + +#include "config.h" +#include "system.h" +#include "coretypes.h" +#include "tm.h" +#include "sreal.h" + +void dump_sreal PARAMS ((FILE *, sreal *)); +static inline void copy PARAMS ((sreal *, sreal *)); +static inline void shift_right PARAMS ((sreal *, int)); +static void normalize PARAMS ((sreal *)); + +/* Print the content of struct sreal. */ + +void +dump_sreal (file, x) + FILE *file; + sreal *x; +{ +#if SREAL_PART_BITS < 32 + fprintf (file, "(("); + fprintf (file, HOST_WIDE_INT_PRINT_UNSIGNED, x->sig_hi); + fprintf (file, " * 2^16 + "); + fprintf (file, HOST_WIDE_INT_PRINT_UNSIGNED, x->sig_lo); + fprintf (file, ") * 2^%d)", x->exp); +#else + fprintf (file, "("); + fprintf (file, HOST_WIDE_INT_PRINT_UNSIGNED, x->sig); + fprintf (file, " * 2^%d)", x->exp); +#endif +} + +/* Copy the sreal number. */ + +static inline void +copy (r, a) + sreal *r; + sreal *a; +{ +#if SREAL_PART_BITS < 32 + r->sig_lo = a->sig_lo; + r->sig_hi = a->sig_hi; +#else + r->sig = a->sig; +#endif + r->exp = a->exp; +} + +/* Shift X right by S bits. Needed: 0 < S <= SREAL_BITS. + When the most significant bit shifted out is 1, add 1 to X (rounding). */ + +static inline void +shift_right (x, s) + sreal *x; + int s; +{ +#ifdef ENABLE_CHECKING + if (s <= 0 || s > SREAL_BITS) + abort (); + if (x->exp + s > SREAL_MAX_EXP) + { + /* Exponent should never be so large because shift_right is used only by + sreal_add and sreal_sub ant thus the number cannot be shifted out from + exponent range. */ + abort (); + } +#endif + + x->exp += s; + +#if SREAL_PART_BITS < 32 + if (s > SREAL_PART_BITS) + { + s -= SREAL_PART_BITS; + x->sig_hi += (uhwi) 1 << (s - 1); + x->sig_lo = x->sig_hi >> s; + x->sig_hi = 0; + } + else + { + x->sig_lo += (uhwi) 1 << (s - 1); + if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS)) + { + x->sig_hi++; + x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS; + } + x->sig_lo >>= s; + x->sig_lo |= (x->sig_hi & (((uhwi) 1 << s) - 1)) << (SREAL_PART_BITS - s); + x->sig_hi >>= s; + } +#else + x->sig += (uhwi) 1 << (s - 1); + x->sig >>= s; +#endif +} + +/* Normalize *X. */ + +static void +normalize (x) + sreal *x; +{ +#if SREAL_PART_BITS < 32 + int shift; + HOST_WIDE_INT mask; + + if (x->sig_lo == 0 && x->sig_hi == 0) + { + x->exp = -SREAL_MAX_EXP; + } + else if (x->sig_hi < SREAL_MIN_SIG) + { + if (x->sig_hi == 0) + { + /* Move lower part of significant to higher part. */ + x->sig_hi = x->sig_lo; + x->sig_lo = 0; + x->exp -= SREAL_PART_BITS; + } + shift = 0; + while (x->sig_hi < SREAL_MIN_SIG) + { + x->sig_hi <<= 1; + x->exp--; + shift++; + } + /* Check underflow. */ + if (x->exp < -SREAL_MAX_EXP) + { + x->exp = -SREAL_MAX_EXP; + x->sig_hi = 0; + x->sig_lo = 0; + } + else if (shift) + { + mask = (1 << SREAL_PART_BITS) - (1 << (SREAL_PART_BITS - shift)); + x->sig_hi |= (x->sig_lo & mask) >> (SREAL_PART_BITS - shift); + x->sig_lo = (x->sig_lo << shift) & (((uhwi) 1 << SREAL_PART_BITS) - 1); + } + } + else if (x->sig_hi > SREAL_MAX_SIG) + { + unsigned HOST_WIDE_INT tmp = x->sig_hi; + + /* Find out how many bits will be shifted. */ + shift = 0; + do + { + tmp >>= 1; + shift++; + } + while (tmp > SREAL_MAX_SIG); + + /* Round the number. */ + x->sig_lo += (uhwi) 1 << (shift - 1); + + x->sig_lo >>= shift; + x->sig_lo += ((x->sig_hi & (((uhwi) 1 << shift) - 1)) + << (SREAL_PART_BITS - shift)); + x->sig_hi >>= shift; + x->exp += shift; + if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS)) + { + x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS; + x->sig_hi++; + if (x->sig_hi > SREAL_MAX_SIG) + { + /* x->sig_hi was SREAL_MAX_SIG before increment + so now last bit is zero. */ + x->sig_hi >>= 1; + x->sig_lo >>= 1; + x->exp++; + } + } + + /* Check overflow. */ + if (x->exp > SREAL_MAX_EXP) + { + x->exp = SREAL_MAX_EXP; + x->sig_hi = SREAL_MAX_SIG; + x->sig_lo = SREAL_MAX_SIG; + } + } +#else + if (x->sig == 0) + { + x->exp = -SREAL_MAX_EXP; + } + else if (x->sig < SREAL_MIN_SIG) + { + do + { + x->sig <<= 1; + x->exp--; + } + while (x->sig < SREAL_MIN_SIG); + + /* Check underflow. */ + if (x->exp < -SREAL_MAX_EXP) + { + x->exp = -SREAL_MAX_EXP; + x->sig = 0; + } + } + else if (x->sig > SREAL_MAX_SIG) + { + int last_bit; + do + { + last_bit = x->sig & 1; + x->sig >>= 1; + x->exp++; + } + while (x->sig > SREAL_MAX_SIG); + + /* Round the number. */ + x->sig += last_bit; + if (x->sig > SREAL_MAX_SIG) + { + x->sig >>= 1; + x->exp++; + } + + /* Check overflow. */ + if (x->exp > SREAL_MAX_EXP) + { + x->exp = SREAL_MAX_EXP; + x->sig = SREAL_MAX_SIG; + } + } +#endif +} + +/* Set *R to SIG * 2 ^ EXP. Return R. */ + +sreal * +sreal_init (r, sig, exp) + sreal *r; + unsigned HOST_WIDE_INT sig; + signed int exp; +{ +#if SREAL_PART_BITS < 32 + r->sig_lo = 0; + r->sig_hi = sig; + r->exp = exp - 16; +#else + r->sig = sig; + r->exp = exp; +#endif + normalize (r); + return r; +} + +/* Return integer value of *R. */ + +HOST_WIDE_INT +sreal_to_int (r) + sreal *r; +{ +#if SREAL_PART_BITS < 32 + if (r->exp <= -SREAL_BITS) + return 0; + if (r->exp >= 0) + return MAX_HOST_WIDE_INT; + return ((r->sig_hi << SREAL_PART_BITS) + r->sig_lo) >> -r->exp; +#else + if (r->exp <= -SREAL_BITS) + return 0; + if (r->exp >= SREAL_PART_BITS) + return MAX_HOST_WIDE_INT; + if (r->exp > 0) + return r->sig << r->exp; + if (r->exp < 0) + return r->sig >> -r->exp; + return r->sig; +#endif +} + +/* Compare *A and *B. Return -1 if *A < *B, 1 if *A > *B and 0 if *A == *B. */ + +int +sreal_compare (a, b) + sreal *a; + sreal *b; +{ + if (a->exp > b->exp) + return 1; + if (a->exp < b->exp) + return -1; +#if SREAL_PART_BITS < 32 + if (a->sig_hi > b->sig_hi) + return 1; + if (a->sig_hi < b->sig_hi) + return -1; + if (a->sig_lo > b->sig_lo) + return 1; + if (a->sig_lo < b->sig_lo) + return -1; +#else + if (a->sig > b->sig) + return 1; + if (a->sig < b->sig) + return -1; +#endif + return 0; +} + +/* *R = *A + *B. Return R. */ + +sreal * +sreal_add (r, a, b) + sreal *r; + sreal *a; + sreal *b; +{ + int dexp; + sreal tmp; + sreal *bb; + + if (sreal_compare (a, b) < 0) + { + sreal *swap; + swap = a; + a = b; + b = swap; + } + + dexp = a->exp - b->exp; + r->exp = a->exp; + if (dexp > SREAL_BITS) + { +#if SREAL_PART_BITS < 32 + r->sig_hi = a->sig_hi; + r->sig_lo = a->sig_lo; +#else + r->sig = a->sig; +#endif + return r; + } + + if (dexp == 0) + bb = b; + else + { + copy (&tmp, b); + shift_right (&tmp, dexp); + bb = &tmp; + } + +#if SREAL_PART_BITS < 32 + r->sig_hi = a->sig_hi + bb->sig_hi; + r->sig_lo = a->sig_lo + bb->sig_lo; + if (r->sig_lo & ((uhwi) 1 << SREAL_PART_BITS)) + { + r->sig_hi++; + r->sig_lo -= (uhwi) 1 << SREAL_PART_BITS; + } +#else + r->sig = a->sig + bb->sig; +#endif + normalize (r); + return r; +} + +/* *R = *A - *B. Return R. */ + +sreal * +sreal_sub (r, a, b) + sreal *r; + sreal *a; + sreal *b; +{ + int dexp; + sreal tmp; + sreal *bb; + + if (sreal_compare (a, b) < 0) + { + abort (); + } + + dexp = a->exp - b->exp; + r->exp = a->exp; + if (dexp > SREAL_BITS) + { +#if SREAL_PART_BITS < 32 + r->sig_hi = a->sig_hi; + r->sig_lo = a->sig_lo; +#else + r->sig = a->sig; +#endif + return r; + } + if (dexp == 0) + bb = b; + else + { + copy (&tmp, b); + shift_right (&tmp, dexp); + bb = &tmp; + } + +#if SREAL_PART_BITS < 32 + if (a->sig_lo < bb->sig_lo) + { + r->sig_hi = a->sig_hi - bb->sig_hi - 1; + r->sig_lo = a->sig_lo + ((uhwi) 1 << SREAL_PART_BITS) - bb->sig_lo; + } + else + { + r->sig_hi = a->sig_hi - bb->sig_hi; + r->sig_lo = a->sig_lo - bb->sig_lo; + } +#else + r->sig = a->sig - bb->sig; +#endif + normalize (r); + return r; +} + +/* *R = *A * *B. Return R. */ + +sreal * +sreal_mul (r, a, b) + sreal *r; + sreal *a; + sreal *b; +{ +#if SREAL_PART_BITS < 32 + if (a->sig_hi < SREAL_MIN_SIG || b->sig_hi < SREAL_MIN_SIG) + { + r->sig_lo = 0; + r->sig_hi = 0; + r->exp = -SREAL_MAX_EXP; + } + else + { + unsigned HOST_WIDE_INT tmp1, tmp2, tmp3; + if (sreal_compare (a, b) < 0) + { + sreal *swap; + swap = a; + a = b; + b = swap; + } + + r->exp = a->exp + b->exp + SREAL_PART_BITS; + + tmp1 = a->sig_lo * b->sig_lo; + tmp2 = a->sig_lo * b->sig_hi; + tmp3 = a->sig_hi * b->sig_lo + (tmp1 >> SREAL_PART_BITS); + + r->sig_hi = a->sig_hi * b->sig_hi; + r->sig_hi += (tmp2 >> SREAL_PART_BITS) + (tmp3 >> SREAL_PART_BITS); + tmp2 &= ((uhwi) 1 << SREAL_PART_BITS) - 1; + tmp3 &= ((uhwi) 1 << SREAL_PART_BITS) - 1; + tmp1 = tmp2 + tmp3; + + r->sig_lo = tmp1 & (((uhwi) 1 << SREAL_PART_BITS) - 1); + r->sig_hi += tmp1 >> SREAL_PART_BITS; + + normalize (r); + } +#else + if (a->sig < SREAL_MIN_SIG || b->sig < SREAL_MIN_SIG) + { + r->sig = 0; + r->exp = -SREAL_MAX_EXP; + } + else + { + r->sig = a->sig * b->sig; + r->exp = a->exp + b->exp; + normalize (r); + } +#endif + return r; +} + +/* *R = *A / *B. Return R. */ + +sreal * +sreal_div (r, a, b) + sreal *r; + sreal *a; + sreal *b; +{ +#if SREAL_PART_BITS < 32 + unsigned HOST_WIDE_INT tmp, tmp1, tmp2; + + if (b->sig_hi < SREAL_MIN_SIG) + { + abort (); + } + else if (a->sig_hi < SREAL_MIN_SIG) + { + r->sig_hi = 0; + r->sig_lo = 0; + r->exp = -SREAL_MAX_EXP; + } + else + { + /* Since division by the whole number is pretty ugly to write + we are dividing by first 3/4 of bits of number. */ + + tmp1 = (a->sig_hi << SREAL_PART_BITS) + a->sig_lo; + tmp2 = ((b->sig_hi << (SREAL_PART_BITS / 2)) + + (b->sig_lo >> (SREAL_PART_BITS / 2))); + if (b->sig_lo & ((uhwi) 1 << ((SREAL_PART_BITS / 2) - 1))) + tmp2++; + + r->sig_lo = 0; + tmp = tmp1 / tmp2; + tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2); + r->sig_hi = tmp << SREAL_PART_BITS; + + tmp = tmp1 / tmp2; + tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2); + r->sig_hi += tmp << (SREAL_PART_BITS / 2); + + tmp = tmp1 / tmp2; + r->sig_hi += tmp; + + r->exp = a->exp - b->exp - SREAL_BITS - SREAL_PART_BITS / 2; + normalize (r); + } +#else + if (b->sig == 0) + { + abort (); + } + else + { + r->sig = (a->sig << SREAL_PART_BITS) / b->sig; + r->exp = a->exp - b->exp - SREAL_PART_BITS; + normalize (r); + } +#endif + return r; +} |