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

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