fltcalc.c 42.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
/*
 * Project:     libFIRM
 * File name:   ir/tv/fltcalc.c
 * Purpose:
 * Author:
 * Modified by:
 * Created:     2003
 * CVS-ID:      $Id$
 * Copyright:   (c) 2003 Universitt Karlsruhe
 * Licence:     This file protected by GPL -  GNU GENERAL PUBLIC LICENSE.
11
12
 */

13

Boris Boesler's avatar
Boris Boesler committed
14
15
16
17
18
#ifdef HAVE_CONFIG_H
# include "config.h"
#endif


19
#include "fltcalc.h"
20
21
22
23
24
25
26
27
#include "strcalc.h"

#include <math.h>    /* need isnan() and isinf() (will be changed)*/
/* undef some reused constants defined by math.h */
#ifdef NAN
#  undef NAN
#endif

28
#ifdef HAVE_INTTYPES_H
Michael Beck's avatar
Michael Beck committed
29
30
31
32
33
34
35
# include <inttypes.h>
#endif
#ifdef HAVE_STRING_H
# include <string.h>
#endif
#ifdef HAVE_STDLIB_H
# include <stdlib.h>
36
#endif
Boris Boesler's avatar
Boris Boesler committed
37
#ifdef HAVE_ALLOCA_H
Boris Boesler's avatar
Boris Boesler committed
38
39
# include <alloca.h>
#endif
40
41
42
#ifdef HAVE_MALLOC_H
# include <malloc.h>
#endif
Michael Beck's avatar
Michael Beck committed
43
44
#include <stdio.h>
#include <assert.h>
45

46
47
#include "xmalloc.h"

48
typedef uint32_t UINT32;
Michael Beck's avatar
Tarval:    
Michael Beck committed
49

50
#ifdef HAVE_LONG_DOUBLE
51
#ifdef WORDS_BIGENDIAN
52
53
54
typedef union {
  struct {
    UINT32 high;
55
56
    UINT32 mid;
    UINT32 low;
57
58
59
60
61
62
63
  } val;
  volatile long double d;
} value_t;
#else
typedef union {
  struct {
    UINT32 low;
64
65
    UINT32 mid;
    UINT32 high;
66
67
68
69
70
  } val;
  volatile long double d;
} value_t;
#endif
#else
71
#ifdef WORDS_BIGENDIAN
72
73
74
typedef union {
  struct {
    UINT32 high;
75
    UINT32 low;
76
77
78
79
80
81
82
  } val;
  volatile double d;
} value_t;
#else
typedef union {
  struct {
    UINT32 low;
83
    UINT32 high;
84
85
86
87
88
  } val;
  volatile double d;
} value_t;
#endif
#endif
89

Michael Beck's avatar
Michael Beck committed
90
91
92
/**
 * possible float states
 */
93
typedef enum {
Michael Beck's avatar
Michael Beck committed
94
95
96
97
98
  NORMAL,       /**< normal representation, implicit 1 */
  ZERO,         /**< +/-0 */
  SUBNORMAL,    /**< denormals, implicit 0 */
  INF,          /**< +/-oo */
  NAN,          /**< Not A Number */
99
100
} value_class_t;

Michael Beck's avatar
Michael Beck committed
101
/** A descriptor for an IEEE float value. */
102
typedef struct {
Michael Beck's avatar
Michael Beck committed
103
104
105
  unsigned char exponent_size;    /**< size of exponent in bits */
  unsigned char mantissa_size;    /**< size of mantissa in bits */
  value_class_t clss;             /**< state of this float */
106
107
} descriptor_t;

108
#define CLEAR_BUFFER(buffer) memset(buffer, 0, calc_buffer_size)
109
110
111
112
113

/* because variable sized structs are impossible, the internal
 * value is represented as a pseudo-struct char array, addressed
 * by macros
 * struct {
114
 *   char sign;             //  0 for positive, 1 for negative
115
116
 *   char exp[value_size];
 *   char mant[value_size];
117
118
119
120
121
122
 *   descriptor_t desc;
 * };
 */
#define _sign(a) (((char*)a)[SIGN_POS])
#define _exp(a) (&((char*)a)[EXPONENT_POS])
#define _mant(a) (&((char*)a)[MANTISSA_POS])
Michael Beck's avatar
Michael Beck committed
123
#define _desc(a) (*(descriptor_t *)&((char *)a)[DESCRIPTOR_POS])
Michael Beck's avatar
Tarval:    
Michael Beck committed
124

125
126
127
#define _save_result(x) memcpy((x), sc_get_buffer(), value_size)
#define _shift_right(x, y, b) sc_shr((x), (y), value_size*4, 0, (b))
#define _shift_left(x, y, b) sc_shl((x), (y), value_size*4, 0, (b))
Michael Beck's avatar
Tarval:    
Michael Beck committed
128

Michael Beck's avatar
Michael Beck committed
129
130
131
132
#define FC_DEFINE1(code) \
char *fc_##code(const void *a, void *result)  {                        \
  return _calc((const char*)a, NULL, FC_##code, (char*)result);        \
}
Michael Beck's avatar
Tarval:    
Michael Beck committed
133

Michael Beck's avatar
Michael Beck committed
134
135
136
137
#define FC_DEFINE2(code) \
char *fc_##code(const void *a, const void *b, void *result) {             \
  return _calc((const char*)a, (const char*)b, FC_##code, (char*)result); \
}
138
139
140
141
142
143
144

#define FUNC_PTR(code) fc_##code

#if FLTCALC_DEBUG
#  define DEBUGPRINTF(x) printf x
#else
#  define DEBUGPRINTF(x) ((void)0)
Michael Beck's avatar
Tarval:    
Michael Beck committed
145
146
#endif

147
148
149
150
151
#if FLTCALC_TRACE_CALC
#  define TRACEPRINTF(x) printf x
#else
#  define TRACEPRINTF(x) ((void)0)
#endif
152

