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