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

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