153
static char *calc_buffer = NULL;
154

155
static fc_rounding_mode_t rounding_mode;
156

157
158
static int calc_buffer_size;
static int value_size;
159
160
161
162
static int SIGN_POS;
static int EXPONENT_POS;
static int MANTISSA_POS;
static int DESCRIPTOR_POS;
163

164
static int max_precision;
165
/********
166
 * private functions
167
 ********/
168
169
#if 0
static void _fail_char(const char *str, unsigned int len, int pos)
170
{
171
172
173
174
175
176
177
178
  if (*(str+pos))
    printf("ERROR: Unexpected character '%c'\n", *(str + pos));
  else
    printf("ERROR: Unexpected end of string\n");
  while (len-- && *str) printf("%c", *str++); printf("\n");
  while (pos--) printf(" "); printf("^\n");
  /* the front end has to to check constant strings */
  exit(-1);
179
}
180
#endif
181

Michael Beck's avatar
Michael Beck committed
182
/** pack machine-like */
183
static char* _pack(const char *int_float, char *packed)
184
{
185
186
187
  char *shift_val;
  char *temp;
  char *val_buffer;
188

189
190
  temp = alloca(value_size);
  shift_val = alloca(value_size);
191

Michael Beck's avatar
Michael Beck committed
192
  switch (_desc(int_float).clss) {
193
    case NAN:
194
      val_buffer = alloca(calc_buffer_size);
195
196
197
198
199
      fc_get_qnan(_desc(int_float).exponent_size, _desc(int_float).mantissa_size, val_buffer);
      int_float = val_buffer;
      break;

    case INF:
200
      val_buffer = alloca(calc_buffer_size);
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
      fc_get_plusinf(_desc(int_float).exponent_size, _desc(int_float).mantissa_size, val_buffer);
      _sign(val_buffer) = _sign(int_float);
      int_float = val_buffer;
      break;

    default:
      break;
  }
  /* pack sign */
  sc_val_from_ulong(_sign(int_float), temp);

  sc_val_from_ulong(_desc(int_float).exponent_size + _desc(int_float).mantissa_size, NULL);
  _shift_left(temp, sc_get_buffer(), packed);

  /* extract exponent */
  sc_val_from_ulong(_desc(int_float).mantissa_size, shift_val);

  _shift_left(_exp(int_float), shift_val, temp);

  sc_or(temp, packed, packed);

  /* extract mantissa */
  /* remove 2 rounding bits */
  sc_val_from_ulong(2, shift_val);
  _shift_right(_mant(int_float), shift_val, temp);

  /* remove leading 1 (or 0 if denormalized) */
  sc_max_from_bits(_desc(int_float).mantissa_size, 0, shift_val); /* all mantissa bits are 1's */
  sc_and(temp, shift_val, temp);

  /* save result */
  sc_or(temp, packed, packed);

  return packed;
235
236
}

