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

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