fltcalc.c 42.7 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
Michael Beck's avatar
Michael Beck committed
37
38
#include <stdio.h>
#include <assert.h>
39

40
41
#include "xmalloc.h"

42
typedef uint32_t UINT32;
Michael Beck's avatar
Tarval:    
Michael Beck committed
43

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

Michael Beck's avatar
Michael Beck committed
84
85
86
/**
 * possible float states
 */
87
typedef enum {
Michael Beck's avatar
Michael Beck committed
88
89
90
91
92
  NORMAL,       /**< normal representation, implicit 1 */
  ZERO,         /**< +/-0 */
  SUBNORMAL,    /**< denormals, implicit 0 */
  INF,          /**< +/-oo */
  NAN,          /**< Not A Number */
93
94
} value_class_t;

Michael Beck's avatar
Michael Beck committed
95
/** A descriptor for an IEEE float value. */
96
typedef struct {
Michael Beck's avatar
Michael Beck committed
97
98
99
  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 */
100
101
} descriptor_t;

102
#define CLEAR_BUFFER(buffer) memset(buffer, 0, calc_buffer_size)
103
104
105
106
107

/* because variable sized structs are impossible, the internal
 * value is represented as a pseudo-struct char array, addressed
 * by macros
 * struct {
108
 *   char sign;             //  0 for positive, 1 for negative
109
110
 *   char exp[value_size];
 *   char mant[value_size];
111
112
113
114
115
116
 *   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
117
#define _desc(a) (*(descriptor_t *)&((char *)a)[DESCRIPTOR_POS])
Michael Beck's avatar
Tarval:    
Michael Beck committed
118

119
120
121
#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
122

Michael Beck's avatar
Michael Beck committed
123
124
125
126
#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
127

Michael Beck's avatar
Michael Beck committed
128
129
130
131
#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); \
}
132
133
134
135
136
137
138

#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
139
140
#endif

141
142
143
144
145
#if FLTCALC_TRACE_CALC
#  define TRACEPRINTF(x) printf x
#else
#  define TRACEPRINTF(x) ((void)0)
#endif
146

147
static char *calc_buffer = NULL;
148

149
static fc_rounding_mode_t rounding_mode;
150

151
152
static int calc_buffer_size;
static int value_size;
153
154
155
156
static int SIGN_POS;
static int EXPONENT_POS;
static int MANTISSA_POS;
static int DESCRIPTOR_POS;
157

158
static int max_precision;
159
/********
160
 * private functions
161
 ********/
162
163
#if 0
static void _fail_char(const char *str, unsigned int len, int pos)
164
{
165
166
167
168
169
170
171
172
  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);
173
}
174
#endif
175

