fltcalc.c

/*
 * Copyright (C) 1995-2007 University of Karlsruhe.  All right reserved.
 *
 * This file is part of libFirm.
 *
 * This file may be distributed and/or modified under the terms of the
 * GNU General Public License version 2 as published by the Free Software
 * Foundation and appearing in the file LICENSE.GPL included in the
 * packaging of this file.
 *
 * Licensees holding valid libFirm Professional Edition licenses may use
 * this file in accordance with the libFirm Commercial License.
 * Agreement provided with the Software.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE.
 */

/*
 * Project:     libFIRM
 * File name:   ir/tv/fltcalc.c
 * Purpose:
 * Author:
 * Modified by:
 * Created:     2003
 * CVS-ID:      $Id$
 * Copyright:   (c) 2003 Universitt Karlsruhe
 */


#ifdef HAVE_CONFIG_H
# include "config.h"
#endif


#include "fltcalc.h"
#include "strcalc.h"

#include <math.h>    /* need isnan() and isinf() (will be changed)*/
/* undef some reused constants defined by math.h */
#ifdef NAN
#  undef NAN
#endif

#ifdef HAVE_INTTYPES_H
# include <inttypes.h>
#endif
#ifdef HAVE_STRING_H
# include <string.h>
#endif
#ifdef HAVE_STDLIB_H
# include <stdlib.h>
#endif
#include <stdio.h>
#include <assert.h>

#include "xmalloc.h"

typedef uint32_t UINT32;

#ifdef HAVE_LONG_DOUBLE
#ifdef WORDS_BIGENDIAN
typedef union {
  struct {
    UINT32 high;
    UINT32 mid;
    UINT32 low;
  } val;
  volatile long double d;
} value_t;
#else
typedef union {
  struct {
    UINT32 low;
    UINT32 mid;
    UINT32 high;
  } val;
  volatile long double d;
} value_t;
#endif
#else
#ifdef WORDS_BIGENDIAN
typedef union {
  struct {
    UINT32 high;
    UINT32 low;
  } val;
  volatile double d;
} value_t;
#else
typedef union {
  struct {
    UINT32 low;
    UINT32 high;
  } val;
  volatile double d;
} value_t;
#endif
#endif

/**
 * possible float states
 */
typedef enum {
  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 {
  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)

/* 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];
 *   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 _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_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

#if FLTCALC_DEBUG
#  define DEBUGPRINTF(x) printf x
#else
#  define DEBUGPRINTF(x) ((void)0)
#endif

#if FLTCALC_TRACE_CALC
#  define TRACEPRINTF(x) printf x
#else
#  define TRACEPRINTF(x) ((void)0)
#endif

static char *calc_buffer = NULL;

static fc_rounding_mode_t rounding_mode;

static int calc_buffer_size;
static int value_size;
static int SIGN_POS;
static int EXPONENT_POS;
static int MANTISSA_POS;
static int DESCRIPTOR_POS;

static int max_precision;
/********
 * private functions
 ********/
#if 0
static void _fail_char(const char *str, unsigned int len, int pos)
{
  if (*(str+pos))
    printf("ERROR: Unexpected character '%c'\n", *(str + pos));
  else
    printf("ERROR: Unexpected end of string\n");
  while (len-- && *str) printf("%c", *str++); printf("\n");
  while (pos--) printf(" "); printf("^\n");
  /* the front end has to to check constant strings */
  exit(-1);
}
#endif

/** 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);

  switch (_desc(int_float).clss) {
    case NAN:
      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);
      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;
      break;

    default:
      break;
  }
  /* pack sign */
  sc_val_from_ulong(_sign(int_float), temp);

  sc_val_from_ulong(_desc(int_float).exponent_size + _desc(int_float).mantissa_size, NULL);
  _shift_left(temp, sc_get_buffer(), packed);

  /* extract exponent */
  sc_val_from_ulong(_desc(int_float).mantissa_size, shift_val);

  _shift_left(_exp(int_float), shift_val, temp);

  sc_or(temp, packed, packed);

  /* extract mantissa */
  /* remove 2 rounding bits */
  sc_val_from_ulong(2, shift_val);
  _shift_right(_mant(int_float), shift_val, temp);

  /* remove leading 1 (or 0 if denormalized) */
  sc_max_from_bits(_desc(int_float).mantissa_size, 0, shift_val); /* all mantissa bits are 1's */
  sc_and(temp, shift_val, temp);

  /* save result */
  sc_or(temp, packed, packed);

  return packed;
}

