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