Michael Beck's avatar
Michael Beck committed
176
/** pack machine-like */
177
static char* _pack(const char *int_float, char *packed)
178
{
179
180
181
  char *shift_val;
  char *temp;
  char *val_buffer;
182

183
184
  temp = alloca(value_size);
  shift_val = alloca(value_size);
185

Michael Beck's avatar
Michael Beck committed
186
  switch (_desc(int_float).clss) {
187
    case NAN:
188
      val_buffer = alloca(calc_buffer_size);
189
190
191
192
193
      fc_get_qnan(_desc(int_float).exponent_size, _desc(int_float).mantissa_size, val_buffer);
      int_float = val_buffer;
      break;

    case INF:
194
      val_buffer = alloca(calc_buffer_size);
195
196
197
198
199
200
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
      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;
229
230
}

231
char* _normalize(const char *in_val, char *out_val, int sticky)
232
{
233
234
235
236
  int hsb;
  char lsb, guard, round, round_dir = 0;
  char *temp;

237
  temp = alloca(value_size);
238
239
240
241
242
243
244
245
246
247

  /* +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
248
  _desc(out_val).clss = NORMAL;
249

250
  /* mantissa all zeros, so zero exponent (because of explicit one)*/
251
252
253
254
255
256
  if (hsb == 2 + _desc(in_val).mantissa_size)
  {
    sc_val_from_ulong(0, _exp(out_val));
    hsb = -1;
  }

257
  /* shift the first 1 into the left of the radix point (i.e. hsb == -1) */
258
259
260
261
262
263
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
  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
293
    _desc(out_val).clss = SUBNORMAL;
294
295
296
297
298
299
300
301
302
  }

  /* 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;

303
  switch (rounding_mode)
304
305
306
307
308
309
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
  {
    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
342
343
  if (sc_is_zero(_mant(out_val)) && (_desc(out_val).clss == SUBNORMAL))
    _desc(out_val).clss = ZERO;
344
345
346

  /* 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
347
  if ((_desc(out_val).clss != SUBNORMAL) && (hsb < -1))
348
349
350
351
352
353
  {
    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
354
  else if ((_desc(out_val).clss == SUBNORMAL) && (hsb == -1))
355
  {
Michael Beck's avatar
Michael Beck committed
356
    /* overflow caused the mantissa to be normal again,
357
358
359
360
     * 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
361
    _desc(out_val).clss = NORMAL;
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
  }
  /* 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 */
391
      switch (rounding_mode) {
392
393
        case FC_TONEAREST:
        case FC_TOPOSITIVE:
Michael Beck's avatar
Michael Beck committed
394
          _desc(out_val).clss = INF;
395
396
397
398
399
400
401
402
          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 */
403
      switch (rounding_mode) {
404
405
        case FC_TONEAREST:
        case FC_TONEGATIVE:
Michael Beck's avatar
Michael Beck committed
406
          _desc(out_val).clss = INF;
407
408
409
410
411
412
413
414
415
416
          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;
417
418
}

Michael Beck's avatar
Michael Beck committed
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
/**
 * 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)


/**
436
437
 * calculate a + b, where a is the value with the bigger exponent
 */
Michael Beck's avatar
Michael Beck committed
438
static char* _fadd(const char* a, const char* b, char* result)
439
{
440
441
442
443
444
445
  char *temp;
  char *exp_diff;

  char sign;
  char sticky;

Michael Beck's avatar
Michael Beck committed
446
  handle_NAN(a, b, result);
447
448
449
450
451
452
453
454

  /* 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
455
  /* produce NaN on inf - inf */
Michael Beck's avatar
Michael Beck committed
456
  if (sign && (_desc(a).clss == INF) && (_desc(b).clss == INF))
457
458
    return fc_get_qnan(_desc(a).exponent_size, _desc(b).mantissa_size, result);

Michael Beck's avatar
Michael Beck committed
459
  temp     = alloca(value_size);
460
  exp_diff = alloca(value_size);
461
462
463
464
465
466
467
468

  /* 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
469
  if (sign && sc_val_to_long(exp_diff) == 0) {
470
471
    switch (sc_comp(_mant(a), _mant(b))) {
      case 1:  /* a > b */
Michael Beck's avatar
Michael Beck committed
472
        _sign(result) = _sign(a);  /* abs(a) is bigger and a is negative */
473
474
        break;
      case 0:  /* a == b */
Michael Beck's avatar
Michael Beck committed
475
        _sign(result) = (rounding_mode == FC_TONEGATIVE);
476
477
        break;
      case -1: /* a < b */
Michael Beck's avatar
Michael Beck committed
478
        _sign(result) = _sign(b); /* abs(b) is bigger and b is negative */
479
480
481
482
483
484
        break;
      default:
        /* can't be reached */
        break;
    }
  }
Michael Beck's avatar
Michael Beck committed
485
486
  else
    _sign(result) = _sign(a);
487
488

  /* sign has been taken care of, check for special cases */
Michael Beck's avatar
Michael Beck committed
489
  if (_desc(a).clss == ZERO) {
490
    if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
491
492
    return result;
  }
Michael Beck's avatar
Michael Beck committed
493
  if (_desc(b).clss == ZERO) {
494
    if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
495
496
497
    return result;
  }

Michael Beck's avatar
Michael Beck committed
498
  if (_desc(a).clss == INF) {
499
    if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
500
501
    return result;
  }
Michael Beck's avatar
Michael Beck committed
502
  if (_desc(b).clss == INF) {
503
    if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
504
505
506
507
508
509
    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
510
  if ((_desc(b).clss == SUBNORMAL) && (_desc(a).clss != SUBNORMAL))
511
512
513
514
515
516
517
518
519
520
521
522
523
  {
    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' */
524
    char *temp1 = alloca(calc_buffer_size);
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
    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
540
  if ((_desc(a).clss == SUBNORMAL) && (_desc(b).clss == SUBNORMAL))
541
542
543
544
545
546
  {
    sc_val_from_ulong(1, NULL);
    _shift_left(_mant(result), sc_get_buffer(), _mant(result));
  }

  /* resulting exponent is the bigger one */
547
  memmove(_exp(result), _exp(a), value_size);
548
549

  return _normalize(result, result, sticky);
550
551
}

Michael Beck's avatar
Michael Beck committed
552
553
554
555
/**
 * calculate a * b
 */
static char* _fmul(const char* a, const char* b, char* result)
556
{
557
  char *temp;
558

Michael Beck's avatar
Michael Beck committed
559
  handle_NAN(a, b, result);
560

561
  temp = alloca(value_size);
562
563
564
565
566
567

  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
568
  /* produce NaN on 0 * inf */
Michael Beck's avatar
Michael Beck committed
569
570
  if (_desc(a).clss == ZERO) {
    if (_desc(b).clss == INF)
571
572
      fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
    else
573
      if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
574
575
    return result;
  }
Michael Beck's avatar
Michael Beck committed
576
577
  if (_desc(b).clss == ZERO) {
    if (_desc(a).clss == INF)
578
579
      fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
    else
580
      if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-1);
581
582
583
    return result;
  }

Michael Beck's avatar
Michael Beck committed
584
  if (_desc(a).clss == INF) {
585
    if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
586
587
    return result;
  }
Michael Beck's avatar
Michael Beck committed
588
  if (_desc(b).clss == INF) {
589
    if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-1);
590
591
592
593
594
595
596
597
598
599
    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
600
  if ((_desc(a).clss == SUBNORMAL) ^ (_desc(b).clss == SUBNORMAL))
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
  {
    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());
618
619
}

Michael Beck's avatar
Michael Beck committed
620
621
622
623
/**
 * calculate a / b
 */
static char* _fdiv(const char* a, const char* b, char* result)
624
{
625
626
  char *temp, *dividend;

Michael Beck's avatar
Michael Beck committed
627
  handle_NAN(a, b, result);
628

629
630
  temp = alloca(value_size);
  dividend = alloca(value_size);
631
632
633
634
635
636
637

  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
638
639
  if (_desc(a).clss == ZERO) {
    if (_desc(b).clss == ZERO)
640
641
642
643
      /* 0/0 -> nan */
      fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
    else
      /* 0/x -> a */
644
      if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
645
646
647
    return result;
  }

Michael Beck's avatar
Michael Beck committed
648
649
  if (_desc(b).clss == INF) {
    if (_desc(a).clss == INF)
650
651
652
653
654
655
656
      /* 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
657
      _desc(result).clss = ZERO;
658
659
660
661
    }
    return result;
  }

Michael Beck's avatar
Michael Beck committed
662
  if (_desc(a).clss == INF) {
663
    /* inf/x -> inf */
664
    if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
665
666
    return result;
  }
Michael Beck's avatar
Michael Beck committed
667
  if (_desc(b).clss == ZERO) {
668
669
670
671
672
673
674
675
676
677
678
679
680
681
    /* 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
682
  if ((_desc(a).clss == SUBNORMAL) ^ (_desc(b).clss == SUBNORMAL))
683
  {
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
    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);

  {
699
    char *divisor = alloca(calc_buffer_size);
700
701
702
    sc_val_from_ulong(1, divisor);
    _shift_right(_mant(b), divisor, divisor);
    sc_div(dividend, divisor, _mant(result));
703
  }
704
705

  return _normalize(result, result, sc_had_carry());
706
707
}

Matthias Braun's avatar
Matthias Braun committed
708
#if 0
Michael Beck's avatar
Michael Beck committed
709
static void _power_of_ten(int exp, descriptor_t *desc, char *result)
710
{
711
712
713
714
715
716
717
718
719
720
  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));

721
722
  build = alloca(value_size);
  temp = alloca(value_size);
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737

  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
738
    /* temp is amount of left shift needed to put the value left of the radix point */
739
    sc_val_from_ulong(_desc(result).mantissa_size + 2, temp);
740

741
742
743
744
    _shift_left(build, temp, _mant(result));

    _normalize(result, result, 0);
  }
745
}
Matthias Braun's avatar
Matthias Braun committed
746
#endif
747

Michael Beck's avatar
Michael Beck committed
748
749
750
751
752
/**
 * Truncate the fractional part away.
 *
 * This does not clip to any integer rang.
 */
753
static char* _trunc(const char *a, char *result)
754
{
Michael Beck's avatar
Michael Beck committed
755
756
  /*
   * When exponent == 0 all bits left of the radix point
757
   * are the integral part of the value. For 15bit exp_size
Michael Beck's avatar
Michael Beck committed
758
   * this would require a left shift of max. 16383 bits which
759
760
761
762
763
   * 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
764
765
   * it does not have fractional parts.
   */
766
767
768
769

  int exp_bias, exp_val;
  char *temp;

770
  temp = alloca(value_size);
771
772
773
774
775

  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
776
  exp_val  = sc_val_to_long(_exp(a)) - exp_bias;
777
778
779
780
781

  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
782
    _desc(result).clss = ZERO;
783
784
785
786
787
788

    return result;
  }

  if (exp_val > _desc(a).mantissa_size) {
    if (a != result)
789
      memcpy(result, a, calc_buffer_size);
790
791
792
793
794
795
796
797
798
799
800
801
802

    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));