237
char* _normalize(const char *in_val, char *out_val, int sticky)
238
{
239
240
241
242
  int hsb;
  char lsb, guard, round, round_dir = 0;
  char *temp;

243
  temp = alloca(value_size);
244
245
246
247
248
249
250
251
252
253

  /* +2: save two rounding bits at the end */
  hsb = 2 + _desc(in_val).mantissa_size - sc_get_highest_set_bit(_mant(in_val)) - 1;

  if (in_val != out_val)
  {
    _sign(out_val) = _sign(in_val);
    memcpy(&_desc(out_val), &_desc(in_val), sizeof(descriptor_t));
  }

Michael Beck's avatar
Michael Beck committed
254
  _desc(out_val).clss = NORMAL;
255

256
  /* mantissa all zeros, so zero exponent (because of explicit one)*/
257
258
259
260
261
262
  if (hsb == 2 + _desc(in_val).mantissa_size)
  {
    sc_val_from_ulong(0, _exp(out_val));
    hsb = -1;
  }

263
  /* shift the first 1 into the left of the radix point (i.e. hsb == -1) */
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
  if (hsb < -1)
  {
    /* shift right */
    sc_val_from_ulong(-hsb-1, temp);

    _shift_right(_mant(in_val), temp, _mant(out_val));

    /* remember if some bits were shifted away */
    if (!sticky) sticky = sc_had_carry();

    sc_add(_exp(in_val), temp, _exp(out_val));
  }
  else if (hsb > -1)
  {
    /* shift left */
    sc_val_from_ulong(hsb+1, temp);

    _shift_left(_mant(in_val), temp, _mant(out_val));

    sc_sub(_exp(in_val), temp, _exp(out_val));
  }

  /* check for exponent underflow */
  if (sc_is_negative(_exp(out_val)) || sc_is_zero(_exp(out_val))) {
    DEBUGPRINTF(("Exponent underflow!\n"));
    /* exponent underflow */
    /* shift the mantissa right to have a zero exponent */
    sc_val_from_ulong(1, temp);
    sc_sub(temp, _exp(out_val), NULL);

    _shift_right(_mant(out_val), sc_get_buffer(), _mant(out_val));
    if (!sticky) sticky = sc_had_carry();
    /* denormalized means exponent of zero */
    sc_val_from_ulong(0, _exp(out_val));

Michael Beck's avatar
Michael Beck committed
299
    _desc(out_val).clss = SUBNORMAL;
300
301
302
303
304
305
306
307
308
  }

  /* perform rounding by adding a value that clears the guard bit and the round bit
   * and either causes a carry to round up or not */
  /* get the last 3 bits of the value */
  lsb = sc_sub_bits(_mant(out_val), _desc(out_val).mantissa_size + 2, 0) & 0x7;
  guard = (lsb&0x2)>>1;
  round = lsb&0x1;

309
  switch (rounding_mode)
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
  {
    case FC_TONEAREST:
      /* round to nearest representable value, if in doubt choose the version
       * with lsb == 0 */
      round_dir = guard && (sticky || round || lsb>>2);
      break;
    case FC_TOPOSITIVE:
      /* if positive: round to one if the exact value is bigger, else to zero */
      round_dir = (!_sign(out_val) && (guard || round || sticky));
      break;
    case FC_TONEGATIVE:
      /* if negative: round to one if the exact value is bigger, else to zero */
      round_dir = (_sign(out_val) && (guard || round || sticky));
      break;
    case FC_TOZERO:
      /* always round to 0 (chopping mode) */
      round_dir = 0;
      break;
  }
  DEBUGPRINTF(("Rounding (s%d, l%d, g%d, r%d, s%d) %s\n", _sign(out_val), lsb>>2, guard, round, sticky, (round_dir)?"up":"down"));

  if (round_dir == 1)
  {
    guard = (round^guard)<<1;
    lsb = !(round || guard)<<2 | guard | round;
  }
  else
  {
    lsb = -((guard<<1) | round);
  }

  /* add the rounded value */
  if (lsb != 0) {
    sc_val_from_long(lsb, temp);
    sc_add(_mant(out_val), temp, _mant(out_val));
  }

  /* could have rounded down to zero */
Michael Beck's avatar
Michael Beck committed
348
349
  if (sc_is_zero(_mant(out_val)) && (_desc(out_val).clss == SUBNORMAL))
    _desc(out_val).clss = ZERO;
350
351
352

  /* check for rounding overflow */
  hsb = 2 + _desc(out_val).mantissa_size - sc_get_highest_set_bit(_mant(out_val)) - 1;
Michael Beck's avatar
Michael Beck committed
353
  if ((_desc(out_val).clss != SUBNORMAL) && (hsb < -1))
354
355
356
357
358
359
  {
    sc_val_from_ulong(1, temp);
    _shift_right(_mant(out_val), temp, _mant(out_val));

    sc_add(_exp(out_val), temp, _exp(out_val));
  }
Michael Beck's avatar
Michael Beck committed
360
  else if ((_desc(out_val).clss == SUBNORMAL) && (hsb == -1))
361
  {
Michael Beck's avatar
Michael Beck committed
362
    /* overflow caused the mantissa to be normal again,
363
364
365
366
     * so adapt the exponent accordingly */
    sc_val_from_ulong(1, temp);
    sc_add(_exp(out_val), temp, _exp(out_val));

Michael Beck's avatar
Michael Beck committed
367
    _desc(out_val).clss = NORMAL;
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
  }
  /* no further rounding is needed, because rounding overflow means
   * the carry of the original rounding was propagated all the way
   * up to the bit left of the radix point. This implies the bits
   * to the right are all zeros (rounding is +1) */

  /* check for exponent overflow */
  sc_val_from_ulong((1 << _desc(out_val).exponent_size) - 1, temp);
  if (sc_comp(_exp(out_val), temp) != -1) {
    DEBUGPRINTF(("Exponent overflow!\n"));
    /* exponent overflow, reaction depends on rounding method:
     *
     * mode        | sign of value |  result
     *--------------------------------------------------------------
     * TO_NEAREST  |      +        |   +inf
     *             |      -        |   -inf
     *--------------------------------------------------------------
     * TO_POSITIVE |      +        |   +inf
     *             |      -        |   smallest representable value
     *--------------------------------------------------------------
     * TO_NEAGTIVE |      +        |   largest representable value
     *             |      -        |   -inf
     *--------------------------------------------------------------
     * TO_ZERO     |      +        |   largest representable value
     *             |      -        |   smallest representable value
     *--------------------------------------------------------------*/
    if (_sign(out_val) == 0)
    {
      /* value is positive */
397
      switch (rounding_mode) {
398
399
        case FC_TONEAREST:
        case FC_TOPOSITIVE:
Michael Beck's avatar
Michael Beck committed
400
          _desc(out_val).clss = INF;
401
402
403
404
405
406
407
408
          break;

        case FC_TONEGATIVE:
        case FC_TOZERO:
          fc_get_max(_desc(out_val).exponent_size, _desc(out_val).mantissa_size, out_val);
      }
    } else {
      /* value is negative */
409
      switch (rounding_mode) {
410
411
        case FC_TONEAREST:
        case FC_TONEGATIVE:
Michael Beck's avatar
Michael Beck committed
412
          _desc(out_val).clss = INF;
413
414
415
416
417
418
419
420
421
422
          break;

        case FC_TOPOSITIVE:
        case FC_TOZERO:
          fc_get_min(_desc(out_val).exponent_size, _desc(out_val).mantissa_size, out_val);
      }
    }
  }

  return out_val;
423
424
}

Michael Beck's avatar
Michael Beck committed
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
/**
 * Operations involving NaN's must return NaN
 */
#define handle_NAN(a, b, result) \
do {                                                      \
  if (_desc(a).clss == NAN) {                             \
    if (a != result) memcpy(result, a, calc_buffer_size); \
    return result;                                        \
  }                                                       \
  if (_desc(b).clss == NAN) {                             \
    if (b != result) memcpy(result, b, calc_buffer_size); \
    return result;                                        \
  }                                                       \
}while (0)


/**
442
443
 * calculate a + b, where a is the value with the bigger exponent
 */
Michael Beck's avatar
Michael Beck committed
444
static char* _fadd(const char* a, const char* b, char* result)
445
{
446
447
448
449
450
451
  char *temp;
  char *exp_diff;

  char sign;
  char sticky;

Michael Beck's avatar
Michael Beck committed
452
  handle_NAN(a, b, result);
453
454
455
456
457
458
459
460

  /* make sure result has a descriptor */
  if (result != a && result != b)
    memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));

  /* determine if this is an addition or subtraction */
  sign = _sign(a) ^ _sign(b);

Michael Beck's avatar
Michael Beck committed
461
  /* produce NaN on inf - inf */
