fltcalc.c 36.6 KB
Newer Older
Christian Würdig's avatar
Christian Würdig committed
1
2
/*
 * This file is part of libFirm.
3
 * Copyright (C) 2012 University of Karlsruhe.
Christian Würdig's avatar
Christian Würdig committed
4
5
 */

Michael Beck's avatar
Michael Beck committed
6
7
8
/**
 * @file
 * @brief    tarval floating point calculations
Matthias Braun's avatar
Matthias Braun committed
9
10
 * @date     2003-2013
 * @author   Mathias Heil, Michael Beck, Matthias Braun
11
12
 */
#include "fltcalc.h"
13
#include "strcalc.h"
14
#include "error.h"
15

16
17
#include <math.h>
#include <inttypes.h>
18
#include <float.h>
19
#include <string.h>
20
#include <stdlib.h>
Michael Beck's avatar
Michael Beck committed
21
22
#include <stdio.h>
#include <assert.h>
23
#include <stdbool.h>
24

25
26
#include "xmalloc.h"

27
static long double string_to_long_double(const char *str)
28
{
29
30
31
32
33
#if __STDC_VERSION__ >= 199901L || _POSIX_C_SOURCE >= 200112L
	return strtold(str, NULL);
#else
	return strtod(str, NULL);
#endif
34
}
35
36
37
38
39
40
41
42

static bool my_isnan(long double val)
{
#if __STDC_VERSION__ >= 199901L
	return isnan(val);
#else
	/* hopefully the compiler does not optimize aggressively (=incorrect) */
	return val != val;
43
#endif
44
45
46
47
48
49
50
51
52
53
54
}

static bool my_isinf(long double val)
{
#if __STDC_VERSION__ >= 199901L
	return isinf(val);
#else
	/* hopefully the compiler does not optimize aggressively (=incorrect) */
	return my_isnan(val-val) && !my_isnan(val);
#endif
}
55

Michael Beck's avatar
Michael Beck committed
56
/** The number of extra precision rounding bits */
57
58
#define ROUNDING_BITS 2

59
typedef union {
Matthias Braun's avatar
Matthias Braun committed
60
	volatile struct {
61
#ifdef WORDS_BIGENDIAN
62
		uint32_t high;
63
#else
64
		uint32_t low;
65
#endif
66
		uint32_t mid;
67
68
69
70
71
72
#ifdef WORDS_BIGENDIAN
		uint32_t low;
#else
		uint32_t high;
#endif
	} val_ld12;
Matthias Braun's avatar
Matthias Braun committed
73
	volatile struct {
74
75
76
77
#ifdef WORDS_BIGENDIAN
		uint32_t high;
#else
		uint32_t low;
78
#endif
79
#ifdef WORDS_BIGENDIAN
80
		uint32_t low;
81
#else
82
		uint32_t high;
83
#endif
84
	} val_ld8;
85
86
	volatile long double d;
} value_t;
87

88
#define CLEAR_BUFFER(buffer) memset(buffer, 0, calc_buffer_size)
89

Michael Beck's avatar
Michael Beck committed
90
/* our floating point value */
91
struct fp_value {
92
	float_descriptor_t desc;
Matthias Braun's avatar
Matthias Braun committed
93
94
	unsigned char      clss;
	char               sign;
Matthias Braun's avatar
Matthias Braun committed
95
96
97
98
	/** exp[value_size] + mant[value_size].
	 * Mantissa has an explizit one at the beginning (contrary to many
	 * floatingpoint formats) */
	char               value[];
Michael Beck's avatar
Michael Beck committed
99
100
101
102
};

#define _exp(a)  &((a)->value[0])
#define _mant(a) &((a)->value[value_size])
Michael Beck's avatar
Tarval:    
Michael Beck committed
103

104
#define _save_result(x) memcpy((x), sc_get_buffer(), value_size)
Michael Beck's avatar
Michael Beck committed
105
106
#define _shift_right(x, y, res) sc_shr((x), (y), value_size*4, 0, (res))
#define _shift_left(x, y, res) sc_shl((x), (y), value_size*4, 0, (res))
Michael Beck's avatar
Tarval:    
Michael Beck committed
107

108
/** A temporary buffer. */
Michael Beck's avatar
Michael Beck committed
109
static fp_value *calc_buffer = NULL;
110

111
/** Current rounding mode.*/
112
static fc_rounding_mode_t rounding_mode;
113

114
115
static int calc_buffer_size;
static int value_size;
116
static int max_precision;
117

118
/** Exact flag. */
Matthias Braun's avatar
Matthias Braun committed
119
static bool fc_exact = true;
120

121
122
static float_descriptor_t long_double_desc;