char* _normalize(const char *in_val, char *out_val, int sticky)
{
  int hsb;
  char lsb, guard, round, round_dir = 0;
  char *temp;

  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;

  if (in_val != out_val)
  {
    _sign(out_val) = _sign(in_val);
    memcpy(&_desc(out_val), &_desc(in_val), sizeof(descriptor_t));
  }

  _desc(out_val).clss = NORMAL;

  /* 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 into the left of the radix point (i.e. hsb == -1) */
  if (hsb < -1)
  {
    /* shift right */
    sc_val_from_ulong(-hsb-1, temp);

    _shift_right(_mant(in_val), temp, _mant(out_val));

    /* remember if some bits were shifted away */
    if (!sticky) sticky = sc_had_carry();

    sc_add(_exp(in_val), temp, _exp(out_val));
  }
  else if (hsb > -1)
  {
    /* shift left */
    sc_val_from_ulong(hsb+1, temp);

    _shift_left(_mant(in_val), temp, _mant(out_val));

    sc_sub(_exp(in_val), temp, _exp(out_val));
  }

  /* check for exponent underflow */
  if (sc_is_negative(_exp(out_val)) || sc_is_zero(_exp(out_val))) {
    DEBUGPRINTF(("Exponent underflow!\n"));
    /* exponent underflow */
    /* shift the mantissa right to have a zero exponent */
    sc_val_from_ulong(1, temp);
    sc_sub(temp, _exp(out_val), NULL);

    _shift_right(_mant(out_val), sc_get_buffer(), _mant(out_val));
    if (!sticky) sticky = sc_had_carry();
    /* denormalized means exponent of zero */
    sc_val_from_ulong(0, _exp(out_val));

    _desc(out_val).clss = SUBNORMAL;
  }

  /* 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), _desc(out_val).mantissa_size + 2, 0) & 0x7;
  guard = (lsb&0x2)>>1;
  round = lsb&0x1;

  switch (rounding_mode)
  {
    case FC_TONEAREST:
      /* round to nearest representable value, if in doubt choose the version
       * with lsb == 0 */
      round_dir = guard && (sticky || round || lsb>>2);
      break;
    case FC_TOPOSITIVE:
      /* if positive: round to one if the exact value is bigger, else to zero */
      round_dir = (!_sign(out_val) && (guard || round || sticky));
      break;
    case FC_TONEGATIVE:
      /* if negative: round to one if the exact value is bigger, else to zero */
      round_dir = (_sign(out_val) && (guard || round || sticky));
      break;
    case FC_TOZERO:
      /* always round to 0 (chopping mode) */
      round_dir = 0;
      break;
  }
  DEBUGPRINTF(("Rounding (s%d, l%d, g%d, r%d, s%d) %s\n", _sign(out_val), lsb>>2, guard, round, sticky, (round_dir)?"up":"down"));

  if (round_dir == 1)
  {
    guard = (round^guard)<<1;
    lsb = !(round || guard)<<2 | guard | round;
  }
  else
  {
    lsb = -((guard<<1) | round);
  }

  /* add the rounded value */
  if (lsb != 0) {
    sc_val_from_long(lsb, temp);
    sc_add(_mant(out_val), temp, _mant(out_val));
  }

  /* could have rounded down to 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).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).clss == SUBNORMAL) && (hsb == -1))
  {
    /* 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).clss = NORMAL;
  }
  /* no further rounding is needed, because rounding overflow means
   * the carry of the original rounding was propagated all the way
   * up to the bit left of the radix point. This implies the bits
   * to the right are all zeros (rounding is +1) */

  /* check for exponent overflow */
  sc_val_from_ulong((1 << _desc(out_val).exponent_size) - 1, temp);
  if (sc_comp(_exp(out_val), temp) != -1) {
    DEBUGPRINTF(("Exponent overflow!\n"));
    /* exponent overflow, reaction depends on rounding method:
     *
     * mode        | sign of value |  result
     *--------------------------------------------------------------
     * TO_NEAREST  |      +        |   +inf
     *             |      -        |   -inf
     *--------------------------------------------------------------
     * TO_POSITIVE |      +        |   +inf
     *             |      -        |   smallest representable value
     *--------------------------------------------------------------
     * TO_NEAGTIVE |      +        |   largest representable value
     *             |      -        |   -inf
     *--------------------------------------------------------------
     * TO_ZERO     |      +        |   largest representable value
     *             |      -        |   smallest representable value
     *--------------------------------------------------------------*/
    if (_sign(out_val) == 0)
    {
      /* value is positive */
      switch (rounding_mode) {
        case FC_TONEAREST:
        case FC_TOPOSITIVE:
          _desc(out_val).clss = INF;
          break;

        case FC_TONEGATIVE:
        case FC_TOZERO:
          fc_get_max(_desc(out_val).exponent_size, _desc(out_val).mantissa_size, out_val);
      }
    } else {
      /* value is negative */
      switch (rounding_mode) {
        case FC_TONEAREST:
        case FC_TONEGATIVE:
          _desc(out_val).clss = INF;
          break;

        case FC_TOPOSITIVE:
        case FC_TOZERO:
          fc_get_min(_desc(out_val).exponent_size, _desc(out_val).mantissa_size, out_val);
      }
    }
  }

  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* _fadd(const char* a, const char* b, char* result)
{
  char *temp;
  char *exp_diff;

  char sign;
  char sticky;

  handle_NAN(a, b, result);

  /* make sure result has a descriptor */
  if (result != a && result != b)
    memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));

  /* determine if this is an addition or subtraction */
  sign = _sign(a) ^ _sign(b);

  /* 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);

  /* get exponent difference */
  sc_sub(_exp(a), _exp(b), exp_diff);

  /* initially set sign to be the sign of a, special treatment of subtraction
   * when exponents are equal is required though.
   * Also special care about the sign is needed when the mantissas are equal
   * (+/- 0 ?) */
  if (sign && sc_val_to_long(exp_diff) == 0) {
    switch (sc_comp(_mant(a), _mant(b))) {
      case 1:  /* a > b */
        _sign(result) = _sign(a);  /* abs(a) is bigger and a is negative */
        break;
      case 0:  /* a == b */
        _sign(result) = (rounding_mode == FC_TONEGATIVE);
        break;
      case -1: /* a < b */
        _sign(result) = _sign(b); /* abs(b) is bigger and b is negative */
        break;
      default:
        /* can't be reached */
        break;
    }
  }
  else
    _sign(result) = _sign(a);

  /* sign has been taken care of, check for special cases */
  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).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).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).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).clss == SUBNORMAL) && (_desc(a).clss != SUBNORMAL))
  {
    sc_val_from_ulong(1, temp);
    sc_sub(exp_diff, temp, exp_diff);
  }

  _shift_right(_mant(b), exp_diff, temp);
  sticky = sc_had_carry();

  if (sticky && sign)
  {
    /* 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);
    sc_val_from_ulong(1, temp1);
    sc_add(temp, temp1, temp);
  }

  if (sign) {
    if (sc_comp(_mant(a), temp) == -1)
      sc_sub(temp, _mant(a), _mant(result));
    else
      sc_sub(_mant(a), temp, _mant(result));
  } else {
    sc_add(_mant(a), temp, _mant(result));
  }

  /* _normalize expects a 'normal' radix point, adding two subnormals
   * results in a subnormal radix point -> shifting before normalizing */
  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);

  return _normalize(result, result, sticky);
}

