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 mantissa 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 left shift 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 left shift 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  */
 811 char* _calc(const char *a, const char *b, int opcode, char *result)
 812 {
 813   char *temp;
 814 #ifdef FLTCALC_TRACE_CALC
 815   char *buffer;
 816
 817   buffer = alloca(100);
 818 #endif
 819
 820   if (result == NULL) result = calc_buffer;
 821
 822   TRACEPRINTF(("%s ", fc_print(a, buffer, 100, FC_PACKED)));
 823   switch (opcode)
 824   {
 825     case FC_add:
 826       /* make the value with the bigger exponent the first one */
 827       TRACEPRINTF(("+ %s ", fc_print(b, buffer, 100, FC_PACKED)));
 828       if (sc_comp(_exp(a), _exp(b)) == -1)
 829         _add(b, a, result);
 830       else
 831         _add(a, b, result);
 832       break;
 833     case FC_sub:
 834       TRACEPRINTF(("- %s ", fc_print(b, buffer, 100, FC_PACKED)));
 835       temp = alloca(calc_buffer_size);
 836       memcpy(temp, b, calc_buffer_size);
 837       _sign(temp) = !_sign(b);
 838       if (sc_comp(_exp(a), _exp(temp)) == -1)
 839         _add(temp, a, result);
 840       else
 841         _add(a, temp, result);
 842       break;
 843     case FC_mul:
 844       TRACEPRINTF(("* %s ", fc_print(b, buffer, 100, FC_PACKED)));
 845       _mul(a, b, result);
 846       break;
 847     case FC_div:
 848       TRACEPRINTF(("/ %s ", fc_print(b, buffer, 100, FC_PACKED)));
 849       _div(a, b, result);
 850       break;
 851     case FC_neg:
 852       TRACEPRINTF(("negated "));
 853       if (a != result) memcpy(result, a, calc_buffer_size);
 854       _sign(result) = !_sign(a);
 855       break;
 856     case FC_int:
 857       _trunc(a, result);
 858       break;
 859     case FC_rnd:
 860       break;
 861   }
 862
 863   TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_PACKED)));
 864   return result;
 865 }
 866
 867 /********
 868  * functions defined in fltcalc.h
 869  ********/
 870 const void *fc_get_buffer(void)
 871 {
 872   return calc_buffer;
 873 }
 874
 875 const int fc_get_buffer_length(void)
 876 {
 877   return calc_buffer_size;
 878 }
 879
 880 char* fc_val_from_str(const char *str, unsigned int len, char exp_size, char mant_size, char *result)
 881 {
 882 #if 0
 883   enum {
 884     START,
 885     LEFT_OF_DOT,
 886     RIGHT_OF_DOT,
 887     EXP_START,
 888     EXPONENT,
 889     END
 890   };
 891
 892   char exp_sign;
 893   int exp_int, hsb, state;
 894
 895   const char *old_str;
 896
 897   int pos;
 898   char *mant_str, *exp_val, *power_val;
 899
 900   if (result == NULL) result = calc_buffer;
 901
 902   exp_val = alloca(value_size);
 903   power_val = alloca(calc_buffer_size);
 904   mant_str = alloca((len)?(len):(strlen(str)));
 905
 906   _desc(result).exponent_size = exp_size;
 907   _desc(result).mantissa_size = mant_size;
 908   _desc(result).class = NORMAL;
 909
 910   old_str = str;
 911   pos = 0;
 912   exp_int = 0;
 913   state = START;
 914
 915   while (len == 0 || str-old_str < len)
 916   {
 917     switch (state) {
 918       case START:
 919         switch (*str) {
 920           case '+':
 921             _sign(result) = 0;
 922             state = LEFT_OF_DOT;
 923             str++;
 924             break;
 925
 926           case '-':
 927             _sign(result) = 1;
 928             state = LEFT_OF_DOT;
 929             str++;
 930             break;
 931
 932           case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
 933             _sign(result) = 0;
 934             state = LEFT_OF_DOT;
 935             break;
 936
 937           case '.':
 938             _sign(result) = 0;
 939             state = RIGHT_OF_DOT;
 940             str++;
 941             break;
 942
 943           case 'n':
 944           case 'N':
 945           case 'i':
 946           case 'I':
 947             break;
 948
 949           default:
 950             _fail_char(old_str, len, str - old_str);
 951         }
 952         break;
 953
 954       case LEFT_OF_DOT:
 955         switch (*str) {
 956           case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
 957             mant_str[pos++] = *(str++);
 958             break;
 959
 960           case '.':
 961             state = RIGHT_OF_DOT;
 962             str++;
 963             break;
 964
 965           case 'e':
 966           case 'E':
 967             state = EXP_START;
 968             str++;
 969             break;
 970
 971           case '\0':
 972             mant_str[pos] = '\0';
 973             goto done;
 974
 975           default:
 976             _fail_char(old_str, len, str - old_str);
 977         }
 978         break;
 979
 980       case RIGHT_OF_DOT:
 981         switch (*str) {
 982           case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
 983             mant_str[pos++] = *(str++);
 984             exp_int++;
 985             break;
 986
 987           case 'e':
 988           case 'E':
 989             state = EXP_START;
 990             str++;
 991             break;
 992
 993           case '\0':
 994             mant_str[pos] = '\0';
 995             goto done;
 996
 997           default:
 998             _fail_char(old_str, len, str - old_str);
 999         }
1000         break;
1001
1002       case EXP_START:
1003         switch (*str) {
1004           case '-':
1005             exp_sign = 1;
1006             /* fall through */
1007           case '+':
1008             if (*(str-1) != 'e' && *(str-1) != 'E') _fail_char(old_str, len, str - old_str);
1009             str++;
1010             break;
1011
1012           case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
1013             mant_str[pos] = '\0';
1014             pos = 1;
1015             str++;
1016             state = EXPONENT;
1017             break;
1018
1019           default:
1020             _fail_char(old_str, len, str - old_str);
1021         }
1022         break;
1023
1024       case EXPONENT:
1025         switch (*str) {
1026           case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
1027             pos++;
1028             str++;
1029             break;
1030
1031           case '\0': goto done;
1032
1033           default:
1034             _fail_char(old_str, len, str - old_str);
1035         }
1036     }
1037   } /*  switch(state) */
1038
1039 done:
1040   sc_val_from_str(mant_str, strlen(mant_str), _mant(result));
1041
1042   /* shift to put value left of radix point */
1043   sc_val_from_ulong(mant_size + 2, exp_val);
1044
1045   _shift_left(_mant(result), exp_val, _mant(result));
1046
1047   sc_val_from_ulong((1 << exp_size)/2-1, _exp(result));
1048
1049   _normalize(result, result, 0);
1050
1051   if (state == EXPONENT) {
1052     exp_int -= atoi(str-pos);
1053   }
1054
1055   _power_of_ten(exp_int, &_desc(result), power_val);
1056
1057   _div(result, power_val, result);
1058
1059   return result;
1060 #else
1061
1062   /* XXX excuse of an implementation to make things work */
1063   LLDBL val;
1064 #ifdef HAVE_LONG_DOUBLE
1065   val = strtold(str, NULL);
1066 #else
1067   val = strtod(str, NULL);
1068 #endif
1069
1070   DEBUGPRINTF(("val_from_str(%s)\n", str));
1071   return fc_val_from_float(val, exp_size, mant_size, result);
1072 #endif
1073 }
1074
1075 char* fc_val_from_float(LLDBL l, char exp_size, char mant_size, char* result)
1076 {
1077   char *temp;
1078   int bias_res, bias_val, mant_val;
1079   value_t srcval;
1080   UINT32 sign, exponent, mantissa0, mantissa1;
1081
1082   srcval.d = l;
1083   bias_res = ((1<<exp_size)/2-1);
1084
1085 #ifdef HAVE_LONG_DOUBLE
1086   mant_val  = 64;
1087   bias_val  = 0x3fff;
1088   sign      = (srcval.val.high & 0x00008000) != 0;
1089   exponent  = (srcval.val.high & 0x00007FFF) ;
1090   mantissa0 = srcval.val.mid;
1091   mantissa1 = srcval.val.low;
1092 #else /* no long double */
1093   mant_val  = 52;
1094   bias_val  = 0x3ff;
1095   sign      = (srcval.val.high & 0x80000000) != 0;
1096   exponent  = (srcval.val.high & 0x7FF00000) >> 20;
1097   mantissa0 = srcval.val.high & 0x000FFFFF;
1098   mantissa1 = srcval.val.low;
1099 #endif
1100
1101 #ifdef HAVE_LONG_DOUBLE
1102   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)); */
1103   DEBUGPRINTF(("(%d-%.4X-%.8X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1104 #else
1105   TRACEPRINTF(("val_from_float(%.8X%.8X)\n", srcval.val.high, srcval.val.low));
1106   DEBUGPRINTF(("(%d-%.3X-%.5X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1107 #endif
1108
1109   if (result == NULL) result = calc_buffer;
1110   temp = alloca(value_size);
1111
1112   _desc(result).exponent_size = exp_size;
1113   _desc(result).mantissa_size = mant_size;
1114
1115   /* extract sign */
1116   _sign(result) = sign;
1117
1118   /* sign and flag suffice to identify nan or inf, no exponent/mantissa
1119    * encoding is needed. the function can return immediately in these cases */
1120   if (isnan(l)) {
1121     _desc(result).class = NAN;
1122     TRACEPRINTF(("val_from_float resulted in NAN\n"));
1123     return result;
1124   }
1125   else if (isinf(l)) {
1126     _desc(result).class = INF;
1127     TRACEPRINTF(("val_from_float resulted in %sINF\n", (_sign(result)==1)?"-":""));
1128     return result;
1129   }
1130
1131   /* build exponent, because input and output exponent and mantissa sizes may differ
1132    * this looks more complicated than it is: unbiased input exponent + output bias,
1133    * minus the mantissa difference which is added again later when the output float
1134    * becomes normalized */
1135 #ifdef HAVE_EXPLICIT_ONE
1136   sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size-1), _exp(result));
1137 #else
1138   sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size), _exp(result));
1139 #endif
1140
1141   /* build mantissa representation */
1142 #ifndef HAVE_EXPLICIT_ONE
1143   if (exponent != 0)
1144   {
1145     /* insert the hidden bit */
1146     sc_val_from_ulong(1, temp);
1147     sc_val_from_ulong(mant_val + 2, NULL);
1148     _shift_left(temp, sc_get_buffer(), NULL);
1149   }
1150   else
1151 #endif
1152   {
1153     sc_val_from_ulong(0, NULL);
1154   }
1155
1156   _save_result(_mant(result));
1157
1158   /* bits from the upper word */
1159   sc_val_from_ulong(mantissa0, temp);
1160   sc_val_from_ulong(34, NULL);
1161   _shift_left(temp, sc_get_buffer(), temp);
1162   sc_or(_mant(result), temp, _mant(result));
1163
1164   /* bits from the lower word */
1165   sc_val_from_ulong(mantissa1, temp);
1166   sc_val_from_ulong(2, NULL);
1167   _shift_left(temp, sc_get_buffer(), temp);
1168   sc_or(_mant(result), temp, _mant(result));
1169
1170   /* _normalize expects the radix point to be normal, so shift mantissa of subnormal
1171    * origin one to the left */
1172   if (exponent == 0)
1173   {
1174     sc_val_from_ulong(1, NULL);
1175     _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1176   }
1177
1178   _normalize(result, result, 0);
1179
1180   TRACEPRINTF(("val_from_float results in %s\n", fc_print(result, temp, calc_buffer_size, FC_PACKED)));
1181
1182   return result;
1183 }
1184
1185 LLDBL fc_val_to_float(const void *val)
1186 {
1187   const char *value;
1188   char *temp = NULL;
1189
1190   int byte_offset;
1191
1192   UINT32 sign;
1193   UINT32 exponent;
1194   UINT32 mantissa0;
1195   UINT32 mantissa1;
1196
1197   value_t buildval;
1198
1199 #ifdef HAVE_LONG_DOUBLE
1200   char result_exponent = 15;
1201   char result_mantissa = 64;
1202 #else
1203   char result_exponent = 11;
1204   char result_mantissa = 52;
1205 #endif
1206
1207   temp = alloca(calc_buffer_size);
1208 #ifdef HAVE_EXPLICIT_ONE
1209   value = fc_cast(val, result_exponent, result_mantissa-1, temp);
1210 #else
1211   value = fc_cast(val, result_exponent, result_mantissa, temp);
1212 #endif
1213
1214   sign = _sign(value);
1215
1216   /* @@@ long double exponent is 15bit, so the use of sc_val_to_long should not
1217    * lead to wrong results */
1218   exponent = sc_val_to_long(_exp(value)) ;
1219
1220   sc_val_from_ulong(2, NULL);
1221   _shift_right(_mant(value), sc_get_buffer(), _mant(value));
1222
1223   mantissa0 = 0;
1224   mantissa1 = 0;
1225
1226   for (byte_offset = 0; byte_offset < 4; byte_offset++)
1227     mantissa1 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << (byte_offset<<3);
1228
1229   for (; (byte_offset<<3) < result_mantissa; byte_offset++)
1230     mantissa0 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << ((byte_offset-4)<<3);
1231
1232 #ifdef HAVE_LONG_DOUBLE
1233   buildval.val.high = sign << 15;
1234   buildval.val.high |= exponent;
1235   buildval.val.mid = mantissa0;
1236   buildval.val.low = mantissa1;
1237 #else /* no long double */
1238   mantissa0 &= 0x000FFFFF;  /* get rid of garbage */
1239   buildval.val.high = sign << 31;
1240   buildval.val.high |= exponent << 20;
1241   buildval.val.high |= mantissa0;
1242   buildval.val.low = mantissa1;
1243 #endif
1244
1245   TRACEPRINTF(("val_to_float: %d-%x-%x%x\n", sign, exponent, mantissa0, mantissa1));
1246   return buildval.d;
1247 }
1248
1249 char* fc_cast(const void *val, char exp_size, char mant_size, char *result)
1250 {
1251   const char *value = (const char*) val;
1252   char *temp;
1253   int exp_offset, val_bias, res_bias;
1254
1255   if (result == NULL) result = calc_buffer;
1256   temp = alloca(value_size);
1257
1258   if (_desc(value).exponent_size == exp_size && _desc(value).mantissa_size == mant_size)
1259   {
1260     if (value != result) memcpy(result, value, calc_buffer_size);
1261     return result;
1262   }
1263
1264   /* set the descriptor of the new value */
1265   _desc(result).exponent_size = exp_size;
1266   _desc(result).mantissa_size = mant_size;
1267   _desc(result).class = _desc(value).class;
1268
1269   _sign(result) = _sign(value);
1270
1271   /* when the mantissa sizes differ normalizing has to shift to align it.
1272    * this would change the exponent, which is unwanted. So calculate this
1273    * offset and add it */
1274   val_bias = (1<<_desc(value).exponent_size)/2-1;
1275   res_bias = (1<<exp_size)/2-1;
1276
1277   exp_offset = (res_bias - val_bias) - (_desc(value).mantissa_size - mant_size);
1278   sc_val_from_long(exp_offset, temp);
1279   sc_add(_exp(value), temp, _exp(result));
1280
1281   /* _normalize expects normalized radix point */
1282   if (_desc(val).class == SUBNORMAL) {
1283     sc_val_from_ulong(1, NULL);
1284     _shift_left(_mant(val), sc_get_buffer(), _mant(result));
1285   } else if (value != result) {
1286     memcpy(_mant(result), _mant(value), value_size);
1287   } else {
1288     memmove(_mant(result), _mant(value), value_size);
1289   }
1290
1291   _normalize(result, result, 0);
1292   TRACEPRINTF(("Cast results in %s\n", fc_print(result, temp, value_size, FC_PACKED)));
1293   return result;
1294 }
1295
1296 char* fc_get_max(unsigned int exponent_size, unsigned int mantissa_size, char* result)
1297 {
1298   if (result == NULL) result = calc_buffer;
1299
1300   _desc(result).exponent_size = exponent_size;
1301   _desc(result).mantissa_size = mantissa_size;
1302   _desc(result).class = NORMAL;
1303
1304   _sign(result) = 0;
1305
1306   sc_val_from_ulong((1<<exponent_size) - 2, _exp(result));
1307
1308   sc_max_from_bits(mantissa_size + 1, 0, _mant(result));
1309   sc_val_from_ulong(2, NULL);
1310   _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1311
1312   return result;
1313 }
1314
1315 char* fc_get_min(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1316 {
1317   if (result == NULL) result = calc_buffer;
1318
1319   fc_get_max(exponent_size, mantissa_size, result);
1320   _sign(result) = 1;
1321
1322   return result;
1323 }
1324
1325 char* fc_get_snan(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1326 {
1327   if (result == NULL) result = calc_buffer;
1328
1329   _desc(result).exponent_size = exponent_size;
1330   _desc(result).mantissa_size = mantissa_size;
1331   _desc(result).class = NAN;
1332
1333   _sign(result) = 0;
1334
1335   sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1336
1337   /* signaling NaN has non-zero mantissa with msb not set */
1338   sc_val_from_ulong(1, _mant(result));
1339
1340   return result;
1341 }
1342
1343 char* fc_get_qnan(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1344 {
1345   if (result == NULL) result = calc_buffer;
1346
1347   _desc(result).exponent_size = exponent_size;
1348   _desc(result).mantissa_size = mantissa_size;
1349   _desc(result).class = NAN;
1350
1351   _sign(result) = 0;
1352
1353   sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1354
1355   /* quiet NaN has the msb of the mantissa set, so shift one there */
1356   sc_val_from_ulong(1, _mant(result));
1357   /* mantissa_size >+< 1 because of two extra rounding bits */
1358   sc_val_from_ulong(mantissa_size + 1, NULL);
1359   _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1360
1361   return result;
1362 }
1363
1364 char* fc_get_plusinf(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1365 {
1366   if (result == NULL) result = calc_buffer;
1367
1368   _desc(result).exponent_size = exponent_size;
1369   _desc(result).mantissa_size = mantissa_size;
1370   _desc(result).class = NORMAL;
1371
1372   _sign(result) = 0;
1373
1374   sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1375
1376   sc_val_from_ulong(0, _mant(result));
1377
1378   return result;
1379 }
1380
1381 char* fc_get_minusinf(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1382 {
1383   if (result == NULL) result = calc_buffer;
1384
1385   fc_get_plusinf(exponent_size, mantissa_size, result);
1386   _sign(result) = 1;
1387
1388   return result;
1389 }
1390
1391 int fc_comp(const void *a, const void *b)
1392 {
1393   const char *val_a = (const char*)a;
1394   const char *val_b = (const char*)b;
1395   int mul = 1;
1396
1397   /*
1398    * shortcut: if both values are identical, they are either
1399    * Unordered if NaN or equal
1400    */
1401   if (a == b)
1402     return _desc(val_a).class == NAN ? 2 : 0;
1403
1404   /* unordered if one is a NaN */
1405   if (_desc(val_a).class == NAN || _desc(val_b).class == NAN)
1406     return 2;
1407
1408   /* zero is equal independent of sign */
1409   if ((_desc(val_a).class == ZERO) && (_desc(val_b).class == ZERO))
1410     return 0;
1411
1412   /* different signs make compare easy */
1413   if (_sign(val_a) != _sign(val_b))
1414     return (_sign(val_a)==0)?(1):(-1);
1415
1416   mul = _sign(a) ? -1 : 1;
1417
1418   /* both infinity means equality */
1419   if ((_desc(val_a).class == INF) && (_desc(val_b).class == INF))
1420     return 0;
1421
1422   /* infinity is bigger than the rest */
1423   if (_desc(val_a).class == INF)
1424     return  1 * mul;
1425   if (_desc(val_b).class == INF)
1426     return -1 * mul;
1427
1428   /* check first exponent, that mantissa if equal */
1429   switch (sc_comp(_exp(val_a), _exp(val_b))) {
1430     case -1:
1431       return -1 * mul;
1432     case  1:
1433       return  1 * mul;
1434     case  0:
1435       return sc_comp(_mant(val_a), _mant(val_b)) * mul;
1436     default:
1437       return 2;
1438   }
1439 }
1440
1441 int fc_is_zero(const void *a)
1442 {
1443   return _desc(a).class == ZERO;
1444 }
1445
1446 int fc_is_negative(const void *a)
1447 {
1448   return _sign(a);
1449 }
1450
1451 int fc_is_inf(const void *a)
1452 {
1453   return _desc(a).class == INF;
1454 }
1455
1456 int fc_is_nan(const void *a)
1457 {
1458   return _desc(a).class == NAN;
1459 }
1460
1461 int fc_is_subnormal(const void *a)
1462 {
1463   return _desc(a).class == SUBNORMAL;
1464 }
1465
1466 char *fc_print(const void *a, char *buf, int buflen, unsigned base)
1467 {
1468   const char *val;
1469   char *mul_1;
1470
1471   val = (const char*)a;
1472
1473   mul_1 = alloca(calc_buffer_size);
1474
1475   switch (base) {
1476     case FC_DEC:
1477       switch (_desc(val).class) {
1478         case INF:
1479           if (buflen >= 8+_sign(val)) sprintf(buf, "%sINFINITY", _sign(val)?"-":"");
1480           else snprintf(buf, buflen, "%sINF", _sign(val)?"-":NULL);
1481           break;
1482         case NAN:
1483           snprintf(buf, buflen, "NAN");
1484           break;
1485         case ZERO:
1486           snprintf(buf, buflen, "0.0");
1487           break;
1488         default:
1489           /* XXX to be implemented */
1490 #ifdef HAVE_LONG_DOUBLE
1491           /* XXX 30 is arbitrary */
1492           snprintf(buf, buflen, "%.30LE", fc_val_to_float(val));
1493 #else
1494           snprintf(buf, buflen, "%.18E", fc_val_to_float(val));
1495 #endif
1496       }
1497       break;
1498
1499     case FC_HEX:
1500       switch (_desc(val).class) {
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 #ifdef HAVE_LONG_DOUBLE
1513           snprintf(buf, buflen, "%LA", fc_val_to_float(val));
1514 #else
1515           snprintf(buf, buflen, "%A", fc_val_to_float(val));
1516 #endif
1517       }
1518       break;
1519
1520     case FC_PACKED:
1521     default:
1522       snprintf(buf, buflen, "%s", sc_print(_pack(val, mul_1), value_size*4, SC_HEX));
1523       break;
1524   }
1525   return buf;
1526 }
1527
1528 unsigned char fc_sub_bits(const void *value, unsigned num_bits, unsigned byte_ofs)
1529 {
1530   /* this is used to cache the packed version of the value */
1531   static char *pack = NULL;
1532
1533   if (pack == NULL) pack = xmalloc(value_size);
1534
1535   if (value != NULL)
1536     _pack((const char*)value, pack);
1537
1538   return sc_sub_bits(pack, num_bits, byte_ofs);
1539 }
1540
1541 fc_rounding_mode_t fc_set_rounding_mode(fc_rounding_mode_t mode)
1542 {
1543   if (mode == FC_TONEAREST || mode == FC_TOPOSITIVE || mode == FC_TONEGATIVE || mode == FC_TOZERO)
1544       rounding_mode = mode;
1545
1546   return rounding_mode;
1547 }
1548
1549 fc_rounding_mode_t fc_get_rounding_mode(void)
1550 {
1551   return rounding_mode;
1552 }
1553
1554 void init_fltcalc(int precision)
1555 {
1556   if (calc_buffer == NULL) {
1557     /* does nothing if already init */
1558     if (precision == 0) precision = FC_DEFAULT_PRECISION;
1559
1560     init_strcalc(precision + 4);
1561
1562     /* needs additionally two bits to round, a bit as explicit 1., and one for
1563      * addition overflow */
1564     max_precision = sc_get_precision() - 4;
1565     if (max_precision < precision)
1566       printf("WARING: not enough precision available, using %d\n", max_precision);
1567
1568     rounding_mode = FC_TONEAREST;
1569     value_size = sc_get_buffer_length();
1570     SIGN_POS = 0;
1571     EXPONENT_POS = SIGN_POS + sizeof(char);
1572     MANTISSA_POS = EXPONENT_POS + value_size;
1573     DESCRIPTOR_POS = MANTISSA_POS + value_size;
1574     calc_buffer_size = DESCRIPTOR_POS + sizeof(descriptor_t);
1575
1576     calc_buffer = xmalloc(calc_buffer_size);
1577     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));
1578 #ifdef HAVE_LONG_DOUBLE
1579     DEBUGPRINTF(("\tUsing long double (1-15-64) interface\n"));
1580 #else
1581     DEBUGPRINTF(("\tUsing double (1-11-52) interface\n"));
1582 #endif
1583 #ifdef WORDS_BIGENDIAN
1584     DEBUGPRINTF(("\tWord order is big endian\n\n"));
1585 #else
1586     DEBUGPRINTF(("\tWord order is little endian\n\n"));
1587 #endif
1588   }
1589 }
1590
1591 void finish_fltcalc (void) {
1592   free(calc_buffer); calc_buffer = NULL;
1593 }
1594
1595 /* definition of interface functions */
1596 FC_DEFINE2(add)
1597 FC_DEFINE2(sub)
1598 FC_DEFINE2(mul)
1599 FC_DEFINE2(div)
1600 FC_DEFINE1(neg)
1601 FC_DEFINE1(int)
1602 FC_DEFINE1(rnd)