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         return;
 822 }
 823
 824 /********
 825  * functions defined in fltcalc.h
 826  ********/
 827 const void *fc_get_buffer(void) {
 828         return calc_buffer;
 829 }
 830
 831 int fc_get_buffer_length(void) {
 832         return calc_buffer_size;
 833 }
 834
 835 void *fc_val_from_str(const char *str, unsigned int len, char exp_size, char mant_size, void *result) {
 836 #if 0
 837         enum {
 838                 START,
 839                 LEFT_OF_DOT,
 840                 RIGHT_OF_DOT,
 841                 EXP_START,
 842                 EXPONENT,
 843                 END
 844         };
 845
 846         char exp_sign;
 847         int exp_int, hsb, state;
 848
 849         const char *old_str;
 850
 851         int pos;
 852         char *mant_str, *exp_val, *power_val;
 853
 854         (void) len;
 855         if (result == NULL) result = calc_buffer;
 856
 857         exp_val = alloca(value_size);
 858         power_val = alloca(calc_buffer_size);
 859         mant_str = alloca((len)?(len):(strlen(str)));
 860
 861         result->desc.exponent_size = exp_size;
 862         result->desc.mantissa_size = mant_size;
 863         result->desc.clss = NORMAL;
 864
 865         old_str = str;
 866         pos = 0;
 867         exp_int = 0;
 868         state = START;
 869
 870         while (len == 0 || str-old_str < len) {
 871                 switch (state) {
 872                 case START:
 873                         switch (*str) {
 874                         case '+':
 875                                 result->sign = 0;
 876                                 state = LEFT_OF_DOT;
 877                                 str++;
 878                                 break;
 879
 880                         case '-':
 881                                 result->sign = 1;
 882                                 state = LEFT_OF_DOT;
 883                                 str++;
 884                                 break;
 885
 886                         case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
 887                                 result->sign = 0;
 888                                 state = LEFT_OF_DOT;
 889                                 break;
 890
 891                         case '.':
 892                                 result->sign = 0;
 893                                 state = RIGHT_OF_DOT;
 894                                 str++;
 895                                 break;
 896
 897                         case 'n':
 898                         case 'N':
 899                         case 'i':
 900                         case 'I':
 901                                 break;
 902
 903                         default:
 904                                 fail_char(old_str, len, str - old_str);
 905                         }
 906                         break;
 907
 908                 case LEFT_OF_DOT:
 909                         switch (*str) {
 910                         case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
 911                                 mant_str[pos++] = *(str++);
 912                                 break;
 913
 914                         case '.':
 915                                 state = RIGHT_OF_DOT;
 916                                 str++;
 917                                 break;
 918
 919                         case 'e':
 920                         case 'E':
 921                                 state = EXP_START;
 922                                 str++;
 923                                 break;
 924
 925                         case '\0':
 926                                 mant_str[pos] = '\0';
 927                                 goto done;
 928
 929                         default:
 930                                 fail_char(old_str, len, str - old_str);
 931                         }
 932                         break;
 933
 934                 case RIGHT_OF_DOT:
 935                         switch (*str) {
 936                         case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
 937                                 mant_str[pos++] = *(str++);
 938                                 exp_int++;
 939                                 break;
 940
 941                         case 'e':
 942                         case 'E':
 943                                 state = EXP_START;
 944                                 str++;
 945                                 break;
 946
 947                         case '\0':
 948                                 mant_str[pos] = '\0';
 949                                 goto done;
 950
 951                         default:
 952                                 fail_char(old_str, len, str - old_str);
 953                         }
 954                         break;
 955
 956                 case EXP_START:
 957                         switch (*str) {
 958                         case '-':
 959                                 exp_sign = 1;
 960                                 /* fall through */
 961                         case '+':
 962                                 if (*(str-1) != 'e' && *(str-1) != 'E') fail_char(old_str, len, str - old_str);
 963                                 str++;
 964                                 break;
 965
 966                         case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
 967                                 mant_str[pos] = '\0';
 968                                 pos = 1;
 969                                 str++;
 970                                 state = EXPONENT;
 971                                 break;
 972
 973                         default:
 974                                 fail_char(old_str, len, str - old_str);
 975                         }
 976                         break;
 977
 978                 case EXPONENT:
 979                         switch (*str) {
 980                         case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
 981                                 pos++;
 982                                 str++;
 983                                 break;
 984
 985                         case '\0': goto done;
 986
 987                         default:
 988                                 fail_char(old_str, len, str - old_str);
 989                         }
 990                 }
 991         } /*  switch(state) */
 992
 993 done:
 994         sc_val_from_str(mant_str, strlen(mant_str), _mant(result));
 995
 996         /* shift to put value left of radix point */
 997         sc_val_from_ulong(mant_size + ROUNDING_BITS, exp_val);
 998
 999         _shift_left(_mant(result), exp_val, _mant(result));
1000
1001         sc_val_from_ulong((1 << (exp_size - 1)) - 1, _exp(result));
1002
1003         _normalize(result, result, 0);
1004
1005         if (state == EXPONENT) {
1006                 exp_int -= atoi(str-pos);
1007         }
1008
1009         _power_of_ten(exp_int, &result->desc, power_val);
1010
1011         _fdiv(result, power_val, result);
1012
1013         return result;
1014 #else
1015         /* XXX excuse of an implementation to make things work */
1016         LLDBL val;
1017         fp_value *tmp = alloca(calc_buffer_size);
1018         (void) len;
1019
1020 #ifdef HAVE_LONG_DOUBLE
1021         val = strtold(str, NULL);
1022         DEBUGPRINTF(("val_from_str(%s)\n", str));
1023         fc_val_from_ieee754(val, 15, 64, tmp);
1024 #else
1025         val = strtod(str, NULL);
1026         DEBUGPRINTF(("val_from_str(%s)\n", str));
1027         fc_val_from_ieee754(val, 11, 52, tmp);
1028 #endif /* HAVE_LONG_DOUBLE */
1029         return fc_cast(tmp, exp_size, mant_size, result);
1030 #endif
1031 }
1032
1033 fp_value *fc_val_from_ieee754(LLDBL l, char exp_size, char mant_size, fp_value *result) {
1034         char *temp;
1035         int bias_res, bias_val, mant_val;
1036         value_t srcval;
1037         UINT32 sign, exponent, mantissa0, mantissa1;
1038
1039         srcval.d = l;
1040         bias_res = ((1 << (exp_size - 1)) - 1);
1041
1042 #ifdef HAVE_LONG_DOUBLE
1043         mant_val  = 64;
1044         bias_val  = 0x3fff;
1045         sign      = (srcval.val.high & 0x00008000) != 0;
1046         exponent  = (srcval.val.high & 0x00007FFF) ;
1047         mantissa0 = srcval.val.mid;
1048         mantissa1 = srcval.val.low;
1049 #else /* no long double */
1050         mant_val  = 52;
1051         bias_val  = 0x3ff;
1052         sign      = (srcval.val.high & 0x80000000) != 0;
1053         exponent  = (srcval.val.high & 0x7FF00000) >> 20;
1054         mantissa0 = srcval.val.high & 0x000FFFFF;
1055         mantissa1 = srcval.val.low;
1056 #endif
1057
1058 #ifdef HAVE_LONG_DOUBLE
1059         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)); */
1060         DEBUGPRINTF(("(%d-%.4X-%.8X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1061 #else
1062         TRACEPRINTF(("val_from_float(%.8X%.8X)\n", srcval.val.high, srcval.val.low));
1063         DEBUGPRINTF(("(%d-%.3X-%.5X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1064 #endif
1065
1066         if (result == NULL) result = calc_buffer;
1067         temp = alloca(value_size);
1068
1069         /* CLEAR the buffer, else some bits might be uninitialised */
1070         memset(result, 0, fc_get_buffer_length());
1071
1072         result->desc.exponent_size = exp_size;
1073         result->desc.mantissa_size = mant_size;
1074
1075         /* extract sign */
1076         result->sign = sign;
1077
1078         /* sign and flag suffice to identify nan or inf, no exponent/mantissa
1079          * encoding is needed. the function can return immediately in these cases */
1080         if (isnan(l)) {
1081                 result->desc.clss = NAN;
1082                 TRACEPRINTF(("val_from_float resulted in NAN\n"));
1083                 return result;
1084         }
1085         else if (isinf(l)) {
1086                 result->desc.clss = INF;
1087                 TRACEPRINTF(("val_from_float resulted in %sINF\n", (result->sign == 1) ? "-" : ""));
1088                 return result;
1089         }
1090
1091         /* build exponent, because input and output exponent and mantissa sizes may differ
1092          * this looks more complicated than it is: unbiased input exponent + output bias,
1093          * minus the mantissa difference which is added again later when the output float
1094          * becomes normalized */
1095 #ifdef HAVE_EXPLICIT_ONE
1096         sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size-1), _exp(result));
1097 #else
1098         sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size), _exp(result));
1099 #endif
1100
1101         /* build mantissa representation */
1102 #ifndef HAVE_EXPLICIT_ONE
1103         if (exponent != 0) {
1104                 /* insert the hidden bit */
1105                 sc_val_from_ulong(1, temp);
1106                 sc_val_from_ulong(mant_val + ROUNDING_BITS, NULL);
1107                 _shift_left(temp, sc_get_buffer(), NULL);
1108         }
1109         else
1110 #endif
1111         {
1112                 sc_val_from_ulong(0, NULL);
1113         }
1114
1115         _save_result(_mant(result));
1116
1117         /* bits from the upper word */
1118         sc_val_from_ulong(mantissa0, temp);
1119         sc_val_from_ulong(34, NULL);
1120         _shift_left(temp, sc_get_buffer(), temp);
1121         sc_or(_mant(result), temp, _mant(result));
1122
1123         /* bits from the lower word */
1124         sc_val_from_ulong(mantissa1, temp);
1125         sc_val_from_ulong(ROUNDING_BITS, NULL);
1126         _shift_left(temp, sc_get_buffer(), temp);
1127         sc_or(_mant(result), temp, _mant(result));
1128
1129         /* _normalize expects the radix point to be normal, so shift mantissa of subnormal
1130          * origin one to the left */
1131         if (exponent == 0) {
1132                 sc_val_from_ulong(1, NULL);
1133                 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1134         }
1135
1136         normalize(result, result, 0);
1137
1138         TRACEPRINTF(("val_from_float results in %s\n", fc_print(result, temp, calc_buffer_size, FC_PACKED)));
1139
1140         return result;
1141 }
1142
1143 LLDBL fc_val_to_ieee754(const fp_value *val) {
1144         fp_value *value;
1145         fp_value *temp = NULL;
1146
1147         int byte_offset;
1148
1149         UINT32 sign;
1150         UINT32 exponent;
1151         UINT32 mantissa0;
1152         UINT32 mantissa1;
1153
1154         value_t buildval;
1155
1156 #ifdef HAVE_LONG_DOUBLE
1157         char result_exponent = 15;
1158         char result_mantissa = 64;
1159 #else
1160         char result_exponent = 11;
1161         char result_mantissa = 52;
1162 #endif
1163
1164         temp = alloca(calc_buffer_size);
1165 #ifdef HAVE_EXPLICIT_ONE
1166         value = fc_cast(val, result_exponent, result_mantissa-1, temp);
1167 #else
1168         value = fc_cast(val, result_exponent, result_mantissa, temp);
1169 #endif
1170
1171         sign = value->sign;
1172
1173         /* @@@ long double exponent is 15bit, so the use of sc_val_to_long should not
1174          * lead to wrong results */
1175         exponent = sc_val_to_long(_exp(value)) ;
1176
1177         sc_val_from_ulong(ROUNDING_BITS, NULL);
1178         _shift_right(_mant(value), sc_get_buffer(), _mant(value));
1179
1180         mantissa0 = 0;
1181         mantissa1 = 0;
1182
1183         for (byte_offset = 0; byte_offset < 4; byte_offset++)
1184                 mantissa1 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << (byte_offset<<3);
1185
1186         for (; (byte_offset<<3) < result_mantissa; byte_offset++)
1187                 mantissa0 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << ((byte_offset-4)<<3);
1188
1189 #ifdef HAVE_LONG_DOUBLE
1190         buildval.val.high = sign << 15;
1191         buildval.val.high |= exponent;
1192         buildval.val.mid = mantissa0;
1193         buildval.val.low = mantissa1;
1194 #else /* no long double */
1195         mantissa0 &= 0x000FFFFF;  /* get rid of garbage */
1196         buildval.val.high = sign << 31;
1197         buildval.val.high |= exponent << 20;
1198         buildval.val.high |= mantissa0;
1199         buildval.val.low = mantissa1;
1200 #endif
1201
1202         TRACEPRINTF(("val_to_float: %d-%x-%x%x\n", sign, exponent, mantissa0, mantissa1));
1203         return buildval.d;
1204 }
1205
1206 fp_value *fc_cast(const fp_value *value, char exp_size, char mant_size, fp_value *result) {
1207         char *temp;
1208         int exp_offset, val_bias, res_bias;
1209
1210         if (result == NULL) result = calc_buffer;
1211         temp = alloca(value_size);
1212
1213         if (value->desc.exponent_size == exp_size && value->desc.mantissa_size == mant_size) {
1214                 if (value != result)
1215                         memcpy(result, value, calc_buffer_size);
1216                 return result;
1217         }
1218
1219         if (value->desc.clss == NAN) {
1220                 if (sc_get_highest_set_bit(_mant(value)) == value->desc.mantissa_size + 1)
1221                         return fc_get_qnan(exp_size, mant_size, result);
1222                 else
1223                         return fc_get_snan(exp_size, mant_size, result);
1224         }
1225
1226         /* set the descriptor of the new value */
1227         result->desc.exponent_size = exp_size;
1228         result->desc.mantissa_size = mant_size;
1229         result->desc.clss = value->desc.clss;
1230
1231         result->sign = value->sign;
1232
1233         /* when the mantissa sizes differ normalizing has to shift to align it.
1234          * this would change the exponent, which is unwanted. So calculate this
1235          * offset and add it */
1236         val_bias = (1 << (value->desc.exponent_size - 1)) - 1;
1237         res_bias = (1 << (exp_size - 1)) - 1;
1238
1239         exp_offset = (res_bias - val_bias) - (value->desc.mantissa_size - mant_size);
1240         sc_val_from_long(exp_offset, temp);
1241         sc_add(_exp(value), temp, _exp(result));
1242
1243         /* _normalize expects normalized radix point */
1244         if (value->desc.clss == SUBNORMAL) {
1245                 sc_val_from_ulong(1, NULL);
1246                 _shift_left(_mant(value), sc_get_buffer(), _mant(result));
1247         } else if (value != result) {
1248                 memcpy(_mant(result), _mant(value), value_size);
1249         } else {
1250                 memmove(_mant(result), _mant(value), value_size);
1251         }
1252
1253         normalize(result, result, 0);
1254         TRACEPRINTF(("Cast results in %s\n", fc_print(result, temp, value_size, FC_PACKED)));
1255         return result;
1256 }
1257
1258 fp_value *fc_get_max(unsigned int exponent_size, unsigned int mantissa_size, fp_value *result) {
1259         if (result == NULL) result = calc_buffer;
1260
1261         result->desc.exponent_size = exponent_size;
1262         result->desc.mantissa_size = mantissa_size;
1263         result->desc.clss = NORMAL;
1264
1265         result->sign = 0;
1266
1267         sc_val_from_ulong((1<<exponent_size) - 2, _exp(result));
1268
1269         sc_max_from_bits(mantissa_size + 1, 0, _mant(result));
1270         sc_val_from_ulong(ROUNDING_BITS, NULL);
1271         _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1272
1273         return result;
1274 }
1275
1276 fp_value *fc_get_min(unsigned int exponent_size, unsigned int mantissa_size, fp_value *result) {
1277         if (result == NULL) result = calc_buffer;
1278
1279         fc_get_max(exponent_size, mantissa_size, result);
1280         result->sign = 1;
1281
1282         return result;
1283 }
1284
1285 fp_value *fc_get_snan(unsigned int exponent_size, unsigned int mantissa_size, fp_value *result) {
1286         if (result == NULL) result = calc_buffer;
1287
1288         result->desc.exponent_size = exponent_size;
1289         result->desc.mantissa_size = mantissa_size;
1290         result->desc.clss = NAN;
1291
1292         result->sign = 0;
1293
1294         sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1295
1296         /* signaling NaN has non-zero mantissa with msb not set */
1297         sc_val_from_ulong(1, _mant(result));
1298
1299         return result;
1300 }
1301
1302 fp_value *fc_get_qnan(unsigned int exponent_size, unsigned int mantissa_size, fp_value *result) {
1303         if (result == NULL) result = calc_buffer;
1304
1305         result->desc.exponent_size = exponent_size;
1306         result->desc.mantissa_size = mantissa_size;
1307         result->desc.clss = NAN;
1308
1309         result->sign = 0;
1310
1311         sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1312
1313         /* quiet NaN has the msb of the mantissa set, so shift one there */
1314         sc_val_from_ulong(1, _mant(result));
1315         /* mantissa_size >+< 1 because of two extra rounding bits */
1316         sc_val_from_ulong(mantissa_size + 1, NULL);
1317         _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1318
1319         return result;
1320 }
1321
1322 fp_value *fc_get_plusinf(unsigned int exponent_size, unsigned int mantissa_size, fp_value *result) {
1323         if (result == NULL) result = calc_buffer;
1324
1325         result->desc.exponent_size = exponent_size;
1326         result->desc.mantissa_size = mantissa_size;
1327         result->desc.clss = NORMAL;
1328
1329         result->sign = 0;
1330
1331         sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1332
1333         sc_val_from_ulong(0, _mant(result));
1334
1335         return result;
1336 }
1337
1338 fp_value *fc_get_minusinf(unsigned int exponent_size, unsigned int mantissa_size, fp_value *result) {
1339         if (result == NULL) result = calc_buffer;
1340
1341         fc_get_plusinf(exponent_size, mantissa_size, result);
1342         result->sign = 1;
1343
1344         return result;
1345 }
1346
1347 int fc_comp(const fp_value *val_a, const fp_value *val_b) {
1348         int mul = 1;
1349
1350         /*
1351          * shortcut: if both values are identical, they are either
1352          * Unordered if NaN or equal
1353          */
1354         if (val_a == val_b)
1355                 return val_a->desc.clss == NAN ? 2 : 0;
1356
1357         /* unordered if one is a NaN */
1358         if (val_a->desc.clss == NAN || val_b->desc.clss == NAN)
1359                 return 2;
1360
1361         /* zero is equal independent of sign */
1362         if ((val_a->desc.clss == ZERO) && (val_b->desc.clss == ZERO))
1363                 return 0;
1364
1365         /* different signs make compare easy */
1366         if (val_a->sign != val_b->sign)
1367                 return (val_a->sign == 0) ? (1) : (-1);
1368
1369         mul = val_a->sign ? -1 : 1;
1370
1371         /* both infinity means equality */
1372         if ((val_a->desc.clss == INF) && (val_b->desc.clss == INF))
1373                 return 0;
1374
1375         /* infinity is bigger than the rest */
1376         if (val_a->desc.clss == INF)
1377                 return  1 * mul;
1378         if (val_b->desc.clss == INF)
1379                 return -1 * mul;
1380
1381         /* check first exponent, that mantissa if equal */
1382         switch (sc_comp(_exp(val_a), _exp(val_b))) {
1383         case -1:
1384                 return -1 * mul;
1385         case  1:
1386                 return  1 * mul;
1387         case  0:
1388                 return sc_comp(_mant(val_a), _mant(val_b)) * mul;
1389         default:
1390                 return 2;
1391         }
1392 }
1393
1394 int fc_is_zero(const fp_value *a) {
1395         return a->desc.clss == ZERO;
1396 }
1397
1398 int fc_is_negative(const fp_value *a) {
1399         return a->sign;
1400 }
1401
1402 int fc_is_inf(const fp_value *a) {
1403         return a->desc.clss == INF;
1404 }
1405
1406 int fc_is_nan(const fp_value *a) {
1407         return a->desc.clss == NAN;
1408 }
1409
1410 int fc_is_subnormal(const fp_value *a) {
1411         return a->desc.clss == SUBNORMAL;
1412 }
1413
1414 char *fc_print(const fp_value *val, char *buf, int buflen, unsigned base) {
1415         char *mul_1;
1416
1417         mul_1 = alloca(calc_buffer_size);
1418
1419         switch (base) {
1420         case FC_DEC:
1421                 switch (val->desc.clss) {
1422                 case INF:
1423                         if (buflen >= 8 + val->sign) sprintf(buf, "%sINFINITY", val->sign ? "-":"");
1424                         else snprintf(buf, buflen, "%sINF", val->sign ? "-":NULL);
1425                         break;
1426                 case NAN:
1427                         snprintf(buf, buflen, "NAN");
1428                         break;
1429                 case ZERO:
1430                         snprintf(buf, buflen, "0.0");
1431                         break;
1432                 default:
1433                         /* XXX to be implemented */
1434 #ifdef HAVE_LONG_DOUBLE
1435                         /* XXX 30 is arbitrary */
1436                         snprintf(buf, buflen, "%.30LE", fc_val_to_ieee754(val));
1437 #else
1438                         snprintf(buf, buflen, "%.18E", fc_val_to_ieee754(val));
1439 #endif
1440                 }
1441                 break;
1442
1443         case FC_HEX:
1444                 switch (val->desc.clss) {
1445                 case INF:
1446                         if (buflen >= 8+val->sign) sprintf(buf, "%sINFINITY", val->sign?"-":"");
1447                         else snprintf(buf, buflen, "%sINF", val->sign?"-":NULL);
1448                         break;
1449                 case NAN:
1450                         snprintf(buf, buflen, "NAN");
1451                         break;
1452                 case ZERO:
1453                         snprintf(buf, buflen, "0.0");
1454                         break;
1455                 default:
1456 #ifdef HAVE_LONG_DOUBLE
1457                         snprintf(buf, buflen, "%LA", fc_val_to_ieee754(val));
1458 #else
1459                         snprintf(buf, buflen, "%A", fc_val_to_ieee754(val));
1460 #endif
1461                 }
1462                 break;
1463
1464         case FC_PACKED:
1465         default:
1466                 snprintf(buf, buflen, "%s", sc_print(pack(val, mul_1), value_size*4, SC_HEX, 0));
1467                 buf[buflen - 1] = '\0';
1468                 break;
1469         }
1470         return buf;
1471 }
1472
1473 unsigned char fc_sub_bits(const fp_value *value, unsigned num_bits, unsigned byte_ofs) {
1474         /* this is used to cache the packed version of the value */
1475         static char *packed_value = NULL;
1476
1477         if (packed_value == NULL) packed_value = xmalloc(value_size);
1478
1479         if (value != NULL)
1480                 pack(value, packed_value);
1481
1482         return sc_sub_bits(packed_value, num_bits, byte_ofs);
1483 }
1484
1485 /* Returns non-zero if the mantissa is zero, i.e. 1.0Exxx */
1486 int fc_zero_mantissa(const fp_value *value) {
1487         return sc_get_lowest_set_bit(_mant(value)) == ROUNDING_BITS + value->desc.mantissa_size;
1488 }
1489
1490 /* Returns the exponent of a value. */
1491 int fc_get_exponent(const fp_value *value) {
1492         int exp_bias = (1 << (value->desc.exponent_size - 1)) - 1;
1493         return sc_val_to_long(_exp(value)) - exp_bias;
1494 }
1495
1496 /* Return non-zero if a given value can be converted lossless into another precision */
1497 int fc_can_lossless_conv_to(const fp_value *value, char exp_size, char mant_size) {
1498         int v;
1499         int exp_bias;
1500
1501         /* handle some special cases first */
1502         switch (value->desc.clss) {
1503         case ZERO:
1504         case INF:
1505         case NAN:
1506                 return 1;
1507         default:
1508                 break;
1509         }
1510
1511         /* check if the exponent can be encoded: note, 0 and all ones are reserved for the exponent */
1512         exp_bias = (1 << (exp_size - 1)) - 1;
1513         v = fc_get_exponent(value) + exp_bias;
1514         if (0 < v && v < (1 << exp_size) - 1) {
1515                 /* check the mantissa */
1516                 v = value->desc.mantissa_size + ROUNDING_BITS - sc_get_lowest_set_bit(_mant(value));
1517                 return v < mant_size;
1518         }
1519         return 0;
1520 }
1521
1522
1523 fc_rounding_mode_t fc_set_rounding_mode(fc_rounding_mode_t mode) {
1524         if (mode == FC_TONEAREST || mode == FC_TOPOSITIVE || mode == FC_TONEGATIVE || mode == FC_TOZERO)
1525                 rounding_mode = mode;
1526
1527         return rounding_mode;
1528 }
1529
1530 fc_rounding_mode_t fc_get_rounding_mode(void) {
1531         return rounding_mode;
1532 }
1533
1534 void init_fltcalc(int precision) {
1535         if (calc_buffer == NULL) {
1536                 /* does nothing if already init */
1537                 if (precision == 0) precision = FC_DEFAULT_PRECISION;
1538
1539                 init_strcalc(precision + 4);
1540
1541                 /* needs additionally two bits to round, a bit as explicit 1., and one for
1542                  * addition overflow */
1543                 max_precision = sc_get_precision() - 4;
1544                 if (max_precision < precision)
1545                         printf("WARNING: not enough precision available, using %d\n", max_precision);
1546
1547                 rounding_mode    = FC_TONEAREST;
1548                 value_size       = sc_get_buffer_length();
1549                 calc_buffer_size = sizeof(fp_value) + 2*value_size - 1;
1550
1551                 calc_buffer = xmalloc(calc_buffer_size);
1552                 memset(calc_buffer, 0, calc_buffer_size);
1553                 DEBUGPRINTF(("init fltcalc:\n\tVALUE_SIZE = %d\ntCALC_BUFFER_SIZE = %d\n\tcalc_buffer = %p\n\n", value_size, calc_buffer_size, calc_buffer));
1554 #ifdef HAVE_LONG_DOUBLE
1555                 DEBUGPRINTF(("\tUsing long double (1-15-64) interface\n"));
1556 #else
1557                 DEBUGPRINTF(("\tUsing double (1-11-52) interface\n"));
1558 #endif
1559 #ifdef WORDS_BIGENDIAN
1560                 DEBUGPRINTF(("\tWord order is big endian\n\n"));
1561 #else
1562                 DEBUGPRINTF(("\tWord order is little endian\n\n"));
1563 #endif
1564         }
1565 }
1566
1567 void finish_fltcalc (void) {
1568         free(calc_buffer); calc_buffer = NULL;
1569 }
1570
1571 #ifdef FLTCALC_TRACE_CALC
1572 static char buffer[100];
1573 #endif
1574
1575 /* definition of interface functions */
1576 fp_value *fc_add(const fp_value *a, const fp_value *b, fp_value *result) {
1577         if (result == NULL) result = calc_buffer;
1578
1579         TRACEPRINTF(("%s ", fc_print(a, buffer, sizeof(buffer), FC_PACKED)));
1580         TRACEPRINTF(("+ %s ", fc_print(b, buffer, sizeof(buffer), FC_PACKED)));
1581
1582         /* make the value with the bigger exponent the first one */
1583         if (sc_comp(_exp(a), _exp(b)) == -1)
1584                 _fadd(b, a, result);
1585         else
1586                 _fadd(a, b, result);
1587
1588         TRACEPRINTF(("= %s\n", fc_print(result, buffer, sizeof(buffer), FC_PACKED)));
1589         return result;
1590 }
1591
1592 fp_value *fc_sub(const fp_value *a, const fp_value *b, fp_value *result) {
1593         fp_value *temp;
1594
1595         if (result == NULL) result = calc_buffer;
1596
1597         TRACEPRINTF(("%s ", fc_print(a, buffer, sizeof(buffer), FC_PACKED)));
1598         TRACEPRINTF(("- %s ", fc_print(b, buffer, sizeof(buffer), FC_PACKED)));
1599
1600         temp = alloca(calc_buffer_size);
1601         memcpy(temp, b, calc_buffer_size);
1602         temp->sign = !b->sign;
1603         if (sc_comp(_exp(a), _exp(temp)) == -1)
1604                 _fadd(temp, a, result);
1605         else
1606                 _fadd(a, temp, result);
1607
1608         TRACEPRINTF(("= %s\n", fc_print(result, buffer, sizeof(buffer), FC_PACKED)));
1609         return result;
1610 }
1611
1612 fp_value *fc_mul(const fp_value *a, const fp_value *b, fp_value *result) {
1613         if (result == NULL) result = calc_buffer;
1614
1615         TRACEPRINTF(("%s ", fc_print(a, buffer, sizeof(buffer), FC_PACKED)));
1616         TRACEPRINTF(("* %s ", fc_print(b, buffer, sizeof(buffer), FC_PACKED)));
1617
1618         _fmul(a, b, result);
1619
1620         TRACEPRINTF(("= %s\n", fc_print(result, buffer, sizeof(buffer), FC_PACKED)));
1621         return result;
1622 }
1623
1624 fp_value *fc_div(const fp_value *a, const fp_value *b, fp_value *result) {
1625         if (result == NULL) result = calc_buffer;
1626
1627         TRACEPRINTF(("%s ", fc_print(a, buffer, sizeof(buffer), FC_PACKED)));
1628         TRACEPRINTF(("/ %s ", fc_print(b, buffer, sizeof(buffer), FC_PACKED)));
1629
1630         _fdiv(a, b, result);
1631
1632         TRACEPRINTF(("= %s\n", fc_print(result, buffer, sizeof(buffer), FC_PACKED)));
1633         return result;
1634 }
1635
1636 fp_value *fc_neg(const fp_value *a, fp_value *result) {
1637         if (result == NULL) result = calc_buffer;
1638
1639         TRACEPRINTF(("- %s ", fc_print(a, buffer, sizeof(buffer), FC_PACKED)));
1640
1641         if (a != result)
1642                 memcpy(result, a, calc_buffer_size);
1643         result->sign = !a->sign;
1644
1645         TRACEPRINTF(("= %s\n", fc_print(result, buffer, sizeof(buffer), FC_PACKED)));
1646         return result;
1647 }
1648
1649 fp_value *fc_int(const fp_value *a, fp_value *result) {
1650         if (result == NULL) result = calc_buffer;
1651
1652         TRACEPRINTF(("%s ", fc_print(a, buffer, sizeof(buffer), FC_PACKED)));
1653         TRACEPRINTF(("truncated to integer "));
1654
1655         _trunc(a, result);
1656
1657         TRACEPRINTF(("= %s\n", fc_print(result, buffer, sizeof(buffer), FC_PACKED)));
1658         return result;
1659 }
1660
1661 fp_value *fc_rnd(const fp_value *a, fp_value *result) {
1662         if (result == NULL) result = calc_buffer;
1663
1664         (void) a;
1665         TRACEPRINTF(("%s ", fc_print(a, buffer, sizeof(buffer), FC_PACKED)));
1666         TRACEPRINTF(("rounded to integer "));
1667
1668         assert(!"fc_rnd() not yet implemented");
1669
1670         TRACEPRINTF(("= %s\n", fc_print(result, buffer, sizeof(buffer), FC_PACKED)));
1671         return result;
1672 }
1673
1674 /*
1675  * convert a floating point value into an sc value ...
1676  */
1677 int fc_flt2int(const fp_value *a, void *result, ir_mode *dst_mode) {
1678         if (a->desc.clss == NORMAL) {
1679                 int exp_bias = (1 << (a->desc.exponent_size - 1)) - 1;
1680                 int exp_val  = sc_val_to_long(_exp(a)) - exp_bias;
1681                 int shift, highest;
1682
1683                 if (a->sign && !mode_is_signed(dst_mode)) {
1684                         /* FIXME: for now we cannot convert this */
1685                         return 0;
1686                 }
1687
1688                 assert(exp_val >= 0 && "floating point value not integral before fc_flt2int() call");
1689                 shift = exp_val - (a->desc.mantissa_size + ROUNDING_BITS);
1690
1691                 if (shift > 0) {
1692                         sc_shlI(_mant(a),  shift, 64, 0, result);
1693                 } else {
1694                         sc_shrI(_mant(a), -shift, 64, 0, result);
1695                 }
1696
1697                 /* check for overflow */
1698                 highest = sc_get_highest_set_bit(result);
1699
1700                 if (mode_is_signed(dst_mode)) {
1701                         if (highest == sc_get_lowest_set_bit(result)) {
1702                                 /* need extra test for MIN_INT */
1703                                 if (highest >= (int) get_mode_size_bits(dst_mode)) {
1704                                         /* FIXME: handle overflow */
1705                                         return 0;
1706                                 }
1707                         } else {
1708                                 if (highest >= (int) get_mode_size_bits(dst_mode) - 1) {
1709                                         /* FIXME: handle overflow */
1710                                         return 0;
1711                                 }
1712                         }
1713                 } else {
1714                         if (highest >= (int) get_mode_size_bits(dst_mode)) {
1715                                 /* FIXME: handle overflow */
1716                                 return 0;
1717                         }
1718                 }
1719
1720                 if (a->sign)
1721                         sc_neg(result, result);
1722
1723                 return 1;
1724         }
1725         else if (a->desc.clss == ZERO) {
1726                 sc_zero(result);
1727                 return 1;
1728         }
1729         return 0;
1730 }
1731
1732
1733 unsigned fc_set_immediate_precision(unsigned bits) {
1734         unsigned old = immediate_prec;
1735
1736         immediate_prec = bits;
1737         return old;
1738 }
1739
1740 int fc_is_exact(void) {
1741         return fc_exact;
1742 }