/**
 * calculate a * b
 */
static char* _fmul(const char* a, const char* b, char* result)
{
  char *temp;

  handle_NAN(a, b, result);

  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).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);
    return result;
  }
  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);
    return result;
  }

  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).clss == INF) {
    if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-1);
    return result;
  }

  /* exp = exp(a) + exp(b) - excess */
  sc_add(_exp(a), _exp(b), _exp(result));

  sc_val_from_ulong((1<<_desc(a).exponent_size)/2-1, temp);
  sc_sub(_exp(result), temp, _exp(result));

  /* mixed normal, subnormal values introduce an error of 1, correct it */
  if ((_desc(a).clss == SUBNORMAL) ^ (_desc(b).clss == SUBNORMAL))
  {
    sc_val_from_ulong(1, temp);
    sc_add(_exp(result), temp, _exp(result));
  }

  sc_mul(_mant(a), _mant(b), _mant(result));

  /* realign result: after a multiplication the digits right of the radix
   * 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 + _desc(result).mantissa_size, temp);

  _shift_right(_mant(result), temp, _mant(result));

  return _normalize(result, result, sc_had_carry());
}

/**
 * calculate a / b
 */
static char* _fdiv(const char* a, const char* b, char* result)
{
  char *temp, *dividend;

  handle_NAN(a, b, result);

  temp = alloca(value_size);
  dividend = 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/0 and inf/inf */
  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);
    return result;
  }

  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 {
      /* x/inf -> 0 */
      sc_val_from_ulong(0, NULL);
      _save_result(_exp(result));
      _save_result(_mant(result));
      _desc(result).clss = ZERO;
    }
    return result;
  }

  if (_desc(a).clss == INF) {
    /* inf/x -> inf */
    if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
    return result;
  }
  if (_desc(b).clss == ZERO) {
    /* division by zero */
    if (_sign(result))
      fc_get_minusinf(_desc(a).exponent_size, _desc(a).mantissa_size, result);
    else
      fc_get_plusinf(_desc(a).exponent_size, _desc(a).mantissa_size, result);
    return result;
  }

  /* exp = exp(a) - exp(b) + excess - 1*/
  sc_sub(_exp(a), _exp(b), _exp(result));
  sc_val_from_ulong((1 << _desc(a).exponent_size)/2-2, temp);
  sc_add(_exp(result), temp, _exp(result));

  /* mixed normal, subnormal values introduce an error of 1, correct it */
  if ((_desc(a).clss == SUBNORMAL) ^ (_desc(b).clss == SUBNORMAL))
  {
    sc_val_from_ulong(1, temp);
    sc_add(_exp(result), temp, _exp(result));
  }

  /* mant(res) = mant(a) / 1/2mant(b) */
  /* to gain more bits of precision in the result the dividend could be
   * shifted left, as this operation does not loose bits. This would not
   * 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 + _desc(result).mantissa_size, temp);

  _shift_left(_mant(a), temp, dividend);

  {
    char *divisor = alloca(calc_buffer_size);
    sc_val_from_ulong(1, divisor);
    _shift_right(_mant(b), divisor, divisor);
    sc_div(dividend, divisor, _mant(result));
  }

  return _normalize(result, result, sc_had_carry());
}

#if 0
static void _power_of_ten(int exp, descriptor_t *desc, char *result)
{
  char *build;
  char *temp;

  /* positive sign */
  _sign(result) = 0;

  /* set new descriptor (else result is supposed to already have one) */
  if (desc != NULL)
    memcpy(&_desc(result), desc, sizeof(descriptor_t));

  build = alloca(value_size);
  temp = alloca(value_size);

  sc_val_from_ulong((1 << _desc(result).exponent_size)/2-1, _exp(result));

  if (exp > 0)
  {
    /* temp is value of ten now */
    sc_val_from_ulong(10, NULL);
    _save_result(temp);

    for (exp--; exp > 0; exp--) {
      _save_result(build);
      sc_mul(build, temp, NULL);
    }
    _save_result(build);

    /* 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));

    _normalize(result, result, 0);
  }
}
#endif

/**
 * 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
   * are the integral part of the value. For 15bit exp_size
   * 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 does not have fractional parts.
   */

  int exp_bias, exp_val;
  char *temp;

  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;

  if (exp_val < 0) {
    sc_val_from_ulong(0, NULL);
    _save_result(_exp(result));
    _save_result(_mant(result));
    _desc(result).clss = ZERO;

    return result;
  }

  if (exp_val > _desc(a).mantissa_size) {
    if (a != result)
      memcpy(result, a, calc_buffer_size);

    return result;
  }

  /* set up a proper mask to delete all bits right of the
   * radix point if the mantissa had been shifted until exp == 0 */
  sc_max_from_bits(1 + exp_val, 0, temp);
  sc_val_from_long(_desc(a).mantissa_size - exp_val + 2, NULL);
  _shift_left(temp, sc_get_buffer(), temp);

  /* and the mask and return the result */
  sc_and(_mant(a), temp, _mant(result));

  if (a != result) memcpy(_exp(result), _exp(a), value_size);

  return 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
 */
char* _calc(const char *a, const char *b, int opcode, char *result)
{
  char *temp;
#ifdef FLTCALC_TRACE_CALC
  char *buffer;

  buffer = alloca(100);
#endif

  if (result == NULL) result = calc_buffer;

  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_PACKED)));
      if (sc_comp(_exp(a), _exp(b)) == -1)
        _fadd(b, a, result);
      else
        _fadd(a, b, result);
      break;
    case FC_sub:
      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)
        _fadd(temp, a, result);
      else
        _fadd(a, temp, result);
      break;
    case FC_mul:
      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_PACKED)));
      _fdiv(a, b, result);
      break;
    case FC_neg:
      TRACEPRINTF(("negated "));
      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_PACKED)));
  return result;
}