Michael Beck's avatar
Michael Beck committed
462
  if (sign && (_desc(a).clss == INF) && (_desc(b).clss == INF))
463
464
    return fc_get_qnan(_desc(a).exponent_size, _desc(b).mantissa_size, result);

Michael Beck's avatar
Michael Beck committed
465
  temp     = alloca(value_size);
466
  exp_diff = alloca(value_size);
467
468
469
470
471
472
473
474

  /* get exponent difference */
  sc_sub(_exp(a), _exp(b), exp_diff);

  /* initially set sign to be the sign of a, special treatment of subtraction
   * when exponents are equal is required though.
   * Also special care about the sign is needed when the mantissas are equal
   * (+/- 0 ?) */
Michael Beck's avatar
Michael Beck committed
475
  if (sign && sc_val_to_long(exp_diff) == 0) {
476
477
    switch (sc_comp(_mant(a), _mant(b))) {
      case 1:  /* a > b */
Michael Beck's avatar
Michael Beck committed
478
        _sign(result) = _sign(a);  /* abs(a) is bigger and a is negative */
479
480
        break;
      case 0:  /* a == b */
Michael Beck's avatar
Michael Beck committed
481
        _sign(result) = (rounding_mode == FC_TONEGATIVE);
482
483
        break;
      case -1: /* a < b */
Michael Beck's avatar
Michael Beck committed
484
        _sign(result) = _sign(b); /* abs(b) is bigger and b is negative */
485
486
487
488
489
490
        break;
      default:
        /* can't be reached */
        break;
    }
  }
Michael Beck's avatar
Michael Beck committed
491
492
  else
    _sign(result) = _sign(a);
493
494

  /* sign has been taken care of, check for special cases */
Michael Beck's avatar
Michael Beck committed
495
  if (_desc(a).clss == ZERO) {
496
    if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
497
498
    return result;
  }
Michael Beck's avatar
Michael Beck committed
499
  if (_desc(b).clss == ZERO) {
500
    if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
501
502
503
    return result;
  }

Michael Beck's avatar
Michael Beck committed
504
  if (_desc(a).clss == INF) {
505
    if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
506
507
    return result;
  }
Michael Beck's avatar
Michael Beck committed
508
  if (_desc(b).clss == INF) {
509
    if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
510
511
512
513
514
515
    return result;
  }

  /* shift the smaller value to the right to align the radix point */
  /* subnormals have their radix point shifted to the right,
   * take care of this first */
Michael Beck's avatar
Michael Beck committed
516
  if ((_desc(b).clss == SUBNORMAL) && (_desc(a).clss != SUBNORMAL))
517
518
519
520
521
522
523
524
525
526
527
528
529
  {
    sc_val_from_ulong(1, temp);
    sc_sub(exp_diff, temp, exp_diff);
  }

  _shift_right(_mant(b), exp_diff, temp);
  sticky = sc_had_carry();

  if (sticky && sign)
  {
    /* if subtracting a little more than the represented value or adding a little
     * more than the represented value to a negative value this, in addition to the
     * still set sticky bit, takes account of the 'little more' */
530
    char *temp1 = alloca(calc_buffer_size);
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
    sc_val_from_ulong(1, temp1);
    sc_add(temp, temp1, temp);
  }

  if (sign) {
    if (sc_comp(_mant(a), temp) == -1)
      sc_sub(temp, _mant(a), _mant(result));
    else
      sc_sub(_mant(a), temp, _mant(result));
  } else {
    sc_add(_mant(a), temp, _mant(result));
  }

  /* _normalize expects a 'normal' radix point, adding two subnormals
   * results in a subnormal radix point -> shifting before normalizing */
Michael Beck's avatar
Michael Beck committed
546
  if ((_desc(a).clss == SUBNORMAL) && (_desc(b).clss == SUBNORMAL))
547
548
549
550
551
552
  {
    sc_val_from_ulong(1, NULL);
    _shift_left(_mant(result), sc_get_buffer(), _mant(result));
  }

  /* resulting exponent is the bigger one */
553
  memmove(_exp(result), _exp(a), value_size);
554
555

  return _normalize(result, result, sticky);
556
557
}

Michael Beck's avatar
Michael Beck committed
558
559
560
561
/**
 * calculate a * b
 */
static char* _fmul(const char* a, const char* b, char* result)
562
{
563
  char *temp;
564

Michael Beck's avatar
Michael Beck committed
565
  handle_NAN(a, b, result);
566

567
  temp = alloca(value_size);
568
569
570
571
572
573

  if (result != a && result != b)
    memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));

  _sign(result) = _sign(a) ^ _sign(b);

Michael Beck's avatar
Michael Beck committed
574
  /* produce NaN on 0 * inf */
Michael Beck's avatar
Michael Beck committed
575
576
  if (_desc(a).clss == ZERO) {
    if (_desc(b).clss == INF)
577
578
      fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
    else
579
      if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
580
581
    return result;
  }
Michael Beck's avatar
Michael Beck committed
582
583
  if (_desc(b).clss == ZERO) {
    if (_desc(a).clss == INF)
584
585
      fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
    else
586
      if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-1);
587
588
589
    return result;
  }

Michael Beck's avatar
Michael Beck committed
590
  if (_desc(a).clss == INF) {
591
    if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
592
593
    return result;
  }
Michael Beck's avatar
Michael Beck committed
594
  if (_desc(b).clss == INF) {
595
    if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-1);
596
597
598
599
600
601
602
603
604
605
    return result;
  }

  /* exp = exp(a) + exp(b) - excess */
  sc_add(_exp(a), _exp(b), _exp(result));

  sc_val_from_ulong((1<<_desc(a).exponent_size)/2-1, temp);
  sc_sub(_exp(result), temp, _exp(result));

  /* mixed normal, subnormal values introduce an error of 1, correct it */