Michael Beck's avatar
Michael Beck committed
123
/** pack machine-like */
124
125
static void *pack(const fp_value *int_float, void *packed)
{
126
	switch ((value_class_t)int_float->clss) {
Matthias Braun's avatar
Matthias Braun committed
127
128
	case FC_NAN: {
		fp_value *val_buffer = (fp_value*) alloca(calc_buffer_size);
Michael Beck's avatar
Michael Beck committed
129
		fc_get_qnan(&int_float->desc, val_buffer);
Michael Beck's avatar
Michael Beck committed
130
131
		int_float = val_buffer;
		break;
Matthias Braun's avatar
Matthias Braun committed
132
	}
Michael Beck's avatar
Michael Beck committed
133

Matthias Braun's avatar
Matthias Braun committed
134
135
	case FC_INF: {
		fp_value *val_buffer = (fp_value*) alloca(calc_buffer_size);
Michael Beck's avatar
Michael Beck committed
136
		fc_get_plusinf(&int_float->desc, val_buffer);
Michael Beck's avatar
Michael Beck committed
137
138
139
		val_buffer->sign = int_float->sign;
		int_float = val_buffer;
		break;
Matthias Braun's avatar
Matthias Braun committed
140
	}
Michael Beck's avatar
Michael Beck committed
141
142
143
144

	default:
		break;
	}
145
146
	assert(int_float->desc.explicit_one <= 1);

Michael Beck's avatar
Michael Beck committed
147
	/* pack sign: move it to the left after exponent AND mantissa */
Matthias Braun's avatar
Matthias Braun committed
148
	char *temp = ALLOCAN(char, value_size);
Michael Beck's avatar
Michael Beck committed
149
150
	sc_val_from_ulong(int_float->sign, temp);

Matthias Braun's avatar
Matthias Braun committed
151
	int pos = int_float->desc.exponent_size + int_float->desc.mantissa_size + int_float->desc.explicit_one;
Michael Beck's avatar
Michael Beck committed
152
	sc_val_from_ulong(pos, NULL);
Michael Beck's avatar
Michael Beck committed
153
154
	_shift_left(temp, sc_get_buffer(), packed);

Michael Beck's avatar
Michael Beck committed
155
156
	/* pack exponent: move it to the left after mantissa */
	pos = int_float->desc.mantissa_size + int_float->desc.explicit_one;
Matthias Braun's avatar
Matthias Braun committed
157
	char *shift_val = ALLOCAN(char, value_size);
Michael Beck's avatar
Michael Beck committed
158
	sc_val_from_ulong(pos, shift_val);
Michael Beck's avatar
Michael Beck committed
159
160
	_shift_left(_exp(int_float), shift_val, temp);

Michael Beck's avatar
Michael Beck committed
161
	/* combine sign|exponent */
Michael Beck's avatar
Michael Beck committed
162
163
164
	sc_or(temp, packed, packed);

	/* extract mantissa */
165
166
	/* remove rounding bits */
	sc_val_from_ulong(ROUNDING_BITS, shift_val);
Michael Beck's avatar
Michael Beck committed
167
168
169
	_shift_right(_mant(int_float), shift_val, temp);

	/* remove leading 1 (or 0 if denormalized) */
Michael Beck's avatar
Michael Beck committed
170
	sc_max_from_bits(pos, 0, shift_val); /* all mantissa bits are 1's */
Michael Beck's avatar
Michael Beck committed
171
172
	sc_and(temp, shift_val, temp);

Michael Beck's avatar
Michael Beck committed
173
	/* combine sign|exponent|mantissa */
Michael Beck's avatar
Michael Beck committed
174
175
176
	sc_or(temp, packed, packed);

	return packed;
177
178
}

179
180
181
/**
 * Normalize a fp_value.
 *
Matthias Braun's avatar
Matthias Braun committed
182
 * @return true if result is exact
183
 */
Matthias Braun's avatar
Matthias Braun committed
184
static bool normalize(const fp_value *in_val, fp_value *out_val, bool sticky)
185
{
186
	/* save rounding bits at the end */
Matthias Braun's avatar
Matthias Braun committed
187
	int hsb = ROUNDING_BITS + in_val->desc.mantissa_size - sc_get_highest_set_bit(_mant(in_val)) - 1;
Michael Beck's avatar
Michael Beck committed
188

Matthias Braun's avatar
Matthias Braun committed
189
	if (in_val != out_val) {
Michael Beck's avatar
Michael Beck committed
190
		out_val->sign = in_val->sign;
191
		out_val->desc = in_val->desc;
Michael Beck's avatar
Michael Beck committed
192
193
	}

194
	out_val->clss = FC_NORMAL;
Michael Beck's avatar
Michael Beck committed
195

196
	/* mantissa all zeros, so zero exponent (because of explicit one) */
Matthias Braun's avatar
Matthias Braun committed
197
	if (hsb == ROUNDING_BITS + in_val->desc.mantissa_size) {
Michael Beck's avatar
Michael Beck committed
198
199
200
201
202
		sc_val_from_ulong(0, _exp(out_val));
		hsb = -1;
	}

	/* shift the first 1 into the left of the radix point (i.e. hsb == -1) */
Matthias Braun's avatar
Matthias Braun committed
203
204
205
	bool  exact = true;
	char *temp  = (char*) alloca(value_size);
	if (hsb < -1) {
Michael Beck's avatar
Michael Beck committed
206
207
208
209
210
211
		/* 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 */
212
		if (sc_had_carry()) {
Matthias Braun's avatar
Matthias Braun committed
213
214
			exact  = false;
			sticky = true;
215
		}
Michael Beck's avatar
Michael Beck committed
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
		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))) {
		/* 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));
234
		if (sc_had_carry()) {
Matthias Braun's avatar
Matthias Braun committed
235
236
			exact  = false;
			sticky = true;
237
		}
Michael Beck's avatar
Michael Beck committed
238
239
240
		/* denormalized means exponent of zero */
		sc_val_from_ulong(0, _exp(out_val));

241
		out_val->clss = FC_SUBNORMAL;
Michael Beck's avatar
Michael Beck committed
242
243
244
245
246
	}

	/* 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 */
