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