Michael Beck's avatar
Michael Beck committed
606
  if ((_desc(a).clss == SUBNORMAL) ^ (_desc(b).clss == SUBNORMAL))
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
  {
    sc_val_from_ulong(1, temp);
    sc_add(_exp(result), temp, _exp(result));
  }

  sc_mul(_mant(a), _mant(b), _mant(result));

  /* realign result: after a multiplication the digits right of the radix
   * point are the sum of the factors' digits after the radix point. As all
   * values are normalized they both have the same amount of these digits,
   * which has to be restored by proper shifting
   * +2 because of the two rounding bits */
  sc_val_from_ulong(2 + _desc(result).mantissa_size, temp);

  _shift_right(_mant(result), temp, _mant(result));

  return _normalize(result, result, sc_had_carry());
624
625
}

Michael Beck's avatar
Michael Beck committed
626
627
628
629
/**
 * calculate a / b
 */
static char* _fdiv(const char* a, const char* b, char* result)
630
{
631
632
  char *temp, *dividend;

Michael Beck's avatar
Michael Beck committed
633
  handle_NAN(a, b, result);
634

635
636
  temp = alloca(value_size);
  dividend = alloca(value_size);
637
638
639
640
641
642
643

  if (result != a && result != b)
    memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));

  _sign(result) = _sign(a) ^ _sign(b);

  /* produce nan on 0/0 and inf/inf */
Michael Beck's avatar
Michael Beck committed
644
645
  if (_desc(a).clss == ZERO) {
    if (_desc(b).clss == ZERO)
646
647
648
649
      /* 0/0 -> nan */
      fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
    else
      /* 0/x -> a */
650
      if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
651
652
653
    return result;
  }

Michael Beck's avatar
Michael Beck committed
654
655
  if (_desc(b).clss == INF) {
    if (_desc(a).clss == INF)
656
657
658
659
660
661
662
      /* inf/inf -> nan */
      fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
    else {
      /* x/inf -> 0 */
      sc_val_from_ulong(0, NULL);
      _save_result(_exp(result));
      _save_result(_mant(result));
Michael Beck's avatar
Michael Beck committed
663
      _desc(result).clss = ZERO;
664
665
666
667
    }
    return result;
  }

Michael Beck's avatar
Michael Beck committed
668
  if (_desc(a).clss == INF) {
669
    /* inf/x -> inf */
670
    if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
671
672
    return result;
  }
Michael Beck's avatar
Michael Beck committed
673
  if (_desc(b).clss == ZERO) {
674
675
676
677
678
679
680
681
682
683
684
685
686
687
    /* division by zero */
    if (_sign(result))
      fc_get_minusinf(_desc(a).exponent_size, _desc(a).mantissa_size, result);
    else
      fc_get_plusinf(_desc(a).exponent_size, _desc(a).mantissa_size, result);
    return result;
  }

  /* exp = exp(a) - exp(b) + excess - 1*/
  sc_sub(_exp(a), _exp(b), _exp(result));
  sc_val_from_ulong((1 << _desc(a).exponent_size)/2-2, temp);
  sc_add(_exp(result), temp, _exp(result));

  /* mixed normal, subnormal values introduce an error of 1, correct it */
Michael Beck's avatar
Michael Beck committed
688
  if ((_desc(a).clss == SUBNORMAL) ^ (_desc(b).clss == SUBNORMAL))
689
  {
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
    sc_val_from_ulong(1, temp);
    sc_add(_exp(result), temp, _exp(result));
  }

  /* mant(res) = mant(a) / 1/2mant(b) */
  /* to gain more bits of precision in the result the dividend could be
   * shifted left, as this operation does not loose bits. This would not
   * fit into the integer precision, but due to the rounding bits (which
   * are always zero because the values are all normalized) the divisor
   * can be shifted right instead to achieve the same result */
  sc_val_from_ulong(2 + _desc(result).mantissa_size, temp);

  _shift_left(_mant(a), temp, dividend);

  {
705
    char *divisor = alloca(calc_buffer_size);
706
707
708
    sc_val_from_ulong(1, divisor);
    _shift_right(_mant(b), divisor, divisor);
    sc_div(dividend, divisor, _mant(result));
709
  }
710
711

  return _normalize(result, result, sc_had_carry());
712
713
}

Michael Beck's avatar
Michael Beck committed
714
static void _power_of_ten(int exp, descriptor_t *desc, char *result)
715
{
716
717
718
719
720
721
722
723
724
725
  char *build;
  char *temp;

  /* positive sign */
  _sign(result) = 0;

  /* set new descriptor (else result is supposed to already have one) */
  if (desc != NULL)
    memcpy(&_desc(result), desc, sizeof(descriptor_t));

726
727
  build = alloca(value_size);
  temp = alloca(value_size);
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742

  sc_val_from_ulong((1 << _desc(result).exponent_size)/2-1, _exp(result));

  if (exp > 0)
  {
    /* temp is value of ten now */
    sc_val_from_ulong(10, NULL);
    _save_result(temp);

    for (exp--; exp > 0; exp--) {
      _save_result(build);
      sc_mul(build, temp, NULL);
    }
    _save_result(build);

Michael Beck's avatar
Michael Beck committed
743
    /* temp is amount of left shift needed to put the value left of the radix point */
744
    sc_val_from_ulong(_desc(result).mantissa_size + 2, temp);
745

746
747
748
749
    _shift_left(build, temp, _mant(result));

    _normalize(result, result, 0);
  }
750
751
}

Michael Beck's avatar
Michael Beck committed
752
753
754
755
756
/**
 * Truncate the fractional part away.
 *
 * This does not clip to any integer rang.
 */
