X-Git-Url: http://nsz.repo.hu/git/?a=blobdiff_plain;f=ir%2Ftv%2Ffltcalc.c;h=7c5d6fdd15a1eed8f34c49b0d1db5674254e7965;hb=e70c22e8aa4cc93e7d72f58453c03cb9bd797e31;hp=60556050c38557e2cf915c9337840dd989ef08e8;hpb=86b5e05a8b444327d8d464453146043fd3b93f76;p=libfirm diff --git a/ir/tv/fltcalc.c b/ir/tv/fltcalc.c index 60556050c..7c5d6fdd1 100644 --- a/ir/tv/fltcalc.c +++ b/ir/tv/fltcalc.c @@ -1,7 +1,21 @@ -/* fltcalc.c - * Authors: Matthias Heil +/* + * Project: libFIRM + * File name: ir/tv/fltcalc.c + * Purpose: + * Author: + * Modified by: + * Created: 2003 + * CVS-ID: $Id$ + * Copyright: (c) 2003 Universität Karlsruhe + * Licence: This file protected by GPL - GNU GENERAL PUBLIC LICENSE. */ + +#ifdef HAVE_CONFIG_H +# include "config.h" +#endif + + #include "fltcalc.h" #include "strcalc.h" @@ -11,110 +25,116 @@ # undef NAN #endif -#include -#include -#include +#ifdef HAVE_INTTYPES_H +# include +#endif +#ifdef HAVE_STRING_H +# include +#endif +#ifdef HAVE_STDLIB_H +# include +#endif +#ifdef HAVE_ALLOCA_H +# include +#endif +#ifdef HAVE_MALLOC_H +# include +#endif #include #include -#include + +#include "xmalloc.h" typedef uint32_t UINT32; -#ifdef WORD_LITTLE_ENDIAN -#undef WORD_LITTLE_ENDIAN -#endif -#ifdef WORD_BIG_ENDIAN -#undef WORD_BIG_ENDIAN -#endif -#define WORD_LITTLE_ENDIAN #ifdef HAVE_LONG_DOUBLE -#ifdef WORD_LITTLE_ENDIAN +#ifdef WORDS_BIGENDIAN typedef union { struct { - UINT32 low; - UINT32 mid; UINT32 high; + UINT32 mid; + UINT32 low; } val; volatile long double d; } value_t; #else typedef union { struct { - UINT32 high; - UINT32 mid; UINT32 low; + UINT32 mid; + UINT32 high; } val; volatile long double d; } value_t; #endif #else -#ifdef WORD_LITTLE_ENDIAN +#ifdef WORDS_BIGENDIAN typedef union { struct { - UINT32 low; UINT32 high; + UINT32 low; } val; volatile double d; } value_t; #else typedef union { struct { - UINT32 high; UINT32 low; + UINT32 high; } val; volatile double d; } value_t; #endif #endif - -/******** - * globals - ********/ +/** + * possible float states + */ typedef enum { - NORMAL, - ZERO, - SUBNORMAL, - INF, - NAN, + NORMAL, /**< normal representation, implicit 1 */ + ZERO, /**< +/-0 */ + SUBNORMAL, /**< denormals, implicit 0 */ + INF, /**< +/-oo */ + NAN, /**< Not A Number */ } value_class_t; +/** A descriptor for an IEEE float value. */ typedef struct { - char exponent_size; - char mantissa_size; - value_class_t class; + unsigned char exponent_size; /**< size of exponent in bits */ + unsigned char mantissa_size; /**< size of mantissa in bits */ + value_class_t clss; /**< state of this float */ } descriptor_t; -#define CLEAR_BUFFER(buffer) memset(buffer, 0, CALC_BUFFER_SIZE) +#define CLEAR_BUFFER(buffer) memset(buffer, 0, calc_buffer_size) /* because variable sized structs are impossible, the internal * value is represented as a pseudo-struct char array, addressed * by macros * struct { - * char sign; // 0 for positive, 1 for negative - * char exp[VALUE_SIZE]; - * char mant[VALUE_SIZE]; + * char sign; // 0 for positive, 1 for negative + * char exp[value_size]; + * char mant[value_size]; * descriptor_t desc; * }; */ #define _sign(a) (((char*)a)[SIGN_POS]) #define _exp(a) (&((char*)a)[EXPONENT_POS]) #define _mant(a) (&((char*)a)[MANTISSA_POS]) -#define _desc(a) (*(descriptor_t *)&((char*)a)[DESCRIPTOR_POS]) +#define _desc(a) (*(descriptor_t *)&((char *)a)[DESCRIPTOR_POS]) -#define _save_result(x) memcpy((x), sc_get_buffer(), VALUE_SIZE) -#define _shift_right(x, y, b) sc_shr((x), (y), VALUE_SIZE*4, 0, (b)) -#define _shift_left(x, y, b) sc_shl((x), (y), VALUE_SIZE*4, 0, (b)) +#define _save_result(x) memcpy((x), sc_get_buffer(), value_size) +#define _shift_right(x, y, b) sc_shr((x), (y), value_size*4, 0, (b)) +#define _shift_left(x, y, b) sc_shl((x), (y), value_size*4, 0, (b)) -#define FC_DEFINE1(code) char* fc_##code(const void *a, void *result) \ - { \ - return _calc((const char*)a, NULL, FC_##code, (char*)result); \ - } +#define FC_DEFINE1(code) \ +char *fc_##code(const void *a, void *result) { \ + return _calc((const char*)a, NULL, FC_##code, (char*)result); \ +} -#define FC_DEFINE2(code) char* fc_##code(const void *a, const void *b, void *result) \ - { \ - return _calc((const char*)a, (const char*)b, FC_##code, (char*)result); \ - } +#define FC_DEFINE2(code) \ +char *fc_##code(const void *a, const void *b, void *result) { \ + return _calc((const char*)a, (const char*)b, FC_##code, (char*)result); \ +} #define FUNC_PTR(code) fc_##code @@ -132,10 +152,10 @@ typedef struct { static char *calc_buffer = NULL; -static fc_rounding_mode_t ROUNDING_MODE; +static fc_rounding_mode_t rounding_mode; -static int CALC_BUFFER_SIZE; -static int VALUE_SIZE; +static int calc_buffer_size; +static int value_size; static int SIGN_POS; static int EXPONENT_POS; static int MANTISSA_POS; @@ -159,25 +179,25 @@ static void _fail_char(const char *str, unsigned int len, int pos) } #endif -/* pack ieee-like */ +/** pack machine-like */ static char* _pack(const char *int_float, char *packed) { char *shift_val; char *temp; char *val_buffer; - temp = alloca(VALUE_SIZE); - shift_val = alloca(VALUE_SIZE); + temp = alloca(value_size); + shift_val = alloca(value_size); - switch (_desc(int_float).class) { + switch (_desc(int_float).clss) { case NAN: - val_buffer = alloca(CALC_BUFFER_SIZE); + val_buffer = alloca(calc_buffer_size); fc_get_qnan(_desc(int_float).exponent_size, _desc(int_float).mantissa_size, val_buffer); int_float = val_buffer; break; case INF: - val_buffer = alloca(CALC_BUFFER_SIZE); + val_buffer = alloca(calc_buffer_size); fc_get_plusinf(_desc(int_float).exponent_size, _desc(int_float).mantissa_size, val_buffer); _sign(val_buffer) = _sign(int_float); int_float = val_buffer; @@ -220,7 +240,7 @@ char* _normalize(const char *in_val, char *out_val, int sticky) char lsb, guard, round, round_dir = 0; char *temp; - temp = alloca(VALUE_SIZE); + temp = alloca(value_size); /* +2: save two rounding bits at the end */ hsb = 2 + _desc(in_val).mantissa_size - sc_get_highest_set_bit(_mant(in_val)) - 1; @@ -231,16 +251,16 @@ char* _normalize(const char *in_val, char *out_val, int sticky) memcpy(&_desc(out_val), &_desc(in_val), sizeof(descriptor_t)); } - _desc(out_val).class = NORMAL; + _desc(out_val).clss = NORMAL; - /* mantissa all zeroes, so zero exponent (because of explicit one)*/ + /* mantissa all zeros, so zero exponent (because of explicit one)*/ if (hsb == 2 + _desc(in_val).mantissa_size) { sc_val_from_ulong(0, _exp(out_val)); hsb = -1; } - /* shift the first 1 ito the left of the radix point (i.e. hsb == -1) */ + /* shift the first 1 into the left of the radix point (i.e. hsb == -1) */ if (hsb < -1) { /* shift right */ @@ -276,7 +296,7 @@ char* _normalize(const char *in_val, char *out_val, int sticky) /* denormalized means exponent of zero */ sc_val_from_ulong(0, _exp(out_val)); - _desc(out_val).class = SUBNORMAL; + _desc(out_val).clss = SUBNORMAL; } /* perform rounding by adding a value that clears the guard bit and the round bit @@ -286,7 +306,7 @@ char* _normalize(const char *in_val, char *out_val, int sticky) guard = (lsb&0x2)>>1; round = lsb&0x1; - switch (ROUNDING_MODE) + switch (rounding_mode) { case FC_TONEAREST: /* round to nearest representable value, if in doubt choose the version @@ -325,26 +345,26 @@ char* _normalize(const char *in_val, char *out_val, int sticky) } /* could have rounded down to zero */ - if (sc_is_zero(_mant(out_val)) && (_desc(out_val).class == SUBNORMAL)) - _desc(out_val).class = ZERO; + if (sc_is_zero(_mant(out_val)) && (_desc(out_val).clss == SUBNORMAL)) + _desc(out_val).clss = ZERO; /* check for rounding overflow */ hsb = 2 + _desc(out_val).mantissa_size - sc_get_highest_set_bit(_mant(out_val)) - 1; - if ((_desc(out_val).class != SUBNORMAL) && (hsb < -1)) + if ((_desc(out_val).clss != SUBNORMAL) && (hsb < -1)) { sc_val_from_ulong(1, temp); _shift_right(_mant(out_val), temp, _mant(out_val)); sc_add(_exp(out_val), temp, _exp(out_val)); } - else if ((_desc(out_val).class == SUBNORMAL) && (hsb == -1)) + else if ((_desc(out_val).clss == SUBNORMAL) && (hsb == -1)) { - /* overflow caused the matissa to be normal again, + /* overflow caused the mantissa to be normal again, * so adapt the exponent accordingly */ sc_val_from_ulong(1, temp); sc_add(_exp(out_val), temp, _exp(out_val)); - _desc(out_val).class = NORMAL; + _desc(out_val).clss = NORMAL; } /* no further rounding is needed, because rounding overflow means * the carry of the original rounding was propagated all the way @@ -374,10 +394,10 @@ char* _normalize(const char *in_val, char *out_val, int sticky) if (_sign(out_val) == 0) { /* value is positive */ - switch (ROUNDING_MODE) { + switch (rounding_mode) { case FC_TONEAREST: case FC_TOPOSITIVE: - _desc(out_val).class = INF; + _desc(out_val).clss = INF; break; case FC_TONEGATIVE: @@ -386,10 +406,10 @@ char* _normalize(const char *in_val, char *out_val, int sticky) } } else { /* value is negative */ - switch (ROUNDING_MODE) { + switch (rounding_mode) { case FC_TONEAREST: case FC_TONEGATIVE: - _desc(out_val).class = INF; + _desc(out_val).clss = INF; break; case FC_TOPOSITIVE: @@ -402,10 +422,26 @@ char* _normalize(const char *in_val, char *out_val, int sticky) return out_val; } -/* +/** + * Operations involving NaN's must return NaN + */ +#define handle_NAN(a, b, result) \ +do { \ + if (_desc(a).clss == NAN) { \ + if (a != result) memcpy(result, a, calc_buffer_size); \ + return result; \ + } \ + if (_desc(b).clss == NAN) { \ + if (b != result) memcpy(result, b, calc_buffer_size); \ + return result; \ + } \ +}while (0) + + +/** * calculate a + b, where a is the value with the bigger exponent */ -static char* _add(const char* a, const char* b, char* result) +static char* _fadd(const char* a, const char* b, char* result) { char *temp; char *exp_diff; @@ -413,14 +449,7 @@ static char* _add(const char* a, const char* b, char* result) char sign; char sticky; - if (_desc(a).class == NAN) { - if (a != result) memcpy(result, a, CALC_BUFFER_SIZE); - return result; - } - if (_desc(b).class == NAN) { - if (b != result) memcpy(result, b, CALC_BUFFER_SIZE); - return result; - } + handle_NAN(a, b, result); /* make sure result has a descriptor */ if (result != a && result != b) @@ -428,14 +457,13 @@ static char* _add(const char* a, const char* b, char* result) /* determine if this is an addition or subtraction */ sign = _sign(a) ^ _sign(b); - DEBUGPRINTF(("sign a: %d, sign b: %d -> %s\n", _sign(a), _sign(b), sign?"sub":"add")); - /* produce nan on inf - inf */ - if (sign && (_desc(a).class == INF) && (_desc(b).class == INF)) + /* produce NaN on inf - inf */ + if (sign && (_desc(a).clss == INF) && (_desc(b).clss == INF)) return fc_get_qnan(_desc(a).exponent_size, _desc(b).mantissa_size, result); - temp = alloca(VALUE_SIZE); - exp_diff = alloca(VALUE_SIZE); + temp = alloca(value_size); + exp_diff = alloca(value_size); /* get exponent difference */ sc_sub(_exp(a), _exp(b), exp_diff); @@ -447,50 +475,45 @@ static char* _add(const char* a, const char* b, char* result) if (sign && sc_val_to_long(exp_diff) == 0) { switch (sc_comp(_mant(a), _mant(b))) { case 1: /* a > b */ - if (_sign(a)) _sign(result) = 1; /* abs(a) is bigger and a is negative */ - else _sign(result) = 0; + _sign(result) = _sign(a); /* abs(a) is bigger and a is negative */ break; case 0: /* a == b */ - if (ROUNDING_MODE == FC_TONEGATIVE) - _sign(result) = 1; - else - _sign(result) = 0; + _sign(result) = (rounding_mode == FC_TONEGATIVE); break; case -1: /* a < b */ - if (_sign(b)) _sign(result) = 1; /* abs(b) is bigger and b is negative */ - else _sign(result) = 0; + _sign(result) = _sign(b); /* abs(b) is bigger and b is negative */ break; default: /* can't be reached */ break; } - } else { - _sign(result) = _sign(a); } + else + _sign(result) = _sign(a); /* sign has been taken care of, check for special cases */ - if (_desc(a).class == ZERO) { - if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, CALC_BUFFER_SIZE-SIGN_POS-1); + if (_desc(a).clss == ZERO) { + if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-SIGN_POS-1); return result; } - if (_desc(b).class == ZERO) { - if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-SIGN_POS-1); + if (_desc(b).clss == ZERO) { + if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-SIGN_POS-1); return result; } - if (_desc(a).class == INF) { - if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-SIGN_POS-1); + if (_desc(a).clss == INF) { + if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-SIGN_POS-1); return result; } - if (_desc(b).class == INF) { - if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, CALC_BUFFER_SIZE-SIGN_POS-1); + if (_desc(b).clss == INF) { + if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-SIGN_POS-1); return result; } /* shift the smaller value to the right to align the radix point */ /* subnormals have their radix point shifted to the right, * take care of this first */ - if ((_desc(b).class == SUBNORMAL) && (_desc(a).class != SUBNORMAL)) + if ((_desc(b).clss == SUBNORMAL) && (_desc(a).clss != SUBNORMAL)) { sc_val_from_ulong(1, temp); sc_sub(exp_diff, temp, exp_diff); @@ -504,7 +527,7 @@ static char* _add(const char* a, const char* b, char* result) /* if subtracting a little more than the represented value or adding a little * more than the represented value to a negative value this, in addition to the * still set sticky bit, takes account of the 'little more' */ - char *temp1 = alloca(CALC_BUFFER_SIZE); + char *temp1 = alloca(calc_buffer_size); sc_val_from_ulong(1, temp1); sc_add(temp, temp1, temp); } @@ -520,60 +543,56 @@ static char* _add(const char* a, const char* b, char* result) /* _normalize expects a 'normal' radix point, adding two subnormals * results in a subnormal radix point -> shifting before normalizing */ - if ((_desc(a).class == SUBNORMAL) && (_desc(b).class == SUBNORMAL)) + if ((_desc(a).clss == SUBNORMAL) && (_desc(b).clss == SUBNORMAL)) { sc_val_from_ulong(1, NULL); _shift_left(_mant(result), sc_get_buffer(), _mant(result)); } /* resulting exponent is the bigger one */ - memmove(_exp(result), _exp(a), VALUE_SIZE); + memmove(_exp(result), _exp(a), value_size); return _normalize(result, result, sticky); } -static char* _mul(const char* a, const char* b, char* result) +/** + * calculate a * b + */ +static char* _fmul(const char* a, const char* b, char* result) { char *temp; - if (_desc(a).class == NAN) { - if (a != result) memcpy(result, a, CALC_BUFFER_SIZE); - return result; - } - if (_desc(b).class == NAN) { - if (b != result) memcpy(result, b, CALC_BUFFER_SIZE); - return result; - } + handle_NAN(a, b, result); - temp = alloca(VALUE_SIZE); + temp = alloca(value_size); if (result != a && result != b) memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t)); _sign(result) = _sign(a) ^ _sign(b); - /* produce nan on 0 * inf */ - if (_desc(a).class == ZERO) { - if (_desc(b).class == INF) + /* produce NaN on 0 * inf */ + if (_desc(a).clss == ZERO) { + if (_desc(b).clss == INF) fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result); else - if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-1); + if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1); return result; } - if (_desc(b).class == ZERO) { - if (_desc(a).class == INF) + if (_desc(b).clss == ZERO) { + if (_desc(a).clss == INF) fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result); else - if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, CALC_BUFFER_SIZE-1); + if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-1); return result; } - if (_desc(a).class == INF) { - if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-1); + if (_desc(a).clss == INF) { + if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1); return result; } - if (_desc(b).class == INF) { - if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, CALC_BUFFER_SIZE-1); + if (_desc(b).clss == INF) { + if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-1); return result; } @@ -584,7 +603,7 @@ static char* _mul(const char* a, const char* b, char* result) sc_sub(_exp(result), temp, _exp(result)); /* mixed normal, subnormal values introduce an error of 1, correct it */ - if ((_desc(a).class == SUBNORMAL) ^ (_desc(b).class == SUBNORMAL)) + if ((_desc(a).clss == SUBNORMAL) ^ (_desc(b).clss == SUBNORMAL)) { sc_val_from_ulong(1, temp); sc_add(_exp(result), temp, _exp(result)); @@ -604,21 +623,17 @@ static char* _mul(const char* a, const char* b, char* result) return _normalize(result, result, sc_had_carry()); } -static char* _div(const char* a, const char* b, char* result) +/** + * calculate a / b + */ +static char* _fdiv(const char* a, const char* b, char* result) { char *temp, *dividend; - if (_desc(a).class == NAN) { - if (a != result) memcpy(result, a, CALC_BUFFER_SIZE); - return result; - } - if (_desc(b).class == NAN) { - if (b != result) memcpy(result, b, CALC_BUFFER_SIZE); - return result; - } + handle_NAN(a, b, result); - temp = alloca(VALUE_SIZE); - dividend = alloca(VALUE_SIZE); + temp = alloca(value_size); + dividend = alloca(value_size); if (result != a && result != b) memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t)); @@ -626,18 +641,18 @@ static char* _div(const char* a, const char* b, char* result) _sign(result) = _sign(a) ^ _sign(b); /* produce nan on 0/0 and inf/inf */ - if (_desc(a).class == ZERO) { - if (_desc(b).class == ZERO) + if (_desc(a).clss == ZERO) { + if (_desc(b).clss == ZERO) /* 0/0 -> nan */ fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result); else /* 0/x -> a */ - if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-1); + if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1); return result; } - if (_desc(b).class == INF) { - if (_desc(a).class == INF) + if (_desc(b).clss == INF) { + if (_desc(a).clss == INF) /* inf/inf -> nan */ fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result); else { @@ -645,17 +660,17 @@ static char* _div(const char* a, const char* b, char* result) sc_val_from_ulong(0, NULL); _save_result(_exp(result)); _save_result(_mant(result)); - _desc(result).class = ZERO; + _desc(result).clss = ZERO; } return result; } - if (_desc(a).class == INF) { + if (_desc(a).clss == INF) { /* inf/x -> inf */ - if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-1); + if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1); return result; } - if (_desc(b).class == ZERO) { + if (_desc(b).clss == ZERO) { /* division by zero */ if (_sign(result)) fc_get_minusinf(_desc(a).exponent_size, _desc(a).mantissa_size, result); @@ -670,7 +685,7 @@ static char* _div(const char* a, const char* b, char* result) sc_add(_exp(result), temp, _exp(result)); /* mixed normal, subnormal values introduce an error of 1, correct it */ - if ((_desc(a).class == SUBNORMAL) ^ (_desc(b).class == SUBNORMAL)) + if ((_desc(a).clss == SUBNORMAL) ^ (_desc(b).clss == SUBNORMAL)) { sc_val_from_ulong(1, temp); sc_add(_exp(result), temp, _exp(result)); @@ -687,7 +702,7 @@ static char* _div(const char* a, const char* b, char* result) _shift_left(_mant(a), temp, dividend); { - char *divisor = alloca(CALC_BUFFER_SIZE); + char *divisor = alloca(calc_buffer_size); sc_val_from_ulong(1, divisor); _shift_right(_mant(b), divisor, divisor); sc_div(dividend, divisor, _mant(result)); @@ -696,7 +711,7 @@ static char* _div(const char* a, const char* b, char* result) return _normalize(result, result, sc_had_carry()); } -void _power_of_ten(int exp, descriptor_t *desc, char *result) +static void _power_of_ten(int exp, descriptor_t *desc, char *result) { char *build; char *temp; @@ -708,8 +723,8 @@ void _power_of_ten(int exp, descriptor_t *desc, char *result) if (desc != NULL) memcpy(&_desc(result), desc, sizeof(descriptor_t)); - build = alloca(VALUE_SIZE); - temp = alloca(VALUE_SIZE); + build = alloca(value_size); + temp = alloca(value_size); sc_val_from_ulong((1 << _desc(result).exponent_size)/2-1, _exp(result)); @@ -725,7 +740,7 @@ void _power_of_ten(int exp, descriptor_t *desc, char *result) } _save_result(build); - /* temp is amount of leftshift needed to put the value left of the radix point */ + /* temp is amount of left shift needed to put the value left of the radix point */ sc_val_from_ulong(_desc(result).mantissa_size + 2, temp); _shift_left(build, temp, _mant(result)); @@ -734,41 +749,48 @@ void _power_of_ten(int exp, descriptor_t *desc, char *result) } } +/** + * Truncate the fractional part away. + * + * This does not clip to any integer rang. + */ static char* _trunc(const char *a, char *result) { - /* when exponent == 0 all bits left of the radix point + /* + * When exponent == 0 all bits left of the radix point * are the integral part of the value. For 15bit exp_size - * this would require a leftshift of max. 16383 bits which + * this would require a left shift of max. 16383 bits which * is too much. * But it is enough to ensure that no bit right of the radix * point remains set. This restricts the interesting * exponents to the interval [0, mant_size-1]. * Outside this interval the truncated value is either 0 or - * it is are already truncated */ + * it does not have fractional parts. + */ int exp_bias, exp_val; char *temp; - temp = alloca(VALUE_SIZE); + temp = alloca(value_size); if (a != result) memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t)); exp_bias = (1<<_desc(a).exponent_size)/2-1; - exp_val = sc_val_to_long(_exp(a)) - exp_bias; + exp_val = sc_val_to_long(_exp(a)) - exp_bias; if (exp_val < 0) { sc_val_from_ulong(0, NULL); _save_result(_exp(result)); _save_result(_mant(result)); - _desc(result).class = ZERO; + _desc(result).clss = ZERO; return result; } if (exp_val > _desc(a).mantissa_size) { if (a != result) - memcpy(result, a, CALC_BUFFER_SIZE); + memcpy(result, a, calc_buffer_size); return result; } @@ -782,7 +804,7 @@ static char* _trunc(const char *a, char *result) /* and the mask and return the result */ sc_and(_mant(a), temp, _mant(result)); - if (a != result) memcpy(_exp(result), _exp(a), VALUE_SIZE); + if (a != result) memcpy(_exp(result), _exp(a), value_size); return result; } @@ -790,7 +812,8 @@ static char* _trunc(const char *a, char *result) /* * This does value sanity checking(or should do it), sets up any prerequisites, * calls the proper internal functions, clears up and returns - * the result */ + * the result + */ char* _calc(const char *a, const char *b, int opcode, char *result) { char *temp; @@ -802,48 +825,51 @@ char* _calc(const char *a, const char *b, int opcode, char *result) if (result == NULL) result = calc_buffer; - TRACEPRINTF(("%s ", fc_print(a, buffer, 100, FC_HEX))); + TRACEPRINTF(("%s ", fc_print(a, buffer, 100, FC_PACKED))); switch (opcode) { case FC_add: /* make the value with the bigger exponent the first one */ - TRACEPRINTF(("+ %s ", fc_print(b, buffer, 100, FC_HEX))); + TRACEPRINTF(("+ %s ", fc_print(b, buffer, 100, FC_PACKED))); if (sc_comp(_exp(a), _exp(b)) == -1) - _add(b, a, result); + _fadd(b, a, result); else - _add(a, b, result); + _fadd(a, b, result); break; case FC_sub: - TRACEPRINTF(("- %s ", fc_print(b, buffer, 100, FC_HEX))); - temp = alloca(CALC_BUFFER_SIZE); - memcpy(temp, b, CALC_BUFFER_SIZE); + TRACEPRINTF(("- %s ", fc_print(b, buffer, 100, FC_PACKED))); + temp = alloca(calc_buffer_size); + memcpy(temp, b, calc_buffer_size); _sign(temp) = !_sign(b); if (sc_comp(_exp(a), _exp(temp)) == -1) - _add(temp, a, result); + _fadd(temp, a, result); else - _add(a, temp, result); + _fadd(a, temp, result); break; case FC_mul: - TRACEPRINTF(("* %s ", fc_print(b, buffer, 100, FC_HEX))); - _mul(a, b, result); + TRACEPRINTF(("* %s ", fc_print(b, buffer, 100, FC_PACKED))); + _fmul(a, b, result); break; case FC_div: - TRACEPRINTF(("/ %s ", fc_print(b, buffer, 100, FC_HEX))); - _div(a, b, result); + TRACEPRINTF(("/ %s ", fc_print(b, buffer, 100, FC_PACKED))); + _fdiv(a, b, result); break; case FC_neg: TRACEPRINTF(("negated ")); - if (a != result) memcpy(result, a, CALC_BUFFER_SIZE); + if (a != result) memcpy(result, a, calc_buffer_size); _sign(result) = !_sign(a); break; case FC_int: + TRACEPRINTF(("truncated to integer ")); _trunc(a, result); break; case FC_rnd: + TRACEPRINTF(("rounded to integer ")); + assert(0); break; } - TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_HEX))); + TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_PACKED))); return result; } @@ -855,9 +881,9 @@ const void *fc_get_buffer(void) return calc_buffer; } -const int fc_get_buffer_length(void) +int fc_get_buffer_length(void) { - return CALC_BUFFER_SIZE; + return calc_buffer_size; } char* fc_val_from_str(const char *str, unsigned int len, char exp_size, char mant_size, char *result) @@ -882,13 +908,13 @@ char* fc_val_from_str(const char *str, unsigned int len, char exp_size, char man if (result == NULL) result = calc_buffer; - exp_val = alloca(VALUE_SIZE); - power_val = alloca(CALC_BUFFER_SIZE); + exp_val = alloca(value_size); + power_val = alloca(calc_buffer_size); mant_str = alloca((len)?(len):(strlen(str))); _desc(result).exponent_size = exp_size; _desc(result).mantissa_size = mant_size; - _desc(result).class = NORMAL; + _desc(result).clss = NORMAL; old_str = str; pos = 0; @@ -1017,7 +1043,7 @@ char* fc_val_from_str(const char *str, unsigned int len, char exp_size, char man _fail_char(old_str, len, str - old_str); } } - } // switch(state) + } /* switch(state) */ done: sc_val_from_str(mant_str, strlen(mant_str), _mant(result)); @@ -1037,7 +1063,7 @@ done: _power_of_ten(exp_int, &_desc(result), power_val); - _div(result, power_val, result); + _fdiv(result, power_val, result); return result; #else @@ -1058,39 +1084,39 @@ done: char* fc_val_from_float(LLDBL l, char exp_size, char mant_size, char* result) { char *temp; - int bias_res = ((1<> 16; - UINT32 mantissa1 = srcval.val.mid; - UINT32 mantissa0 = srcval.val.low; + mant_val = 64; + bias_val = 0x3fff; + sign = (srcval.val.high & 0x00008000) != 0; + exponent = (srcval.val.high & 0x00007FFF) ; + mantissa0 = srcval.val.mid; + mantissa1 = srcval.val.low; #else /* no long double */ - mant_val = 52; - bias_val = 0x3ff; - UINT32 sign = (srcval.val.high & 0x80000000) != 0; - UINT32 exponent = (srcval.val.high & 0x7FF00000) >> 20; - UINT32 mantissa0 = srcval.val.high & 0x000FFFFF; - UINT32 mantissa1 = srcval.val.low; + mant_val = 52; + bias_val = 0x3ff; + sign = (srcval.val.high & 0x80000000) != 0; + exponent = (srcval.val.high & 0x7FF00000) >> 20; + mantissa0 = srcval.val.high & 0x000FFFFF; + mantissa1 = srcval.val.low; #endif #ifdef HAVE_LONG_DOUBLE - TRACEPRINTF(("val_from_float(%.8X%.8X%.8X)\n", ((int*)&l)[2], ((int*)&l)[1], ((int*)&l)[0])); + TRACEPRINTF(("val_from_float(%.8X%.8X%.8X)\n", ((int*)&l)[2], ((int*)&l)[1], ((int*)&l)[0]));/* srcval.val.high, srcval.val.mid, srcval.val.low)); */ DEBUGPRINTF(("(%d-%.4X-%.8X%.8X)\n", sign, exponent, mantissa0, mantissa1)); #else - TRACEPRINTF(("val_from_float(%.8X%.8X)\n", ((int*)&l)[1], ((int*)&l)[0])); + TRACEPRINTF(("val_from_float(%.8X%.8X)\n", srcval.val.high, srcval.val.low)); DEBUGPRINTF(("(%d-%.3X-%.5X%.8X)\n", sign, exponent, mantissa0, mantissa1)); #endif if (result == NULL) result = calc_buffer; - temp = alloca(VALUE_SIZE); + temp = alloca(value_size); _desc(result).exponent_size = exp_size; _desc(result).mantissa_size = mant_size; @@ -1101,20 +1127,28 @@ char* fc_val_from_float(LLDBL l, char exp_size, char mant_size, char* result) /* sign and flag suffice to identify nan or inf, no exponent/mantissa * encoding is needed. the function can return immediately in these cases */ if (isnan(l)) { - _desc(result).class = NAN; + _desc(result).clss = NAN; TRACEPRINTF(("val_from_float resulted in NAN\n")); return result; } else if (isinf(l)) { - _desc(result).class = INF; + _desc(result).clss = INF; TRACEPRINTF(("val_from_float resulted in %sINF\n", (_sign(result)==1)?"-":"")); return result; } - /* extract exponent */ + /* build exponent, because input and output exponent and mantissa sizes may differ + * this looks more complicated than it is: unbiased input exponent + output bias, + * minus the mantissa difference which is added again later when the output float + * becomes normalized */ +#ifdef HAVE_EXPLICIT_ONE + sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size-1), _exp(result)); +#else sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size), _exp(result)); +#endif - /* extract mantissa */ + /* build mantissa representation */ +#ifndef HAVE_EXPLICIT_ONE if (exponent != 0) { /* insert the hidden bit */ @@ -1123,6 +1157,7 @@ char* fc_val_from_float(LLDBL l, char exp_size, char mant_size, char* result) _shift_left(temp, sc_get_buffer(), NULL); } else +#endif { sc_val_from_ulong(0, NULL); } @@ -1151,7 +1186,7 @@ char* fc_val_from_float(LLDBL l, char exp_size, char mant_size, char* result) _normalize(result, result, 0); - TRACEPRINTF(("val_from_float results in %s\n", fc_print(result, temp, CALC_BUFFER_SIZE, FC_HEX))); + TRACEPRINTF(("val_from_float results in %s\n", fc_print(result, temp, calc_buffer_size, FC_PACKED))); return result; } @@ -1160,7 +1195,6 @@ LLDBL fc_val_to_float(const void *val) { const char *value; char *temp = NULL; - char *pack = NULL; int byte_offset; @@ -1172,43 +1206,45 @@ LLDBL fc_val_to_float(const void *val) value_t buildval; #ifdef HAVE_LONG_DOUBLE - char result_mantissa = 64; char result_exponent = 15; + char result_mantissa = 64; #else - char result_mantissa = 52; char result_exponent = 11; + char result_mantissa = 52; #endif - temp = alloca(CALC_BUFFER_SIZE); - pack = alloca(VALUE_SIZE); + temp = alloca(calc_buffer_size); +#ifdef HAVE_EXPLICIT_ONE + value = fc_cast(val, result_exponent, result_mantissa-1, temp); +#else value = fc_cast(val, result_exponent, result_mantissa, temp); +#endif sign = _sign(value); - /* long double exponent is 15bit, so the use of sc_val_to_long should not + + /* @@@ long double exponent is 15bit, so the use of sc_val_to_long should not * lead to wrong results */ exponent = sc_val_to_long(_exp(value)) ; - _pack(value, pack); + sc_val_from_ulong(2, NULL); + _shift_right(_mant(value), sc_get_buffer(), _mant(value)); mantissa0 = 0; mantissa1 = 0; for (byte_offset = 0; byte_offset < 4; byte_offset++) - mantissa1 |= sc_sub_bits(pack, result_mantissa, byte_offset) << (byte_offset<<3); + mantissa1 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << (byte_offset<<3); for (; (byte_offset<<3) < result_mantissa; byte_offset++) - mantissa0 |= sc_sub_bits(pack, result_mantissa, byte_offset) << ((byte_offset-4)<<3); - -#ifndef HAVE_LONG_DOUBLE - mantissa0 &= 0x000FFFFF; -#endif + mantissa0 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << ((byte_offset-4)<<3); #ifdef HAVE_LONG_DOUBLE - buildval.val.high = sign << 31; - buildval.val.high |= exponent << 16; - buildval.val.mid = mantissa1; - buildval.val.low = mantissa0; + buildval.val.high = sign << 15; + buildval.val.high |= exponent; + buildval.val.mid = mantissa0; + buildval.val.low = mantissa1; #else /* no long double */ + mantissa0 &= 0x000FFFFF; /* get rid of garbage */ buildval.val.high = sign << 31; buildval.val.high |= exponent << 20; buildval.val.high |= mantissa0; @@ -1216,7 +1252,6 @@ LLDBL fc_val_to_float(const void *val) #endif TRACEPRINTF(("val_to_float: %d-%x-%x%x\n", sign, exponent, mantissa0, mantissa1)); - //printf("val_to_float: %d-%x-%x%x\n", sign, exponent, mantissa0, mantissa1); return buildval.d; } @@ -1227,18 +1262,18 @@ char* fc_cast(const void *val, char exp_size, char mant_size, char *result) int exp_offset, val_bias, res_bias; if (result == NULL) result = calc_buffer; - temp = alloca(CALC_BUFFER_SIZE); + temp = alloca(value_size); if (_desc(value).exponent_size == exp_size && _desc(value).mantissa_size == mant_size) { - if (value != result) memcpy(result, value, CALC_BUFFER_SIZE); + if (value != result) memcpy(result, value, calc_buffer_size); return result; } /* set the descriptor of the new value */ _desc(result).exponent_size = exp_size; _desc(result).mantissa_size = mant_size; - _desc(result).class = _desc(value).class; + _desc(result).clss = _desc(value).clss; _sign(result) = _sign(value); @@ -1247,15 +1282,24 @@ char* fc_cast(const void *val, char exp_size, char mant_size, char *result) * offset and add it */ val_bias = (1<<_desc(value).exponent_size)/2-1; res_bias = (1<+< 1 because of two extra rounding bits */ sc_val_from_ulong(mantissa_size + 1, NULL); @@ -1332,7 +1376,7 @@ char* fc_get_plusinf(unsigned int exponent_size, unsigned int mantissa_size, cha _desc(result).exponent_size = exponent_size; _desc(result).mantissa_size = mantissa_size; - _desc(result).class = NORMAL; + _desc(result).clss = NORMAL; _sign(result) = 0; @@ -1357,26 +1401,47 @@ int fc_comp(const void *a, const void *b) { const char *val_a = (const char*)a; const char *val_b = (const char*)b; + int mul = 1; + + /* + * shortcut: if both values are identical, they are either + * Unordered if NaN or equal + */ + if (a == b) + return _desc(val_a).clss == NAN ? 2 : 0; + + /* unordered if one is a NaN */ + if (_desc(val_a).clss == NAN || _desc(val_b).clss == NAN) + return 2; - /* unordered */ - if (_desc(val_a).class == NAN || _desc(val_b).class == NAN) return 2; /* zero is equal independent of sign */ - if ((_desc(val_a).class == ZERO) && (_desc(val_b).class == ZERO)) return 0; + if ((_desc(val_a).clss == ZERO) && (_desc(val_b).clss == ZERO)) + return 0; + /* different signs make compare easy */ - if (_sign(val_a) != _sign(val_b)) return (_sign(val_a)==0)?(1):(-1); + if (_sign(val_a) != _sign(val_b)) + return (_sign(val_a)==0)?(1):(-1); + + mul = _sign(a) ? -1 : 1; + /* both infinity means equality */ - if ((_desc(val_a).class == INF) && (_desc(val_b).class == INF)) return 0; + if ((_desc(val_a).clss == INF) && (_desc(val_b).clss == INF)) + return 0; + /* infinity is bigger than the rest */ - if (_desc(val_a).class == INF) return _sign(val_a)?(-1):(1); - if (_desc(val_b).class == INF) return _sign(val_b)?(1):(-1); + if (_desc(val_a).clss == INF) + return 1 * mul; + if (_desc(val_b).clss == INF) + return -1 * mul; + /* check first exponent, that mantissa if equal */ switch (sc_comp(_exp(val_a), _exp(val_b))) { case -1: - return -1; + return -1 * mul; case 1: - return 1; + return 1 * mul; case 0: - return sc_comp(_mant(val_a), _mant(val_b)); + return sc_comp(_mant(val_a), _mant(val_b)) * mul; default: return 2; } @@ -1384,27 +1449,27 @@ int fc_comp(const void *a, const void *b) int fc_is_zero(const void *a) { - return _desc((const char*)a).class == ZERO; + return _desc(a).clss == ZERO; } int fc_is_negative(const void *a) { - return _sign((const char*)a); + return _sign(a); } int fc_is_inf(const void *a) { - return _desc(a).class == INF; + return _desc(a).clss == INF; } int fc_is_nan(const void *a) { - return _desc(a).class == NAN; + return _desc(a).clss == NAN; } int fc_is_subnormal(const void *a) { - return _desc(a).class == SUBNORMAL; + return _desc(a).clss == SUBNORMAL; } char *fc_print(const void *a, char *buf, int buflen, unsigned base) @@ -1414,11 +1479,11 @@ char *fc_print(const void *a, char *buf, int buflen, unsigned base) val = (const char*)a; - mul_1 = alloca(CALC_BUFFER_SIZE); + mul_1 = alloca(calc_buffer_size); switch (base) { case FC_DEC: - switch (_desc(val).class) { + switch (_desc(val).clss) { case INF: if (buflen >= 8+_sign(val)) sprintf(buf, "%sINFINITY", _sign(val)?"-":""); else snprintf(buf, buflen, "%sINF", _sign(val)?"-":NULL); @@ -1439,8 +1504,31 @@ char *fc_print(const void *a, char *buf, int buflen, unsigned base) #endif } break; + case FC_HEX: - snprintf(buf, buflen, "%s", sc_print(_pack(val, mul_1), 0, SC_HEX)); + switch (_desc(val).clss) { + case INF: + if (buflen >= 8+_sign(val)) sprintf(buf, "%sINFINITY", _sign(val)?"-":""); + else snprintf(buf, buflen, "%sINF", _sign(val)?"-":NULL); + break; + case NAN: + snprintf(buf, buflen, "NAN"); + break; + case ZERO: + snprintf(buf, buflen, "0.0"); + break; + default: +#ifdef HAVE_LONG_DOUBLE + snprintf(buf, buflen, "%LA", fc_val_to_float(val)); +#else + snprintf(buf, buflen, "%A", fc_val_to_float(val)); +#endif + } + break; + + case FC_PACKED: + default: + snprintf(buf, buflen, "%s", sc_print(_pack(val, mul_1), value_size*4, SC_HEX)); break; } return buf; @@ -1451,7 +1539,7 @@ unsigned char fc_sub_bits(const void *value, unsigned num_bits, unsigned byte_of /* this is used to cache the packed version of the value */ static char *pack = NULL; - if (pack == NULL) pack = malloc(VALUE_SIZE); + if (pack == NULL) pack = xmalloc(value_size); if (value != NULL) _pack((const char*)value, pack); @@ -1462,14 +1550,14 @@ unsigned char fc_sub_bits(const void *value, unsigned num_bits, unsigned byte_of fc_rounding_mode_t fc_set_rounding_mode(fc_rounding_mode_t mode) { if (mode == FC_TONEAREST || mode == FC_TOPOSITIVE || mode == FC_TONEGATIVE || mode == FC_TOZERO) - ROUNDING_MODE = mode; + rounding_mode = mode; - return ROUNDING_MODE; + return rounding_mode; } fc_rounding_mode_t fc_get_rounding_mode(void) { - return ROUNDING_MODE; + return rounding_mode; } void init_fltcalc(int precision) @@ -1486,19 +1574,33 @@ void init_fltcalc(int precision) if (max_precision < precision) printf("WARING: not enough precision available, using %d\n", max_precision); - ROUNDING_MODE = FC_TONEAREST; - VALUE_SIZE = sc_get_buffer_length(); - SIGN_POS = 0; - EXPONENT_POS = SIGN_POS + sizeof(char); - MANTISSA_POS = EXPONENT_POS + VALUE_SIZE; - DESCRIPTOR_POS = MANTISSA_POS + VALUE_SIZE; - CALC_BUFFER_SIZE = DESCRIPTOR_POS + sizeof(descriptor_t); + rounding_mode = FC_TONEAREST; + value_size = sc_get_buffer_length(); + SIGN_POS = 0; + EXPONENT_POS = SIGN_POS + sizeof(char); + MANTISSA_POS = EXPONENT_POS + value_size; + DESCRIPTOR_POS = MANTISSA_POS + value_size; + calc_buffer_size = DESCRIPTOR_POS + sizeof(descriptor_t); - calc_buffer = malloc(CALC_BUFFER_SIZE); - DEBUGPRINTF(("init fltcalc:\n\tVALUE_SIZE = %d\n\tSIGN_POS = %d\n\tEXPONENT_POS = %d\n\tMANTISSA_POS = %d\n\tDESCRIPTOR_POS = %d\n\tCALC_BUFFER_SIZE = %d\n\tcalc_buffer = %p\n\n", VALUE_SIZE, SIGN_POS, EXPONENT_POS, MANTISSA_POS, DESCRIPTOR_POS, CALC_BUFFER_SIZE, calc_buffer)); + calc_buffer = xmalloc(calc_buffer_size); + DEBUGPRINTF(("init fltcalc:\n\tVALUE_SIZE = %d\n\tSIGN_POS = %d\n\tEXPONENT_POS = %d\n\tMANTISSA_POS = %d\n\tDESCRIPTOR_POS = %d\n\tCALC_BUFFER_SIZE = %d\n\tcalc_buffer = %p\n\n", value_size, SIGN_POS, EXPONENT_POS, MANTISSA_POS, DESCRIPTOR_POS, calc_buffer_size, calc_buffer)); +#ifdef HAVE_LONG_DOUBLE + DEBUGPRINTF(("\tUsing long double (1-15-64) interface\n")); +#else + DEBUGPRINTF(("\tUsing double (1-11-52) interface\n")); +#endif +#ifdef WORDS_BIGENDIAN + DEBUGPRINTF(("\tWord order is big endian\n\n")); +#else + DEBUGPRINTF(("\tWord order is little endian\n\n")); +#endif } } +void finish_fltcalc (void) { + free(calc_buffer); calc_buffer = NULL; +} + /* definition of interface functions */ FC_DEFINE2(add) FC_DEFINE2(sub)