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

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