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

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