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

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