Matthias Braun's avatar
Matthias Braun committed
247
248
249
250
	char lsb       = sc_sub_bits(_mant(out_val), out_val->desc.mantissa_size + ROUNDING_BITS, 0) & 0x7;
	char guard     = (lsb&0x2)>>1;
	char round     = lsb&0x1;
	char round_dir = 0;
Michael Beck's avatar
Michael Beck committed
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
	switch (rounding_mode) {
	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 = (!out_val->sign && (guard || round || sticky));
		break;
	case FC_TONEGATIVE:
		/* if negative: round to one if the exact value is bigger, else to zero */
		round_dir = (out_val->sign && (guard || round || sticky));
		break;
	case FC_TOZERO:
		/* always round to 0 (chopping mode) */
		round_dir = 0;
		break;
	}

	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));
Matthias Braun's avatar
Matthias Braun committed
282
		exact = false;
Michael Beck's avatar
Michael Beck committed
283
284
285
	}

	/* could have rounded down to zero */
286
287
	if (sc_is_zero(_mant(out_val)) && (out_val->clss == FC_SUBNORMAL))
		out_val->clss = FC_ZERO;
Michael Beck's avatar
Michael Beck committed
288
289

	/* check for rounding overflow */
290
	hsb = ROUNDING_BITS + out_val->desc.mantissa_size - sc_get_highest_set_bit(_mant(out_val)) - 1;
291
	if ((out_val->clss != FC_SUBNORMAL) && (hsb < -1)) {
Michael Beck's avatar
Michael Beck committed
292
293
		sc_val_from_ulong(1, temp);
		_shift_right(_mant(out_val), temp, _mant(out_val));
294
		if (exact && sc_had_carry())
Matthias Braun's avatar
Matthias Braun committed
295
			exact = false;
Michael Beck's avatar
Michael Beck committed
296
		sc_add(_exp(out_val), temp, _exp(out_val));
297
	} else if ((out_val->clss == FC_SUBNORMAL) && (hsb == -1)) {
Michael Beck's avatar
Michael Beck committed
298
299
300
301
302
		/* overflow caused the mantissa to be normal again,
		 * so adapt the exponent accordingly */
		sc_val_from_ulong(1, temp);
		sc_add(_exp(out_val), temp, _exp(out_val));

303
		out_val->clss = FC_NORMAL;
Michael Beck's avatar
Michael Beck committed
304
305
306
307
308
309
310
311
	}
	/* 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 << out_val->desc.exponent_size) - 1, temp);
312
	if (sc_comp(_exp(out_val), temp) != ir_relation_less) {
Michael Beck's avatar
Michael Beck committed
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
		/* 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 (out_val->sign == 0) {
			/* value is positive */
			switch (rounding_mode) {
			case FC_TONEAREST:
			case FC_TOPOSITIVE:
334
				out_val->clss = FC_INF;
Michael Beck's avatar
Michael Beck committed
335
336
337
338
				break;

			case FC_TONEGATIVE:
			case FC_TOZERO:
Michael Beck's avatar
Michael Beck committed
339
				fc_get_max(&out_val->desc, out_val);
Michael Beck's avatar
Michael Beck committed
340
341
342
343
344
345
			}
		} else {
			/* value is negative */
			switch (rounding_mode) {
			case FC_TONEAREST:
			case FC_TONEGATIVE:
346
				out_val->clss = FC_INF;
Michael Beck's avatar
Michael Beck committed
347
348
349
350
				break;

			case FC_TOPOSITIVE:
			case FC_TOZERO:
Michael Beck's avatar
Michael Beck committed
351
				fc_get_min(&out_val->desc, out_val);
Michael Beck's avatar
Michael Beck committed
352
353
354
			}
		}
	}
355
	return exact;
356
357
}

Michael Beck's avatar
Michael Beck committed
358
/**
359
360
 * Operations involving NaN's must return NaN.
 * They are NOT exact.
Michael Beck's avatar
Michael Beck committed
361
362
363
 */
#define handle_NAN(a, b, result) \
do {                                                      \
364
  if (a->clss == FC_NAN) {                                \
Michael Beck's avatar
Michael Beck committed
365
    if (a != result) memcpy(result, a, calc_buffer_size); \
Matthias Braun's avatar
Matthias Braun committed
366
    fc_exact = false;                                     \
Michael Beck's avatar
Michael Beck committed
367
    return;                                               \
Michael Beck's avatar
Michael Beck committed
368
  }                                                       \
369
  if (b->clss == FC_NAN) {                                \
Michael Beck's avatar
Michael Beck committed
370
    if (b != result) memcpy(result, b, calc_buffer_size); \
Matthias Braun's avatar
Matthias Braun committed
371
    fc_exact = false;                                     \
Michael Beck's avatar
Michael Beck committed
372
    return;                                               \
Michael Beck's avatar
Michael Beck committed
373
  }                                                       \
Matthias Braun's avatar
Matthias Braun committed
374
} while (0)
Michael Beck's avatar
Michael Beck committed
375
376
377