/********
 * functions defined in fltcalc.h
 ********/
const void *fc_get_buffer(void)
{
  return calc_buffer;
}

int fc_get_buffer_length(void)
{
  return calc_buffer_size;
}

char* fc_val_from_str(const char *str, unsigned int len, char exp_size, char mant_size, char *result)
{
#if 0
  enum {
    START,
    LEFT_OF_DOT,
    RIGHT_OF_DOT,
    EXP_START,
    EXPONENT,
    END
  };

  char exp_sign;
  int exp_int, hsb, state;

  const char *old_str;

  int pos;
  char *mant_str, *exp_val, *power_val;

  if (result == NULL) result = calc_buffer;

  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).clss = NORMAL;

  old_str = str;
  pos = 0;
  exp_int = 0;
  state = START;

  while (len == 0 || str-old_str < len)
  {
    switch (state) {
      case START:
        switch (*str) {
          case '+':
            _sign(result) = 0;
            state = LEFT_OF_DOT;
            str++;
            break;

          case '-':
            _sign(result) = 1;
            state = LEFT_OF_DOT;
            str++;
            break;

          case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
            _sign(result) = 0;
            state = LEFT_OF_DOT;
            break;

          case '.':
            _sign(result) = 0;
            state = RIGHT_OF_DOT;
            str++;
            break;

          case 'n':
          case 'N':
          case 'i':
          case 'I':
            break;

          default:
            _fail_char(old_str, len, str - old_str);
        }
        break;

      case LEFT_OF_DOT:
        switch (*str) {
          case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
            mant_str[pos++] = *(str++);
            break;

          case '.':
            state = RIGHT_OF_DOT;
            str++;
            break;

          case 'e':
          case 'E':
            state = EXP_START;
            str++;
            break;

          case '\0':
            mant_str[pos] = '\0';
            goto done;

          default:
            _fail_char(old_str, len, str - old_str);
        }
        break;

      case RIGHT_OF_DOT:
        switch (*str) {
          case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
            mant_str[pos++] = *(str++);
            exp_int++;
            break;

          case 'e':
          case 'E':
            state = EXP_START;
            str++;
            break;

          case '\0':
            mant_str[pos] = '\0';
            goto done;

          default:
            _fail_char(old_str, len, str - old_str);
        }
        break;

      case EXP_START:
        switch (*str) {
          case '-':
            exp_sign = 1;
            /* fall through */
          case '+':
            if (*(str-1) != 'e' && *(str-1) != 'E') _fail_char(old_str, len, str - old_str);
            str++;
            break;

          case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
            mant_str[pos] = '\0';
            pos = 1;
            str++;
            state = EXPONENT;
            break;

          default:
            _fail_char(old_str, len, str - old_str);
        }
        break;

      case EXPONENT:
        switch (*str) {
          case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
            pos++;
            str++;
            break;

          case '\0': goto done;

          default:
            _fail_char(old_str, len, str - old_str);
        }
    }
  } /*  switch(state) */

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);

  _shift_left(_mant(result), exp_val, _mant(result));

  sc_val_from_ulong((1 << exp_size)/2-1, _exp(result));

  _normalize(result, result, 0);

  if (state == EXPONENT) {
    exp_int -= atoi(str-pos);
  }

  _power_of_ten(exp_int, &_desc(result), power_val);

  _fdiv(result, power_val, result);

  return result;
