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