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