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

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