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

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