803
  if (a != result) memcpy(_exp(result), _exp(a), value_size);
804
805

  return result;
806
807
}

808
809
810
/*
 * 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
811
812
 * the result
 */
813
char* _calc(const char *a, const char *b, int opcode, char *result)
814
{
815
816
817
818
819
820
821
822
823
  char *temp;
#ifdef FLTCALC_TRACE_CALC
  char *buffer;

  buffer = alloca(100);
#endif

  if (result == NULL) result = calc_buffer;

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

868
  TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_PACKED)));
869
  return result;
870
871
}

872
873
874
875
/********
 * functions defined in fltcalc.h
 ********/
const void *fc_get_buffer(void)
876
{
877
878
879
  return calc_buffer;
}

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

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;

907
908
  exp_val = alloca(value_size);
  power_val = alloca(calc_buffer_size);
909
910
911
912
  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
913
  _desc(result).clss = NORMAL;
914
915
916
917
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

  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);
        }
    }
1042
  } /*  switch(state) */
1043
1044
1045

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

1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
  /* 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
1062
  _fdiv(result, power_val, result);
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082

  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;
1083
  int bias_res, bias_val, mant_val;
1084
  value_t srcval;
Michael Beck's avatar
Michael Beck committed
1085
  UINT32 sign, exponent, mantissa0, mantissa1;
1086
1087

  srcval.d = l;
1088
  bias_res = ((1<<exp_size)/2-1);
1089
1090

#ifdef HAVE_LONG_DOUBLE
Michael Beck's avatar
Michael Beck committed
1091
1092
1093
1094
1095
1096
  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;
1097
#else /* no long double */
Michael Beck's avatar
Michael Beck committed
1098
1099
1100
1101
1102
1103
  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;
