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
1056   srcval.d = l;
1057   bias_res = ((1<<exp_size)/2-1);
1058
1059 #ifdef HAVE_LONG_DOUBLE
1060   mant_val = 64;
1061   bias_val = 0x3fff;
1062   UINT32 sign = (srcval.val.high & 0x00008000) != 0;
1063   UINT32 exponent = (srcval.val.high & 0x00007FFF) ;
1064   UINT32 mantissa0 = srcval.val.mid;
1065   UINT32 mantissa1 = srcval.val.low;
1066 #else /* no long double */
1067   mant_val = 52;
1068   bias_val = 0x3ff;
1069   UINT32 sign = (srcval.val.high & 0x80000000) != 0;
1070   UINT32 exponent = (srcval.val.high & 0x7FF00000) >> 20;
1071   UINT32 mantissa0 = srcval.val.high & 0x000FFFFF;
1072   UINT32 mantissa1 = srcval.val.low;
1073 #endif
1074
1075 #ifdef HAVE_LONG_DOUBLE
1076   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));
1077   DEBUGPRINTF(("(%d-%.4X-%.8X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1078 #else
1079   TRACEPRINTF(("val_from_float(%.8X%.8X)\n", srcval.val.high, srcval.val.low));
1080   DEBUGPRINTF(("(%d-%.3X-%.5X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1081 #endif
1082
1083   if (result == NULL) result = calc_buffer;
1084   temp = alloca(VALUE_SIZE);
1085
1086   _desc(result).exponent_size = exp_size;
1087   _desc(result).mantissa_size = mant_size;
1088
1089   /* extract sign */
1090   _sign(result) = sign;
1091
1092   /* sign and flag suffice to identify nan or inf, no exponent/mantissa
1093    * encoding is needed. the function can return immediately in these cases */
1094   if (isnan(l)) {
1095     _desc(result).class = NAN;
1096     TRACEPRINTF(("val_from_float resulted in NAN\n"));
1097     return result;
1098   }
1099   else if (isinf(l)) {
1100     _desc(result).class = INF;
1101     TRACEPRINTF(("val_from_float resulted in %sINF\n", (_sign(result)==1)?"-":""));
1102     return result;
1103   }
1104
1105   /* build exponent, because input and output exponent and mantissa sizes may differ
1106    * this looks more complicated than it is: unbiased input exponent + output bias,
1107    * minus the mantissa difference which is added again later when the output float
1108    * becomes normalized */
1109 #ifdef HAVE_EXPLICIT_ONE
1110   sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size-1), _exp(result));
1111 #else
1112   sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size), _exp(result));
1113 #endif
1114
1115   /* build mantissa representation */
1116 #ifndef HAVE_EXPLICIT_ONE
1117   if (exponent != 0)
1118   {
1119     /* insert the hidden bit */
1120     sc_val_from_ulong(1, temp);
1121     sc_val_from_ulong(mant_val + 2, NULL);
1122     _shift_left(temp, sc_get_buffer(), NULL);
1123   }
1124   else
1125 #endif
1126   {
1127     sc_val_from_ulong(0, NULL);
1128   }
1129
1130   _save_result(_mant(result));
1131
1132   /* bits from the upper word */
1133   sc_val_from_ulong(mantissa0, temp);
1134   sc_val_from_ulong(34, NULL);
1135   _shift_left(temp, sc_get_buffer(), temp);
1136   sc_or(_mant(result), temp, _mant(result));
1137
1138   /* bits from the lower word */
1139   sc_val_from_ulong(mantissa1, temp);
1140   sc_val_from_ulong(2, NULL);
1141   _shift_left(temp, sc_get_buffer(), temp);
1142   sc_or(_mant(result), temp, _mant(result));
1143
1144   /* _normalize expects the radix point to be normal, so shift mantissa of subnormal
1145    * origin one to the left */
1146   if (exponent == 0)
1147   {
1148     sc_val_from_ulong(1, NULL);
1149     _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1150   }
1151
1152   _normalize(result, result, 0);
1153
1154   TRACEPRINTF(("val_from_float results in %s\n", fc_print(result, temp, CALC_BUFFER_SIZE, FC_PACKED)));
1155
1156   return result;
1157 }
1158
1159 LLDBL fc_val_to_float(const void *val)
1160 {
1161   const char *value;
1162   char *temp = NULL;
1163
1164   int byte_offset;
1165
1166   UINT32 sign;
1167   UINT32 exponent;
1168   UINT32 mantissa0;
1169   UINT32 mantissa1;
1170
1171   value_t buildval;
1172
1173 #ifdef HAVE_LONG_DOUBLE
1174   char result_exponent = 15;
1175   char result_mantissa = 64;
1176 #else
1177   char result_exponent = 11;
1178   char result_mantissa = 52;
1179 #endif
1180
1181   temp = alloca(CALC_BUFFER_SIZE);
1182 #ifdef HAVE_EXPLICIT_ONE
1183   value = fc_cast(val, result_exponent, result_mantissa-1, temp);
1184 #else
1185   value = fc_cast(val, result_exponent, result_mantissa, temp);
1186 #endif
1187
1188   sign = _sign(value);
1189
1190   /* @@@ long double exponent is 15bit, so the use of sc_val_to_long should not
1191    * lead to wrong results */
1192   exponent = sc_val_to_long(_exp(value)) ;
1193
1194   sc_val_from_ulong(2, NULL);
1195   _shift_right(_mant(value), sc_get_buffer(), _mant(value));
1196
1197   mantissa0 = 0;
1198   mantissa1 = 0;
1199
1200   for (byte_offset = 0; byte_offset < 4; byte_offset++)
1201     mantissa1 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << (byte_offset<<3);
1202
1203   for (; (byte_offset<<3) < result_mantissa; byte_offset++)
1204     mantissa0 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << ((byte_offset-4)<<3);
1205
1206 #ifndef HAVE_LONG_DOUBLE
1207   mantissa0 &= 0x000FFFFF;  /* get rid of garbage */
1208 #endif
1209
1210 #ifdef HAVE_LONG_DOUBLE
1211   buildval.val.high = sign << 15;
1212   buildval.val.high |= exponent;
1213   buildval.val.mid = mantissa0;
1214   buildval.val.low = mantissa1;
1215 #else /* no long double */
1216   buildval.val.high = sign << 31;
1217   buildval.val.high |= exponent << 20;
1218   buildval.val.high |= mantissa0;
1219   buildval.val.low = mantissa1;
1220 #endif
1221
1222   TRACEPRINTF(("val_to_float: %d-%x-%x%x\n", sign, exponent, mantissa0, mantissa1));
1223   return buildval.d;
1224 }
1225
1226 char* fc_cast(const void *val, char exp_size, char mant_size, char *result)
1227 {
1228   const char *value = (const char*) val;
1229   char *temp;
1230   int exp_offset, val_bias, res_bias;
1231
1232   if (result == NULL) result = calc_buffer;
1233   temp = alloca(VALUE_SIZE);
1234
1235   if (_desc(value).exponent_size == exp_size && _desc(value).mantissa_size == mant_size)
1236   {
1237     if (value != result) memcpy(result, value, CALC_BUFFER_SIZE);
1238     return result;
1239   }
1240
1241   /* set the descriptor of the new value */
1242   _desc(result).exponent_size = exp_size;
1243   _desc(result).mantissa_size = mant_size;
1244   _desc(result).class = _desc(value).class;
1245
1246   _sign(result) = _sign(value);
1247
1248   /* when the mantissa sizes differ normalizing has to shift to align it.
1249    * this would change the exponent, which is unwanted. So calculate this
1250    * offset and add it */
1251   val_bias = (1<<_desc(value).exponent_size)/2-1;
1252   res_bias = (1<<exp_size)/2-1;
1253
1254   exp_offset = (res_bias - val_bias) - (_desc(value).mantissa_size - mant_size);
1255   sc_val_from_long(exp_offset, temp);
1256   sc_add(_exp(value), temp, _exp(result));
1257
1258   /* _normalize expects normalized radix point */
1259   if (_desc(val).class == SUBNORMAL) {
1260     sc_val_from_ulong(1, NULL);
1261     _shift_left(_mant(val), sc_get_buffer(), _mant(result));
1262   } else if (value != result) {
1263     memcpy(_mant(result), _mant(value), VALUE_SIZE);
1264   } else {
1265     memmove(_mant(result), _mant(value), VALUE_SIZE);
1266   }
1267
1268   _normalize(result, result, 0);
1269   TRACEPRINTF(("Cast results in %s\n", fc_print(result, temp, VALUE_SIZE, FC_PACKED)));
1270   return result;
1271 }
1272
1273 char* fc_get_max(unsigned int exponent_size, unsigned int mantissa_size, char* result)
1274 {
1275   if (result == NULL) result = calc_buffer;
1276
1277   _desc(result).exponent_size = exponent_size;
1278   _desc(result).mantissa_size = mantissa_size;
1279   _desc(result).class = NORMAL;
1280
1281   _sign(result) = 0;
1282
1283   sc_val_from_ulong((1<<exponent_size) - 2, _exp(result));
1284
1285   sc_max_from_bits(mantissa_size + 1, 0, _mant(result));
1286   sc_val_from_ulong(2, NULL);
1287   _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1288
1289   return result;
1290 }
1291
1292 char* fc_get_min(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1293 {
1294   if (result == NULL) result = calc_buffer;
1295
1296   fc_get_max(exponent_size, mantissa_size, result);
1297   _sign(result) = 1;
1298
1299   return result;
1300 }
1301
1302 char* fc_get_snan(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).class = NAN;
1309
1310   _sign(result) = 0;
1311
1312   sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1313
1314   /* signalling nan has non-zero mantissa with msb not set */
1315   sc_val_from_ulong(1, _mant(result));
1316
1317   return result;
1318 }
1319
1320 char* fc_get_qnan(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1321 {
1322   if (result == NULL) result = calc_buffer;
1323
1324   _desc(result).exponent_size = exponent_size;
1325   _desc(result).mantissa_size = mantissa_size;
1326   _desc(result).class = NAN;
1327
1328   _sign(result) = 0;
1329
1330   sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1331
1332   /* quiet nan has the msb of the mantissa set, so shift one there */
1333   sc_val_from_ulong(1, _mant(result));
1334   /* mantissa_size >+< 1 because of two extra rounding bits */
1335   sc_val_from_ulong(mantissa_size + 1, NULL);
1336   _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1337
1338   return result;
1339 }
1340
1341 char* fc_get_plusinf(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1342 {
1343   if (result == NULL) result = calc_buffer;
1344
1345   _desc(result).exponent_size = exponent_size;
1346   _desc(result).mantissa_size = mantissa_size;
1347   _desc(result).class = NORMAL;
1348
1349   _sign(result) = 0;
1350
1351   sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1352
1353   sc_val_from_ulong(0, _mant(result));
1354
1355   return result;
1356 }
1357
1358 char* fc_get_minusinf(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1359 {
1360   if (result == NULL) result = calc_buffer;
1361
1362   fc_get_plusinf(exponent_size, mantissa_size, result);
1363   _sign(result) = 1;
1364
1365   return result;
1366 }
1367
1368 int fc_comp(const void *a, const void *b)
1369 {
1370   const char *val_a = (const char*)a;
1371   const char *val_b = (const char*)b;
1372
1373   /* unordered */
1374   if (_desc(val_a).class == NAN || _desc(val_b).class == NAN) return 2;
1375   /* zero is equal independent of sign */
1376   if ((_desc(val_a).class == ZERO) && (_desc(val_b).class == ZERO)) return 0;
1377   /* different signs make compare easy */
1378   if (_sign(val_a) != _sign(val_b)) return (_sign(val_a)==0)?(1):(-1);
1379   /* both infinity means equality */
1380   if ((_desc(val_a).class == INF) && (_desc(val_b).class == INF)) return 0;
1381   /* infinity is bigger than the rest */
1382   if (_desc(val_a).class == INF) return _sign(val_a)?(-1):(1);
1383   if (_desc(val_b).class == INF) return _sign(val_b)?(1):(-1);
1384
1385   switch (sc_comp(_exp(val_a), _exp(val_b))) {
1386     case -1:
1387       return -1;
1388     case  1:
1389       return  1;
1390     case  0:
1391       return sc_comp(_mant(val_a), _mant(val_b));
1392     default:
1393       return 2;
1394   }
1395 }
1396
1397 int fc_is_zero(const void *a)
1398 {
1399   return _desc((const char*)a).class == ZERO;
1400 }
1401
1402 int fc_is_negative(const void *a)
1403 {
1404   return _sign((const char*)a);
1405 }
1406
1407 int fc_is_inf(const void *a)
1408 {
1409   return _desc(a).class == INF;
1410 }
1411
1412 int fc_is_nan(const void *a)
1413 {
1414   return _desc(a).class == NAN;
1415 }
1416
1417 int fc_is_subnormal(const void *a)
1418 {
1419   return _desc(a).class == SUBNORMAL;
1420 }
1421
1422 char *fc_print(const void *a, char *buf, int buflen, unsigned base)
1423 {
1424   const char *val;
1425   char *mul_1;
1426
1427   val = (const char*)a;
1428
1429   mul_1 = alloca(CALC_BUFFER_SIZE);
1430
1431   switch (base) {
1432     case FC_DEC:
1433       switch (_desc(val).class) {
1434         case INF:
1435           if (buflen >= 8+_sign(val)) sprintf(buf, "%sINFINITY", _sign(val)?"-":"");
1436           else snprintf(buf, buflen, "%sINF", _sign(val)?"-":NULL);
1437           break;
1438         case NAN:
1439           snprintf(buf, buflen, "NAN");
1440           break;
1441         case ZERO:
1442           snprintf(buf, buflen, "0.0");
1443           break;
1444         default:
1445           /* XXX to be implemented */
1446 #ifdef HAVE_LONG_DOUBLE
1447           /* XXX 30 is arbitrary */
1448           snprintf(buf, buflen, "%.30LE", fc_val_to_float(val));
1449 #else
1450           snprintf(buf, buflen, "%.18E", fc_val_to_float(val));
1451 #endif
1452       }
1453       break;
1454
1455     case FC_HEX:
1456       switch (_desc(val).class) {
1457         case INF:
1458           if (buflen >= 8+_sign(val)) sprintf(buf, "%sINFINITY", _sign(val)?"-":"");
1459           else snprintf(buf, buflen, "%sINF", _sign(val)?"-":NULL);
1460           break;
1461         case NAN:
1462           snprintf(buf, buflen, "NAN");
1463           break;
1464         case ZERO:
1465           snprintf(buf, buflen, "0.0");
1466           break;
1467         default:
1468 #ifdef HAVE_LONG_DOUBLE
1469           snprintf(buf, buflen, "%LA", fc_val_to_float(val));
1470 #else
1471           snprintf(buf, buflen, "%A", fc_val_to_float(val));
1472 #endif
1473       }
1474       break;
1475
1476     case FC_PACKED:
1477     default:
1478       snprintf(buf, buflen, "%s", sc_print(_pack(val, mul_1), VALUE_SIZE*4, SC_HEX));
1479       break;
1480   }
1481   return buf;
1482 }
1483
1484 unsigned char fc_sub_bits(const void *value, unsigned num_bits, unsigned byte_ofs)
1485 {
1486   /* this is used to cache the packed version of the value */
1487   static char *pack = NULL;
1488
1489   if (pack == NULL) pack = malloc(VALUE_SIZE);
1490
1491   if (value != NULL)
1492     _pack((const char*)value, pack);
1493
1494   return sc_sub_bits(pack, num_bits, byte_ofs);
1495 }
1496
1497 fc_rounding_mode_t fc_set_rounding_mode(fc_rounding_mode_t mode)
1498 {
1499   if (mode == FC_TONEAREST || mode == FC_TOPOSITIVE || mode == FC_TONEGATIVE || mode == FC_TOZERO)
1500       ROUNDING_MODE = mode;
1501
1502   return ROUNDING_MODE;
1503 }
1504
1505 fc_rounding_mode_t fc_get_rounding_mode(void)
1506 {
1507   return ROUNDING_MODE;
1508 }
1509
1510 void init_fltcalc(int precision)
1511 {
1512   if (calc_buffer == NULL) {
1513     /* does nothing if already init */
1514     if (precision == 0) precision = FC_DEFAULT_PRECISION;
1515
1516     init_strcalc(precision + 4);
1517
1518     /* needs additionally two bits to round, a bit as explicit 1., and one for
1519      * addition overflow */
1520     max_precision = sc_get_precision() - 4;
1521     if (max_precision < precision)
1522       printf("WARING: not enough precision available, using %d\n", max_precision);
1523
1524     ROUNDING_MODE = FC_TONEAREST;
1525     VALUE_SIZE = sc_get_buffer_length();
1526     SIGN_POS = 0;
1527     EXPONENT_POS = SIGN_POS + sizeof(char);
1528     MANTISSA_POS = EXPONENT_POS + VALUE_SIZE;
1529     DESCRIPTOR_POS = MANTISSA_POS + VALUE_SIZE;
1530     CALC_BUFFER_SIZE = DESCRIPTOR_POS + sizeof(descriptor_t);
1531
1532     calc_buffer = malloc(CALC_BUFFER_SIZE);
1533     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));
1534 #ifdef HAVE_LONG_DOUBLE
1535     DEBUGPRINTF(("\tUsing long double (1-15-64) interface\n"));
1536 #else
1537     DEBUGPRINTF(("\tUsing double (1-11-52) interface\n"));
1538 #endif
1539 #ifdef WORDS_BIGENDIAN
1540     DEBUGPRINTF(("\tWord order is big endian\n\n"));
1541 #else
1542     DEBUGPRINTF(("\tWord order is little endian\n\n"));
1543 #endif
1544   }
1545 }
1546
1547 /* definition of interface functions */
1548 FC_DEFINE2(add)
1549 FC_DEFINE2(sub)
1550 FC_DEFINE2(mul)
1551 FC_DEFINE2(div)
1552 FC_DEFINE1(neg)
1553 FC_DEFINE1(int)
1554 FC_DEFINE1(rnd)