#else

  /* XXX excuse of an implementation to make things work */
  LLDBL val;
#ifdef HAVE_LONG_DOUBLE
  val = strtold(str, NULL);
#else
  val = strtod(str, NULL);
#endif

  DEBUGPRINTF(("val_from_str(%s)\n", str));
  return fc_val_from_float(val, exp_size, mant_size, result);
#endif
}

char* fc_val_from_float(LLDBL l, char exp_size, char mant_size, char* result)
{
  char *temp;
  int bias_res, bias_val, mant_val;
  value_t srcval;
  UINT32 sign, exponent, mantissa0, mantissa1;

  srcval.d = l;
  bias_res = ((1<<exp_size)/2-1);

#ifdef HAVE_LONG_DOUBLE
  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;
  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]));/* 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", 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);

  _desc(result).exponent_size = exp_size;
  _desc(result).mantissa_size = mant_size;

  /* extract sign */
  _sign(result) = sign;

  /* 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).clss = NAN;
    TRACEPRINTF(("val_from_float resulted in NAN\n"));
    return result;
  }
  else if (isinf(l)) {
    _desc(result).clss = INF;
    TRACEPRINTF(("val_from_float resulted in %sINF\n", (_sign(result)==1)?"-":""));
    return result;
  }

  /* 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

  /* build mantissa representation */