757
static char* _trunc(const char *a, char *result)
758
{
Michael Beck's avatar
Michael Beck committed
759
760
  /*
   * When exponent == 0 all bits left of the radix point
761
   * are the integral part of the value. For 15bit exp_size
Michael Beck's avatar
Michael Beck committed
762
   * this would require a left shift of max. 16383 bits which
763
764
765
766
767
   * is too much.
   * But it is enough to ensure that no bit right of the radix
   * point remains set. This restricts the interesting
   * exponents to the interval [0, mant_size-1].
   * Outside this interval the truncated value is either 0 or
Michael Beck's avatar
Michael Beck committed
768
769
   * it does not have fractional parts.
   */
770
771
772
773

  int exp_bias, exp_val;
  char *temp;

774
  temp = alloca(value_size);
775
776
777
778
779

  if (a != result)
    memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));

  exp_bias = (1<<_desc(a).exponent_size)/2-1;
Michael Beck's avatar
Michael Beck committed
780
  exp_val  = sc_val_to_long(_exp(a)) - exp_bias;
781
782
783
784
785

  if (exp_val < 0) {
    sc_val_from_ulong(0, NULL);
    _save_result(_exp(result));
    _save_result(_mant(result));
Michael Beck's avatar
Michael Beck committed
786
    _desc(result).clss = ZERO;
787
788
789
790
791
792

    return result;
  }

  if (exp_val > _desc(a).mantissa_size) {
    if (a != result)
793
      memcpy(result, a, calc_buffer_size);
794
795
796
797
798
799
800
801
802
803
804
805
806

    return result;
  }

  /* set up a proper mask to delete all bits right of the
   * radix point if the mantissa had been shifted until exp == 0 */
  sc_max_from_bits(1 + exp_val, 0, temp);
  sc_val_from_long(_desc(a).mantissa_size - exp_val + 2, NULL);
  _shift_left(temp, sc_get_buffer(), temp);

  /* and the mask and return the result */
  sc_and(_mant(a), temp, _mant(result));

807
  if (a != result) memcpy(_exp(result), _exp(a), value_size);
808
809

  return result;
810
811
}

812
813
814
/*
 * This does value sanity checking(or should do it), sets up any prerequisites,
 * calls the proper internal functions, clears up and returns
Michael Beck's avatar
Michael Beck committed
815
816
 * the result
 */
817
char* _calc(const char *a, const char *b, int opcode, char *result)
818
{
819
820
821
822
823
824
825
826
827
  char *temp;
#ifdef FLTCALC_TRACE_CALC
  char *buffer;

  buffer = alloca(100);
#endif

  if (result == NULL) result = calc_buffer;

828
  TRACEPRINTF(("%s ", fc_print(a, buffer, 100, FC_PACKED)));
829
830
  switch (opcode)
  {
831
832
    case FC_add:
      /* make the value with the bigger exponent the first one */
833
      TRACEPRINTF(("+ %s ", fc_print(b, buffer, 100, FC_PACKED)));
834
      if (sc_comp(_exp(a), _exp(b)) == -1)
Michael Beck's avatar
Michael Beck committed
835
        _fadd(b, a, result);
836
      else
Michael Beck's avatar
Michael Beck committed
837
        _fadd(a, b, result);
838
839
      break;
    case FC_sub:
840
      TRACEPRINTF(("- %s ", fc_print(b, buffer, 100, FC_PACKED)));
841
842
      temp = alloca(calc_buffer_size);
      memcpy(temp, b, calc_buffer_size);
843
844
      _sign(temp) = !_sign(b);
      if (sc_comp(_exp(a), _exp(temp)) == -1)
Michael Beck's avatar
Michael Beck committed
845
        _fadd(temp, a, result);
846
      else
Michael Beck's avatar
Michael Beck committed
847
        _fadd(a, temp, result);
848
      break;
849
    case FC_mul:
850
      TRACEPRINTF(("* %s ", fc_print(b, buffer, 100, FC_PACKED)));
Michael Beck's avatar
Michael Beck committed
851
      _fmul(a, b, result);
852
      break;
853
    case FC_div:
854
      TRACEPRINTF(("/ %s ", fc_print(b, buffer, 100, FC_PACKED)));
Michael Beck's avatar
Michael Beck committed
855
      _fdiv(a, b, result);
856
      break;
857
858
    case FC_neg:
      TRACEPRINTF(("negated "));
859
      if (a != result) memcpy(result, a, calc_buffer_size);
860
      _sign(result) = !_sign(a);
861
      break;
862
    case FC_int:
Michael Beck's avatar
Michael Beck committed
863
      TRACEPRINTF(("truncated to integer "));
864
865
866
      _trunc(a, result);
      break;
    case FC_rnd:
Michael Beck's avatar
Michael Beck committed
867
868
      TRACEPRINTF(("rounded to integer "));
      assert(0);
869
      break;
870
  }
871

872
  TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_PACKED)));
873
  return result;
874
875
}

876
877
878
879
/********
 * functions defined in fltcalc.h
 ********/
const void *fc_get_buffer(void)
880
{
881
882
883
  return calc_buffer;
}

Michael Beck's avatar
Michael Beck committed
884
int fc_get_buffer_length(void)
885
{
886
  return calc_buffer_size;
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
}

