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

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