- fixed warnings
[libfirm] / ir / tv / fltcalc.c
index 4e498ea..821183d 100644 (file)
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 1995-2007 University of Karlsruhe.  All right reserved.
+ * Copyright (C) 1995-2008 University of Karlsruhe.  All right reserved.
  *
  * This file is part of libFirm.
  *
@@ -52,6 +52,9 @@
 
 #include "xmalloc.h"
 
+/** The number of extra precesion rounding bits */
+#define ROUNDING_BITS 2
+
 typedef uint32_t UINT32;
 
 #ifdef HAVE_LONG_DOUBLE
@@ -154,6 +157,9 @@ static int calc_buffer_size;
 static int value_size;
 static int max_precision;
 
+/** Exact flag. */
+static int fc_exact = 1;
+
 #if 0
 static void fail_char(const char *str, unsigned int len, int pos) {
        if (*(str+pos))
@@ -207,8 +213,8 @@ static void *pack(const fp_value *int_float, void *packed) {
        sc_or(temp, packed, packed);
 
        /* extract mantissa */
-       /* remove rounding bits */
-       sc_val_from_ulong(2, shift_val);
+       /* remove rounding bits */
+       sc_val_from_ulong(ROUNDING_BITS, shift_val);
        _shift_right(_mant(int_float), shift_val, temp);
 
        /* remove leading 1 (or 0 if denormalized) */
@@ -221,13 +227,19 @@ static void *pack(const fp_value *int_float, void *packed) {
        return packed;
 }
 
-static void normalize(const fp_value *in_val, fp_value *out_val, int sticky) {
+/**
+ * Normalize a fp_value.
+ *
+ * @return non-zero if result is exact
+ */
+static int normalize(const fp_value *in_val, fp_value *out_val, int sticky) {
+       int exact = 1;
        int hsb;
        char lsb, guard, round, round_dir = 0;
        char *temp = alloca(value_size);
 
-       /* +2: save two rounding bits at the end */
-       hsb = 2 + in_val->desc.mantissa_size - sc_get_highest_set_bit(_mant(in_val)) - 1;
+       /* save rounding bits at the end */
+       hsb = ROUNDING_BITS + in_val->desc.mantissa_size - sc_get_highest_set_bit(_mant(in_val)) - 1;
 
        if (in_val != out_val)   {
                out_val->sign = in_val->sign;
@@ -236,8 +248,8 @@ static void normalize(const fp_value *in_val, fp_value *out_val, int sticky) {
 
        out_val->desc.clss = NORMAL;
 
-       /* mantissa all zeros, so zero exponent (because of explicit one)*/
-       if (hsb == 2 + in_val->desc.mantissa_size)   {
+       /* mantissa all zeros, so zero exponent (because of explicit one) */
+       if (hsb == ROUNDING_BITS + in_val->desc.mantissa_size)   {
                sc_val_from_ulong(0, _exp(out_val));
                hsb = -1;
        }
@@ -250,8 +262,10 @@ static void normalize(const fp_value *in_val, fp_value *out_val, int sticky) {
                _shift_right(_mant(in_val), temp, _mant(out_val));
 
                /* remember if some bits were shifted away */
-               if (!sticky) sticky = sc_had_carry();
-
+               if (sc_had_carry()) {
+                       exact = 0;
+                       sticky = 1;
+               }
                sc_add(_exp(in_val), temp, _exp(out_val));
        } else if (hsb > -1) {
                /* shift left */
@@ -271,7 +285,10 @@ static void normalize(const fp_value *in_val, fp_value *out_val, int sticky) {
                sc_sub(temp, _exp(out_val), NULL);
 
                _shift_right(_mant(out_val), sc_get_buffer(), _mant(out_val));
-               if (!sticky) sticky = sc_had_carry();
+               if (sc_had_carry()) {
+                       exact  = 0;
+                       sticky = 1;
+               }
                /* denormalized means exponent of zero */
                sc_val_from_ulong(0, _exp(out_val));
 
@@ -281,7 +298,7 @@ static void normalize(const fp_value *in_val, fp_value *out_val, int sticky) {
        /* perform rounding by adding a value that clears the guard bit and the round bit
         * and either causes a carry to round up or not */
        /* get the last 3 bits of the value */
-       lsb = sc_sub_bits(_mant(out_val), out_val->desc.mantissa_size + 2, 0) & 0x7;
+       lsb = sc_sub_bits(_mant(out_val), out_val->desc.mantissa_size + ROUNDING_BITS, 0) & 0x7;
        guard = (lsb&0x2)>>1;
        round = lsb&0x1;
 
@@ -317,6 +334,7 @@ static void normalize(const fp_value *in_val, fp_value *out_val, int sticky) {
        if (lsb != 0) {
                sc_val_from_long(lsb, temp);
                sc_add(_mant(out_val), temp, _mant(out_val));
+               exact = 0;
        }
 
        /* could have rounded down to zero */
@@ -324,11 +342,12 @@ static void normalize(const fp_value *in_val, fp_value *out_val, int sticky) {
                out_val->desc.clss = ZERO;
 
        /* check for rounding overflow */
-       hsb = 2 + out_val->desc.mantissa_size - sc_get_highest_set_bit(_mant(out_val)) - 1;
+       hsb = ROUNDING_BITS + out_val->desc.mantissa_size - sc_get_highest_set_bit(_mant(out_val)) - 1;
        if ((out_val->desc.clss != SUBNORMAL) && (hsb < -1)) {
                sc_val_from_ulong(1, temp);
                _shift_right(_mant(out_val), temp, _mant(out_val));
-
+               if (exact && sc_had_carry())
+                       exact = 0;
                sc_add(_exp(out_val), temp, _exp(out_val));
        } else if ((out_val->desc.clss == SUBNORMAL) && (hsb == -1)) {
                /* overflow caused the mantissa to be normal again,
@@ -389,6 +408,7 @@ static void normalize(const fp_value *in_val, fp_value *out_val, int sticky) {
                        }
                }
        }
+       return exact;
 }
 
 /**
@@ -417,6 +437,8 @@ static void _fadd(const fp_value *a, const fp_value *b, fp_value *result) {
        char sign, res_sign;
        char sticky;
 
+       fc_exact = 1;
+
        handle_NAN(a, b, result);
 
        /* make sure result has a descriptor */
@@ -487,6 +509,7 @@ static void _fadd(const fp_value *a, const fp_value *b, fp_value *result) {
 
        _shift_right(_mant(b), exp_diff, temp);
        sticky = sc_had_carry();
+       fc_exact &= !sticky;
 
        if (sticky && sign) {
                /* if subtracting a little more than the represented value or adding a little
@@ -516,16 +539,19 @@ static void _fadd(const fp_value *a, const fp_value *b, fp_value *result) {
        /* resulting exponent is the bigger one */
        memmove(_exp(result), _exp(a), value_size);
 
-       normalize(result, result, sticky);
+       fc_exact &= normalize(result, result, sticky);
 }
 
 /**
  * calculate a * b
  */
 static void _fmul(const fp_value *a, const fp_value *b, fp_value *result) {
+       int sticky;
        char *temp;
        char res_sign;
 
+       fc_exact = 1;
+
        handle_NAN(a, b, result);
 
        temp = alloca(value_size);
@@ -588,21 +614,26 @@ static void _fmul(const fp_value *a, const fp_value *b, fp_value *result) {
         * point are the sum of the factors' digits after the radix point. As all
         * values are normalized they both have the same amount of these digits,
         * which has to be restored by proper shifting
-        * +2 because of the two rounding bits */
-       sc_val_from_ulong(2 + result->desc.mantissa_size, temp);
+        * because of the rounding bits */
+       sc_val_from_ulong(ROUNDING_BITS + result->desc.mantissa_size, temp);
 
        _shift_right(_mant(result), temp, _mant(result));
+       sticky = sc_had_carry();
+       fc_exact &= !sticky;
 
-       normalize(result, result, sc_had_carry());
+       fc_exact &= normalize(result, result, sticky);
 }
 
 /**
  * calculate a / b
  */
 static void _fdiv(const fp_value *a, const fp_value *b, fp_value *result) {
+       int sticky;
        char *temp, *dividend;
        char res_sign;
 
+       fc_exact = 1;
+
        handle_NAN(a, b, result);
 
        temp = alloca(value_size);
@@ -674,7 +705,7 @@ static void _fdiv(const fp_value *a, const fp_value *b, fp_value *result) {
         * fit into the integer precision, but due to the rounding bits (which
         * are always zero because the values are all normalized) the divisor
         * can be shifted right instead to achieve the same result */
-       sc_val_from_ulong(2 + result->desc.mantissa_size, temp);
+       sc_val_from_ulong(ROUNDING_BITS + result->desc.mantissa_size, temp);
 
        _shift_left(_mant(a), temp, dividend);
 
@@ -683,9 +714,11 @@ static void _fdiv(const fp_value *a, const fp_value *b, fp_value *result) {
                sc_val_from_ulong(1, divisor);
                _shift_right(_mant(b), divisor, divisor);
                sc_div(dividend, divisor, _mant(result));
+               sticky = sc_had_carry();
+               fc_exact &= !sticky;
        }
 
-       normalize(result, result, sc_had_carry());
+       fc_exact &= normalize(result, result, sticky);
 }
 
 #if 0
@@ -703,7 +736,7 @@ static void _power_of_ten(int exp, descriptor_t *desc, char *result) {
        build = alloca(value_size);
        temp = alloca(value_size);
 
-       sc_val_from_ulong((1 << result->desc.exponent_size)/2-1, _exp(result));
+       sc_val_from_ulong((1 << (result->desc.exponent_size - 1)) - 1, _exp(result));
 
        if (exp > 0) {
                /* temp is value of ten now */
@@ -717,7 +750,7 @@ static void _power_of_ten(int exp, descriptor_t *desc, char *result) {
                _save_result(build);
 
                /* temp is amount of left shift needed to put the value left of the radix point */
-               sc_val_from_ulong(result->desc.mantissa_size + 2, temp);
+               sc_val_from_ulong(result->desc.mantissa_size + ROUNDING_BITS, temp);
 
                _shift_left(build, temp, _mant(result));
 
@@ -747,6 +780,9 @@ static void _trunc(const fp_value *a, fp_value *result) {
        int exp_bias, exp_val;
        char *temp;
 
+       /* fixme: can be exact */
+       fc_exact = 0;
+
        temp = alloca(value_size);
 
        if (a != result)
@@ -958,11 +994,11 @@ done:
        sc_val_from_str(mant_str, strlen(mant_str), _mant(result));
 
        /* shift to put value left of radix point */
-       sc_val_from_ulong(mant_size + 2, exp_val);
+       sc_val_from_ulong(mant_size + ROUNDING_BITS, exp_val);
 
        _shift_left(_mant(result), exp_val, _mant(result));
 
-       sc_val_from_ulong((1 << exp_size)/2-1, _exp(result));
+       sc_val_from_ulong((1 << (exp_size - 1)) - 1, _exp(result));
 
        _normalize(result, result, 0);
 
@@ -1001,7 +1037,7 @@ fp_value *fc_val_from_ieee754(LLDBL l, char exp_size, char mant_size, fp_value *
        UINT32 sign, exponent, mantissa0, mantissa1;
 
        srcval.d = l;
-       bias_res = ((1<<exp_size)/2-1);
+       bias_res = ((1 << (exp_size - 1)) - 1);
 
 #ifdef HAVE_LONG_DOUBLE
        mant_val  = 64;
@@ -1030,6 +1066,9 @@ fp_value *fc_val_from_ieee754(LLDBL l, char exp_size, char mant_size, fp_value *
        if (result == NULL) result = calc_buffer;
        temp = alloca(value_size);
 
+       /* CLEAR the buffer, else some bits might be uninitialised */
+       memset(result, 0, fc_get_buffer_length());
+
        result->desc.exponent_size = exp_size;
        result->desc.mantissa_size = mant_size;
 
@@ -1064,7 +1103,7 @@ fp_value *fc_val_from_ieee754(LLDBL l, char exp_size, char mant_size, fp_value *
        if (exponent != 0) {
                /* insert the hidden bit */
                sc_val_from_ulong(1, temp);
-               sc_val_from_ulong(mant_val + 2, NULL);
+               sc_val_from_ulong(mant_val + ROUNDING_BITS, NULL);
                _shift_left(temp, sc_get_buffer(), NULL);
        }
        else
@@ -1083,7 +1122,7 @@ fp_value *fc_val_from_ieee754(LLDBL l, char exp_size, char mant_size, fp_value *
 
        /* bits from the lower word */
        sc_val_from_ulong(mantissa1, temp);
-       sc_val_from_ulong(2, NULL);
+       sc_val_from_ulong(ROUNDING_BITS, NULL);
        _shift_left(temp, sc_get_buffer(), temp);
        sc_or(_mant(result), temp, _mant(result));
 
@@ -1135,7 +1174,7 @@ LLDBL fc_val_to_ieee754(const fp_value *val) {
         * lead to wrong results */
        exponent = sc_val_to_long(_exp(value)) ;
 
-       sc_val_from_ulong(2, NULL);
+       sc_val_from_ulong(ROUNDING_BITS, NULL);
        _shift_right(_mant(value), sc_get_buffer(), _mant(value));
 
        mantissa0 = 0;
@@ -1177,6 +1216,13 @@ fp_value *fc_cast(const fp_value *value, char exp_size, char mant_size, fp_value
                return result;
        }
 
+       if (value->desc.clss == NAN) {
+               if (sc_get_highest_set_bit(_mant(value)) == value->desc.mantissa_size + 1)
+                       return fc_get_qnan(exp_size, mant_size, result);
+               else
+                       return fc_get_snan(exp_size, mant_size, result);
+       }
+
        /* set the descriptor of the new value */
        result->desc.exponent_size = exp_size;
        result->desc.mantissa_size = mant_size;
@@ -1221,7 +1267,7 @@ fp_value *fc_get_max(unsigned int exponent_size, unsigned int mantissa_size, fp_
        sc_val_from_ulong((1<<exponent_size) - 2, _exp(result));
 
        sc_max_from_bits(mantissa_size + 1, 0, _mant(result));
-       sc_val_from_ulong(2, NULL);
+       sc_val_from_ulong(ROUNDING_BITS, NULL);
        _shift_left(_mant(result), sc_get_buffer(), _mant(result));
 
        return result;
@@ -1436,15 +1482,43 @@ unsigned char fc_sub_bits(const fp_value *value, unsigned num_bits, unsigned byt
        return sc_sub_bits(packed_value, num_bits, byte_ofs);
 }
 
+/* Returns non-zero if the mantissa is zero, i.e. 1.0Exxx */
 int fc_zero_mantissa(const fp_value *value) {
-       return sc_get_lowest_set_bit(_mant(value)) == 2 + value->desc.mantissa_size;
+       return sc_get_lowest_set_bit(_mant(value)) == ROUNDING_BITS + value->desc.mantissa_size;
 }
 
+/* Returns the exponent of a value. */
 int fc_get_exponent(const fp_value *value) {
        int exp_bias = (1 << (value->desc.exponent_size - 1)) - 1;
        return sc_val_to_long(_exp(value)) - exp_bias;
 }
 
+/* Return non-zero if a given value can be converted lossless into another precision */
+int fc_can_lossless_conv_to(const fp_value *value, char exp_size, char mant_size) {
+       int v;
+       int exp_bias;
+
+       /* handle some special cases first */
+       switch (value->desc.clss) {
+       case ZERO:
+       case INF:
+       case NAN:
+               return 1;
+       default:
+               break;
+       }
+
+       /* check if the exponent can be encoded: note, 0 and all ones are reserved for the exponent */
+       exp_bias = (1 << (exp_size - 1)) - 1;
+       v = fc_get_exponent(value) + exp_bias;
+       if (0 < v && v < (1 << exp_size) - 1) {
+               /* check the mantissa */
+               v = value->desc.mantissa_size + ROUNDING_BITS - sc_get_lowest_set_bit(_mant(value));
+               return v < mant_size;
+       }
+       return 0;
+}
+
 
 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)
@@ -1494,7 +1568,9 @@ void finish_fltcalc (void) {
        free(calc_buffer); calc_buffer = NULL;
 }
 
+#ifdef FLTCALC_TRACE_CALC
 static char buffer[100];
+#endif
 
 /* definition of interface functions */
 fp_value *fc_add(const fp_value *a, const fp_value *b, fp_value *result) {
@@ -1585,6 +1661,7 @@ fp_value *fc_int(const fp_value *a, fp_value *result) {
 fp_value *fc_rnd(const fp_value *a, fp_value *result) {
        if (result == NULL) result = calc_buffer;
 
+       (void) a;
        TRACEPRINTF(("%s ", fc_print(a, buffer, sizeof(buffer), FC_PACKED)));
        TRACEPRINTF(("rounded to integer "));
 
@@ -1594,9 +1671,72 @@ fp_value *fc_rnd(const fp_value *a, fp_value *result) {
        return result;
 }
 
+/*
+ * convert a floating point value into an sc value ...
+ */
+int fc_flt2int(const fp_value *a, void *result, ir_mode *dst_mode) {
+       if (a->desc.clss == NORMAL) {
+               int exp_bias = (1 << (a->desc.exponent_size - 1)) - 1;
+               int exp_val  = sc_val_to_long(_exp(a)) - exp_bias;
+               int shift, highest;
+
+               if (a->sign && !mode_is_signed(dst_mode)) {
+                       /* FIXME: for now we cannot convert this */
+                       return 0;
+               }
+
+               assert(exp_val >= 0 && "floating point value not integral before fc_flt2int() call");
+               shift = exp_val - (a->desc.mantissa_size + ROUNDING_BITS);
+
+               if (shift > 0) {
+                       sc_shlI(_mant(a),  shift, 64, 0, result);
+               } else {
+                       sc_shrI(_mant(a), -shift, 64, 0, result);
+               }
+
+               /* check for overflow */
+               highest = sc_get_highest_set_bit(result);
+
+               if (mode_is_signed(dst_mode)) {
+                       if (highest == sc_get_lowest_set_bit(result)) {
+                               /* need extra test for MIN_INT */
+                               if (highest >= (int) get_mode_size_bits(dst_mode)) {
+                                       /* FIXME: handle overflow */
+                                       return 0;
+                               }
+                       } else {
+                               if (highest >= (int) get_mode_size_bits(dst_mode) - 1) {
+                                       /* FIXME: handle overflow */
+                                       return 0;
+                               }
+                       }
+               } else {
+                       if (highest >= (int) get_mode_size_bits(dst_mode)) {
+                               /* FIXME: handle overflow */
+                               return 0;
+                       }
+               }
+
+               if (a->sign)
+                       sc_neg(result, result);
+
+               return 1;
+       }
+       else if (a->desc.clss == ZERO) {
+               sc_zero(result);
+               return 1;
+       }
+       return 0;
+}
+
+
 unsigned fc_set_immediate_precision(unsigned bits) {
        unsigned old = immediate_prec;
 
        immediate_prec = bits;
        return old;
 }
+
+int fc_is_exact(void) {
+       return fc_exact;
+}