/**
378
379
 * calculate a + b, where a is the value with the bigger exponent
 */
380
381
static void _fadd(const fp_value *a, const fp_value *b, fp_value *result)
{
Matthias Braun's avatar
Matthias Braun committed
382
	fc_exact = true;
383

Michael Beck's avatar
Michael Beck committed
384
385
386
387
	handle_NAN(a, b, result);

	/* make sure result has a descriptor */
	if (result != a && result != b)
388
		result->desc = a->desc;
Michael Beck's avatar
Michael Beck committed
389
390

	/* determine if this is an addition or subtraction */
Matthias Braun's avatar
Matthias Braun committed
391
	char sign = a->sign ^ b->sign;
Michael Beck's avatar
Michael Beck committed
392
393

	/* produce NaN on inf - inf */
394
	if (sign && (a->clss == FC_INF) && (b->clss == FC_INF)) {
Matthias Braun's avatar
Matthias Braun committed
395
		fc_exact = false;
Michael Beck's avatar
Michael Beck committed
396
		fc_get_qnan(&a->desc, result);
Michael Beck's avatar
Michael Beck committed
397
398
399
		return;
	}

Matthias Braun's avatar
Matthias Braun committed
400
401
	char *temp     = ALLOCAN(char, value_size);
	char *exp_diff = ALLOCAN(char, value_size);
Michael Beck's avatar
Michael Beck committed
402
403
404
405
406
407
408
409

	/* 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 ?) */
Matthias Braun's avatar
Matthias Braun committed
410
	char res_sign;
Michael Beck's avatar
Michael Beck committed
411
412
	if (sign && sc_val_to_long(exp_diff) == 0) {
		switch (sc_comp(_mant(a), _mant(b))) {
413
		case ir_relation_greater:  /* a > b */
Michael Beck's avatar
Michael Beck committed
414
415
			res_sign = a->sign;  /* abs(a) is bigger and a is negative */
			break;
416
		case ir_relation_equal:  /* a == b */
Michael Beck's avatar
Michael Beck committed
417
418
			res_sign = (rounding_mode == FC_TONEGATIVE);
			break;
419
		case ir_relation_less: /* a < b */
Michael Beck's avatar
Michael Beck committed
420
421
422
			res_sign = b->sign; /* abs(b) is bigger and b is negative */
			break;
		default:
423
			panic("invalid comparison result");
Michael Beck's avatar
Michael Beck committed
424
		}
Matthias Braun's avatar
Matthias Braun committed
425
	} else {
Michael Beck's avatar
Michael Beck committed
426
		res_sign = a->sign;
Matthias Braun's avatar
Matthias Braun committed
427
	}
Michael Beck's avatar
Michael Beck committed
428
429
430
	result->sign = res_sign;

	/* sign has been taken care of, check for special cases */
431
	if (a->clss == FC_ZERO || b->clss == FC_INF) {
Michael Beck's avatar
Michael Beck committed
432
433
		if (b != result)
			memcpy(result, b, calc_buffer_size);
434
		fc_exact = b->clss == FC_NORMAL;
Michael Beck's avatar
Michael Beck committed
435
436
437
		result->sign = res_sign;
		return;
	}
438
	if (b->clss == FC_ZERO || a->clss == FC_INF) {
Michael Beck's avatar
Michael Beck committed
439
440
		if (a != result)
			memcpy(result, a, calc_buffer_size);
441
		fc_exact = a->clss == FC_NORMAL;
Michael Beck's avatar
Michael Beck committed
442
443
444
445
446
447
448
		result->sign = res_sign;
		return;
	}

	/* 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 */
449
	if ((b->clss == FC_SUBNORMAL) && (a->clss != FC_SUBNORMAL)) {
Michael Beck's avatar
Michael Beck committed
450
451
452
453
454
		sc_val_from_ulong(1, temp);
		sc_sub(exp_diff, temp, exp_diff);
	}

	_shift_right(_mant(b), exp_diff, temp);
Matthias Braun's avatar
Matthias Braun committed
455
	bool sticky = sc_had_carry();
456
	fc_exact &= !sticky;
Michael Beck's avatar
Michael Beck committed
457
458

	if (sticky && sign) {
Matthias Braun's avatar
Matthias Braun committed
459
460
461
462
463
		/* 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' */
		char *temp1 = ALLOCAN(char, value_size);
Michael Beck's avatar
Michael Beck committed
464
465
466
467
468
		sc_val_from_ulong(1, temp1);
		sc_add(temp, temp1, temp);
	}

	if (sign) {
469
		if (sc_comp(_mant(a), temp) == ir_relation_less)
Michael Beck's avatar
Michael Beck committed
470
471
472
473
474
475
476
477
478
			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 */
479
	if ((a->clss == FC_SUBNORMAL) && (b->clss == FC_SUBNORMAL)) {
Michael Beck's avatar
Michael Beck committed
480
481
482
483
484
485
486
		sc_val_from_ulong(1, NULL);
		_shift_left(_mant(result), sc_get_buffer(), _mant(result));
	}

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

487
	fc_exact &= normalize(result, result, sticky);
488
489
}

Michael Beck's avatar
Michael Beck committed
490
491
492
/**
 * calculate a * b
 */
493
494
static void _fmul(const fp_value *a, const fp_value *b, fp_value *result)
{
Matthias Braun's avatar
Matthias Braun committed
495
	fc_exact = true;
496

Michael Beck's avatar
Michael Beck committed
497
498
499
	handle_NAN(a, b, result);

	if (result != a && result != b)
500
		result->desc = a->desc;
Michael Beck's avatar
Michael Beck committed
501

Matthias Braun's avatar
Matthias Braun committed
502
	char res_sign;
Michael Beck's avatar
Michael Beck committed
503
504
505
	result->sign = res_sign = a->sign ^ b->sign;

	/* produce NaN on 0 * inf */
506
507
	if (a->clss == FC_ZERO) {
		if (b->clss == FC_INF) {
Michael Beck's avatar
Michael Beck committed
508
			fc_get_qnan(&a->desc, result);
Matthias Braun's avatar
Matthias Braun committed
509
			fc_exact = false;
510
		} else {
Michael Beck's avatar
Michael Beck committed
511
512
513
514
515
516
			if (a != result)
				memcpy(result, a, calc_buffer_size);
			result->sign = res_sign;
		}
		return;
	}
517
518
	if (b->clss == FC_ZERO) {
		if (a->clss == FC_INF) {
Michael Beck's avatar
Michael Beck committed
519
			fc_get_qnan(&a->desc, result);
Matthias Braun's avatar
Matthias Braun committed
520
			fc_exact = false;
521
		} else {
Michael Beck's avatar
Michael Beck committed
522
523
524
525
526
527
528
			if (b != result)
				memcpy(result, b, calc_buffer_size);
			result->sign = res_sign;
		}
		return;
	}

529
	if (a->clss == FC_INF) {
Matthias Braun's avatar
Matthias Braun committed
530
		fc_exact = false;
Michael Beck's avatar
Michael Beck committed
531
532
533
534
535
		if (a != result)
			memcpy(result, a, calc_buffer_size);
		result->sign = res_sign;
		return;
	}
536
	if (b->clss == FC_INF) {
Matthias Braun's avatar
Matthias Braun committed
537
		fc_exact = false;
Michael Beck's avatar
Michael Beck committed
538
539
540
541
542
543
544
545
546
		if (b != result)
			memcpy(result, b, calc_buffer_size);
		result->sign = res_sign;
		return;
	}

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

Matthias Braun's avatar
Matthias Braun committed
547
	char *temp = ALLOCAN(char, value_size);
Michael Beck's avatar
Michael Beck committed
548
549
550
551
	sc_val_from_ulong((1 << (a->desc.exponent_size - 1)) - 1, temp);
	sc_sub(_exp(result), temp, _exp(result));

	/* mixed normal, subnormal values introduce an error of 1, correct it */
552
	if ((a->clss == FC_SUBNORMAL) ^ (b->clss == FC_SUBNORMAL)) {
Michael Beck's avatar
Michael Beck committed
553
554
555
556
557
558
559
560
561
562
		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
563
564
	 * because of the rounding bits */
	sc_val_from_ulong(ROUNDING_BITS + result->desc.mantissa_size, temp);
Michael Beck's avatar
Michael Beck committed
565
566

	_shift_right(_mant(result), temp, _mant(result));
Matthias Braun's avatar
Matthias Braun committed
567
	bool sticky = sc_had_carry();
568
	fc_exact &= !sticky;
Michael Beck's avatar
Michael Beck committed
569

570
	fc_exact &= normalize(result, result, sticky);
571
572
}

Michael Beck's avatar
Michael Beck committed
573
574
575
/**
 * calculate a / b
 */
576
577
static void _fdiv(const fp_value *a, const fp_value *b, fp_value *result)
{
Matthias Braun's avatar
Matthias Braun committed
578
	fc_exact = true;
579

Michael Beck's avatar
Michael Beck committed
580
581
582
	handle_NAN(a, b, result);

	if (result != a && result != b)
583
		result->desc = a->desc;
Michael Beck's avatar
Michael Beck committed
584

Matthias Braun's avatar
Matthias Braun committed
585
	char res_sign;
Michael Beck's avatar
Michael Beck committed
586
587
	result->sign = res_sign = a->sign ^ b->sign;

588
	/* produce FC_NAN on 0/0 and inf/inf */
589
590
	if (a->clss == FC_ZERO) {
		if (b->clss == FC_ZERO) {
Michael Beck's avatar
Michael Beck committed
591
592
			/* 0/0 -> NaN */
			fc_get_qnan(&a->desc, result);
Matthias Braun's avatar
Matthias Braun committed
593
			fc_exact = false;
594
		} else {
Michael Beck's avatar
Michael Beck committed
595
596
597
598
599
600
601
602
			/* 0/x -> a */
			if (a != result)
				memcpy(result, a, calc_buffer_size);
			result->sign = res_sign;
		}
		return;
	}

603
	if (b->clss == FC_INF) {
Matthias Braun's avatar
Matthias Braun committed
604
		fc_exact = false;
605
		if (a->clss == FC_INF) {
Michael Beck's avatar
Michael Beck committed
606
607
			/* inf/inf -> NaN */
			fc_get_qnan(&a->desc, result);
608
		} else {
Michael Beck's avatar
Michael Beck committed
609
610
611
612
			/* x/inf -> 0 */
			sc_val_from_ulong(0, NULL);
			_save_result(_exp(result));
			_save_result(_mant(result));
613
			result->clss = FC_ZERO;
Michael Beck's avatar
Michael Beck committed
614
615
616
617
		}
		return;
	}

618
	if (a->clss == FC_INF) {
Matthias Braun's avatar
Matthias Braun committed
619
		fc_exact = false;
Michael Beck's avatar
Michael Beck committed
620
621
622
623
624
625
		/* inf/x -> inf */
		if (a != result)
			memcpy(result, a, calc_buffer_size);
		result->sign = res_sign;
		return;
	}
626
	if (b->clss == FC_ZERO) {
Matthias Braun's avatar
Matthias Braun committed
627
		fc_exact = false;
Michael Beck's avatar
Michael Beck committed
628
629
		/* division by zero */
		if (result->sign)
Michael Beck's avatar
Michael Beck committed
630
			fc_get_minusinf(&a->desc, result);
Michael Beck's avatar
Michael Beck committed
631
		else
Michael Beck's avatar
Michael Beck committed
632
			fc_get_plusinf(&a->desc, result);
Michael Beck's avatar
Michael Beck committed
633
634
635
636
		return;
	}

	/* exp = exp(a) - exp(b) + excess - 1*/
Matthias Braun's avatar
Matthias Braun committed
637
	char *temp = ALLOCAN(char, value_size);
Michael Beck's avatar
Michael Beck committed
638
639
640
641
642
	sc_sub(_exp(a), _exp(b), _exp(result));
	sc_val_from_ulong((1 << (a->desc.exponent_size - 1)) - 2, temp);
	sc_add(_exp(result), temp, _exp(result));

	/* mixed normal, subnormal values introduce an error of 1, correct it */
643
	if ((a->clss == FC_SUBNORMAL) ^ (b->clss == FC_SUBNORMAL)) {
Michael Beck's avatar
Michael Beck committed
644
645
646
647
648
649
650
651
652
653
		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 */
654
	sc_val_from_ulong(ROUNDING_BITS + result->desc.mantissa_size, temp);
Michael Beck's avatar
Michael Beck committed
655

Matthias Braun's avatar
Matthias Braun committed
656
	char *dividend = ALLOCAN(char, value_size);
Michael Beck's avatar
Michael Beck committed
657
658
	_shift_left(_mant(a), temp, dividend);

Matthias Braun's avatar
Matthias Braun committed
659
660
661
662
663
664
	char *divisor = ALLOCAN(char, value_size);
	sc_val_from_ulong(1, divisor);
	_shift_right(_mant(b), divisor, divisor);
	sc_div(dividend, divisor, _mant(result));
	bool sticky = sc_had_carry();
	fc_exact &= !sticky;
Michael Beck's avatar
Michael Beck committed
665

666
	fc_exact &= normalize(result, result, sticky);
667
668
}

Michael Beck's avatar
Michael Beck committed
669
670
671
/**
 * Truncate the fractional part away.
 *
672
 * This does not clip to any integer range.
Michael Beck's avatar
Michael Beck committed
673
 */
674
675
static void _trunc(const fp_value *a, fp_value *result)
{
Matthias Braun's avatar
Matthias Braun committed
676
	/* When exponent == 0 all bits left of the radix point
Michael Beck's avatar
Michael Beck committed
677
678
679
680
681
682
683
	 * are the integral part of the value. For 15bit exp_size
	 * this would require a left shift of max. 16383 bits which
	 * 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
Matthias Braun's avatar
Matthias Braun committed
684
	 * it does not have fractional parts. */
Michael Beck's avatar
Michael Beck committed
685

686
	/* fixme: can be exact */
Matthias Braun's avatar
Matthias Braun committed
687
	fc_exact = false;
Michael Beck's avatar
Michael Beck committed
688

689
	if (a != result) {
690
		result->desc = a->desc;
691
692
		result->clss = a->clss;
	}
Michael Beck's avatar
Michael Beck committed
693

Matthias Braun's avatar
Matthias Braun committed
694
695
	int exp_bias = (1 << (a->desc.exponent_size - 1)) - 1;
	int exp_val  = sc_val_to_long(_exp(a)) - exp_bias;
Michael Beck's avatar
Michael Beck committed
696
697
698
699
	if (exp_val < 0) {
		sc_val_from_ulong(0, NULL);
		_save_result(_exp(result));
		_save_result(_mant(result));
700
		result->clss = FC_ZERO;
Michael Beck's avatar
Michael Beck committed
701
702
703
		return;
	}

704
	if (exp_val > (long)a->desc.mantissa_size) {
Michael Beck's avatar
Michael Beck committed
705
706
707
708
709
710
711
		if (a != result)
			memcpy(result, a, calc_buffer_size);
		return;
	}

	/* set up a proper mask to delete all bits right of the
	 * radix point if the mantissa had been shifted until exp == 0 */
Matthias Braun's avatar
Matthias Braun committed
712
	char *temp = ALLOCAN(char, value_size);
Michael Beck's avatar
Michael Beck committed
713
714
715
716
717
718
719
	sc_max_from_bits(1 + exp_val, 0, temp);
	sc_val_from_long(a->desc.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));

720
721
722
723
	if (a != result) {
		memcpy(_exp(result), _exp(a), value_size);
		result->sign = a->sign;
	}
724
725
}

726
727
728
/********
 * functions defined in fltcalc.h
 ********/
729
730
const void *fc_get_buffer(void)
{
Michael Beck's avatar
Michael Beck committed
731
	return calc_buffer;
732
733
}

734
735
int fc_get_buffer_length(void)
{
Michael Beck's avatar
Michael Beck committed
736
	return calc_buffer_size;
737
738
}

739
void *fc_val_from_str(const char *str, size_t len, void *result)
740
{
741
	char *buffer = alloca(len + 1);
742
743
	memcpy(buffer, str, len);
	buffer[len] = '\0';
744
745
	long double val = string_to_long_double(buffer);
	return fc_val_from_ieee754(val, result);
746
747
}

748
fp_value *fc_val_from_ieee754(long double l, fp_value *result)
749
{
750
	value_t  srcval;
Michael Beck's avatar
Michael Beck committed
751
	srcval.d = l;
Matthias Braun's avatar
Matthias Braun committed
752
753
754
755
756
757
758
	int  bias_res = ((1 << (long_double_desc.exponent_size - 1)) - 1);
	int  bias_val;
	int  mant_val;
	char sign;
	uint32_t exponent;
	uint32_t mantissa0;
	uint32_t mantissa1;
759
760
	if (long_double_desc.exponent_size == 11 && long_double_desc.mantissa_size == 52) {
		assert(sizeof(long double) == 8);
761
762
763
764
765
766
		mant_val  = 52;
		bias_val  = 0x3ff;
		sign      = (srcval.val_ld8.high & 0x80000000) != 0;
		exponent  = (srcval.val_ld8.high & 0x7FF00000) >> 20;
		mantissa0 = srcval.val_ld8.high & 0x000FFFFF;
		mantissa1 = srcval.val_ld8.low;
767
	} else if (long_double_desc.exponent_size == 15 && long_double_desc.mantissa_size == 63) {
768
		/* we assume an x86-like 80bit representation of the value... */
769
		assert(sizeof(long double) == 12 || sizeof(long double) == 16);
770
771
772
773
774
775
		mant_val  = 63;
		bias_val  = 0x3fff;
		sign      = (srcval.val_ld12.high & 0x00008000) != 0;
		exponent  = (srcval.val_ld12.high & 0x00007FFF) ;
		mantissa0 = srcval.val_ld12.mid;
		mantissa1 = srcval.val_ld12.low;
776
777
	} else {
		panic("unsupported long double format");
778
	}
779

780
781
	if (result == NULL)
		result = calc_buffer;
Michael Beck's avatar
Michael Beck committed
782

Michael Beck's avatar
Michael Beck committed
783
	/* CLEAR the buffer, else some bits might be uninitialized */
784
	memset(result, 0, fc_get_buffer_length());
785

786
	result->desc = long_double_desc;
787
	result->clss = FC_NORMAL;
Michael Beck's avatar
Michael Beck committed
788
789
	result->sign = sign;

Michael Beck's avatar
Michael Beck committed
790
	/* sign and flag suffice to identify NaN or inf, no exponent/mantissa
Michael Beck's avatar
Michael Beck committed
791
	 * encoding is needed. the function can return immediately in these cases */
792
	if (my_isnan(l)) {
793
		result->clss = FC_NAN;
Michael Beck's avatar
Michael Beck committed
794
		return result;
795
	} else if (my_isinf(l)) {
796
		result->clss = FC_INF;
Michael Beck's avatar
Michael Beck committed
797
798
799
		return result;
	}

Matthias Braun's avatar
Matthias Braun committed
800
801
802
803
	/* 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 */
804
	sc_val_from_long((exponent - bias_val + bias_res) - (mant_val - long_double_desc.mantissa_size), _exp(result));
805

Michael Beck's avatar
Michael Beck committed
806
	/* build mantissa representation */
Matthias Braun's avatar
Matthias Braun committed
807
	char *temp = ALLOCAN(char, value_size);
Michael Beck's avatar
Michael Beck committed
808
809
810
	if (exponent != 0) {
		/* insert the hidden bit */
		sc_val_from_ulong(1, temp);
811
		sc_val_from_ulong(mant_val + ROUNDING_BITS, NULL);
Michael Beck's avatar
Michael Beck committed
812
		_shift_left(temp, sc_get_buffer(), NULL);
Matthias Braun's avatar
Matthias Braun committed
813
	} else {
Michael Beck's avatar
Michael Beck committed
814
815
		sc_val_from_ulong(0, NULL);
	}
816

Michael Beck's avatar
Michael Beck committed
817
	_save_result(_mant(result));
818

Michael Beck's avatar
Michael Beck committed
819
820
821
822
823
	/* 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));
824

Michael Beck's avatar
Michael Beck committed
825
826
	/* bits from the lower word */
	sc_val_from_ulong(mantissa1, temp);
827
	sc_val_from_ulong(ROUNDING_BITS, NULL);
Michael Beck's avatar
Michael Beck committed
828
829
	_shift_left(temp, sc_get_buffer(), temp);
	sc_or(_mant(result), temp, _mant(result));
830

Matthias Braun's avatar
Matthias Braun committed
831
832
	/* _normalize expects the radix point to be normal, so shift mantissa of
	 * subnormal origin one to the left */
Michael Beck's avatar
Michael Beck committed
833
834
835
836
	if (exponent == 0) {
		sc_val_from_ulong(1, NULL);
		_shift_left(_mant(result), sc_get_buffer(), _mant(result));
	}
837

838
	normalize(result, result, 0);
839

Michael Beck's avatar
Michael Beck committed
840
	return result;
841
842
}

843
long double fc_val_to_ieee754(const fp_value *val)
844
{
Matthias Braun's avatar
Matthias Braun committed
845
846
	unsigned mantissa_size
		= long_double_desc.mantissa_size + long_double_desc.explicit_one;
847

Matthias Braun's avatar
Matthias Braun committed
848
849
	fp_value *temp  = (fp_value*) alloca(calc_buffer_size);
	fp_value *value = fc_cast(val, &long_double_desc, temp);
850

Matthias Braun's avatar
Matthias Braun committed
851
	uint32_t sign = value->sign;
852

Matthias Braun's avatar
Matthias Braun committed
853
854
855
	/* @@@ long double exponent is 15bit, so the use of sc_val_to_long should
	 * not lead to wrong results */
	uint32_t exponent = sc_val_to_long(_exp(value)) ;
856

857
	sc_val_from_ulong(ROUNDING_BITS, NULL);
Michael Beck's avatar
Michael Beck committed
858
	_shift_right(_mant(value), sc_get_buffer(), _mant(value));
859

Matthias Braun's avatar
Matthias Braun committed
860
861
	uint32_t mantissa0 = 0;
	uint32_t mantissa1 = 0;
862

Matthias Braun's avatar
Matthias Braun committed
863
	unsigned byte_offset;
Michael Beck's avatar
Michael Beck committed
864
	for (byte_offset = 0; byte_offset < 4; byte_offset++)
865
		mantissa1 |= sc_sub_bits(_mant(value), mantissa_size, byte_offset) << (byte_offset << 3);
866

867
	for (; (byte_offset<<3) < long_double_desc.mantissa_size; byte_offset++)
868
		mantissa0 |= sc_sub_bits(_mant(value), mantissa_size, byte_offset) << ((byte_offset - 4) << 3);
869

Matthias Braun's avatar
Matthias Braun committed
870
	value_t buildval;
871
872
	if (long_double_desc.exponent_size == 11 && long_double_desc.mantissa_size == 52) {
		assert(sizeof(long double) == 8);
873
874
875
876
877
		mantissa0 &= 0x000FFFFF;  /* get rid of garbage */
		buildval.val_ld8.high = sign << 31;
		buildval.val_ld8.high |= exponent << 20;
		buildval.val_ld8.high |= mantissa0;
		buildval.val_ld8.low = mantissa1;
878
879
880
	} else if (long_double_desc.exponent_size == 15 && long_double_desc.mantissa_size == 63) {
		/* we assume an x86-like 80bit representation of the value... */
		assert(sizeof(long double) == 12 || sizeof(long double) == 16);
881
882
883
884
		buildval.val_ld12.high = sign << 15;
		buildval.val_ld12.high |= exponent;
		buildval.val_ld12.mid = mantissa0;
		buildval.val_ld12.low = mantissa1;
885
886
	} else {
		panic("unsupported long double format");
887
	}
888

Michael Beck's avatar
Michael Beck committed
889
	return buildval.d;
890
891
}

892
fp_value *fc_cast(const fp_value *value, const float_descriptor_t *desc,
893
                  fp_value *result)
894
{
895
896
897
	if (result == NULL)
		result = calc_buffer;
	assert(value != result);
Michael Beck's avatar
Michael Beck committed
898

Michael Beck's avatar
Michael Beck committed
899
900
901
	if (value->desc.exponent_size == desc->exponent_size &&
		value->desc.mantissa_size == desc->mantissa_size &&
		value->desc.explicit_one  == desc->explicit_one) {
Michael Beck's avatar
Michael Beck committed
902
903
904
905
906
		if (value != result)
			memcpy(result, value, calc_buffer_size);
		return result;
	}

907
	if (value->clss == FC_NAN) {
Michael Beck's avatar
Michael Beck committed
908
		if (sc_get_highest_set_bit(_mant(value)) == value->desc.mantissa_size + 1)
Michael Beck's avatar
Michael Beck committed
909
			return fc_get_qnan(desc, result);
Michael Beck's avatar
Michael Beck committed
910
		else
Michael Beck's avatar
Michael Beck committed
911
			return fc_get_snan(desc, result);
Michael Beck's avatar
Michael Beck committed
912
	}
913
	else if (value->clss == FC_INF) {
914
915
916
917
918
		if (value->sign == 0)
			return fc_get_plusinf(desc, result);
		else
			return fc_get_minusinf(desc, result);
	}
Michael Beck's avatar
Michael Beck committed
919

Michael Beck's avatar
Michael Beck committed
920
	/* set the descriptor of the new value */
921
922
	result->desc = *desc;
	result->clss = value->clss;