char* fc_val_from_str(const char *str, unsigned int len, char exp_size, char mant_size, char *result)
{
#if 0
  enum {
    START,
    LEFT_OF_DOT,
    RIGHT_OF_DOT,
    EXP_START,
    EXPONENT,
    END
  };

  char exp_sign;
  int exp_int, hsb, state;

  const char *old_str;

  int pos;
  char *mant_str, *exp_val, *power_val;

  if (result == NULL) result = calc_buffer;

911
912
  exp_val = alloca(value_size);
  power_val = alloca(calc_buffer_size);
913
914
915
916
  mant_str = alloca((len)?(len):(strlen(str)));

  _desc(result).exponent_size = exp_size;
  _desc(result).mantissa_size = mant_size;
Michael Beck's avatar
Michael Beck committed
917
  _desc(result).clss = NORMAL;
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045

  old_str = str;
  pos = 0;
  exp_int = 0;
  state = START;

  while (len == 0 || str-old_str < len)
  {
    switch (state) {
      case START:
        switch (*str) {
          case '+':
            _sign(result) = 0;
            state = LEFT_OF_DOT;
            str++;
            break;

          case '-':
            _sign(result) = 1;
            state = LEFT_OF_DOT;
            str++;
            break;

          case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
            _sign(result) = 0;
            state = LEFT_OF_DOT;
            break;

          case '.':
            _sign(result) = 0;
            state = RIGHT_OF_DOT;
            str++;
            break;

          case 'n':
          case 'N':
          case 'i':
          case 'I':
            break;

          default:
            _fail_char(old_str, len, str - old_str);
        }
        break;

      case LEFT_OF_DOT:
        switch (*str) {
          case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
            mant_str[pos++] = *(str++);
            break;

          case '.':
            state = RIGHT_OF_DOT;
            str++;
            break;

          case 'e':
          case 'E':
            state = EXP_START;
            str++;
            break;

          case '\0':
            mant_str[pos] = '\0';
            goto done;

          default:
            _fail_char(old_str, len, str - old_str);
        }
        break;

      case RIGHT_OF_DOT:
        switch (*str) {
          case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
            mant_str[pos++] = *(str++);
            exp_int++;
            break;

          case 'e':
          case 'E':
            state = EXP_START;
            str++;
            break;

          case '\0':
            mant_str[pos] = '\0';
            goto done;

          default:
            _fail_char(old_str, len, str - old_str);
        }
        break;

      case EXP_START:
        switch (*str) {
          case '-':
            exp_sign = 1;
            /* fall through */
          case '+':
            if (*(str-1) != 'e' && *(str-1) != 'E') _fail_char(old_str, len, str - old_str);
            str++;
            break;

          case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
            mant_str[pos] = '\0';
            pos = 1;
            str++;
            state = EXPONENT;
            break;

          default:
            _fail_char(old_str, len, str - old_str);
        }
        break;

      case EXPONENT:
        switch (*str) {
          case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
            pos++;
            str++;
            break;

          case '\0': goto done;

          default:
            _fail_char(old_str, len, str - old_str);
        }
    }
1046
  } /*  switch(state) */
1047
1048
1049

done:
  sc_val_from_str(mant_str, strlen(mant_str), _mant(result));
1050

1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
  /* shift to put value left of radix point */
  sc_val_from_ulong(mant_size + 2, exp_val);

  _shift_left(_mant(result), exp_val, _mant(result));

  sc_val_from_ulong((1 << exp_size)/2-1, _exp(result));

  _normalize(result, result, 0);

  if (state == EXPONENT) {
    exp_int -= atoi(str-pos);
  }

  _power_of_ten(exp_int, &_desc(result), power_val);

Michael Beck's avatar
Michael Beck committed
1066
  _fdiv(result, power_val, result);
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086

  return result;
#else

  /* XXX excuse of an implementation to make things work */
  LLDBL val;
#ifdef HAVE_LONG_DOUBLE
  val = strtold(str, NULL);
#else
  val = strtod(str, NULL);
#endif

  DEBUGPRINTF(("val_from_str(%s)\n", str));
  return fc_val_from_float(val, exp_size, mant_size, result);
#endif
}

char* fc_val_from_float(LLDBL l, char exp_size, char mant_size, char* result)
{
  char *temp;
1087
  int bias_res, bias_val, mant_val;
1088
  value_t srcval;
Michael Beck's avatar
Michael Beck committed
1089
  UINT32 sign, exponent, mantissa0, mantissa1;
1090
1091

  srcval.d = l;
1092
  bias_res = ((1<<exp_size)/2-1);
1093
1094

#ifdef HAVE_LONG_DOUBLE
Michael Beck's avatar
Michael Beck committed
1095
1096
1097
1098
1099
1100
  mant_val  = 64;
  bias_val  = 0x3fff;
  sign      = (srcval.val.high & 0x00008000) != 0;
  exponent  = (srcval.val.high & 0x00007FFF) ;
  mantissa0 = srcval.val.mid;
  mantissa1 = srcval.val.low;
1101
#else /* no long double */
Michael Beck's avatar
Michael Beck committed
1102
1103
1104
1105
1106
1107
  mant_val  = 52;
  bias_val  = 0x3ff;
  sign      = (srcval.val.high & 0x80000000) != 0;
  exponent  = (srcval.val.high & 0x7FF00000) >> 20;
  mantissa0 = srcval.val.high & 0x000FFFFF;
  mantissa1 = srcval.val.low;
1108
1109
1110
#endif

#ifdef HAVE_LONG_DOUBLE
1111
  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)); */
1112
1113
  DEBUGPRINTF(("(%d-%.4X-%.8X%.8X)\n", sign, exponent, mantissa0, mantissa1));
#else
1114
  TRACEPRINTF(("val_from_float(%.8X%.8X)\n", srcval.val.high, srcval.val.low));
1115
1116
1117
1118
  DEBUGPRINTF(("(%d-%.3X-%.5X%.8X)\n", sign, exponent, mantissa0, mantissa1));
#endif

  if (result == NULL) result = calc_buffer;
1119
  temp = alloca(value_size);
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129

  _desc(result).exponent_size = exp_size;
  _desc(result).mantissa_size = mant_size;

  /* extract sign */
  _sign(result) = sign;

  /* sign and flag suffice to identify nan or inf, no exponent/mantissa
   * encoding is needed. the function can return immediately in these cases */
  if (isnan(l)) {
Michael Beck's avatar
Michael Beck committed
1130
    _desc(result).clss = NAN;
1131
1132
    TRACEPRINTF(("val_from_float resulted in NAN\n"));
    return result;
1133
  }
1134
  else if (isinf(l)) {
Michael Beck's avatar
Michael Beck committed
1135
    _desc(result).clss = INF;
1136
1137
    TRACEPRINTF(("val_from_float resulted in %sINF\n", (_sign(result)==1)?"-":""));
    return result;
1138
  }