1104
1105
1106
#endif

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

  if (result == NULL) result = calc_buffer;
1115
  temp = alloca(value_size);
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125

  _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
1126
    _desc(result).clss = NAN;
1127
1128
    TRACEPRINTF(("val_from_float resulted in NAN\n"));
    return result;
1129
  }
1130
  else if (isinf(l)) {
Michael Beck's avatar
Michael Beck committed
1131
    _desc(result).clss = INF;
1132
1133
    TRACEPRINTF(("val_from_float resulted in %sINF\n", (_sign(result)==1)?"-":""));
    return result;
1134
  }
1135

1136
1137
1138
1139
1140
1141
1142
  /* 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
1143
  sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size), _exp(result));
1144
#endif
1145

1146
1147
  /* build mantissa representation */
#ifndef HAVE_EXPLICIT_ONE
1148
1149
1150
1151
1152
1153
  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);
1154
  }
1155
  else
1156
#endif
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
  {
    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);

1185
  TRACEPRINTF(("val_from_float results in %s\n", fc_print(result, temp, calc_buffer_size, FC_PACKED)));
1186
1187

  return result;
1188
1189
}

1190
LLDBL fc_val_to_float(const void *val)
1191
{
1192
1193
1194
1195
  const char *value;
  char *temp = NULL;

  int byte_offset;
1196

1197
1198
1199
1200
1201
1202
1203
1204
1205
  UINT32 sign;
  UINT32 exponent;
  UINT32 mantissa0;
  UINT32 mantissa1;

  value_t buildval;

#ifdef HAVE_LONG_DOUBLE
  char result_exponent = 15;
1206
  char result_mantissa = 64;
Michael Beck's avatar
Tarval:    
Michael Beck committed
1207
#else
1208
  char result_exponent = 11;
1209
  char result_mantissa = 52;
Michael Beck's avatar
Tarval:    
Michael Beck committed
1210
#endif
1211

1212
  temp = alloca(calc_buffer_size);
1213
1214
1215
#ifdef HAVE_EXPLICIT_ONE
  value = fc_cast(val, result_exponent, result_mantissa-1, temp);
#else
1216
  value = fc_cast(val, result_exponent, result_mantissa, temp);
1217
#endif
1218
1219

  sign = _sign(value);
1220
1221

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

1225
1226
  sc_val_from_ulong(2, NULL);
  _shift_right(_mant(value), sc_get_buffer(), _mant(value));
1227
1228
1229
1230
1231

  mantissa0 = 0;
  mantissa1 = 0;

  for (byte_offset = 0; byte_offset < 4; byte_offset++)
1232
    mantissa1 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << (byte_offset<<3);
1233
1234

  for (; (byte_offset<<3) < result_mantissa; byte_offset++)
1235
    mantissa0 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << ((byte_offset-4)<<3);
1236
1237

#ifdef HAVE_LONG_DOUBLE
1238
1239
1240
1241
  buildval.val.high = sign << 15;
  buildval.val.high |= exponent;
  buildval.val.mid = mantissa0;
  buildval.val.low = mantissa1;
1242
#else /* no long double */
Michael Beck's avatar
Michael Beck committed
1243
  mantissa0 &= 0x000FFFFF;  /* get rid of garbage */
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
  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;
1261
  temp = alloca(value_size);
1262
1263
1264

  if (_desc(value).exponent_size == exp_size && _desc(value).mantissa_size == mant_size)
  {
1265
    if (value != result) memcpy(result, value, calc_buffer_size);
1266
1267
1268
1269
1270
1271
    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
1272
  _desc(result).clss = _desc(value).clss;
1273
1274
1275
1276
1277
1278
1279
1280
1281

  _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;

1282
  exp_offset = (res_bias - val_bias) - (_desc(value).mantissa_size - mant_size);
1283
1284
1285
  sc_val_from_long(exp_offset, temp);
  sc_add(_exp(value), temp, _exp(result));

1286
  /* _normalize expects normalized radix point */