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