1139

1140
1141
1142
1143
1144
1145
1146
  /* build exponent, because input and output exponent and mantissa sizes may differ
   * this looks more complicated than it is: unbiased input exponent + output bias,
   * minus the mantissa difference which is added again later when the output float
   * becomes normalized */
#ifdef HAVE_EXPLICIT_ONE
  sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size-1), _exp(result));
#else
1147
  sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size), _exp(result));
1148
#endif
1149

1150
1151
  /* build mantissa representation */
#ifndef HAVE_EXPLICIT_ONE
1152
1153
1154
1155
1156
1157
  if (exponent != 0)
  {
    /* insert the hidden bit */
    sc_val_from_ulong(1, temp);
    sc_val_from_ulong(mant_val + 2, NULL);
    _shift_left(temp, sc_get_buffer(), NULL);
1158
  }
1159
  else
1160
#endif
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
  {
    sc_val_from_ulong(0, NULL);
  }

  _save_result(_mant(result));

  /* bits from the upper word */
  sc_val_from_ulong(mantissa0, temp);
  sc_val_from_ulong(34, NULL);
  _shift_left(temp, sc_get_buffer(), temp);
  sc_or(_mant(result), temp, _mant(result));

  /* bits from the lower word */
  sc_val_from_ulong(mantissa1, temp);
  sc_val_from_ulong(2, NULL);
  _shift_left(temp, sc_get_buffer(), temp);
  sc_or(_mant(result), temp, _mant(result));

  /* _normalize expects the radix point to be normal, so shift mantissa of subnormal
   * origin one to the left */
  if (exponent == 0)
  {
    sc_val_from_ulong(1, NULL);
    _shift_left(_mant(result), sc_get_buffer(), _mant(result));
  }

  _normalize(result, result, 0);

1189
  TRACEPRINTF(("val_from_float results in %s\n", fc_print(result, temp, calc_buffer_size, FC_PACKED)));
1190
1191

  return result;
1192
1193
}

1194
LLDBL fc_val_to_float(const void *val)
1195
{
1196
1197
1198
1199
  const char *value;
  char *temp = NULL;

  int byte_offset;
1200

1201
1202
1203
1204
1205
1206
1207
1208
1209
  UINT32 sign;
  UINT32 exponent;
  UINT32 mantissa0;
  UINT32 mantissa1;

  value_t buildval;

#ifdef HAVE_LONG_DOUBLE
  char result_exponent = 15;
1210
  char result_mantissa = 64;
Michael Beck's avatar
Tarval:    
Michael Beck committed
1211
#else
1212
  char result_exponent = 11;
1213
  char result_mantissa = 52;
Michael Beck's avatar
Tarval:    
Michael Beck committed
1214
#endif
1215

1216
  temp = alloca(calc_buffer_size);
1217
1218
1219
#ifdef HAVE_EXPLICIT_ONE
  value = fc_cast(val, result_exponent, result_mantissa-1, temp);
#else
1220
  value = fc_cast(val, result_exponent, result_mantissa, temp);
1221
#endif
1222
1223

  sign = _sign(value);
1224
1225

  /* @@@ long double exponent is 15bit, so the use of sc_val_to_long should not
1226
   * lead to wrong results */
Michael Beck's avatar
Michael Beck committed
1227
  exponent = sc_val_to_long(_exp(value)) ;
1228

1229
1230
  sc_val_from_ulong(2, NULL);
  _shift_right(_mant(value), sc_get_buffer(), _mant(value));
1231
1232
1233
1234
1235

  mantissa0 = 0;
  mantissa1 = 0;

  for (byte_offset = 0; byte_offset < 4; byte_offset++)
1236
    mantissa1 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << (byte_offset<<3);
1237
1238

  for (; (byte_offset<<3) < result_mantissa; byte_offset++)
1239
    mantissa0 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << ((byte_offset-4)<<3);
1240
1241

#ifdef HAVE_LONG_DOUBLE
1242
1243
1244
1245
  buildval.val.high = sign << 15;
  buildval.val.high |= exponent;
  buildval.val.mid = mantissa0;
  buildval.val.low = mantissa1;
1246
#else /* no long double */
Michael Beck's avatar
Michael Beck committed
1247
  mantissa0 &= 0x000FFFFF;  /* get rid of garbage */
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
  buildval.val.high = sign << 31;
  buildval.val.high |= exponent << 20;
  buildval.val.high |= mantissa0;
  buildval.val.low = mantissa1;
#endif

  TRACEPRINTF(("val_to_float: %d-%x-%x%x\n", sign, exponent, mantissa0, mantissa1));
  return buildval.d;
}

char* fc_cast(const void *val, char exp_size, char mant_size, char *result)
{
  const char *value = (const char*) val;
  char *temp;
  int exp_offset, val_bias, res_bias;

  if (result == NULL) result = calc_buffer;
1265
  temp = alloca(value_size);
1266
1267
1268

  if (_desc(value).exponent_size == exp_size && _desc(value).mantissa_size == mant_size)
  {
1269
    if (value != result) memcpy(result, value, calc_buffer_size);
1270
1271
1272
1273
1274
1275
    return result;
  }

  /* set the descriptor of the new value */
  _desc(result).exponent_size = exp_size;
  _desc(result).mantissa_size = mant_size;
Michael Beck's avatar
Michael Beck committed
1276
  _desc(result).clss = _desc(value).clss;
1277
1278
1279
1280
1281
1282
1283
1284
1285

  _sign(result) = _sign(value);

  /* when the mantissa sizes differ normalizing has to shift to align it.
   * this would change the exponent, which is unwanted. So calculate this
   * offset and add it */
  val_bias = (1<<_desc(value).exponent_size)/2-1;
  res_bias = (1<<exp_size)/2-1;

1286
  exp_offset = (res_bias - val_bias) - (_desc(value).mantissa_size - mant_size);