nsz Git - libfirm/blob - ir/tv/fltcalc.c

   1 /*
   2  * Project:     libFIRM
   3  * File name:   ir/tv/fltcalc.c
   4  * Purpose:
   5  * Author:
   6  * Modified by:
   7  * Created:     2003
   8  * CVS-ID:      $Id$
   9  * Copyright:   (c) 2003 Universität Karlsruhe
  10  * Licence:     This file protected by GPL -  GNU GENERAL PUBLIC LICENSE.
  11  */
  12
  13
  14 #ifdef HAVE_CONFIG_H
  15 # include "config.h"
  16 #endif
  17
  18
  19 #include "fltcalc.h"
  20 #include "strcalc.h"
  21
  22 #include <math.h>    /* need isnan() and isinf() (will be changed)*/
  23 /* undef some reused constants defined by math.h */
  24 #ifdef NAN
  25 #  undef NAN
  26 #endif
  27
  28 #ifdef HAVE_INTTYPES_H
  29 # include <inttypes.h>
  30 #endif
  31 #ifdef HAVE_STRING_H
  32 # include <string.h>
  33 #endif
  34 #ifdef HAVE_STDLIB_H
  35 # include <stdlib.h>
  36 #endif
  37 #ifdef HAVE_ALLOCA_H
  38 # include <alloca.h>
  39 #endif
  40 #ifdef HAVE_MALLOC_H
  41 # include <malloc.h>
  42 #endif
  43 #include <stdio.h>
  44 #include <assert.h>
  45
  46 #include "xmalloc.h"
  47
  48 typedef uint32_t UINT32;
  49
  50 #ifdef HAVE_LONG_DOUBLE
  51 #ifdef WORDS_BIGENDIAN
  52 typedef union {
  53   struct {
  54     UINT32 high;
  55     UINT32 mid;
  56     UINT32 low;
  57   } val;
  58   volatile long double d;
  59 } value_t;
  60 #else
  61 typedef union {
  62   struct {
  63     UINT32 low;
  64     UINT32 mid;
  65     UINT32 high;
  66   } val;
  67   volatile long double d;
  68 } value_t;
  69 #endif
  70 #else
  71 #ifdef WORDS_BIGENDIAN
  72 typedef union {
  73   struct {
  74     UINT32 high;
  75     UINT32 low;
  76   } val;
  77   volatile double d;
  78 } value_t;
  79 #else
  80 typedef union {
  81   struct {
  82     UINT32 low;
  83     UINT32 high;
  84   } val;
  85   volatile double d;
  86 } value_t;
  87 #endif
  88 #endif
  89
  90 /**
  91  * possible float states
  92  */
  93 typedef enum {
  94   NORMAL,       /**< normal representation, implicit 1 */
  95   ZERO,         /**< +/-0 */
  96   SUBNORMAL,    /**< denormals, implicit 0 */
  97   INF,          /**< +/-oo */
  98   NAN,          /**< Not A Number */
  99 } value_class_t;
 100
 101 /** A descriptor for an IEEE float value. */
 102 typedef struct {
 103   unsigned char exponent_size;    /**< size of exponent in bits */
 104   unsigned char mantissa_size;    /**< size of mantissa in bits */
 105   value_class_t clss;             /**< state of this float */
 106 } descriptor_t;
 107
 108 #define CLEAR_BUFFER(buffer) memset(buffer, 0, calc_buffer_size)
 109
 110 /* because variable sized structs are impossible, the internal
 111  * value is represented as a pseudo-struct char array, addressed
 112  * by macros
 113  * struct {
 114  *   char sign;             //  0 for positive, 1 for negative
 115  *   char exp[value_size];
 116  *   char mant[value_size];
 117  *   descriptor_t desc;
 118  * };
 119  */
 120 #define _sign(a) (((char*)a)[SIGN_POS])
 121 #define _exp(a) (&((char*)a)[EXPONENT_POS])
 122 #define _mant(a) (&((char*)a)[MANTISSA_POS])
 123 #define _desc(a) (*(descriptor_t *)&((char *)a)[DESCRIPTOR_POS])
 124
 125 #define _save_result(x) memcpy((x), sc_get_buffer(), value_size)
 126 #define _shift_right(x, y, b) sc_shr((x), (y), value_size*4, 0, (b))
 127 #define _shift_left(x, y, b) sc_shl((x), (y), value_size*4, 0, (b))
 128
 129 #define FC_DEFINE1(code) \
 130 char *fc_##code(const void *a, void *result)  {                        \
 131   return _calc((const char*)a, NULL, FC_##code, (char*)result);        \
 132 }
 133
 134 #define FC_DEFINE2(code) \
 135 char *fc_##code(const void *a, const void *b, void *result) {             \
 136   return _calc((const char*)a, (const char*)b, FC_##code, (char*)result); \
 137 }
 138
 139 #define FUNC_PTR(code) fc_##code
 140
 141 #if FLTCALC_DEBUG
 142 #  define DEBUGPRINTF(x) printf x
 143 #else
 144 #  define DEBUGPRINTF(x) ((void)0)
 145 #endif
 146
 147 #if FLTCALC_TRACE_CALC
 148 #  define TRACEPRINTF(x) printf x
 149 #else
 150 #  define TRACEPRINTF(x) ((void)0)
 151 #endif
 152
 153 static char *calc_buffer = NULL;
 154
 155 static fc_rounding_mode_t rounding_mode;
 156
 157 static int calc_buffer_size;
 158 static int value_size;
 159 static int SIGN_POS;
 160 static int EXPONENT_POS;
 161 static int MANTISSA_POS;
 162 static int DESCRIPTOR_POS;
 163
 164 static int max_precision;
 165 /********
 166  * private functions
 167  ********/
 168 #if 0
 169 static void _fail_char(const char *str, unsigned int len, int pos)
 170 {
 171   if (*(str+pos))
 172     printf("ERROR: Unexpected character '%c'\n", *(str + pos));
 173   else
 174     printf("ERROR: Unexpected end of string\n");
 175   while (len-- && *str) printf("%c", *str++); printf("\n");
 176   while (pos--) printf(" "); printf("^\n");
 177   /* the front end has to to check constant strings */
 178   exit(-1);
 179 }
 180 #endif
 181
 182 /** pack machine-like */
 183 static char* _pack(const char *int_float, char *packed)
 184 {
 185   char *shift_val;
 186   char *temp;
 187   char *val_buffer;
 188
 189   temp = alloca(value_size);
 190   shift_val = alloca(value_size);
 191
 192   switch (_desc(int_float).clss) {
 193     case NAN:
 194       val_buffer = alloca(calc_buffer_size);
 195       fc_get_qnan(_desc(int_float).exponent_size, _desc(int_float).mantissa_size, val_buffer);
 196       int_float = val_buffer;
 197       break;
 198
 199     case INF:
 200       val_buffer = alloca(calc_buffer_size);
 201       fc_get_plusinf(_desc(int_float).exponent_size, _desc(int_float).mantissa_size, val_buffer);
 202       _sign(val_buffer) = _sign(int_float);
 203       int_float = val_buffer;
 204       break;
 205
 206     default:
 207       break;
 208   }
 209   /* pack sign */
 210   sc_val_from_ulong(_sign(int_float), temp);
 211
 212   sc_val_from_ulong(_desc(int_float).exponent_size + _desc(int_float).mantissa_size, NULL);
 213   _shift_left(temp, sc_get_buffer(), packed);
 214
 215   /* extract exponent */
 216   sc_val_from_ulong(_desc(int_float).mantissa_size, shift_val);
 217
 218   _shift_left(_exp(int_float), shift_val, temp);
 219
 220   sc_or(temp, packed, packed);
 221
 222   /* extract mantissa */
 223   /* remove 2 rounding bits */
 224   sc_val_from_ulong(2, shift_val);
 225   _shift_right(_mant(int_float), shift_val, temp);
 226
 227   /* remove leading 1 (or 0 if denormalized) */
 228   sc_max_from_bits(_desc(int_float).mantissa_size, 0, shift_val); /* all mantissa bits are 1's */
 229   sc_and(temp, shift_val, temp);
 230
 231   /* save result */
 232   sc_or(temp, packed, packed);
 233
 234   return packed;
 235 }
 236
 237 char* _normalize(const char *in_val, char *out_val, int sticky)
 238 {
 239   int hsb;
 240   char lsb, guard, round, round_dir = 0;
 241   char *temp;
 242
 243   temp = alloca(value_size);
 244
 245   /* +2: save two rounding bits at the end */
 246   hsb = 2 + _desc(in_val).mantissa_size - sc_get_highest_set_bit(_mant(in_val)) - 1;
 247
 248   if (in_val != out_val)
 249   {
 250     _sign(out_val) = _sign(in_val);
 251     memcpy(&_desc(out_val), &_desc(in_val), sizeof(descriptor_t));
 252   }
 253
 254   _desc(out_val).clss = NORMAL;
 255
 256   /* mantissa all zeros, so zero exponent (because of explicit one)*/
 257   if (hsb == 2 + _desc(in_val).mantissa_size)
 258   {
 259     sc_val_from_ulong(0, _exp(out_val));
 260     hsb = -1;
 261   }
 262
 263   /* shift the first 1 into the left of the radix point (i.e. hsb == -1) */
 264   if (hsb < -1)
 265   {
 266     /* shift right */
 267     sc_val_from_ulong(-hsb-1, temp);
 268
 269     _shift_right(_mant(in_val), temp, _mant(out_val));
 270
 271     /* remember if some bits were shifted away */
 272     if (!sticky) sticky = sc_had_carry();
 273
 274     sc_add(_exp(in_val), temp, _exp(out_val));
 275   }
 276   else if (hsb > -1)
 277   {
 278     /* shift left */
 279     sc_val_from_ulong(hsb+1, temp);
 280
 281     _shift_left(_mant(in_val), temp, _mant(out_val));
 282
 283     sc_sub(_exp(in_val), temp, _exp(out_val));
 284   }
 285
 286   /* check for exponent underflow */
 287   if (sc_is_negative(_exp(out_val)) || sc_is_zero(_exp(out_val))) {
 288     DEBUGPRINTF(("Exponent underflow!\n"));
 289     /* exponent underflow */
 290     /* shift the mantissa right to have a zero exponent */
 291     sc_val_from_ulong(1, temp);
 292     sc_sub(temp, _exp(out_val), NULL);
 293
 294     _shift_right(_mant(out_val), sc_get_buffer(), _mant(out_val));
 295     if (!sticky) sticky = sc_had_carry();
 296     /* denormalized means exponent of zero */
 297     sc_val_from_ulong(0, _exp(out_val));
 298
 299     _desc(out_val).clss = SUBNORMAL;
 300   }
 301
 302   /* perform rounding by adding a value that clears the guard bit and the round bit
 303    * and either causes a carry to round up or not */
 304   /* get the last 3 bits of the value */
 305   lsb = sc_sub_bits(_mant(out_val), _desc(out_val).mantissa_size + 2, 0) & 0x7;
 306   guard = (lsb&0x2)>>1;
 307   round = lsb&0x1;
 308
 309   switch (rounding_mode)
 310   {
 311     case FC_TONEAREST:
 312       /* round to nearest representable value, if in doubt choose the version
 313        * with lsb == 0 */
 314       round_dir = guard && (sticky || round || lsb>>2);
 315       break;
 316     case FC_TOPOSITIVE:
 317       /* if positive: round to one if the exact value is bigger, else to zero */
 318       round_dir = (!_sign(out_val) && (guard || round || sticky));
 319       break;
 320     case FC_TONEGATIVE:
 321       /* if negative: round to one if the exact value is bigger, else to zero */
 322       round_dir = (_sign(out_val) && (guard || round || sticky));
 323       break;
 324     case FC_TOZERO:
 325       /* always round to 0 (chopping mode) */
 326       round_dir = 0;
 327       break;
 328   }
 329   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"));
 330
 331   if (round_dir == 1)
 332   {
 333     guard = (round^guard)<<1;
 334     lsb = !(round || guard)<<2 | guard | round;
 335   }
 336   else
 337   {
 338     lsb = -((guard<<1) | round);
 339   }
 340
 341   /* add the rounded value */
 342   if (lsb != 0) {
 343     sc_val_from_long(lsb, temp);
 344     sc_add(_mant(out_val), temp, _mant(out_val));
 345   }
 346
 347   /* could have rounded down to zero */
 348   if (sc_is_zero(_mant(out_val)) && (_desc(out_val).clss == SUBNORMAL))
 349     _desc(out_val).clss = ZERO;
 350
 351   /* check for rounding overflow */
 352   hsb = 2 + _desc(out_val).mantissa_size - sc_get_highest_set_bit(_mant(out_val)) - 1;
 353   if ((_desc(out_val).clss != SUBNORMAL) && (hsb < -1))
 354   {
 355     sc_val_from_ulong(1, temp);
 356     _shift_right(_mant(out_val), temp, _mant(out_val));
 357
 358     sc_add(_exp(out_val), temp, _exp(out_val));
 359   }
 360   else if ((_desc(out_val).clss == SUBNORMAL) && (hsb == -1))
 361   {
 362     /* overflow caused the mantissa to be normal again,
 363      * so adapt the exponent accordingly */
 364     sc_val_from_ulong(1, temp);
 365     sc_add(_exp(out_val), temp, _exp(out_val));
 366
 367     _desc(out_val).clss = NORMAL;
 368   }
 369   /* no further rounding is needed, because rounding overflow means
 370    * the carry of the original rounding was propagated all the way
 371    * up to the bit left of the radix point. This implies the bits
 372    * to the right are all zeros (rounding is +1) */
 373
 374   /* check for exponent overflow */
 375   sc_val_from_ulong((1 << _desc(out_val).exponent_size) - 1, temp);
 376   if (sc_comp(_exp(out_val), temp) != -1) {
 377     DEBUGPRINTF(("Exponent overflow!\n"));
 378     /* exponent overflow, reaction depends on rounding method:
 379      *
 380      * mode        | sign of value |  result
 381      *--------------------------------------------------------------
 382      * TO_NEAREST  |      +        |   +inf
 383      *             |      -        |   -inf
 384      *--------------------------------------------------------------
 385      * TO_POSITIVE |      +        |   +inf
 386      *             |      -        |   smallest representable value
 387      *--------------------------------------------------------------
 388      * TO_NEAGTIVE |      +        |   largest representable value
 389      *             |      -        |   -inf
 390      *--------------------------------------------------------------
 391      * TO_ZERO     |      +        |   largest representable value
 392      *             |      -        |   smallest representable value
 393      *--------------------------------------------------------------*/
 394     if (_sign(out_val) == 0)
 395     {
 396       /* value is positive */
 397       switch (rounding_mode) {
 398         case FC_TONEAREST:
 399         case FC_TOPOSITIVE:
 400           _desc(out_val).clss = INF;
 401           break;
 402
 403         case FC_TONEGATIVE:
 404         case FC_TOZERO:
 405           fc_get_max(_desc(out_val).exponent_size, _desc(out_val).mantissa_size, out_val);
 406       }
 407     } else {
 408       /* value is negative */
 409       switch (rounding_mode) {
 410         case FC_TONEAREST:
 411         case FC_TONEGATIVE:
 412           _desc(out_val).clss = INF;
 413           break;
 414
 415         case FC_TOPOSITIVE:
 416         case FC_TOZERO:
 417           fc_get_min(_desc(out_val).exponent_size, _desc(out_val).mantissa_size, out_val);
 418       }
 419     }
 420   }
 421
 422   return out_val;
 423 }
 424
 425 /**
 426  * Operations involving NaN's must return NaN
 427  */
 428 #define handle_NAN(a, b, result) \
 429 do {                                                      \
 430   if (_desc(a).clss == NAN) {                             \
 431     if (a != result) memcpy(result, a, calc_buffer_size); \
 432     return result;                                        \
 433   }                                                       \
 434   if (_desc(b).clss == NAN) {                             \
 435     if (b != result) memcpy(result, b, calc_buffer_size); \
 436     return result;                                        \
 437   }                                                       \
 438 }while (0)
 439
 440
 441 /**
 442  * calculate a + b, where a is the value with the bigger exponent
 443  */
 444 static char* _fadd(const char* a, const char* b, char* result)
 445 {
 446   char *temp;
 447   char *exp_diff;
 448
 449   char sign;
 450   char sticky;
 451
 452   handle_NAN(a, b, result);
 453
 454   /* make sure result has a descriptor */
 455   if (result != a && result != b)
 456     memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
 457
 458   /* determine if this is an addition or subtraction */
 459   sign = _sign(a) ^ _sign(b);
 460
 461   /* produce NaN on inf - inf */
 462   if (sign && (_desc(a).clss == INF) && (_desc(b).clss == INF))
 463     return fc_get_qnan(_desc(a).exponent_size, _desc(b).mantissa_size, result);
 464
 465   temp     = alloca(value_size);
 466   exp_diff = alloca(value_size);
 467
 468   /* get exponent difference */
 469   sc_sub(_exp(a), _exp(b), exp_diff);
 470
 471   /* initially set sign to be the sign of a, special treatment of subtraction
 472    * when exponents are equal is required though.
 473    * Also special care about the sign is needed when the mantissas are equal
 474    * (+/- 0 ?) */
 475   if (sign && sc_val_to_long(exp_diff) == 0) {
 476     switch (sc_comp(_mant(a), _mant(b))) {
 477       case 1:  /* a > b */
 478         _sign(result) = _sign(a);  /* abs(a) is bigger and a is negative */
 479         break;
 480       case 0:  /* a == b */
 481         _sign(result) = (rounding_mode == FC_TONEGATIVE);
 482         break;
 483       case -1: /* a < b */
 484         _sign(result) = _sign(b); /* abs(b) is bigger and b is negative */
 485         break;
 486       default:
 487         /* can't be reached */
 488         break;
 489     }
 490   }
 491   else
 492     _sign(result) = _sign(a);
 493
 494   /* sign has been taken care of, check for special cases */
 495   if (_desc(a).clss == ZERO) {
 496     if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
 497     return result;
 498   }
 499   if (_desc(b).clss == ZERO) {
 500     if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
 501     return result;
 502   }
 503
 504   if (_desc(a).clss == INF) {
 505     if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
 506     return result;
 507   }
 508   if (_desc(b).clss == INF) {
 509     if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
 510     return result;
 511   }
 512
 513   /* shift the smaller value to the right to align the radix point */
 514   /* subnormals have their radix point shifted to the right,
 515    * take care of this first */
 516   if ((_desc(b).clss == SUBNORMAL) && (_desc(a).clss != SUBNORMAL))
 517   {
 518     sc_val_from_ulong(1, temp);
 519     sc_sub(exp_diff, temp, exp_diff);
 520   }
 521
 522   _shift_right(_mant(b), exp_diff, temp);
 523   sticky = sc_had_carry();
 524
 525   if (sticky && sign)
 526   {
 527     /* if subtracting a little more than the represented value or adding a little
 528      * more than the represented value to a negative value this, in addition to the
 529      * still set sticky bit, takes account of the 'little more' */
 530     char *temp1 = alloca(calc_buffer_size);
 531     sc_val_from_ulong(1, temp1);
 532     sc_add(temp, temp1, temp);
 533   }
 534
 535   if (sign) {
 536     if (sc_comp(_mant(a), temp) == -1)
 537       sc_sub(temp, _mant(a), _mant(result));
 538     else
 539       sc_sub(_mant(a), temp, _mant(result));
 540   } else {
 541     sc_add(_mant(a), temp, _mant(result));
 542   }
 543
 544   /* _normalize expects a 'normal' radix point, adding two subnormals
 545    * results in a subnormal radix point -> shifting before normalizing */
 546   if ((_desc(a).clss == SUBNORMAL) && (_desc(b).clss == SUBNORMAL))
 547   {
 548     sc_val_from_ulong(1, NULL);
 549     _shift_left(_mant(result), sc_get_buffer(), _mant(result));
 550   }
 551
 552   /* resulting exponent is the bigger one */
 553   memmove(_exp(result), _exp(a), value_size);
 554
 555   return _normalize(result, result, sticky);
 556 }
 557
 558 /**
 559  * calculate a * b
 560  */
 561 static char* _fmul(const char* a, const char* b, char* result)
 562 {
 563   char *temp;
 564
 565   handle_NAN(a, b, result);
 566
 567   temp = alloca(value_size);
 568
 569   if (result != a && result != b)
 570     memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
 571
 572   _sign(result) = _sign(a) ^ _sign(b);
 573
 574   /* produce NaN on 0 * inf */
 575   if (_desc(a).clss == ZERO) {
 576     if (_desc(b).clss == INF)
 577       fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
 578     else
 579       if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
 580     return result;
 581   }
 582   if (_desc(b).clss == ZERO) {
 583     if (_desc(a).clss == INF)
 584       fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
 585     else
 586       if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-1);
 587     return result;
 588   }
 589
 590   if (_desc(a).clss == INF) {
 591     if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
 592     return result;
 593   }
 594   if (_desc(b).clss == INF) {
 595     if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-1);
 596     return result;
 597   }
 598
 599   /* exp = exp(a) + exp(b) - excess */
 600   sc_add(_exp(a), _exp(b), _exp(result));
 601
 602   sc_val_from_ulong((1<<_desc(a).exponent_size)/2-1, temp);
 603   sc_sub(_exp(result), temp, _exp(result));
 604
 605   /* mixed normal, subnormal values introduce an error of 1, correct it */
 606   if ((_desc(a).clss == SUBNORMAL) ^ (_desc(b).clss == SUBNORMAL))
 607   {
 608     sc_val_from_ulong(1, temp);
 609     sc_add(_exp(result), temp, _exp(result));
 610   }
 611
 612   sc_mul(_mant(a), _mant(b), _mant(result));
 613
 614   /* realign result: after a multiplication the digits right of the radix
 615    * point are the sum of the factors' digits after the radix point. As all
 616    * values are normalized they both have the same amount of these digits,
 617    * which has to be restored by proper shifting
 618    * +2 because of the two rounding bits */
 619   sc_val_from_ulong(2 + _desc(result).mantissa_size, temp);
 620
 621   _shift_right(_mant(result), temp, _mant(result));
 622
 623   return _normalize(result, result, sc_had_carry());
 624 }
 625
 626 /**
 627  * calculate a / b
 628  */
 629 static char* _fdiv(const char* a, const char* b, char* result)
 630 {
 631   char *temp, *dividend;
 632
 633   handle_NAN(a, b, result);
 634
 635   temp = alloca(value_size);
 636   dividend = alloca(value_size);
 637
 638   if (result != a && result != b)
 639     memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
 640
 641   _sign(result) = _sign(a) ^ _sign(b);
 642
 643   /* produce nan on 0/0 and inf/inf */
 644   if (_desc(a).clss == ZERO) {
 645     if (_desc(b).clss == ZERO)
 646       /* 0/0 -> nan */
 647       fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
 648     else
 649       /* 0/x -> a */
 650       if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
 651     return result;
 652   }
 653
 654   if (_desc(b).clss == INF) {
 655     if (_desc(a).clss == INF)
 656       /* inf/inf -> nan */
 657       fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
 658     else {
 659       /* x/inf -> 0 */
 660       sc_val_from_ulong(0, NULL);
 661       _save_result(_exp(result));
 662       _save_result(_mant(result));
 663       _desc(result).clss = ZERO;
 664     }
 665     return result;
 666   }
 667
 668   if (_desc(a).clss == INF) {
 669     /* inf/x -> inf */
 670     if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
 671     return result;
 672   }
 673   if (_desc(b).clss == ZERO) {
 674     /* division by zero */
 675     if (_sign(result))
 676       fc_get_minusinf(_desc(a).exponent_size, _desc(a).mantissa_size, result);
 677     else
 678       fc_get_plusinf(_desc(a).exponent_size, _desc(a).mantissa_size, result);
 679     return result;
 680   }
 681
 682   /* exp = exp(a) - exp(b) + excess - 1*/
 683   sc_sub(_exp(a), _exp(b), _exp(result));
 684   sc_val_from_ulong((1 << _desc(a).exponent_size)/2-2, temp);
 685   sc_add(_exp(result), temp, _exp(result));
 686
 687   /* mixed normal, subnormal values introduce an error of 1, correct it */
 688   if ((_desc(a).clss == SUBNORMAL) ^ (_desc(b).clss == SUBNORMAL))
 689   {
 690     sc_val_from_ulong(1, temp);
 691     sc_add(_exp(result), temp, _exp(result));
 692   }
 693
 694   /* mant(res) = mant(a) / 1/2mant(b) */
 695   /* to gain more bits of precision in the result the dividend could be
 696    * shifted left, as this operation does not loose bits. This would not
 697    * fit into the integer precision, but due to the rounding bits (which
 698    * are always zero because the values are all normalized) the divisor
 699    * can be shifted right instead to achieve the same result */
 700   sc_val_from_ulong(2 + _desc(result).mantissa_size, temp);
 701
 702   _shift_left(_mant(a), temp, dividend);
 703
 704   {
 705     char *divisor = alloca(calc_buffer_size);
 706     sc_val_from_ulong(1, divisor);
 707     _shift_right(_mant(b), divisor, divisor);
 708     sc_div(dividend, divisor, _mant(result));
 709   }
 710
 711   return _normalize(result, result, sc_had_carry());
 712 }
 713
 714 #if 0
 715 static void _power_of_ten(int exp, descriptor_t *desc, char *result)
 716 {
 717   char *build;
 718   char *temp;
 719
 720   /* positive sign */
 721   _sign(result) = 0;
 722
 723   /* set new descriptor (else result is supposed to already have one) */
 724   if (desc != NULL)
 725     memcpy(&_desc(result), desc, sizeof(descriptor_t));
 726
 727   build = alloca(value_size);
 728   temp = alloca(value_size);
 729
 730   sc_val_from_ulong((1 << _desc(result).exponent_size)/2-1, _exp(result));
 731
 732   if (exp > 0)
 733   {
 734     /* temp is value of ten now */
 735     sc_val_from_ulong(10, NULL);
 736     _save_result(temp);
 737
 738     for (exp--; exp > 0; exp--) {
 739       _save_result(build);
 740       sc_mul(build, temp, NULL);
 741     }
 742     _save_result(build);
 743
 744     /* temp is amount of left shift needed to put the value left of the radix point */
 745     sc_val_from_ulong(_desc(result).mantissa_size + 2, temp);
 746
 747     _shift_left(build, temp, _mant(result));
 748
 749     _normalize(result, result, 0);
 750   }
 751 }
 752 #endif
 753
 754 /**
 755  * Truncate the fractional part away.
 756  *
 757  * This does not clip to any integer rang.
 758  */
 759 static char* _trunc(const char *a, char *result)
 760 {
 761   /*
 762    * When exponent == 0 all bits left of the radix point
 763    * are the integral part of the value. For 15bit exp_size
 764    * this would require a left shift of max. 16383 bits which
 765    * is too much.
 766    * But it is enough to ensure that no bit right of the radix
 767    * point remains set. This restricts the interesting
 768    * exponents to the interval [0, mant_size-1].
 769    * Outside this interval the truncated value is either 0 or
 770    * it does not have fractional parts.
 771    */
 772
 773   int exp_bias, exp_val;
 774   char *temp;
 775
 776   temp = alloca(value_size);
 777
 778   if (a != result)
 779     memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
 780
 781   exp_bias = (1<<_desc(a).exponent_size)/2-1;
 782   exp_val  = sc_val_to_long(_exp(a)) - exp_bias;
 783
 784   if (exp_val < 0) {
 785     sc_val_from_ulong(0, NULL);
 786     _save_result(_exp(result));
 787     _save_result(_mant(result));
 788     _desc(result).clss = ZERO;
 789
 790     return result;
 791   }
 792
 793   if (exp_val > _desc(a).mantissa_size) {
 794     if (a != result)
 795       memcpy(result, a, calc_buffer_size);
 796
 797     return result;
 798   }
 799
 800   /* set up a proper mask to delete all bits right of the
 801    * radix point if the mantissa had been shifted until exp == 0 */
 802   sc_max_from_bits(1 + exp_val, 0, temp);
 803   sc_val_from_long(_desc(a).mantissa_size - exp_val + 2, NULL);
 804   _shift_left(temp, sc_get_buffer(), temp);
 805
 806   /* and the mask and return the result */
 807   sc_and(_mant(a), temp, _mant(result));
 808
 809   if (a != result) memcpy(_exp(result), _exp(a), value_size);
 810
 811   return result;
 812 }
 813
 814 /*
 815  * This does value sanity checking(or should do it), sets up any prerequisites,
 816  * calls the proper internal functions, clears up and returns
 817  * the result
 818  */
 819 char* _calc(const char *a, const char *b, int opcode, char *result)
 820 {
 821   char *temp;
 822 #ifdef FLTCALC_TRACE_CALC
 823   char *buffer;
 824
 825   buffer = alloca(100);
 826 #endif
 827
 828   if (result == NULL) result = calc_buffer;
 829
 830   TRACEPRINTF(("%s ", fc_print(a, buffer, 100, FC_PACKED)));
 831   switch (opcode)
 832   {
 833     case FC_add:
 834       /* make the value with the bigger exponent the first one */
 835       TRACEPRINTF(("+ %s ", fc_print(b, buffer, 100, FC_PACKED)));
 836       if (sc_comp(_exp(a), _exp(b)) == -1)
 837         _fadd(b, a, result);
 838       else
 839         _fadd(a, b, result);
 840       break;
 841     case FC_sub:
 842       TRACEPRINTF(("- %s ", fc_print(b, buffer, 100, FC_PACKED)));
 843       temp = alloca(calc_buffer_size);
 844       memcpy(temp, b, calc_buffer_size);
 845       _sign(temp) = !_sign(b);
 846       if (sc_comp(_exp(a), _exp(temp)) == -1)
 847         _fadd(temp, a, result);
 848       else
 849         _fadd(a, temp, result);
 850       break;
 851     case FC_mul:
 852       TRACEPRINTF(("* %s ", fc_print(b, buffer, 100, FC_PACKED)));
 853       _fmul(a, b, result);
 854       break;
 855     case FC_div:
 856       TRACEPRINTF(("/ %s ", fc_print(b, buffer, 100, FC_PACKED)));
 857       _fdiv(a, b, result);
 858       break;
 859     case FC_neg:
 860       TRACEPRINTF(("negated "));
 861       if (a != result) memcpy(result, a, calc_buffer_size);
 862       _sign(result) = !_sign(a);
 863       break;
 864     case FC_int:
 865       TRACEPRINTF(("truncated to integer "));
 866       _trunc(a, result);
 867       break;
 868     case FC_rnd:
 869       TRACEPRINTF(("rounded to integer "));
 870       assert(0);
 871       break;
 872   }
 873
 874   TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_PACKED)));
 875   return result;
 876 }
 877
 878 /********
 879  * functions defined in fltcalc.h
 880  ********/
 881 const void *fc_get_buffer(void)
 882 {
 883   return calc_buffer;
 884 }
 885
 886 int fc_get_buffer_length(void)
 887 {
 888   return calc_buffer_size;
 889 }
 890
 891 char* fc_val_from_str(const char *str, unsigned int len, char exp_size, char mant_size, char *result)
 892 {
 893 #if 0
 894   enum {
 895     START,
 896     LEFT_OF_DOT,
 897     RIGHT_OF_DOT,
 898     EXP_START,
 899     EXPONENT,
 900     END
 901   };
 902
 903   char exp_sign;
 904   int exp_int, hsb, state;
 905
 906   const char *old_str;
 907
 908   int pos;
 909   char *mant_str, *exp_val, *power_val;
 910
 911   if (result == NULL) result = calc_buffer;
 912
 913   exp_val = alloca(value_size);
 914   power_val = alloca(calc_buffer_size);
 915   mant_str = alloca((len)?(len):(strlen(str)));
 916
 917   _desc(result).exponent_size = exp_size;
 918   _desc(result).mantissa_size = mant_size;
 919   _desc(result).clss = NORMAL;
 920
 921   old_str = str;
 922   pos = 0;
 923   exp_int = 0;
 924   state = START;
 925
 926   while (len == 0 || str-old_str < len)
 927   {
 928     switch (state) {
 929       case START:
 930         switch (*str) {
 931           case '+':
 932             _sign(result) = 0;
 933             state = LEFT_OF_DOT;
 934             str++;
 935             break;
 936
 937           case '-':
 938             _sign(result) = 1;
 939             state = LEFT_OF_DOT;
 940             str++;
 941             break;
 942
 943           case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
 944             _sign(result) = 0;
 945             state = LEFT_OF_DOT;
 946             break;
 947
 948           case '.':
 949             _sign(result) = 0;
 950             state = RIGHT_OF_DOT;
 951             str++;
 952             break;
 953
 954           case 'n':
 955           case 'N':
 956           case 'i':
 957           case 'I':
 958             break;
 959
 960           default:
 961             _fail_char(old_str, len, str - old_str);
 962         }
 963         break;
 964
 965       case LEFT_OF_DOT:
 966         switch (*str) {
 967           case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
 968             mant_str[pos++] = *(str++);
 969             break;
 970
 971           case '.':
 972             state = RIGHT_OF_DOT;
 973             str++;
 974             break;
 975
 976           case 'e':
 977           case 'E':
 978             state = EXP_START;
 979             str++;
 980             break;
 981
 982           case '\0':
 983             mant_str[pos] = '\0';
 984             goto done;
 985
 986           default:
 987             _fail_char(old_str, len, str - old_str);
 988         }
 989         break;
 990
 991       case RIGHT_OF_DOT:
 992         switch (*str) {
 993           case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
 994             mant_str[pos++] = *(str++);
 995             exp_int++;
 996             break;
 997
 998           case 'e':
 999           case 'E':
1000             state = EXP_START;
1001             str++;
1002             break;
1003
1004           case '\0':
1005             mant_str[pos] = '\0';
1006             goto done;
1007
1008           default:
1009             _fail_char(old_str, len, str - old_str);
1010         }
1011         break;
1012
1013       case EXP_START:
1014         switch (*str) {
1015           case '-':
1016             exp_sign = 1;
1017             /* fall through */
1018           case '+':
1019             if (*(str-1) != 'e' && *(str-1) != 'E') _fail_char(old_str, len, str - old_str);
1020             str++;
1021             break;
1022
1023           case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
1024             mant_str[pos] = '\0';
1025             pos = 1;
1026             str++;
1027             state = EXPONENT;
1028             break;
1029
1030           default:
1031             _fail_char(old_str, len, str - old_str);
1032         }
1033         break;
1034
1035       case EXPONENT:
1036         switch (*str) {
1037           case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
1038             pos++;
1039             str++;
1040             break;
1041
1042           case '\0': goto done;
1043
1044           default:
1045             _fail_char(old_str, len, str - old_str);
1046         }
1047     }
1048   } /*  switch(state) */
1049
1050 done:
1051   sc_val_from_str(mant_str, strlen(mant_str), _mant(result));
1052
1053   /* shift to put value left of radix point */
1054   sc_val_from_ulong(mant_size + 2, exp_val);
1055
1056   _shift_left(_mant(result), exp_val, _mant(result));
1057
1058   sc_val_from_ulong((1 << exp_size)/2-1, _exp(result));
1059
1060   _normalize(result, result, 0);
1061
1062   if (state == EXPONENT) {
1063     exp_int -= atoi(str-pos);
1064   }
1065
1066   _power_of_ten(exp_int, &_desc(result), power_val);
1067
1068   _fdiv(result, power_val, result);
1069
1070   return result;
1071 #else
1072
1073   /* XXX excuse of an implementation to make things work */
1074   LLDBL val;
1075 #ifdef HAVE_LONG_DOUBLE
1076   val = strtold(str, NULL);
1077 #else
1078   val = strtod(str, NULL);
1079 #endif
1080
1081   DEBUGPRINTF(("val_from_str(%s)\n", str));
1082   return fc_val_from_float(val, exp_size, mant_size, result);
1083 #endif
1084 }
1085
1086 char* fc_val_from_float(LLDBL l, char exp_size, char mant_size, char* result)
1087 {
1088   char *temp;
1089   int bias_res, bias_val, mant_val;
1090   value_t srcval;
1091   UINT32 sign, exponent, mantissa0, mantissa1;
1092
1093   srcval.d = l;
1094   bias_res = ((1<<exp_size)/2-1);
1095
1096 #ifdef HAVE_LONG_DOUBLE
1097   mant_val  = 64;
1098   bias_val  = 0x3fff;
1099   sign      = (srcval.val.high & 0x00008000) != 0;
1100   exponent  = (srcval.val.high & 0x00007FFF) ;
1101   mantissa0 = srcval.val.mid;
1102   mantissa1 = srcval.val.low;
1103 #else /* no long double */
1104   mant_val  = 52;
1105   bias_val  = 0x3ff;
1106   sign      = (srcval.val.high & 0x80000000) != 0;
1107   exponent  = (srcval.val.high & 0x7FF00000) >> 20;
1108   mantissa0 = srcval.val.high & 0x000FFFFF;
1109   mantissa1 = srcval.val.low;
1110 #endif
1111
1112 #ifdef HAVE_LONG_DOUBLE
1113   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)); */
1114   DEBUGPRINTF(("(%d-%.4X-%.8X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1115 #else
1116   TRACEPRINTF(("val_from_float(%.8X%.8X)\n", srcval.val.high, srcval.val.low));
1117   DEBUGPRINTF(("(%d-%.3X-%.5X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1118 #endif
1119
1120   if (result == NULL) result = calc_buffer;
1121   temp = alloca(value_size);
1122
1123   _desc(result).exponent_size = exp_size;
1124   _desc(result).mantissa_size = mant_size;
1125
1126   /* extract sign */
1127   _sign(result) = sign;
1128
1129   /* sign and flag suffice to identify nan or inf, no exponent/mantissa
1130    * encoding is needed. the function can return immediately in these cases */
1131   if (isnan(l)) {
1132     _desc(result).clss = NAN;
1133     TRACEPRINTF(("val_from_float resulted in NAN\n"));
1134     return result;
1135   }
1136   else if (isinf(l)) {
1137     _desc(result).clss = INF;
1138     TRACEPRINTF(("val_from_float resulted in %sINF\n", (_sign(result)==1)?"-":""));
1139     return result;
1140   }
1141
1142   /* build exponent, because input and output exponent and mantissa sizes may differ
1143    * this looks more complicated than it is: unbiased input exponent + output bias,
1144    * minus the mantissa difference which is added again later when the output float
1145    * becomes normalized */
1146 #ifdef HAVE_EXPLICIT_ONE
1147   sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size-1), _exp(result));
1148 #else
1149   sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size), _exp(result));
1150 #endif
1151
1152   /* build mantissa representation */
1153 #ifndef HAVE_EXPLICIT_ONE
1154   if (exponent != 0)
1155   {
1156     /* insert the hidden bit */
1157     sc_val_from_ulong(1, temp);
1158     sc_val_from_ulong(mant_val + 2, NULL);
1159     _shift_left(temp, sc_get_buffer(), NULL);
1160   }
1161   else
1162 #endif
1163   {
1164     sc_val_from_ulong(0, NULL);
1165   }
1166
1167   _save_result(_mant(result));
1168
1169   /* bits from the upper word */
1170   sc_val_from_ulong(mantissa0, temp);
1171   sc_val_from_ulong(34, NULL);
1172   _shift_left(temp, sc_get_buffer(), temp);
1173   sc_or(_mant(result), temp, _mant(result));
1174
1175   /* bits from the lower word */
1176   sc_val_from_ulong(mantissa1, temp);
1177   sc_val_from_ulong(2, NULL);
1178   _shift_left(temp, sc_get_buffer(), temp);
1179   sc_or(_mant(result), temp, _mant(result));
1180
1181   /* _normalize expects the radix point to be normal, so shift mantissa of subnormal
1182    * origin one to the left */
1183   if (exponent == 0)
1184   {
1185     sc_val_from_ulong(1, NULL);
1186     _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1187   }
1188
1189   _normalize(result, result, 0);
1190
1191   TRACEPRINTF(("val_from_float results in %s\n", fc_print(result, temp, calc_buffer_size, FC_PACKED)));
1192
1193   return result;
1194 }
1195
1196 LLDBL fc_val_to_float(const void *val)
1197 {
1198   const char *value;
1199   char *temp = NULL;
1200
1201   int byte_offset;
1202
1203   UINT32 sign;
1204   UINT32 exponent;
1205   UINT32 mantissa0;
1206   UINT32 mantissa1;
1207
1208   value_t buildval;
1209
1210 #ifdef HAVE_LONG_DOUBLE
1211   char result_exponent = 15;
1212   char result_mantissa = 64;
1213 #else
1214   char result_exponent = 11;
1215   char result_mantissa = 52;
1216 #endif
1217
1218   temp = alloca(calc_buffer_size);
1219 #ifdef HAVE_EXPLICIT_ONE
1220   value = fc_cast(val, result_exponent, result_mantissa-1, temp);
1221 #else
1222   value = fc_cast(val, result_exponent, result_mantissa, temp);
1223 #endif
1224
1225   sign = _sign(value);
1226
1227   /* @@@ long double exponent is 15bit, so the use of sc_val_to_long should not
1228    * lead to wrong results */
1229   exponent = sc_val_to_long(_exp(value)) ;
1230
1231   sc_val_from_ulong(2, NULL);
1232   _shift_right(_mant(value), sc_get_buffer(), _mant(value));
1233
1234   mantissa0 = 0;
1235   mantissa1 = 0;
1236
1237   for (byte_offset = 0; byte_offset < 4; byte_offset++)
1238     mantissa1 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << (byte_offset<<3);
1239
1240   for (; (byte_offset<<3) < result_mantissa; byte_offset++)
1241     mantissa0 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << ((byte_offset-4)<<3);
1242
1243 #ifdef HAVE_LONG_DOUBLE
1244   buildval.val.high = sign << 15;
1245   buildval.val.high |= exponent;
1246   buildval.val.mid = mantissa0;
1247   buildval.val.low = mantissa1;
1248 #else /* no long double */
1249   mantissa0 &= 0x000FFFFF;  /* get rid of garbage */
1250   buildval.val.high = sign << 31;
1251   buildval.val.high |= exponent << 20;
1252   buildval.val.high |= mantissa0;
1253   buildval.val.low = mantissa1;
1254 #endif
1255
1256   TRACEPRINTF(("val_to_float: %d-%x-%x%x\n", sign, exponent, mantissa0, mantissa1));
1257   return buildval.d;
1258 }
1259
1260 char* fc_cast(const void *val, char exp_size, char mant_size, char *result)
1261 {
1262   const char *value = (const char*) val;
1263   char *temp;
1264   int exp_offset, val_bias, res_bias;
1265
1266   if (result == NULL) result = calc_buffer;
1267   temp = alloca(value_size);
1268
1269   if (_desc(value).exponent_size == exp_size && _desc(value).mantissa_size == mant_size)
1270   {
1271     if (value != result) memcpy(result, value, calc_buffer_size);
1272     return result;
1273   }
1274
1275   /* set the descriptor of the new value */
1276   _desc(result).exponent_size = exp_size;
1277   _desc(result).mantissa_size = mant_size;
1278   _desc(result).clss = _desc(value).clss;
1279
1280   _sign(result) = _sign(value);
1281
1282   /* when the mantissa sizes differ normalizing has to shift to align it.
1283    * this would change the exponent, which is unwanted. So calculate this
1284    * offset and add it */
1285   val_bias = (1<<_desc(value).exponent_size)/2-1;
1286   res_bias = (1<<exp_size)/2-1;
1287
1288   exp_offset = (res_bias - val_bias) - (_desc(value).mantissa_size - mant_size);
1289   sc_val_from_long(exp_offset, temp);
1290   sc_add(_exp(value), temp, _exp(result));
1291
1292   /* _normalize expects normalized radix point */
1293   if (_desc(val).clss == SUBNORMAL) {
1294     sc_val_from_ulong(1, NULL);
1295     _shift_left(_mant(val), sc_get_buffer(), _mant(result));
1296   } else if (value != result) {
1297     memcpy(_mant(result), _mant(value), value_size);
1298   } else {
1299     memmove(_mant(result), _mant(value), value_size);
1300   }
1301
1302   _normalize(result, result, 0);
1303   TRACEPRINTF(("Cast results in %s\n", fc_print(result, temp, value_size, FC_PACKED)));
1304   return result;
1305 }
1306
1307 char* fc_get_max(unsigned int exponent_size, unsigned int mantissa_size, char* result)
1308 {
1309   if (result == NULL) result = calc_buffer;
1310
1311   _desc(result).exponent_size = exponent_size;
1312   _desc(result).mantissa_size = mantissa_size;
1313   _desc(result).clss = NORMAL;
1314
1315   _sign(result) = 0;
1316
1317   sc_val_from_ulong((1<<exponent_size) - 2, _exp(result));
1318
1319   sc_max_from_bits(mantissa_size + 1, 0, _mant(result));
1320   sc_val_from_ulong(2, NULL);
1321   _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1322
1323   return result;
1324 }
1325
1326 char* fc_get_min(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1327 {
1328   if (result == NULL) result = calc_buffer;
1329
1330   fc_get_max(exponent_size, mantissa_size, result);
1331   _sign(result) = 1;
1332
1333   return result;
1334 }
1335
1336 char* fc_get_snan(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1337 {
1338   if (result == NULL) result = calc_buffer;
1339
1340   _desc(result).exponent_size = exponent_size;
1341   _desc(result).mantissa_size = mantissa_size;
1342   _desc(result).clss = NAN;
1343
1344   _sign(result) = 0;
1345
1346   sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1347
1348   /* signaling NaN has non-zero mantissa with msb not set */
1349   sc_val_from_ulong(1, _mant(result));
1350
1351   return result;
1352 }
1353
1354 char* fc_get_qnan(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1355 {
1356   if (result == NULL) result = calc_buffer;
1357
1358   _desc(result).exponent_size = exponent_size;
1359   _desc(result).mantissa_size = mantissa_size;
1360   _desc(result).clss = NAN;
1361
1362   _sign(result) = 0;
1363
1364   sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1365
1366   /* quiet NaN has the msb of the mantissa set, so shift one there */
1367   sc_val_from_ulong(1, _mant(result));
1368   /* mantissa_size >+< 1 because of two extra rounding bits */
1369   sc_val_from_ulong(mantissa_size + 1, NULL);
1370   _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1371
1372   return result;
1373 }
1374
1375 char* fc_get_plusinf(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1376 {
1377   if (result == NULL) result = calc_buffer;
1378
1379   _desc(result).exponent_size = exponent_size;
1380   _desc(result).mantissa_size = mantissa_size;
1381   _desc(result).clss = NORMAL;
1382
1383   _sign(result) = 0;
1384
1385   sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1386
1387   sc_val_from_ulong(0, _mant(result));
1388
1389   return result;
1390 }
1391
1392 char* fc_get_minusinf(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1393 {
1394   if (result == NULL) result = calc_buffer;
1395
1396   fc_get_plusinf(exponent_size, mantissa_size, result);
1397   _sign(result) = 1;
1398
1399   return result;
1400 }
1401
1402 int fc_comp(const void *a, const void *b)
1403 {
1404   const char *val_a = (const char*)a;
1405   const char *val_b = (const char*)b;
1406   int mul = 1;
1407
1408   /*
1409    * shortcut: if both values are identical, they are either
1410    * Unordered if NaN or equal
1411    */
1412   if (a == b)
1413     return _desc(val_a).clss == NAN ? 2 : 0;
1414
1415   /* unordered if one is a NaN */
1416   if (_desc(val_a).clss == NAN || _desc(val_b).clss == NAN)
1417     return 2;
1418
1419   /* zero is equal independent of sign */
1420   if ((_desc(val_a).clss == ZERO) && (_desc(val_b).clss == ZERO))
1421     return 0;
1422
1423   /* different signs make compare easy */
1424   if (_sign(val_a) != _sign(val_b))
1425     return (_sign(val_a)==0)?(1):(-1);
1426
1427   mul = _sign(a) ? -1 : 1;
1428
1429   /* both infinity means equality */
1430   if ((_desc(val_a).clss == INF) && (_desc(val_b).clss == INF))
1431     return 0;
1432
1433   /* infinity is bigger than the rest */
1434   if (_desc(val_a).clss == INF)
1435     return  1 * mul;
1436   if (_desc(val_b).clss == INF)
1437     return -1 * mul;
1438
1439   /* check first exponent, that mantissa if equal */
1440   switch (sc_comp(_exp(val_a), _exp(val_b))) {
1441     case -1:
1442       return -1 * mul;
1443     case  1:
1444       return  1 * mul;
1445     case  0:
1446       return sc_comp(_mant(val_a), _mant(val_b)) * mul;
1447     default:
1448       return 2;
1449   }
1450 }
1451
1452 int fc_is_zero(const void *a)
1453 {
1454   return _desc(a).clss == ZERO;
1455 }
1456
1457 int fc_is_negative(const void *a)
1458 {
1459   return _sign(a);
1460 }
1461
1462 int fc_is_inf(const void *a)
1463 {
1464   return _desc(a).clss == INF;
1465 }
1466
1467 int fc_is_nan(const void *a)
1468 {
1469   return _desc(a).clss == NAN;
1470 }
1471
1472 int fc_is_subnormal(const void *a)
1473 {
1474   return _desc(a).clss == SUBNORMAL;
1475 }
1476
1477 char *fc_print(const void *a, char *buf, int buflen, unsigned base)
1478 {
1479   const char *val;
1480   char *mul_1;
1481
1482   val = (const char*)a;
1483
1484   mul_1 = alloca(calc_buffer_size);
1485
1486   switch (base) {
1487     case FC_DEC:
1488       switch (_desc(val).clss) {
1489         case INF:
1490           if (buflen >= 8+_sign(val)) sprintf(buf, "%sINFINITY", _sign(val)?"-":"");
1491           else snprintf(buf, buflen, "%sINF", _sign(val)?"-":NULL);
1492           break;
1493         case NAN:
1494           snprintf(buf, buflen, "NAN");
1495           break;
1496         case ZERO:
1497           snprintf(buf, buflen, "0.0");
1498           break;
1499         default:
1500           /* XXX to be implemented */
1501 #ifdef HAVE_LONG_DOUBLE
1502           /* XXX 30 is arbitrary */
1503           snprintf(buf, buflen, "%.30LE", fc_val_to_float(val));
1504 #else
1505           snprintf(buf, buflen, "%.18E", fc_val_to_float(val));
1506 #endif
1507       }
1508       break;
1509
1510     case FC_HEX:
1511       switch (_desc(val).clss) {
1512         case INF:
1513           if (buflen >= 8+_sign(val)) sprintf(buf, "%sINFINITY", _sign(val)?"-":"");
1514           else snprintf(buf, buflen, "%sINF", _sign(val)?"-":NULL);
1515           break;
1516         case NAN:
1517           snprintf(buf, buflen, "NAN");
1518           break;
1519         case ZERO:
1520           snprintf(buf, buflen, "0.0");
1521           break;
1522         default:
1523 #ifdef HAVE_LONG_DOUBLE
1524           snprintf(buf, buflen, "%LA", fc_val_to_float(val));
1525 #else
1526           snprintf(buf, buflen, "%A", fc_val_to_float(val));
1527 #endif
1528       }
1529       break;
1530
1531     case FC_PACKED:
1532     default:
1533       snprintf(buf, buflen, "%s", sc_print(_pack(val, mul_1), value_size*4, SC_HEX, 0));
1534       buf[buflen - 1] = '\0';
1535       break;
1536   }
1537   return buf;
1538 }
1539
1540 unsigned char fc_sub_bits(const void *value, unsigned num_bits, unsigned byte_ofs)
1541 {
1542   /* this is used to cache the packed version of the value */
1543   static char *pack = NULL;
1544
1545   if (pack == NULL) pack = xmalloc(value_size);
1546
1547   if (value != NULL)
1548     _pack((const char*)value, pack);
1549
1550   return sc_sub_bits(pack, num_bits, byte_ofs);
1551 }
1552
1553 fc_rounding_mode_t fc_set_rounding_mode(fc_rounding_mode_t mode)
1554 {
1555   if (mode == FC_TONEAREST || mode == FC_TOPOSITIVE || mode == FC_TONEGATIVE || mode == FC_TOZERO)
1556       rounding_mode = mode;
1557
1558   return rounding_mode;
1559 }
1560
1561 fc_rounding_mode_t fc_get_rounding_mode(void)
1562 {
1563   return rounding_mode;
1564 }
1565
1566 void init_fltcalc(int precision)
1567 {
1568   if (calc_buffer == NULL) {
1569     /* does nothing if already init */
1570     if (precision == 0) precision = FC_DEFAULT_PRECISION;
1571
1572     init_strcalc(precision + 4);
1573
1574     /* needs additionally two bits to round, a bit as explicit 1., and one for
1575      * addition overflow */
1576     max_precision = sc_get_precision() - 4;
1577     if (max_precision < precision)
1578       printf("WARING: not enough precision available, using %d\n", max_precision);
1579
1580     rounding_mode    = FC_TONEAREST;
1581     value_size       = sc_get_buffer_length();
1582     SIGN_POS         = 0;
1583     EXPONENT_POS     = SIGN_POS + sizeof(char);
1584     MANTISSA_POS     = EXPONENT_POS + value_size;
1585     DESCRIPTOR_POS   = MANTISSA_POS + value_size;
1586     calc_buffer_size = DESCRIPTOR_POS + sizeof(descriptor_t);
1587
1588     calc_buffer = xmalloc(calc_buffer_size);
1589     memset(calc_buffer, 0, calc_buffer_size);
1590     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));
1591 #ifdef HAVE_LONG_DOUBLE
1592     DEBUGPRINTF(("\tUsing long double (1-15-64) interface\n"));
1593 #else
1594     DEBUGPRINTF(("\tUsing double (1-11-52) interface\n"));
1595 #endif
1596 #ifdef WORDS_BIGENDIAN
1597     DEBUGPRINTF(("\tWord order is big endian\n\n"));
1598 #else
1599     DEBUGPRINTF(("\tWord order is little endian\n\n"));
1600 #endif
1601   }
1602 }
1603
1604 void finish_fltcalc (void) {
1605   free(calc_buffer); calc_buffer = NULL;
1606 }
1607
1608 /* definition of interface functions */
1609 FC_DEFINE2(add)
1610 FC_DEFINE2(sub)
1611 FC_DEFINE2(mul)
1612 FC_DEFINE2(div)
1613 FC_DEFINE1(neg)
1614 FC_DEFINE1(int)
1615 FC_DEFINE1(rnd)