#ifndef HAVE_EXPLICIT_ONE
  if (exponent != 0)
  {
    /* insert the hidden bit */
    sc_val_from_ulong(1, temp);
    sc_val_from_ulong(mant_val + 2, NULL);
    _shift_left(temp, sc_get_buffer(), NULL);
  }
  else
#endif
  {
    sc_val_from_ulong(0, NULL);
  }

  _save_result(_mant(result));

  /* bits from the upper word */
  sc_val_from_ulong(mantissa0, temp);
  sc_val_from_ulong(34, NULL);
  _shift_left(temp, sc_get_buffer(), temp);
  sc_or(_mant(result), temp, _mant(result));

  /* bits from the lower word */
  sc_val_from_ulong(mantissa1, temp);
  sc_val_from_ulong(2, NULL);
  _shift_left(temp, sc_get_buffer(), temp);
  sc_or(_mant(result), temp, _mant(result));

  /* _normalize expects the radix point to be normal, so shift mantissa of subnormal
   * origin one to the left */
  if (exponent == 0)
  {
    sc_val_from_ulong(1, NULL);
    _shift_left(_mant(result), sc_get_buffer(), _mant(result));
  }

  _normalize(result, result, 0);

  TRACEPRINTF(("val_from_float results in %s\n", fc_print(result, temp, calc_buffer_size, FC_PACKED)));

  return result;
}

LLDBL fc_val_to_float(const void *val)
{
  const char *value;
  char *temp = NULL;

  int byte_offset;

  UINT32 sign;
  UINT32 exponent;
  UINT32 mantissa0;
  UINT32 mantissa1;

  value_t buildval;

#ifdef HAVE_LONG_DOUBLE
  char result_exponent = 15;
  char result_mantissa = 64;
#else
  char result_exponent = 11;
  char result_mantissa = 52;
#endif

  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
   * lead to wrong results */
  exponent = sc_val_to_long(_exp(value)) ;

  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(_mant(value), result_mantissa, byte_offset) << (byte_offset<<3);

  for (; (byte_offset<<3) < result_mantissa; byte_offset++)
    mantissa0 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << ((byte_offset-4)<<3);

#ifdef HAVE_LONG_DOUBLE
  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;
  buildval.val.low = mantissa1;
#endif

  TRACEPRINTF(("val_to_float: %d-%x-%x%x\n", sign, exponent, mantissa0, mantissa1));
  return buildval.d;
}

char* fc_cast(const void *val, char exp_size, char mant_size, char *result)
{
  const char *value = (const char*) val;
  char *temp;
  int exp_offset, val_bias, res_bias;

  if (result == NULL) result = calc_buffer;
  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);
    return result;
  }

  /* set the descriptor of the new value */
  _desc(result).exponent_size = exp_size;
  _desc(result).mantissa_size = mant_size;
  _desc(result).clss = _desc(value).clss;

  _sign(result) = _sign(value);

  /* when the mantissa sizes differ normalizing has to shift to align it.
   * this would change the exponent, which is unwanted. So calculate this
   * offset and add it */
  val_bias = (1<<_desc(value).exponent_size)/2-1;
  res_bias = (1<<exp_size)/2-1;