diff options
Diffstat (limited to 'sim/common/sim-fpu.c')
| -rw-r--r-- | sim/common/sim-fpu.c | 2578 |
1 files changed, 2578 insertions, 0 deletions
diff --git a/sim/common/sim-fpu.c b/sim/common/sim-fpu.c new file mode 100644 index 0000000..abf746a --- /dev/null +++ b/sim/common/sim-fpu.c @@ -0,0 +1,2578 @@ +/* This is a software floating point library which can be used instead + of the floating point routines in libgcc1.c for targets without + hardware floating point. */ + +/* Copyright (C) 1994,1997-1998 Free Software Foundation, Inc. + +This file is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +In addition to the permissions in the GNU General Public License, the +Free Software Foundation gives you unlimited permission to link the +compiled version of this file with other programs, and to distribute +those programs without any restriction coming from the use of this +file. (The General Public License restrictions do apply in other +respects; for example, they cover modification of the file, and +distribution when not linked into another program.) + +This file is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; see the file COPYING. If not, write to +the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */ + +/* As a special exception, if you link this library with other files, + some of which are compiled with GCC, to produce an executable, + this library does not by itself cause the resulting executable + to be covered by the GNU General Public License. + This exception does not however invalidate any other reasons why + the executable file might be covered by the GNU General Public License. */ + +/* This implements IEEE 754 format arithmetic, but does not provide a + mechanism for setting the rounding mode, or for generating or handling + exceptions. + + The original code by Steve Chamberlain, hacked by Mark Eichin and Jim + Wilson, all of Cygnus Support. */ + + +#ifndef SIM_FPU_C +#define SIM_FPU_C + +#include "sim-basics.h" +#include "sim-fpu.h" + +#include "sim-io.h" +#include "sim-assert.h" + + +/* Debugging support. */ + +static void +print_bits (unsigned64 x, + int msbit, + sim_fpu_print_func print, + void *arg) +{ + unsigned64 bit = LSBIT64 (msbit); + int i = 4; + while (bit) + { + if (i == 0) + print (arg, ","); + if ((x & bit)) + print (arg, "1"); + else + print (arg, "0"); + bit >>= 1; + i = (i + 1) % 4; + } +} + + + +/* Quick and dirty conversion between a host double and host 64bit int */ + +typedef union { + double d; + unsigned64 i; +} sim_fpu_map; + + +/* A packed IEEE floating point number. + + Form is <SIGN:1><BIASEDEXP:NR_EXPBITS><FRAC:NR_FRACBITS> for both + 32 and 64 bit numbers. This number is interpreted as: + + Normalized (0 < BIASEDEXP && BIASEDEXP < EXPMAX): + (sign ? '-' : '+') 1.<FRAC> x 2 ^ (BIASEDEXP - EXPBIAS) + + Denormalized (0 == BIASEDEXP && FRAC != 0): + (sign ? "-" : "+") 0.<FRAC> x 2 ^ (- EXPBIAS) + + Zero (0 == BIASEDEXP && FRAC == 0): + (sign ? "-" : "+") 0.0 + + Infinity (BIASEDEXP == EXPMAX && FRAC == 0): + (sign ? "-" : "+") "infinity" + + SignalingNaN (BIASEDEXP == EXPMAX && FRAC > 0 && FRAC < QUIET_NAN): + SNaN.FRAC + + QuietNaN (BIASEDEXP == EXPMAX && FRAC > 0 && FRAC > QUIET_NAN): + QNaN.FRAC + + */ + +#define NR_EXPBITS (is_double ? 11 : 8) +#define NR_FRACBITS (is_double ? 52 : 23) +#define SIGNBIT (is_double ? MSBIT64 (0) : MSBIT64 (32)) + +#define EXPMAX32 (255) +#define EXMPAX64 (2047) +#define EXPMAX ((unsigned) (is_double ? EXMPAX64 : EXPMAX32)) + +#define EXPBIAS32 (127) +#define EXPBIAS64 (1023) +#define EXPBIAS (is_double ? EXPBIAS64 : EXPBIAS32) + +#define QUIET_NAN LSBIT64 (NR_FRACBITS - 1) + + + +/* An unpacked floating point number. + + When unpacked, the fraction of both a 32 and 64 bit floating point + number is stored using the same format: + + 64 bit - <IMPLICIT_1:1><FRACBITS:52><GUARDS:8><PAD:00> + 32 bit - <IMPLICIT_1:1><FRACBITS:23><GUARDS:7><PAD:30> */ + +#define NR_PAD32 (30) +#define NR_PAD64 (0) +#define NR_PAD (is_double ? NR_PAD64 : NR_PAD32) +#define PADMASK (is_double ? 0 : LSMASK64 (NR_PAD32 - 1, 0)) + +#define NR_GUARDS32 (7 + NR_PAD32) +#define NR_GUARDS64 (8 + NR_PAD64) +#define NR_GUARDS (is_double ? NR_GUARDS64 : NR_GUARDS32) +#define GUARDMASK LSMASK64 (NR_GUARDS - 1, 0) + +#define GUARDMSB LSBIT64 (NR_GUARDS - 1) +#define GUARDLSB LSBIT64 (NR_PAD) +#define GUARDROUND LSMASK64 (NR_GUARDS - 2, 0) + +#define NR_FRAC_GUARD (60) +#define IMPLICIT_1 LSBIT64 (NR_FRAC_GUARD) +#define IMPLICIT_2 LSBIT64 (NR_FRAC_GUARD + 1) +#define IMPLICIT_4 LSBIT64 (NR_FRAC_GUARD + 2) +#define NR_SPARE 2 + +#define FRAC32MASK LSMASK64 (63, NR_FRAC_GUARD - 32 + 1) + +#define NORMAL_EXPMIN (-(EXPBIAS)+1) + +#define NORMAL_EXPMAX32 (EXPBIAS32) +#define NORMAL_EXPMAX64 (EXPBIAS64) +#define NORMAL_EXPMAX (EXPBIAS) + + +/* Integer constants */ + +#define MAX_INT32 ((signed64) LSMASK64 (30, 0)) +#define MAX_UINT32 LSMASK64 (31, 0) +#define MIN_INT32 ((signed64) LSMASK64 (63, 31)) + +#define MAX_INT64 ((signed64) LSMASK64 (62, 0)) +#define MAX_UINT64 LSMASK64 (63, 0) +#define MIN_INT64 ((signed64) LSMASK64 (63, 63)) + +#define MAX_INT (is_64bit ? MAX_INT64 : MAX_INT32) +#define MIN_INT (is_64bit ? MIN_INT64 : MIN_INT32) +#define MAX_UINT (is_64bit ? MAX_UINT64 : MAX_UINT32) +#define NR_INTBITS (is_64bit ? 64 : 32) + +/* Squeese an unpacked sim_fpu struct into a 32/64 bit integer */ +STATIC_INLINE_SIM_FPU (unsigned64) +pack_fpu (const sim_fpu *src, + int is_double) +{ + int sign; + unsigned64 exp; + unsigned64 fraction; + unsigned64 packed; + + switch (src->class) + { + /* create a NaN */ + case sim_fpu_class_qnan: + sign = src->sign; + exp = EXPMAX; + /* force fraction to correct class */ + fraction = src->fraction; + fraction >>= NR_GUARDS; + fraction |= QUIET_NAN; + break; + case sim_fpu_class_snan: + sign = src->sign; + exp = EXPMAX; + /* force fraction to correct class */ + fraction = src->fraction; + fraction >>= NR_GUARDS; + fraction &= ~QUIET_NAN; + break; + case sim_fpu_class_infinity: + sign = src->sign; + exp = EXPMAX; + fraction = 0; + break; + case sim_fpu_class_zero: + sign = src->sign; + exp = 0; + fraction = 0; + break; + case sim_fpu_class_number: + case sim_fpu_class_denorm: + ASSERT (src->fraction >= IMPLICIT_1); + ASSERT (src->fraction < IMPLICIT_2); + if (src->normal_exp < NORMAL_EXPMIN) + { + /* This number's exponent is too low to fit into the bits + available in the number We'll denormalize the number by + storing zero in the exponent and shift the fraction to + the right to make up for it. */ + int nr_shift = NORMAL_EXPMIN - src->normal_exp; + if (nr_shift > NR_FRACBITS) + { + /* underflow, just make the number zero */ + sign = src->sign; + exp = 0; + fraction = 0; + } + else + { + sign = src->sign; + exp = 0; + /* Shift by the value */ + fraction = src->fraction; + fraction >>= NR_GUARDS; + fraction >>= nr_shift; + } + } + else if (src->normal_exp > NORMAL_EXPMAX) + { + /* Infinity */ + sign = src->sign; + exp = EXPMAX; + fraction = 0; + } + else + { + exp = (src->normal_exp + EXPBIAS); + sign = src->sign; + fraction = src->fraction; + /* FIXME: Need to round according to WITH_SIM_FPU_ROUNDING + or some such */ + /* Round to nearest: If the guard bits are the all zero, but + the first, then we're half way between two numbers, + choose the one which makes the lsb of the answer 0. */ + if ((fraction & GUARDMASK) == GUARDMSB) + { + if ((fraction & (GUARDMSB << 1))) + fraction += (GUARDMSB << 1); + } + else + { + /* Add a one to the guards to force round to nearest */ + fraction += GUARDROUND; + } + if ((fraction & IMPLICIT_2)) /* rounding resulted in carry */ + { + exp += 1; + fraction >>= 1; + } + fraction >>= NR_GUARDS; + /* When exp == EXPMAX (overflow from carry) fraction must + have been made zero */ + ASSERT ((exp == EXPMAX) <= ((fraction & ~IMPLICIT_1) == 0)); + } + break; + default: + abort (); + } + + packed = ((sign ? SIGNBIT : 0) + | (exp << NR_FRACBITS) + | LSMASKED64 (fraction, NR_FRACBITS - 1, 0)); + + /* trace operation */ +#if 0 + if (is_double) + { + } + else + { + printf ("pack_fpu: "); + printf ("-> %c%0lX.%06lX\n", + LSMASKED32 (packed, 31, 31) ? '8' : '0', + (long) LSEXTRACTED32 (packed, 30, 23), + (long) LSEXTRACTED32 (packed, 23 - 1, 0)); + } +#endif + + return packed; +} + + +/* Unpack a 32/64 bit integer into a sim_fpu structure */ +STATIC_INLINE_SIM_FPU (void) +unpack_fpu (sim_fpu *dst, unsigned64 packed, int is_double) +{ + unsigned64 fraction = LSMASKED64 (packed, NR_FRACBITS - 1, 0); + unsigned exp = LSEXTRACTED64 (packed, NR_EXPBITS + NR_FRACBITS - 1, NR_FRACBITS); + int sign = (packed & SIGNBIT) != 0; + + if (exp == 0) + { + /* Hmm. Looks like 0 */ + if (fraction == 0) + { + /* tastes like zero */ + dst->class = sim_fpu_class_zero; + dst->sign = sign; + } + else + { + /* Zero exponent with non zero fraction - it's denormalized, + so there isn't a leading implicit one - we'll shift it so + it gets one. */ + dst->normal_exp = exp - EXPBIAS + 1; + dst->class = sim_fpu_class_denorm; + dst->sign = sign; + fraction <<= NR_GUARDS; + while (fraction < IMPLICIT_1) + { + fraction <<= 1; + dst->normal_exp--; + } + dst->fraction = fraction; + } + } + else if (exp == EXPMAX) + { + /* Huge exponent*/ + if (fraction == 0) + { + /* Attached to a zero fraction - means infinity */ + dst->class = sim_fpu_class_infinity; + dst->sign = sign; + /* dst->normal_exp = EXPBIAS; */ + /* dst->fraction = 0; */ + } + else + { + /* Non zero fraction, means NaN */ + dst->sign = sign; + dst->fraction = (fraction << NR_GUARDS); + if (fraction >= QUIET_NAN) + dst->class = sim_fpu_class_qnan; + else + dst->class = sim_fpu_class_snan; + } + } + else + { + /* Nothing strange about this number */ + dst->class = sim_fpu_class_number; + dst->sign = sign; + dst->fraction = ((fraction << NR_GUARDS) | IMPLICIT_1); + dst->normal_exp = exp - EXPBIAS; + } + + /* trace operation */ +#if 0 + if (is_double) + { + } + else + { + printf ("unpack_fpu: %c%02lX.%06lX ->\n", + LSMASKED32 (packed, 31, 31) ? '8' : '0', + (long) LSEXTRACTED32 (packed, 30, 23), + (long) LSEXTRACTED32 (packed, 23 - 1, 0)); + } +#endif + + /* sanity checks */ + { + sim_fpu_map val; + val.i = pack_fpu (dst, 1); + if (is_double) + { + ASSERT (val.i == packed); + } + else + { + unsigned32 val = pack_fpu (dst, 0); + unsigned32 org = packed; + ASSERT (val == org); + } + } +} + + +/* Convert a floating point into an integer */ +STATIC_INLINE_SIM_FPU (int) +fpu2i (signed64 *i, + const sim_fpu *s, + int is_64bit, + sim_fpu_round round) +{ + unsigned64 tmp; + int shift; + int status = 0; + if (sim_fpu_is_zero (s)) + { + *i = 0; + return 0; + } + if (sim_fpu_is_snan (s)) + { + *i = MIN_INT; /* FIXME */ + return sim_fpu_status_invalid_cvi; + } + if (sim_fpu_is_qnan (s)) + { + *i = MIN_INT; /* FIXME */ + return sim_fpu_status_invalid_cvi; + } + /* map infinity onto MAX_INT... */ + if (sim_fpu_is_infinity (s)) + { + *i = s->sign ? MIN_INT : MAX_INT; + return sim_fpu_status_invalid_cvi; + } + /* it is a number, but a small one */ + if (s->normal_exp < 0) + { + *i = 0; + return sim_fpu_status_inexact; + } + /* Is the floating point MIN_INT or just close? */ + if (s->sign && s->normal_exp == (NR_INTBITS - 1)) + { + *i = MIN_INT; + ASSERT (s->fraction >= IMPLICIT_1); + if (s->fraction == IMPLICIT_1) + return 0; /* exact */ + if (is_64bit) /* can't round */ + return sim_fpu_status_invalid_cvi; /* must be overflow */ + /* For a 32bit with MAX_INT, rounding is possible */ + switch (round) + { + case sim_fpu_round_default: + abort (); + case sim_fpu_round_zero: + if ((s->fraction & FRAC32MASK) != IMPLICIT_1) + return sim_fpu_status_invalid_cvi; + else + return sim_fpu_status_inexact; + break; + case sim_fpu_round_near: + { + if ((s->fraction & FRAC32MASK) != IMPLICIT_1) + return sim_fpu_status_invalid_cvi; + else if ((s->fraction & !FRAC32MASK) >= (~FRAC32MASK >> 1)) + return sim_fpu_status_invalid_cvi; + else + return sim_fpu_status_inexact; + } + case sim_fpu_round_up: + if ((s->fraction & FRAC32MASK) == IMPLICIT_1) + return sim_fpu_status_inexact; + else + return sim_fpu_status_invalid_cvi; + case sim_fpu_round_down: + return sim_fpu_status_invalid_cvi; + } + } + /* Would right shifting result in the FRAC being shifted into + (through) the integer's sign bit? */ + if (s->normal_exp > (NR_INTBITS - 2)) + { + *i = s->sign ? MIN_INT : MAX_INT; + return sim_fpu_status_invalid_cvi; + } + /* normal number shift it into place */ + tmp = s->fraction; + shift = (s->normal_exp - (NR_FRAC_GUARD)); + if (shift > 0) + { + tmp <<= shift; + } + else + { + shift = -shift; + if (tmp & ((SIGNED64 (1) << shift) - 1)) + status |= sim_fpu_status_inexact; + tmp >>= shift; + } + *i = s->sign ? (-tmp) : (tmp); + return status; +} + +/* convert an integer into a floating point */ +STATIC_INLINE_SIM_FPU (int) +i2fpu (sim_fpu *f, signed64 i, int is_64bit) +{ + int status = 0; + if (i == 0) + { + f->class = sim_fpu_class_zero; + f->sign = 0; + } + else + { + f->class = sim_fpu_class_number; + f->sign = (i < 0); + f->normal_exp = NR_FRAC_GUARD; + + if (f->sign) + { + /* Special case for minint, since there is no corresponding + +ve integer representation for it */ + if (i == MIN_INT) + { + f->fraction = IMPLICIT_1; + f->normal_exp = NR_INTBITS - 1; + } + else + f->fraction = (-i); + } + else + f->fraction = i; + + if (f->fraction >= IMPLICIT_2) + { + do + { + f->fraction >>= 1; + f->normal_exp += 1; + } + while (f->fraction >= IMPLICIT_2); + } + else if (f->fraction < IMPLICIT_1) + { + do + { + f->fraction <<= 1; + f->normal_exp -= 1; + } + while (f->fraction < IMPLICIT_1); + } + } + + /* trace operation */ +#if 0 + { + printf ("i2fpu: 0x%08lX ->\n", (long) i); + } +#endif + + /* sanity check */ + { + signed64 val; + fpu2i (&val, f, is_64bit, sim_fpu_round_zero); + if (i >= MIN_INT32 && i <= MAX_INT32) + { + ASSERT (val == i); + } + } + + return status; +} + + +/* Convert a floating point into an integer */ +STATIC_INLINE_SIM_FPU (int) +fpu2u (unsigned64 *u, const sim_fpu *s, int is_64bit) +{ + const int is_double = 1; + unsigned64 tmp; + int shift; + if (sim_fpu_is_zero (s)) + { + *u = 0; + return 0; + } + if (sim_fpu_is_nan (s)) + { + *u = 0; + return 0; + } + /* it is a negative number */ + if (s->sign) + { + *u = 0; + return 0; + } + /* get reasonable MAX_USI_INT... */ + if (sim_fpu_is_infinity (s)) + { + *u = MAX_UINT; + return 0; + } + /* it is a number, but a small one */ + if (s->normal_exp < 0) + { + *u = 0; + return 0; + } + /* overflow */ + if (s->normal_exp > (NR_INTBITS - 1)) + { + *u = MAX_UINT; + return 0; + } + /* normal number */ + tmp = (s->fraction & ~PADMASK); + shift = (s->normal_exp - (NR_FRACBITS + NR_GUARDS)); + if (shift > 0) + { + tmp <<= shift; + } + else + { + shift = -shift; + tmp >>= shift; + } + *u = tmp; + return 0; +} + +/* Convert an unsigned integer into a floating point */ +STATIC_INLINE_SIM_FPU (int) +u2fpu (sim_fpu *f, unsigned64 u, int is_64bit) +{ + if (u == 0) + { + f->class = sim_fpu_class_zero; + f->sign = 0; + } + else + { + f->class = sim_fpu_class_number; + f->sign = 0; + f->normal_exp = NR_FRAC_GUARD; + f->fraction = u; + + while (f->fraction < IMPLICIT_1) + { + f->fraction <<= 1; + f->normal_exp -= 1; + } + } + return 0; +} + + +/* register <-> sim_fpu */ + +INLINE_SIM_FPU (void) +sim_fpu_32to (sim_fpu *f, unsigned32 s) +{ + unpack_fpu (f, s, 0); +} + + +INLINE_SIM_FPU (void) +sim_fpu_232to (sim_fpu *f, unsigned32 h, unsigned32 l) +{ + unsigned64 s = h; + s = (s << 32) | l; + unpack_fpu (f, s, 1); +} + + +INLINE_SIM_FPU (void) +sim_fpu_64to (sim_fpu *f, unsigned64 s) +{ + unpack_fpu (f, s, 1); +} + + +INLINE_SIM_FPU (void) +sim_fpu_to32 (unsigned32 *s, + const sim_fpu *f) +{ + *s = pack_fpu (f, 0); +} + + +INLINE_SIM_FPU (void) +sim_fpu_to232 (unsigned32 *h, unsigned32 *l, + const sim_fpu *f) +{ + unsigned64 s = pack_fpu (f, 1); + *l = s; + *h = (s >> 32); +} + + +INLINE_SIM_FPU (void) +sim_fpu_to64 (unsigned64 *u, + const sim_fpu *f) +{ + *u = pack_fpu (f, 1); +} + + +INLINE_SIM_FPU (void) +sim_fpu_fractionto (sim_fpu *f, + int sign, + int normal_exp, + unsigned64 fraction, + int precision) +{ + int shift = (NR_FRAC_GUARD - precision); + f->class = sim_fpu_class_number; + f->sign = sign; + f->normal_exp = normal_exp; + /* shift the fraction to where sim-fpu expects it */ + if (shift >= 0) + f->fraction = (fraction << shift); + else + f->fraction = (fraction >> -shift); + f->fraction |= IMPLICIT_1; +} + + +INLINE_SIM_FPU (unsigned64) +sim_fpu_tofraction (const sim_fpu *d, + int precision) +{ + /* we have NR_FRAC_GUARD bits, we want only PRECISION bits */ + int shift = (NR_FRAC_GUARD - precision); + unsigned64 fraction = (d->fraction & ~IMPLICIT_1); + if (shift >= 0) + return fraction >> shift; + else + return fraction << -shift; +} + + +/* Rounding */ + +STATIC_INLINE_SIM_FPU (int) +do_normal_overflow (sim_fpu *f, + int is_double, + sim_fpu_round round) +{ + switch (round) + { + case sim_fpu_round_default: + return 0; + case sim_fpu_round_near: + f->class = sim_fpu_class_infinity; + break; + case sim_fpu_round_up: + if (!f->sign) + f->class = sim_fpu_class_infinity; + break; + case sim_fpu_round_down: + if (f->sign) + f->class = sim_fpu_class_infinity; + break; + case sim_fpu_round_zero: + break; + } + f->normal_exp = NORMAL_EXPMAX; + f->fraction = LSMASK64 (NR_FRAC_GUARD, NR_GUARDS); + return (sim_fpu_status_overflow | sim_fpu_status_inexact); +} + +STATIC_INLINE_SIM_FPU (int) +do_normal_underflow (sim_fpu *f, + int is_double, + sim_fpu_round round) +{ + switch (round) + { + case sim_fpu_round_default: + return 0; + case sim_fpu_round_near: + f->class = sim_fpu_class_zero; + break; + case sim_fpu_round_up: + if (f->sign) + f->class = sim_fpu_class_zero; + break; + case sim_fpu_round_down: + if (!f->sign) + f->class = sim_fpu_class_zero; + break; + case sim_fpu_round_zero: + f->class = sim_fpu_class_zero; + break; + } + f->normal_exp = NORMAL_EXPMIN - NR_FRACBITS; + f->fraction = IMPLICIT_1; + return (sim_fpu_status_inexact | sim_fpu_status_underflow); +} + + + +/* Round a number using NR_GUARDS. + Will return the rounded number or F->FRACTION == 0 when underflow */ + +STATIC_INLINE_SIM_FPU (int) +do_normal_round (sim_fpu *f, + int nr_guards, + sim_fpu_round round) +{ + unsigned64 guardmask = LSMASK64 (nr_guards - 1, 0); + unsigned64 guardmsb = LSBIT64 (nr_guards - 1); + unsigned64 fraclsb = guardmsb << 1; + if ((f->fraction & guardmask)) + { + int status = sim_fpu_status_inexact; + switch (round) + { + case sim_fpu_round_default: + return 0; + case sim_fpu_round_near: + if ((f->fraction & guardmsb)) + { + if ((f->fraction & fraclsb)) + { + status |= sim_fpu_status_rounded; + } + else if ((f->fraction & (guardmask >> 1))) + { + status |= sim_fpu_status_rounded; + } + } + break; + case sim_fpu_round_up: + if (!f->sign) + status |= sim_fpu_status_rounded; + break; + case sim_fpu_round_down: + if (f->sign) + status |= sim_fpu_status_rounded; + break; + case sim_fpu_round_zero: + break; + } + f->fraction &= ~guardmask; + /* round if needed, handle resulting overflow */ + if ((status & sim_fpu_status_rounded)) + { + f->fraction += fraclsb; + if ((f->fraction & IMPLICIT_2)) + { + f->fraction >>= 1; + f->normal_exp += 1; + } + } + return status; + } + else + return 0; +} + + +STATIC_INLINE_SIM_FPU (int) +do_round (sim_fpu *f, + int is_double, + sim_fpu_round round, + sim_fpu_denorm denorm) +{ + switch (f->class) + { + case sim_fpu_class_qnan: + case sim_fpu_class_zero: + case sim_fpu_class_infinity: + return 0; + break; + case sim_fpu_class_snan: + /* Quieten a SignalingNaN */ + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + break; + case sim_fpu_class_number: + case sim_fpu_class_denorm: + { + int status; + ASSERT (f->fraction < IMPLICIT_2); + ASSERT (f->fraction >= IMPLICIT_1); + if (f->normal_exp < NORMAL_EXPMIN) + { + /* This number's exponent is too low to fit into the bits + available in the number. Round off any bits that will be + discarded as a result of denormalization. Edge case is + the implicit bit shifted to GUARD0 and then rounded + up. */ + int shift = NORMAL_EXPMIN - f->normal_exp; + if (shift + NR_GUARDS <= NR_FRAC_GUARD + 1 + && !(denorm & sim_fpu_denorm_zero)) + { + status = do_normal_round (f, shift + NR_GUARDS, round); + if (f->fraction == 0) /* rounding underflowed */ + { + status |= do_normal_underflow (f, is_double, round); + } + else if (f->normal_exp < NORMAL_EXPMIN) /* still underflow? */ + { + status |= sim_fpu_status_denorm; + /* Any loss of precision when denormalizing is + underflow. Some processors check for underflow + before rounding, some after! */ + if (status & sim_fpu_status_inexact) + status |= sim_fpu_status_underflow; + /* Flag that resultant value has been denormalized */ + f->class = sim_fpu_class_denorm; + } + else if ((denorm & sim_fpu_denorm_underflow_inexact)) + { + if ((status & sim_fpu_status_inexact)) + status |= sim_fpu_status_underflow; + } + } + else + { + status = do_normal_underflow (f, is_double, round); + } + } + else if (f->normal_exp > NORMAL_EXPMAX) + { + /* Infinity */ + status = do_normal_overflow (f, is_double, round); + } + else + { + status = do_normal_round (f, NR_GUARDS, round); + if (f->fraction == 0) + /* f->class = sim_fpu_class_zero; */ + status |= do_normal_underflow (f, is_double, round); + else if (f->normal_exp > NORMAL_EXPMAX) + /* oops! rounding caused overflow */ + status |= do_normal_overflow (f, is_double, round); + } + ASSERT ((f->class == sim_fpu_class_number + || f->class == sim_fpu_class_denorm) + <= (f->fraction < IMPLICIT_2 && f->fraction >= IMPLICIT_1)); + return status; + } + } + return 0; +} + +INLINE_SIM_FPU (int) +sim_fpu_round_32 (sim_fpu *f, + sim_fpu_round round, + sim_fpu_denorm denorm) +{ + return do_round (f, 0, round, denorm); +} + +INLINE_SIM_FPU (int) +sim_fpu_round_64 (sim_fpu *f, + sim_fpu_round round, + sim_fpu_denorm denorm) +{ + return do_round (f, 1, round, denorm); +} + + + +/* Arithmetic ops */ + +INLINE_SIM_FPU (int) +sim_fpu_add (sim_fpu *f, + const sim_fpu *l, + const sim_fpu *r) +{ + if (sim_fpu_is_snan (l)) + { + *f = *l; + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_snan (r)) + { + *f = *r; + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_qnan (l)) + { + *f = *l; + return 0; + } + if (sim_fpu_is_qnan (r)) + { + *f = *r; + return 0; + } + if (sim_fpu_is_infinity (l)) + { + if (sim_fpu_is_infinity (r) + && l->sign != r->sign) + { + *f = sim_fpu_qnan; + return sim_fpu_status_invalid_isi; + } + *f = *l; + return 0; + } + if (sim_fpu_is_infinity (r)) + { + *f = *r; + return 0; + } + if (sim_fpu_is_zero (l)) + { + if (sim_fpu_is_zero (r)) + { + *f = sim_fpu_zero; + f->sign = l->sign & r->sign; + } + else + *f = *r; + return 0; + } + if (sim_fpu_is_zero (r)) + { + *f = *l; + return 0; + } + { + int status = 0; + int shift = l->normal_exp - r->normal_exp; + unsigned64 lfraction; + unsigned64 rfraction; + /* use exp of larger */ + if (shift >= NR_FRAC_GUARD) + { + /* left has much bigger magnitute */ + *f = *l; + return sim_fpu_status_inexact; + } + if (shift <= - NR_FRAC_GUARD) + { + /* right has much bigger magnitute */ + *f = *r; + return sim_fpu_status_inexact; + } + lfraction = l->fraction; + rfraction = r->fraction; + if (shift > 0) + { + f->normal_exp = l->normal_exp; + if (rfraction & LSMASK64 (shift - 1, 0)) + { + status |= sim_fpu_status_inexact; + rfraction |= LSBIT64 (shift); /* stick LSBit */ + } + rfraction >>= shift; + } + else if (shift < 0) + { + f->normal_exp = r->normal_exp; + if (lfraction & LSMASK64 (- shift - 1, 0)) + { + status |= sim_fpu_status_inexact; + lfraction |= LSBIT64 (- shift); /* stick LSBit */ + } + lfraction >>= -shift; + } + else + { + f->normal_exp = r->normal_exp; + } + + /* perform the addition */ + if (l->sign) + lfraction = - lfraction; + if (r->sign) + rfraction = - rfraction; + f->fraction = lfraction + rfraction; + + /* zero? */ + if (f->fraction == 0) + { + *f = sim_fpu_zero; + return 0; + } + + /* sign? */ + f->class = sim_fpu_class_number; + if ((signed64) f->fraction >= 0) + f->sign = 0; + else + { + f->sign = 1; + f->fraction = - f->fraction; + } + + /* normalize it */ + if ((f->fraction & IMPLICIT_2)) + { + f->fraction = (f->fraction >> 1) | (f->fraction & 1); + f->normal_exp ++; + } + else if (f->fraction < IMPLICIT_1) + { + do + { + f->fraction <<= 1; + f->normal_exp --; + } + while (f->fraction < IMPLICIT_1); + } + ASSERT (f->fraction >= IMPLICIT_1 && f->fraction < IMPLICIT_2); + return status; + } +} + + +INLINE_SIM_FPU (int) +sim_fpu_sub (sim_fpu *f, + const sim_fpu *l, + const sim_fpu *r) +{ + if (sim_fpu_is_snan (l)) + { + *f = *l; + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_snan (r)) + { + *f = *r; + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_qnan (l)) + { + *f = *l; + return 0; + } + if (sim_fpu_is_qnan (r)) + { + *f = *r; + return 0; + } + if (sim_fpu_is_infinity (l)) + { + if (sim_fpu_is_infinity (r) + && l->sign == r->sign) + { + *f = sim_fpu_qnan; + return sim_fpu_status_invalid_isi; + } + *f = *l; + return 0; + } + if (sim_fpu_is_infinity (r)) + { + *f = *r; + f->sign = !r->sign; + return 0; + } + if (sim_fpu_is_zero (l)) + { + if (sim_fpu_is_zero (r)) + { + *f = sim_fpu_zero; + f->sign = l->sign & !r->sign; + } + else + { + *f = *r; + f->sign = !r->sign; + } + return 0; + } + if (sim_fpu_is_zero (r)) + { + *f = *l; + return 0; + } + { + int status = 0; + int shift = l->normal_exp - r->normal_exp; + unsigned64 lfraction; + unsigned64 rfraction; + /* use exp of larger */ + if (shift >= NR_FRAC_GUARD) + { + /* left has much bigger magnitute */ + *f = *l; + return sim_fpu_status_inexact; + } + if (shift <= - NR_FRAC_GUARD) + { + /* right has much bigger magnitute */ + *f = *r; + f->sign = !r->sign; + return sim_fpu_status_inexact; + } + lfraction = l->fraction; + rfraction = r->fraction; + if (shift > 0) + { + f->normal_exp = l->normal_exp; + if (rfraction & LSMASK64 (shift - 1, 0)) + { + status |= sim_fpu_status_inexact; + rfraction |= LSBIT64 (shift); /* stick LSBit */ + } + rfraction >>= shift; + } + else if (shift < 0) + { + f->normal_exp = r->normal_exp; + if (lfraction & LSMASK64 (- shift - 1, 0)) + { + status |= sim_fpu_status_inexact; + lfraction |= LSBIT64 (- shift); /* stick LSBit */ + } + lfraction >>= -shift; + } + else + { + f->normal_exp = r->normal_exp; + } + + /* perform the subtraction */ + if (l->sign) + lfraction = - lfraction; + if (!r->sign) + rfraction = - rfraction; + f->fraction = lfraction + rfraction; + + /* zero? */ + if (f->fraction == 0) + { + *f = sim_fpu_zero; + return 0; + } + + /* sign? */ + f->class = sim_fpu_class_number; + if ((signed64) f->fraction >= 0) + f->sign = 0; + else + { + f->sign = 1; + f->fraction = - f->fraction; + } + + /* normalize it */ + if ((f->fraction & IMPLICIT_2)) + { + f->fraction = (f->fraction >> 1) | (f->fraction & 1); + f->normal_exp ++; + } + else if (f->fraction < IMPLICIT_1) + { + do + { + f->fraction <<= 1; + f->normal_exp --; + } + while (f->fraction < IMPLICIT_1); + } + ASSERT (f->fraction >= IMPLICIT_1 && f->fraction < IMPLICIT_2); + return status; + } +} + + +INLINE_SIM_FPU (int) +sim_fpu_mul (sim_fpu *f, + const sim_fpu *l, + const sim_fpu *r) +{ + if (sim_fpu_is_snan (l)) + { + *f = *l; + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_snan (r)) + { + *f = *r; + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_qnan (l)) + { + *f = *l; + return 0; + } + if (sim_fpu_is_qnan (r)) + { + *f = *r; + return 0; + } + if (sim_fpu_is_infinity (l)) + { + if (sim_fpu_is_zero (r)) + { + *f = sim_fpu_qnan; + return sim_fpu_status_invalid_imz; + } + *f = *l; + f->sign = l->sign ^ r->sign; + return 0; + } + if (sim_fpu_is_infinity (r)) + { + if (sim_fpu_is_zero (l)) + { + *f = sim_fpu_qnan; + return sim_fpu_status_invalid_imz; + } + *f = *r; + f->sign = l->sign ^ r->sign; + return 0; + } + if (sim_fpu_is_zero (l) || sim_fpu_is_zero (r)) + { + *f = sim_fpu_zero; + f->sign = l->sign ^ r->sign; + return 0; + } + /* Calculate the mantissa by multiplying both 64bit numbers to get a + 128 bit number */ + { + unsigned64 low; + unsigned64 high; + unsigned64 nl = l->fraction & 0xffffffff; + unsigned64 nh = l->fraction >> 32; + unsigned64 ml = r->fraction & 0xffffffff; + unsigned64 mh = r->fraction >>32; + unsigned64 pp_ll = ml * nl; + unsigned64 pp_hl = mh * nl; + unsigned64 pp_lh = ml * nh; + unsigned64 pp_hh = mh * nh; + unsigned64 res2 = 0; + unsigned64 res0 = 0; + unsigned64 ps_hh__ = pp_hl + pp_lh; + if (ps_hh__ < pp_hl) + res2 += UNSIGNED64 (0x100000000); + pp_hl = (ps_hh__ << 32) & UNSIGNED64 (0xffffffff00000000); + res0 = pp_ll + pp_hl; + if (res0 < pp_ll) + res2++; + res2 += ((ps_hh__ >> 32) & 0xffffffff) + pp_hh; + high = res2; + low = res0; + + f->normal_exp = l->normal_exp + r->normal_exp; + f->sign = l->sign ^ r->sign; + f->class = sim_fpu_class_number; + + /* Input is bounded by [1,2) ; [2^60,2^61) + Output is bounded by [1,4) ; [2^120,2^122) */ + + /* Adjust the exponent according to where the decimal point ended + up in the high 64 bit word. In the source the decimal point + was at NR_FRAC_GUARD. */ + f->normal_exp += NR_FRAC_GUARD + 64 - (NR_FRAC_GUARD * 2); + + /* The high word is bounded according to the above. Consequently + it has never overflowed into IMPLICIT_2. */ + ASSERT (high < LSBIT64 (((NR_FRAC_GUARD + 1) * 2) - 64)); + ASSERT (high >= LSBIT64 ((NR_FRAC_GUARD * 2) - 64)); + ASSERT (LSBIT64 (((NR_FRAC_GUARD + 1) * 2) - 64) < IMPLICIT_1); + +#if 0 + printf ("\n"); + print_bits (high, 63, (sim_fpu_print_func*)fprintf, stdout); + printf (";"); + print_bits (low, 63, (sim_fpu_print_func*)fprintf, stdout); + printf ("\n"); +#endif + + /* normalize */ + do + { + f->normal_exp--; + high <<= 1; + if (low & LSBIT64 (63)) + high |= 1; + low <<= 1; + } + while (high < IMPLICIT_1); + +#if 0 + print_bits (high, 63, (sim_fpu_print_func*)fprintf, stdout); + printf (";"); + print_bits (low, 63, (sim_fpu_print_func*)fprintf, stdout); + printf ("\n"); +#endif + + ASSERT (high >= IMPLICIT_1 && high < IMPLICIT_2); + if (low != 0) + { + f->fraction = (high | 1); /* sticky */ + return sim_fpu_status_inexact; + } + else + { + f->fraction = high; + return 0; + } + return 0; + } +} + +INLINE_SIM_FPU (int) +sim_fpu_div (sim_fpu *f, + const sim_fpu *l, + const sim_fpu *r) +{ + if (sim_fpu_is_snan (l)) + { + *f = *l; + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_snan (r)) + { + *f = *r; + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_qnan (l)) + { + *f = *l; + f->class = sim_fpu_class_qnan; + return 0; + } + if (sim_fpu_is_qnan (r)) + { + *f = *r; + f->class = sim_fpu_class_qnan; + return 0; + } + if (sim_fpu_is_infinity (l)) + { + if (sim_fpu_is_infinity (r)) + { + *f = sim_fpu_qnan; + return sim_fpu_status_invalid_idi; + } + else + { + *f = *l; + f->sign = l->sign ^ r->sign; + return 0; + } + } + if (sim_fpu_is_zero (l)) + { + if (sim_fpu_is_zero (r)) + { + *f = sim_fpu_qnan; + return sim_fpu_status_invalid_zdz; + } + else + { + *f = *l; + f->sign = l->sign ^ r->sign; + return 0; + } + } + if (sim_fpu_is_infinity (r)) + { + *f = sim_fpu_zero; + f->sign = l->sign ^ r->sign; + return 0; + } + if (sim_fpu_is_zero (r)) + { + f->class = sim_fpu_class_infinity; + f->sign = l->sign ^ r->sign; + return sim_fpu_status_invalid_div0; + } + + /* Calculate the mantissa by multiplying both 64bit numbers to get a + 128 bit number */ + { + /* quotient = ( ( numerator / denominator) + x 2^(numerator exponent - denominator exponent) + */ + unsigned64 numerator; + unsigned64 denominator; + unsigned64 quotient; + unsigned64 bit; + + f->class = sim_fpu_class_number; + f->sign = l->sign ^ r->sign; + f->normal_exp = l->normal_exp - r->normal_exp; + + numerator = l->fraction; + denominator = r->fraction; + + /* Fraction will be less than 1.0 */ + if (numerator < denominator) + { + numerator <<= 1; + f->normal_exp--; + } + ASSERT (numerator >= denominator); + + /* Gain extra precision, already used one spare bit */ + numerator <<= NR_SPARE; + denominator <<= NR_SPARE; + + /* Does divide one bit at a time. Optimize??? */ + quotient = 0; + bit = (IMPLICIT_1 << NR_SPARE); + while (bit) + { + if (numerator >= denominator) + { + quotient |= bit; + numerator -= denominator; + } + bit >>= 1; + numerator <<= 1; + } + +#if 0 + printf ("\n"); + print_bits (quotient, 63, (sim_fpu_print_func*)fprintf, stdout); + printf ("\n"); + print_bits (numerator, 63, (sim_fpu_print_func*)fprintf, stdout); + printf ("\n"); + print_bits (denominator, 63, (sim_fpu_print_func*)fprintf, stdout); + printf ("\n"); +#endif + + /* discard (but save) the extra bits */ + if ((quotient & LSMASK64 (NR_SPARE -1, 0))) + quotient = (quotient >> NR_SPARE) | 1; + else + quotient = (quotient >> NR_SPARE); + + f->fraction = quotient; + ASSERT (f->fraction >= IMPLICIT_1 && f->fraction < IMPLICIT_2); + if (numerator != 0) + { + f->fraction |= 1; /* stick remaining bits */ + return sim_fpu_status_inexact; + } + else + return 0; + } +} + + +INLINE_SIM_FPU (int) +sim_fpu_max (sim_fpu *f, + const sim_fpu *l, + const sim_fpu *r) +{ + if (sim_fpu_is_snan (l)) + { + *f = *l; + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_snan (r)) + { + *f = *r; + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_qnan (l)) + { + *f = *l; + return 0; + } + if (sim_fpu_is_qnan (r)) + { + *f = *r; + return 0; + } + if (sim_fpu_is_infinity (l)) + { + if (sim_fpu_is_infinity (r) + && l->sign == r->sign) + { + *f = sim_fpu_qnan; + return sim_fpu_status_invalid_isi; + } + if (l->sign) + *f = *r; /* -inf < anything */ + else + *f = *l; /* +inf > anthing */ + return 0; + } + if (sim_fpu_is_infinity (r)) + { + if (r->sign) + *f = *l; /* anything > -inf */ + else + *f = *r; /* anthing < +inf */ + return 0; + } + if (l->sign > r->sign) + { + *f = *r; /* -ve < +ve */ + return 0; + } + if (l->sign < r->sign) + { + *f = *l; /* +ve > -ve */ + return 0; + } + ASSERT (l->sign == r->sign); + if (l->normal_exp > r->normal_exp + || (l->normal_exp == r->normal_exp && + l->fraction > r->fraction)) + { + /* |l| > |r| */ + if (l->sign) + *f = *r; /* -ve < -ve */ + else + *f = *l; /* +ve > +ve */ + return 0; + } + else + { + /* |l| <= |r| */ + if (l->sign) + *f = *l; /* -ve > -ve */ + else + *f = *r; /* +ve < +ve */ + return 0; + } +} + + +INLINE_SIM_FPU (int) +sim_fpu_min (sim_fpu *f, + const sim_fpu *l, + const sim_fpu *r) +{ + if (sim_fpu_is_snan (l)) + { + *f = *l; + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_snan (r)) + { + *f = *r; + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_qnan (l)) + { + *f = *l; + return 0; + } + if (sim_fpu_is_qnan (r)) + { + *f = *r; + return 0; + } + if (sim_fpu_is_infinity (l)) + { + if (sim_fpu_is_infinity (r) + && l->sign == r->sign) + { + *f = sim_fpu_qnan; + return sim_fpu_status_invalid_isi; + } + if (l->sign) + *f = *l; /* -inf < anything */ + else + *f = *r; /* +inf > anthing */ + return 0; + } + if (sim_fpu_is_infinity (r)) + { + if (r->sign) + *f = *r; /* anything > -inf */ + else + *f = *l; /* anything < +inf */ + return 0; + } + if (l->sign > r->sign) + { + *f = *l; /* -ve < +ve */ + return 0; + } + if (l->sign < r->sign) + { + *f = *r; /* +ve > -ve */ + return 0; + } + ASSERT (l->sign == r->sign); + if (l->normal_exp > r->normal_exp + || (l->normal_exp == r->normal_exp && + l->fraction > r->fraction)) + { + /* |l| > |r| */ + if (l->sign) + *f = *l; /* -ve < -ve */ + else + *f = *r; /* +ve > +ve */ + return 0; + } + else + { + /* |l| <= |r| */ + if (l->sign) + *f = *r; /* -ve > -ve */ + else + *f = *l; /* +ve < +ve */ + return 0; + } +} + + +INLINE_SIM_FPU (int) +sim_fpu_neg (sim_fpu *f, + const sim_fpu *r) +{ + if (sim_fpu_is_snan (r)) + { + *f = *r; + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_qnan (r)) + { + *f = *r; + return 0; + } + *f = *r; + f->sign = !r->sign; + return 0; +} + + +INLINE_SIM_FPU (int) +sim_fpu_abs (sim_fpu *f, + const sim_fpu *r) +{ + if (sim_fpu_is_snan (r)) + { + *f = *r; + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_qnan (r)) + { + *f = *r; + return 0; + } + *f = *r; + f->sign = 0; + return 0; +} + + +INLINE_SIM_FPU (int) +sim_fpu_inv (sim_fpu *f, + const sim_fpu *r) +{ + if (sim_fpu_is_snan (r)) + { + *f = *r; + f->class = sim_fpu_class_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_qnan (r)) + { + *f = *r; + f->class = sim_fpu_class_qnan; + return 0; + } + if (sim_fpu_is_infinity (r)) + { + *f = sim_fpu_zero; + f->sign = r->sign; + return 0; + } + if (sim_fpu_is_zero (r)) + { + f->class = sim_fpu_class_infinity; + f->sign = r->sign; + return sim_fpu_status_invalid_div0; + } + *f = *r; + f->normal_exp = - r->normal_exp; + return 0; +} + + +INLINE_SIM_FPU (int) +sim_fpu_sqrt (sim_fpu *f, + const sim_fpu *r) +{ + if (sim_fpu_is_snan (r)) + { + *f = sim_fpu_qnan; + return sim_fpu_status_invalid_snan; + } + if (sim_fpu_is_qnan (r)) + { + *f = sim_fpu_qnan; + return 0; + } + if (sim_fpu_is_zero (r)) + { + f->class = sim_fpu_class_zero; + f->sign = r->sign; + return 0; + } + if (sim_fpu_is_infinity (r)) + { + if (r->sign) + { + *f = sim_fpu_qnan; + return sim_fpu_status_invalid_sqrt; + } + else + { + f->class = sim_fpu_class_infinity; + f->sign = 0; + f->sign = 0; + return 0; + } + } + if (r->sign) + { + *f = sim_fpu_qnan; + return sim_fpu_status_invalid_sqrt; + } + + /* @(#)e_sqrt.c 5.1 93/09/24 */ + /* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + + /* __ieee754_sqrt(x) + * Return correctly rounded sqrt. + * ------------------------------------------ + * | Use the hardware sqrt if you have one | + * ------------------------------------------ + * Method: + * Bit by bit method using integer arithmetic. (Slow, but portable) + * 1. Normalization + * Scale x to y in [1,4) with even powers of 2: + * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then + * sqrt(x) = 2^k * sqrt(y) + - + - Since: + - sqrt ( x*2^(2m) ) = sqrt(x).2^m ; m even + - sqrt ( x*2^(2m + 1) ) = sqrt(2.x).2^m ; m odd + - Define: + - y = ((m even) ? x : 2.x) + - Then: + - y in [1, 4) ; [IMPLICIT_1,IMPLICIT_4) + - And: + - sqrt (y) in [1, 2) ; [IMPLICIT_1,IMPLICIT_2) + - + * 2. Bit by bit computation + * Let q = sqrt(y) truncated to i bit after binary point (q = 1), + * i 0 + * i+1 2 + * s = 2*q , and y = 2 * ( y - q ). (1) + * i i i i + * + * To compute q from q , one checks whether + * i+1 i + * + * -(i+1) 2 + * (q + 2 ) <= y. (2) + * i + * -(i+1) + * If (2) is false, then q = q ; otherwise q = q + 2 . + * i+1 i i+1 i + * + * With some algebric manipulation, it is not difficult to see + * that (2) is equivalent to + * -(i+1) + * s + 2 <= y (3) + * i i + * + * The advantage of (3) is that s and y can be computed by + * i i + * the following recurrence formula: + * if (3) is false + * + * s = s , y = y ; (4) + * i+1 i i+1 i + * + - + - NOTE: y = 2*y + - i+1 i + - + * otherwise, + * -i -(i+1) + * s = s + 2 , y = y - s - 2 (5) + * i+1 i i+1 i i + * + - + - -(i+1) + - NOTE: y = 2 (y - s - 2 ) + - i+1 i i + - + * One may easily use induction to prove (4) and (5). + * Note. Since the left hand side of (3) contain only i+2 bits, + * it does not necessary to do a full (53-bit) comparison + * in (3). + * 3. Final rounding + * After generating the 53 bits result, we compute one more bit. + * Together with the remainder, we can decide whether the + * result is exact, bigger than 1/2ulp, or less than 1/2ulp + * (it will never equal to 1/2ulp). + * The rounding mode can be detected by checking whether + * huge + tiny is equal to huge, and whether huge - tiny is + * equal to huge for some floating point number "huge" and "tiny". + * + * Special cases: + * sqrt(+-0) = +-0 ... exact + * sqrt(inf) = inf + * sqrt(-ve) = NaN ... with invalid signal + * sqrt(NaN) = NaN ... with invalid signal for signaling NaN + * + * Other methods : see the appended file at the end of the program below. + *--------------- + */ + + { + /* generate sqrt(x) bit by bit */ + unsigned64 y; + unsigned64 q; + unsigned64 s; + unsigned64 b; + + f->class = sim_fpu_class_number; + f->sign = 0; + y = r->fraction; + f->normal_exp = (r->normal_exp >> 1); /* exp = [exp/2] */ + + /* odd exp, double x to make it even */ + ASSERT (y >= IMPLICIT_1 && y < IMPLICIT_4); + if ((r->normal_exp & 1)) + { + y += y; + } + ASSERT (y >= IMPLICIT_1 && y < (IMPLICIT_2 << 1)); + + /* Let loop determine first value of s (either 1 or 2) */ + b = IMPLICIT_1; + q = 0; + s = 0; + + while (b) + { + unsigned64 t = s + b; + if (t <= y) + { + s |= (b << 1); + y -= t; + q |= b; + } + y <<= 1; + b >>= 1; + } + + ASSERT (q >= IMPLICIT_1 && q < IMPLICIT_2); + f->fraction = q; + if (y != 0) + { + f->fraction |= 1; /* stick remaining bits */ + return sim_fpu_status_inexact; + } + else + return 0; + } +} + + +/* int/long <-> sim_fpu */ + +INLINE_SIM_FPU (int) +sim_fpu_i32to (sim_fpu *f, + signed32 i, + sim_fpu_round round) +{ + i2fpu (f, i, 0); + return 0; +} + +INLINE_SIM_FPU (int) +sim_fpu_u32to (sim_fpu *f, + unsigned32 u, + sim_fpu_round round) +{ + u2fpu (f, u, 0); + return 0; +} + +INLINE_SIM_FPU (int) +sim_fpu_i64to (sim_fpu *f, + signed64 i, + sim_fpu_round round) +{ + i2fpu (f, i, 1); + return 0; +} + +INLINE_SIM_FPU (int) +sim_fpu_u64to (sim_fpu *f, + unsigned64 u, + sim_fpu_round round) +{ + u2fpu (f, u, 1); + return 0; +} + + +INLINE_SIM_FPU (int) +sim_fpu_to32i (signed32 *i, + const sim_fpu *f, + sim_fpu_round round) +{ + signed64 i64; + int status = fpu2i (&i64, f, 0, round); + *i = i64; + return status; +} + +INLINE_SIM_FPU (int) +sim_fpu_to32u (unsigned32 *u, + const sim_fpu *f, + sim_fpu_round round) +{ + unsigned64 u64; + int status = fpu2u (&u64, f, 0); + *u = u64; + return status; +} + +INLINE_SIM_FPU (int) +sim_fpu_to64i (signed64 *i, + const sim_fpu *f, + sim_fpu_round round) +{ + return fpu2i (i, f, 1, round); +} + + +INLINE_SIM_FPU (int) +sim_fpu_to64u (unsigned64 *u, + const sim_fpu *f, + sim_fpu_round round) +{ + return fpu2u (u, f, 1); +} + + + +/* sim_fpu -> host format */ + +#if 0 +INLINE_SIM_FPU (float) +sim_fpu_2f (const sim_fpu *f) +{ + return fval.d; +} +#endif + + +INLINE_SIM_FPU (double) +sim_fpu_2d (const sim_fpu *s) +{ + sim_fpu_map val; + if (sim_fpu_is_snan (s)) + { + /* gag SNaN's */ + sim_fpu n = *s; + n.class = sim_fpu_class_qnan; + val.i = pack_fpu (&n, 1); + } + else + { + val.i = pack_fpu (s, 1); + } + return val.d; +} + + +#if 0 +INLINE_SIM_FPU (void) +sim_fpu_f2 (sim_fpu *f, + float s) +{ + sim_fpu_map val; + val.d = s; + unpack_fpu (f, val.i, 1); +} +#endif + + +INLINE_SIM_FPU (void) +sim_fpu_d2 (sim_fpu *f, + double d) +{ + sim_fpu_map val; + val.d = d; + unpack_fpu (f, val.i, 1); +} + + +/* General */ + +INLINE_SIM_FPU (int) +sim_fpu_is_nan (const sim_fpu *d) +{ + switch (d->class) + { + case sim_fpu_class_qnan: + case sim_fpu_class_snan: + return 1; + default: + return 0; + } +} + +INLINE_SIM_FPU (int) +sim_fpu_is_qnan (const sim_fpu *d) +{ + switch (d->class) + { + case sim_fpu_class_qnan: + return 1; + default: + return 0; + } +} + +INLINE_SIM_FPU (int) +sim_fpu_is_snan (const sim_fpu *d) +{ + switch (d->class) + { + case sim_fpu_class_snan: + return 1; + default: + return 0; + } +} + +INLINE_SIM_FPU (int) +sim_fpu_is_zero (const sim_fpu *d) +{ + switch (d->class) + { + case sim_fpu_class_zero: + return 1; + default: + return 0; + } +} + +INLINE_SIM_FPU (int) +sim_fpu_is_infinity (const sim_fpu *d) +{ + switch (d->class) + { + case sim_fpu_class_infinity: + return 1; + default: + return 0; + } +} + +INLINE_SIM_FPU (int) +sim_fpu_is_number (const sim_fpu *d) +{ + switch (d->class) + { + case sim_fpu_class_denorm: + case sim_fpu_class_number: + return 1; + default: + return 0; + } +} + +INLINE_SIM_FPU (int) +sim_fpu_is_denorm (const sim_fpu *d) +{ + switch (d->class) + { + case sim_fpu_class_denorm: + return 1; + default: + return 0; + } +} + + +INLINE_SIM_FPU (int) +sim_fpu_sign (const sim_fpu *d) +{ + return d->sign; +} + + +INLINE_SIM_FPU (int) +sim_fpu_exp (const sim_fpu *d) +{ + return d->normal_exp; +} + + + +INLINE_SIM_FPU (int) +sim_fpu_is (const sim_fpu *d) +{ + switch (d->class) + { + case sim_fpu_class_qnan: + return SIM_FPU_IS_QNAN; + case sim_fpu_class_snan: + return SIM_FPU_IS_SNAN; + case sim_fpu_class_infinity: + if (d->sign) + return SIM_FPU_IS_NINF; + else + return SIM_FPU_IS_PINF; + case sim_fpu_class_number: + if (d->sign) + return SIM_FPU_IS_NNUMBER; + else + return SIM_FPU_IS_PNUMBER; + case sim_fpu_class_denorm: + if (d->sign) + return SIM_FPU_IS_NDENORM; + else + return SIM_FPU_IS_PDENORM; + case sim_fpu_class_zero: + if (d->sign) + return SIM_FPU_IS_NZERO; + else + return SIM_FPU_IS_PZERO; + default: + return -1; + abort (); + } +} + +INLINE_SIM_FPU (int) +sim_fpu_cmp (const sim_fpu *l, const sim_fpu *r) +{ + sim_fpu res; + sim_fpu_sub (&res, l, r); + return sim_fpu_is (&res); +} + +INLINE_SIM_FPU (int) +sim_fpu_is_lt (const sim_fpu *l, const sim_fpu *r) +{ + int status; + sim_fpu_lt (&status, l, r); + return status; +} + +INLINE_SIM_FPU (int) +sim_fpu_is_le (const sim_fpu *l, const sim_fpu *r) +{ + int is; + sim_fpu_le (&is, l, r); + return is; +} + +INLINE_SIM_FPU (int) +sim_fpu_is_eq (const sim_fpu *l, const sim_fpu *r) +{ + int is; + sim_fpu_eq (&is, l, r); + return is; +} + +INLINE_SIM_FPU (int) +sim_fpu_is_ne (const sim_fpu *l, const sim_fpu *r) +{ + int is; + sim_fpu_ne (&is, l, r); + return is; +} + +INLINE_SIM_FPU (int) +sim_fpu_is_ge (const sim_fpu *l, const sim_fpu *r) +{ + int is; + sim_fpu_ge (&is, l, r); + return is; +} + +INLINE_SIM_FPU (int) +sim_fpu_is_gt (const sim_fpu *l, const sim_fpu *r) +{ + int is; + sim_fpu_gt (&is, l, r); + return is; +} + + +/* Compare operators */ + +INLINE_SIM_FPU (int) +sim_fpu_lt (int *is, + const sim_fpu *l, + const sim_fpu *r) +{ + if (!sim_fpu_is_nan (l) && !sim_fpu_is_nan (r)) + { + sim_fpu_map lval; + sim_fpu_map rval; + lval.i = pack_fpu (l, 1); + rval.i = pack_fpu (r, 1); + (*is) = (lval.d < rval.d); + return 0; + } + else if (sim_fpu_is_snan (l) || sim_fpu_is_snan (r)) + { + *is = 0; + return sim_fpu_status_invalid_snan; + } + else + { + *is = 0; + return sim_fpu_status_invalid_qnan; + } +} + +INLINE_SIM_FPU (int) +sim_fpu_le (int *is, + const sim_fpu *l, + const sim_fpu *r) +{ + if (!sim_fpu_is_nan (l) && !sim_fpu_is_nan (r)) + { + sim_fpu_map lval; + sim_fpu_map rval; + lval.i = pack_fpu (l, 1); + rval.i = pack_fpu (r, 1); + *is = (lval.d <= rval.d); + return 0; + } + else if (sim_fpu_is_snan (l) || sim_fpu_is_snan (r)) + { + *is = 0; + return sim_fpu_status_invalid_snan; + } + else + { + *is = 0; + return sim_fpu_status_invalid_qnan; + } +} + +INLINE_SIM_FPU (int) +sim_fpu_eq (int *is, + const sim_fpu *l, + const sim_fpu *r) +{ + if (!sim_fpu_is_nan (l) && !sim_fpu_is_nan (r)) + { + sim_fpu_map lval; + sim_fpu_map rval; + lval.i = pack_fpu (l, 1); + rval.i = pack_fpu (r, 1); + (*is) = (lval.d == rval.d); + return 0; + } + else if (sim_fpu_is_snan (l) || sim_fpu_is_snan (r)) + { + *is = 0; + return sim_fpu_status_invalid_snan; + } + else + { + *is = 0; + return sim_fpu_status_invalid_qnan; + } +} + +INLINE_SIM_FPU (int) +sim_fpu_ne (int *is, + const sim_fpu *l, + const sim_fpu *r) +{ + if (!sim_fpu_is_nan (l) && !sim_fpu_is_nan (r)) + { + sim_fpu_map lval; + sim_fpu_map rval; + lval.i = pack_fpu (l, 1); + rval.i = pack_fpu (r, 1); + (*is) = (lval.d != rval.d); + return 0; + } + else if (sim_fpu_is_snan (l) || sim_fpu_is_snan (r)) + { + *is = 0; + return sim_fpu_status_invalid_snan; + } + else + { + *is = 0; + return sim_fpu_status_invalid_qnan; + } +} + +INLINE_SIM_FPU (int) +sim_fpu_ge (int *is, + const sim_fpu *l, + const sim_fpu *r) +{ + return sim_fpu_le (is, r, l); +} + +INLINE_SIM_FPU (int) +sim_fpu_gt (int *is, + const sim_fpu *l, + const sim_fpu *r) +{ + return sim_fpu_lt (is, r, l); +} + + +/* A number of useful constants */ + +#if EXTERN_SIM_FPU_P +const sim_fpu sim_fpu_zero = { + sim_fpu_class_zero, +}; +const sim_fpu sim_fpu_qnan = { + sim_fpu_class_qnan, +}; +const sim_fpu sim_fpu_one = { + sim_fpu_class_number, 0, IMPLICIT_1, 1 +}; +const sim_fpu sim_fpu_two = { + sim_fpu_class_number, 0, IMPLICIT_1, 2 +}; +const sim_fpu sim_fpu_max32 = { + sim_fpu_class_number, 0, LSMASK64 (NR_FRAC_GUARD, NR_GUARDS32), NORMAL_EXPMAX32 +}; +const sim_fpu sim_fpu_max64 = { + sim_fpu_class_number, 0, LSMASK64 (NR_FRAC_GUARD, NR_GUARDS64), NORMAL_EXPMAX64 +}; +#endif + + +/* For debugging */ + +INLINE_SIM_FPU (void) +sim_fpu_print_fpu (const sim_fpu *f, + sim_fpu_print_func *print, + void *arg) +{ + print (arg, "%s", f->sign ? "-" : "+"); + switch (f->class) + { + case sim_fpu_class_qnan: + print (arg, "0."); + print_bits (f->fraction, NR_FRAC_GUARD - 1, print, arg); + print (arg, "*QuietNaN"); + break; + case sim_fpu_class_snan: + print (arg, "0."); + print_bits (f->fraction, NR_FRAC_GUARD - 1, print, arg); + print (arg, "*SignalNaN"); + break; + case sim_fpu_class_zero: + print (arg, "0.0"); + break; + case sim_fpu_class_infinity: + print (arg, "INF"); + break; + case sim_fpu_class_number: + case sim_fpu_class_denorm: + print (arg, "1."); + print_bits (f->fraction, NR_FRAC_GUARD - 1, print, arg); + print (arg, "*2^%+-5d", f->normal_exp); + ASSERT (f->fraction >= IMPLICIT_1); + ASSERT (f->fraction < IMPLICIT_2); + } +} + + +INLINE_SIM_FPU (void) +sim_fpu_print_status (int status, + sim_fpu_print_func *print, + void *arg) +{ + int i = 1; + char *prefix = ""; + while (status >= i) + { + switch ((sim_fpu_status) (status & i)) + { + case sim_fpu_status_denorm: + print (arg, "%sD", prefix); + break; + case sim_fpu_status_invalid_snan: + print (arg, "%sSNaN", prefix); + break; + case sim_fpu_status_invalid_qnan: + print (arg, "%sQNaN", prefix); + break; + case sim_fpu_status_invalid_isi: + print (arg, "%sISI", prefix); + break; + case sim_fpu_status_invalid_idi: + print (arg, "%sIDI", prefix); + break; + case sim_fpu_status_invalid_zdz: + print (arg, "%sZDZ", prefix); + break; + case sim_fpu_status_invalid_imz: + print (arg, "%sIMZ", prefix); + break; + case sim_fpu_status_invalid_cvi: + print (arg, "%sCVI", prefix); + break; + case sim_fpu_status_invalid_cmp: + print (arg, "%sCMP", prefix); + break; + case sim_fpu_status_invalid_sqrt: + print (arg, "%sSQRT", prefix); + break; + break; + case sim_fpu_status_inexact: + print (arg, "%sX", prefix); + break; + break; + case sim_fpu_status_overflow: + print (arg, "%sO", prefix); + break; + break; + case sim_fpu_status_underflow: + print (arg, "%sU", prefix); + break; + break; + case sim_fpu_status_invalid_div0: + print (arg, "%s/", prefix); + break; + break; + case sim_fpu_status_rounded: + print (arg, "%sR", prefix); + break; + break; + } + i <<= 1; + prefix = ","; + } +} + +#endif |