1523 lines
31 KiB
C
1523 lines
31 KiB
C
/*-
|
|
* Copyright (c) 1998 Doug Rabson
|
|
* All rights reserved.
|
|
*
|
|
* Redistribution and use in source and binary forms, with or without
|
|
* modification, are permitted provided that the following conditions
|
|
* are met:
|
|
* 1. Redistributions of source code must retain the above copyright
|
|
* notice, this list of conditions and the following disclaimer.
|
|
* 2. Redistributions in binary form must reproduce the above copyright
|
|
* notice, this list of conditions and the following disclaimer in the
|
|
* documentation and/or other materials provided with the distribution.
|
|
*
|
|
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
|
|
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
|
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
|
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
|
|
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
|
|
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
|
|
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
|
|
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
|
|
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
|
|
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
|
|
* SUCH DAMAGE.
|
|
*
|
|
*/
|
|
|
|
/*
|
|
* An implementation of IEEE 754 floating point arithmetic supporting
|
|
* multiply, divide, addition, subtraction and conversion to and from
|
|
* integer. Probably not the fastest floating point code in the world
|
|
* but it should be pretty accurate.
|
|
*
|
|
* A special thanks to John Polstra for pointing out some problems
|
|
* with an earlier version of this code and for educating me as to the
|
|
* correct use of sticky bits.
|
|
*/
|
|
|
|
#include <sys/cdefs.h>
|
|
__FBSDID("$FreeBSD$");
|
|
|
|
#include <sys/types.h>
|
|
#ifdef TEST
|
|
#include "../include/fpu.h"
|
|
#include "ieee_float.h"
|
|
#else
|
|
#include <sys/param.h>
|
|
#include <sys/systm.h>
|
|
#include <sys/sysproto.h>
|
|
#include <sys/sysent.h>
|
|
#include <sys/proc.h>
|
|
#include <machine/fpu.h>
|
|
#include <alpha/alpha/ieee_float.h>
|
|
#endif
|
|
|
|
/*
|
|
* The number of fraction bits in a T format float.
|
|
*/
|
|
#define T_FRACBITS 52
|
|
|
|
/*
|
|
* The number of fraction bits in a S format float.
|
|
*/
|
|
#define S_FRACBITS 23
|
|
|
|
/*
|
|
* Mask the fraction part of a float to contain only those bits which
|
|
* should be in single precision number.
|
|
*/
|
|
#define S_FRACMASK ((1ULL << 52) - (1ULL << 29))
|
|
|
|
/*
|
|
* The number of extra zero bits we shift into the fraction part
|
|
* to gain accuracy. Two guard bits and one sticky bit are required
|
|
* to ensure accurate rounding.
|
|
*/
|
|
#define FRAC_SHIFT 3
|
|
|
|
/*
|
|
* Values for 1.0 and 2.0 fractions (including the extra FRAC_SHIFT
|
|
* bits).
|
|
*/
|
|
#define ONE (1ULL << (T_FRACBITS + FRAC_SHIFT))
|
|
#define TWO (ONE + ONE)
|
|
|
|
/*
|
|
* The maximum and minimum values for S and T format exponents.
|
|
*/
|
|
#define T_MAXEXP 0x3ff
|
|
#define T_MINEXP -0x3fe
|
|
#define S_MAXEXP 0x7f
|
|
#define S_MINEXP -0x7e
|
|
|
|
/*
|
|
* Exponent values in registers are biased by adding this value.
|
|
*/
|
|
#define BIAS_EXP 0x3ff
|
|
|
|
/*
|
|
* Exponent value for INF and NaN.
|
|
*/
|
|
#define NAN_EXP 0x7ff
|
|
|
|
/*
|
|
* If this bit is set in the fraction part of a NaN, then the number
|
|
* is a quiet NaN, i.e. no traps are generated.
|
|
*/
|
|
#define QNAN_BIT (1ULL << 51)
|
|
|
|
/*
|
|
* Return true if the number is any kind of NaN.
|
|
*/
|
|
static __inline int
|
|
isNaN(fp_register_t f)
|
|
{
|
|
return f.t.exponent == NAN_EXP && f.t.fraction != 0;
|
|
}
|
|
|
|
/*
|
|
* Return true if the number is a quiet NaN.
|
|
*/
|
|
static __inline int
|
|
isQNaN(fp_register_t f)
|
|
{
|
|
return f.t.exponent == NAN_EXP && (f.t.fraction & QNAN_BIT);
|
|
}
|
|
|
|
/*
|
|
* Return true if the number is a signalling NaN.
|
|
*/
|
|
static __inline int
|
|
isSNaN(fp_register_t f)
|
|
{
|
|
return isNaN(f) && !isQNaN(f);
|
|
}
|
|
|
|
/*
|
|
* Return true if the number is +/- INF.
|
|
*/
|
|
static __inline int
|
|
isINF(fp_register_t f)
|
|
{
|
|
return f.t.exponent == NAN_EXP && f.t.fraction == 0;
|
|
}
|
|
|
|
/*
|
|
* Return true if the number is +/- 0.
|
|
*/
|
|
static __inline int
|
|
isZERO(fp_register_t f)
|
|
{
|
|
return f.t.exponent == 0 && f.t.fraction == 0;
|
|
}
|
|
|
|
/*
|
|
* Return true if the number is denormalised.
|
|
*/
|
|
static __inline int
|
|
isDENORM(fp_register_t f)
|
|
{
|
|
return f.t.exponent == 0 && f.t.fraction != 0;
|
|
}
|
|
|
|
/*
|
|
* Extract the exponent part of a float register. If the exponent is
|
|
* zero, the number may be denormalised (if the fraction is nonzero).
|
|
* If so, return the minimum exponent for the source datatype.
|
|
*/
|
|
static __inline int
|
|
getexp(fp_register_t f, int src)
|
|
{
|
|
int minexp[] = { S_MINEXP, 0, T_MINEXP, 0 };
|
|
if (f.t.exponent == 0) {
|
|
if (f.t.fraction)
|
|
return minexp[src];
|
|
else
|
|
return 0;
|
|
}
|
|
return f.t.exponent - BIAS_EXP;
|
|
}
|
|
|
|
/*
|
|
* Extract the fraction part of a float register, shift it up a bit
|
|
* to give extra accuracy and add in the implicit 1 bit. Must be
|
|
* careful to handle denormalised numbers and zero correctly.
|
|
*/
|
|
static __inline u_int64_t
|
|
getfrac(fp_register_t f)
|
|
{
|
|
if (f.t.exponent == 0)
|
|
return f.t.fraction << FRAC_SHIFT;
|
|
else
|
|
return (f.t.fraction << FRAC_SHIFT) | ONE;
|
|
}
|
|
|
|
/*
|
|
* Make a float (in register format) from a sign, exponent and
|
|
* fraction, normalising and rounding as necessary.
|
|
* Return the float and set *status if any traps are generated.
|
|
*/
|
|
static fp_register_t
|
|
makefloat(int sign, int exp, u_int64_t frac,
|
|
int src, int rnd,
|
|
u_int64_t control, u_int64_t *status)
|
|
{
|
|
fp_register_t f;
|
|
int minexp = 0, maxexp = 0, alpha = 0;
|
|
u_int64_t epsilon = 0, max = 0;
|
|
|
|
if (frac == 0) {
|
|
f.t.sign = sign;
|
|
f.t.exponent = 0;
|
|
f.t.fraction = 0;
|
|
return f;
|
|
}
|
|
|
|
if (frac >= TWO) {
|
|
/*
|
|
* Fraction is >= 2.0.
|
|
* Shift the fraction down, preserving the 'sticky'
|
|
* bit.
|
|
*/
|
|
while (frac >= TWO) {
|
|
frac = (frac >> 1) | (frac & 1);
|
|
exp++;
|
|
}
|
|
} else if (frac < ONE) {
|
|
/*
|
|
* Fraction is < 1.0. Shift it up.
|
|
*/
|
|
while (frac < ONE) {
|
|
frac = (frac << 1) | (frac & 1);
|
|
exp--;
|
|
}
|
|
}
|
|
|
|
switch (src) {
|
|
case S_FORMAT:
|
|
minexp = S_MINEXP;
|
|
maxexp = S_MAXEXP;
|
|
alpha = 0xc0;
|
|
epsilon = (1ULL << (T_FRACBITS - S_FRACBITS + FRAC_SHIFT));
|
|
max = TWO - epsilon;
|
|
break;
|
|
|
|
case T_FORMAT:
|
|
minexp = T_MINEXP;
|
|
maxexp = T_MAXEXP;
|
|
alpha = 0x600;
|
|
epsilon = (1ULL << FRAC_SHIFT);
|
|
max = TWO - epsilon;
|
|
break;
|
|
}
|
|
|
|
/*
|
|
* Handle underflow before rounding so that denormalised
|
|
* numbers are rounded correctly.
|
|
*/
|
|
if (exp < minexp) {
|
|
*status |= FPCR_INE;
|
|
if (control & IEEE_TRAP_ENABLE_UNF) {
|
|
*status |= FPCR_UNF;
|
|
exp += alpha;
|
|
} else {
|
|
/* denormalise */
|
|
while (exp < minexp) {
|
|
exp++;
|
|
frac = (frac >> 1) | (frac & 1);
|
|
}
|
|
exp = minexp - 1;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Round the fraction according to the rounding mode.
|
|
*/
|
|
if (frac & (epsilon - 1)) {
|
|
u_int64_t fraclo, frachi;
|
|
u_int64_t difflo, diffhi;
|
|
|
|
fraclo = frac & max;
|
|
frachi = fraclo + epsilon;
|
|
switch (rnd) {
|
|
case ROUND_CHOP:
|
|
frac = fraclo;
|
|
break;
|
|
case ROUND_MINUS_INF:
|
|
if (f.t.sign)
|
|
frac = frachi;
|
|
else
|
|
frac = fraclo;
|
|
break;
|
|
case ROUND_NORMAL:
|
|
difflo = frac - fraclo;
|
|
diffhi = frachi - frac;
|
|
if (difflo < diffhi)
|
|
frac = fraclo;
|
|
else if (diffhi < difflo)
|
|
frac = frachi;
|
|
else
|
|
/* round to even */
|
|
if (fraclo & epsilon)
|
|
frac = frachi;
|
|
else
|
|
frac = fraclo;
|
|
break;
|
|
case ROUND_PLUS_INF:
|
|
if (f.t.sign)
|
|
frac = fraclo;
|
|
else
|
|
frac = frachi;
|
|
break;
|
|
}
|
|
|
|
/*
|
|
* Rounding up may take us to TWO if
|
|
* fraclo == (TWO - epsilon). Also If fraclo has been
|
|
* denormalised to (ONE - epsilon) then there is a
|
|
* possibility that we will round to ONE exactly.
|
|
*/
|
|
if (frac >= TWO) {
|
|
frac = (frac >> 1) & ~(epsilon - 1);
|
|
exp++;
|
|
} else if (exp == minexp - 1 && frac == ONE) {
|
|
/* Renormalise to ONE * 2^minexp */
|
|
exp = minexp;
|
|
}
|
|
|
|
*status |= FPCR_INE;
|
|
}
|
|
|
|
/*
|
|
* Check for overflow and round to the correct INF as needed.
|
|
*/
|
|
if (exp > maxexp) {
|
|
*status |= FPCR_OVF | FPCR_INE;
|
|
if (control & IEEE_TRAP_ENABLE_OVF) {
|
|
exp -= alpha;
|
|
} else {
|
|
switch (rnd) {
|
|
case ROUND_CHOP:
|
|
exp = maxexp;
|
|
frac = max;
|
|
break;
|
|
case ROUND_MINUS_INF:
|
|
if (sign) {
|
|
exp = maxexp + 1; /* INF */
|
|
frac = 0;
|
|
} else {
|
|
exp = maxexp;
|
|
frac = max;
|
|
}
|
|
break;
|
|
case ROUND_NORMAL:
|
|
exp = maxexp + 1; /* INF */
|
|
frac = 0;
|
|
break;
|
|
case ROUND_PLUS_INF:
|
|
if (sign) {
|
|
exp = maxexp;
|
|
frac = max;
|
|
} else {
|
|
exp = maxexp + 1; /* INF */
|
|
frac = 0;
|
|
}
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
f.t.sign = sign;
|
|
if (exp > maxexp) /* NaN, INF */
|
|
f.t.exponent = NAN_EXP;
|
|
else if (exp < minexp) /* denorm, zero */
|
|
f.t.exponent = 0;
|
|
else
|
|
f.t.exponent = exp + BIAS_EXP;
|
|
f.t.fraction = (frac & ~ONE) >> FRAC_SHIFT;
|
|
return f;
|
|
}
|
|
|
|
/*
|
|
* Return the canonical quiet NaN in register format.
|
|
*/
|
|
static fp_register_t
|
|
makeQNaN(void)
|
|
{
|
|
fp_register_t f;
|
|
f.t.sign = 0;
|
|
f.t.exponent = NAN_EXP;
|
|
f.t.fraction = QNAN_BIT;
|
|
return f;
|
|
}
|
|
|
|
/*
|
|
* Return +/- INF.
|
|
*/
|
|
static fp_register_t
|
|
makeINF(int sign)
|
|
{
|
|
fp_register_t f;
|
|
f.t.sign = sign;
|
|
f.t.exponent = NAN_EXP;
|
|
f.t.fraction = 0;
|
|
return f;
|
|
}
|
|
|
|
/*
|
|
* Return +/- 0.
|
|
*/
|
|
static fp_register_t
|
|
makeZERO(int sign)
|
|
{
|
|
fp_register_t f;
|
|
f.t.sign = sign;
|
|
f.t.exponent = 0;
|
|
f.t.fraction = 0;
|
|
return f;
|
|
}
|
|
|
|
fp_register_t
|
|
ieee_add(fp_register_t fa, fp_register_t fb,
|
|
int src, int rnd,
|
|
u_int64_t control, u_int64_t *status)
|
|
{
|
|
int shift;
|
|
int expa, expb, exp;
|
|
u_int64_t fraca, fracb, frac;
|
|
int sign, sticky;
|
|
|
|
/* First handle NaNs */
|
|
if (isNaN(fa) || isNaN(fb)) {
|
|
fp_register_t result;
|
|
|
|
/* Instructions Descriptions (I) section 4.7.10.4 */
|
|
if (isQNaN(fb))
|
|
result = fb;
|
|
else if (isSNaN(fb)) {
|
|
result = fb;
|
|
result.t.fraction |= QNAN_BIT;
|
|
} else if (isQNaN(fa))
|
|
result = fa;
|
|
else if (isSNaN(fa)) {
|
|
result = fa;
|
|
result.t.fraction |= QNAN_BIT;
|
|
}
|
|
|
|
/* If either operand is a signalling NaN, trap. */
|
|
if (isSNaN(fa) || isSNaN(fb))
|
|
*status |= FPCR_INV;
|
|
|
|
return result;
|
|
}
|
|
|
|
/* Handle +/- INF */
|
|
if (isINF(fa))
|
|
if (isINF(fb))
|
|
if (fa.t.sign != fb.t.sign) {
|
|
/* If adding -INF to +INF, generate a trap. */
|
|
*status |= FPCR_INV;
|
|
return makeQNaN();
|
|
} else
|
|
return fa;
|
|
else
|
|
return fa;
|
|
else if (isINF(fb))
|
|
return fb;
|
|
|
|
/*
|
|
* Unpack the registers.
|
|
*/
|
|
expa = getexp(fa, src);
|
|
expb = getexp(fb, src);
|
|
fraca = getfrac(fa);
|
|
fracb = getfrac(fb);
|
|
shift = expa - expb;
|
|
if (shift < 0) {
|
|
shift = -shift;
|
|
exp = expb;
|
|
sticky = (fraca & ((1ULL << shift) - 1)) != 0;
|
|
if (shift >= 64)
|
|
fraca = sticky;
|
|
else
|
|
fraca = (fraca >> shift) | sticky;
|
|
} else if (shift > 0) {
|
|
exp = expa;
|
|
sticky = (fracb & ((1ULL << shift) - 1)) != 0;
|
|
if (shift >= 64)
|
|
fracb = sticky;
|
|
else
|
|
fracb = (fracb >> shift) | sticky;
|
|
} else
|
|
exp = expa;
|
|
if (fa.t.sign) fraca = -fraca;
|
|
if (fb.t.sign) fracb = -fracb;
|
|
frac = fraca + fracb;
|
|
if (frac >> 63) {
|
|
sign = 1;
|
|
frac = -frac;
|
|
} else
|
|
sign = 0;
|
|
|
|
/* -0 + -0 = -0 */
|
|
if (fa.t.exponent == 0 && fa.t.fraction == 0
|
|
&& fb.t.exponent == 0 && fb.t.fraction == 0)
|
|
sign = fa.t.sign && fb.t.sign;
|
|
|
|
return makefloat(sign, exp, frac, src, rnd, control, status);
|
|
}
|
|
|
|
fp_register_t
|
|
ieee_sub(fp_register_t fa, fp_register_t fb,
|
|
int src, int rnd,
|
|
u_int64_t control, u_int64_t *status)
|
|
{
|
|
fb.t.sign = !fb.t.sign;
|
|
return ieee_add(fa, fb, src, rnd, control, status);
|
|
}
|
|
|
|
typedef struct {
|
|
u_int64_t lo;
|
|
u_int64_t hi;
|
|
} u_int128_t;
|
|
|
|
#define SRL128(x, b) \
|
|
do { \
|
|
x.lo >>= b; \
|
|
x.lo |= x.hi << (64 - b); \
|
|
x.hi >>= b; \
|
|
} while (0)
|
|
|
|
#define SLL128(x, b) \
|
|
do { \
|
|
if (b >= 64) { \
|
|
x.hi = x.lo << (b - 64); \
|
|
x.lo = 0; \
|
|
} else { \
|
|
x.hi <<= b; \
|
|
x.hi |= x.lo >> (64 - b); \
|
|
x.lo <<= b; \
|
|
} \
|
|
} while (0)
|
|
|
|
#define SUB128(a, b) \
|
|
do { \
|
|
int borrow = a.lo < b.lo; \
|
|
a.lo = a.lo - b.lo; \
|
|
a.hi = a.hi - b.hi - borrow; \
|
|
} while (0)
|
|
|
|
#define LESS128(a, b) (a.hi < b.hi || (a.hi == b.hi && a.lo < b.lo))
|
|
|
|
fp_register_t
|
|
ieee_mul(fp_register_t fa, fp_register_t fb,
|
|
int src, int rnd,
|
|
u_int64_t control, u_int64_t *status)
|
|
{
|
|
int expa, expb, exp;
|
|
u_int64_t fraca, fracb, tmp;
|
|
u_int128_t frac;
|
|
int sign;
|
|
|
|
/* First handle NaNs */
|
|
if (isNaN(fa) || isNaN(fb)) {
|
|
fp_register_t result;
|
|
|
|
/* Instructions Descriptions (I) section 4.7.10.4 */
|
|
if (isQNaN(fb))
|
|
result = fb;
|
|
else if (isSNaN(fb)) {
|
|
result = fb;
|
|
result.t.fraction |= QNAN_BIT;
|
|
} else if (isQNaN(fa))
|
|
result = fa;
|
|
else if (isSNaN(fa)) {
|
|
result = fa;
|
|
result.t.fraction |= QNAN_BIT;
|
|
}
|
|
|
|
/* If either operand is a signalling NaN, trap. */
|
|
if (isSNaN(fa) || isSNaN(fb))
|
|
*status |= FPCR_INV;
|
|
|
|
return result;
|
|
}
|
|
|
|
/* Handle INF and 0 */
|
|
if ((isINF(fa) && isZERO(fb)) || (isZERO(fa) && isINF(fb))) {
|
|
/* INF * 0 = NaN */
|
|
*status |= FPCR_INV;
|
|
return makeQNaN();
|
|
} else
|
|
/* If either is INF or zero, get the sign right */
|
|
if (isINF(fa) || isINF(fb))
|
|
return makeINF(fa.t.sign ^ fb.t.sign);
|
|
else if (isZERO(fa) || isZERO(fb))
|
|
return makeZERO(fa.t.sign ^ fb.t.sign);
|
|
|
|
/*
|
|
* Unpack the registers.
|
|
*/
|
|
expa = getexp(fa, src);
|
|
expb = getexp(fb, src);
|
|
fraca = getfrac(fa);
|
|
fracb = getfrac(fb);
|
|
sign = fa.t.sign ^ fb.t.sign;
|
|
|
|
#define LO32(x) ((x) & ((1ULL << 32) - 1))
|
|
#define HI32(x) ((x) >> 32)
|
|
|
|
/*
|
|
* Calculate the 128bit result of multiplying fraca and fracb.
|
|
*/
|
|
frac.lo = fraca * fracb;
|
|
#ifdef __alpha__
|
|
/*
|
|
* The alpha has a handy instruction to find the high word.
|
|
*/
|
|
__asm__ __volatile__ ("umulh %1,%2,%0"
|
|
: "=r"(tmp)
|
|
: "r"(fraca), "r"(fracb));
|
|
frac.hi = tmp;
|
|
#else
|
|
/*
|
|
* Do the multiply longhand otherwise.
|
|
*/
|
|
frac.hi = HI32(LO32(fraca) * HI32(fracb)
|
|
+ HI32(fraca) * LO32(fracb)
|
|
+ HI32(LO32(fraca) * LO32(fracb)))
|
|
+ HI32(fraca) * HI32(fracb);
|
|
#endif
|
|
exp = expa + expb - (T_FRACBITS + FRAC_SHIFT);
|
|
|
|
while (frac.hi > 0) {
|
|
int sticky;
|
|
exp++;
|
|
sticky = frac.lo & 1;
|
|
SRL128(frac, 1);
|
|
frac.lo |= sticky;
|
|
}
|
|
|
|
return makefloat(sign, exp, frac.lo, src, rnd, control, status);
|
|
}
|
|
|
|
static u_int128_t
|
|
divide_128(u_int128_t a, u_int128_t b)
|
|
{
|
|
u_int128_t result;
|
|
u_int64_t bit;
|
|
int i;
|
|
|
|
/*
|
|
* Make a couple of assumptions on the numbers passed in. The
|
|
* value in 'a' will have bits set in the upper 64 bits only
|
|
* and the number in 'b' will have zeros in the upper 64 bits.
|
|
* Also, 'b' will not be zero.
|
|
*/
|
|
#ifdef TEST
|
|
if (a.hi == 0 || b.hi != 0 || b.lo == 0)
|
|
abort();
|
|
#endif
|
|
|
|
/*
|
|
* Find out how many bits of zeros are at the beginning of the divisor.
|
|
*/
|
|
i = 64;
|
|
bit = 1ULL << 63;
|
|
while (i < 127) {
|
|
if (b.lo & bit)
|
|
break;
|
|
i++;
|
|
bit >>= 1;
|
|
}
|
|
|
|
/*
|
|
* Find out how much to shift the divisor so that its msb
|
|
* matches the msb of the dividend.
|
|
*/
|
|
bit = 1ULL << 63;
|
|
while (i) {
|
|
if (a.hi & bit)
|
|
break;
|
|
i--;
|
|
bit >>= 1;
|
|
}
|
|
|
|
result.lo = 0;
|
|
result.hi = 0;
|
|
SLL128(b, i);
|
|
|
|
/*
|
|
* Calculate the result in two parts to avoid keeping a 128bit
|
|
* value for the result bit.
|
|
*/
|
|
if (i >= 64) {
|
|
bit = 1ULL << (i - 64);
|
|
while (bit) {
|
|
if (!LESS128(a, b)) {
|
|
result.hi |= bit;
|
|
SUB128(a, b);
|
|
if (!a.lo && !a.hi)
|
|
return result;
|
|
}
|
|
bit >>= 1;
|
|
SRL128(b, 1);
|
|
}
|
|
i = 63;
|
|
}
|
|
bit = 1ULL << i;
|
|
while (bit) {
|
|
if (!LESS128(a, b)) {
|
|
result.lo |= bit;
|
|
SUB128(a, b);
|
|
if (!a.lo && !a.hi)
|
|
return result;
|
|
}
|
|
bit >>= 1;
|
|
SRL128(b, 1);
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
fp_register_t
|
|
ieee_div(fp_register_t fa, fp_register_t fb,
|
|
int src, int rnd,
|
|
u_int64_t control, u_int64_t *status)
|
|
{
|
|
int expa, expb, exp;
|
|
u_int128_t fraca, fracb, frac;
|
|
int sign;
|
|
|
|
/* First handle NaNs, INFs and ZEROs */
|
|
if (isNaN(fa) || isNaN(fb)) {
|
|
fp_register_t result;
|
|
|
|
/* Instructions Descriptions (I) section 4.7.10.4 */
|
|
if (isQNaN(fb))
|
|
result = fb;
|
|
else if (isSNaN(fb)) {
|
|
result = fb;
|
|
result.t.fraction |= QNAN_BIT;
|
|
} else if (isQNaN(fa))
|
|
result = fa;
|
|
else if (isSNaN(fa)) {
|
|
result = fa;
|
|
result.t.fraction |= QNAN_BIT;
|
|
}
|
|
|
|
/* If either operand is a signalling NaN, trap. */
|
|
if (isSNaN(fa) || isSNaN(fb))
|
|
*status |= FPCR_INV;
|
|
|
|
return result;
|
|
}
|
|
|
|
/* Handle INF and 0 */
|
|
if (isINF(fa) && isINF(fb)) {
|
|
*status |= FPCR_INV;
|
|
return makeQNaN();
|
|
} else if (isZERO(fb))
|
|
if (isZERO(fa)) {
|
|
*status |= FPCR_INV;
|
|
return makeQNaN();
|
|
} else {
|
|
*status |= FPCR_DZE;
|
|
return makeINF(fa.t.sign ^ fb.t.sign);
|
|
}
|
|
else if (isZERO(fa))
|
|
return makeZERO(fa.t.sign ^ fb.t.sign);
|
|
|
|
/*
|
|
* Unpack the registers.
|
|
*/
|
|
expa = getexp(fa, src);
|
|
expb = getexp(fb, src);
|
|
fraca.hi = getfrac(fa);
|
|
fraca.lo = 0;
|
|
fracb.lo = getfrac(fb);
|
|
fracb.hi = 0;
|
|
sign = fa.t.sign ^ fb.t.sign;
|
|
|
|
frac = divide_128(fraca, fracb);
|
|
|
|
exp = expa - expb - (64 - T_FRACBITS - FRAC_SHIFT);
|
|
while (frac.hi > 0) {
|
|
int sticky;
|
|
exp++;
|
|
sticky = frac.lo & 1;
|
|
SRL128(frac, 1);
|
|
frac.lo |= sticky;
|
|
}
|
|
frac.lo |= 1;
|
|
|
|
return makefloat(sign, exp, frac.lo, src, rnd, control, status);
|
|
}
|
|
|
|
#define IEEE_TRUE 0x4000000000000000ULL
|
|
#define IEEE_FALSE 0
|
|
|
|
fp_register_t
|
|
ieee_cmpun(fp_register_t fa, fp_register_t fb, u_int64_t *status)
|
|
{
|
|
fp_register_t result;
|
|
if (isNaN(fa) || isNaN(fb)) {
|
|
if (isSNaN(fa) || isSNaN(fb))
|
|
*status |= FPCR_INV;
|
|
result.q = IEEE_TRUE;
|
|
} else
|
|
result.q = IEEE_FALSE;
|
|
|
|
return result;
|
|
}
|
|
|
|
fp_register_t
|
|
ieee_cmpeq(fp_register_t fa, fp_register_t fb, u_int64_t *status)
|
|
{
|
|
fp_register_t result;
|
|
if (isNaN(fa) || isNaN(fb)) {
|
|
if (isSNaN(fa) || isSNaN(fb))
|
|
*status |= FPCR_INV;
|
|
result.q = IEEE_FALSE;
|
|
} else {
|
|
if (isZERO(fa) && isZERO(fb))
|
|
result.q = IEEE_TRUE;
|
|
else if (fa.q == fb.q)
|
|
result.q = IEEE_TRUE;
|
|
else
|
|
result.q = IEEE_FALSE;
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
fp_register_t
|
|
ieee_cmplt(fp_register_t fa, fp_register_t fb, u_int64_t *status)
|
|
{
|
|
fp_register_t result;
|
|
if (isNaN(fa) || isNaN(fb)) {
|
|
if (isSNaN(fa) || isSNaN(fb))
|
|
*status |= FPCR_INV;
|
|
result.q = IEEE_FALSE;
|
|
} else {
|
|
if (isZERO(fa) && isZERO(fb))
|
|
result.q = IEEE_FALSE;
|
|
else if (fa.t.sign) {
|
|
/* fa is negative */
|
|
if (!fb.t.sign)
|
|
/* fb is positive, return true */
|
|
result.q = IEEE_TRUE;
|
|
else if (fa.t.exponent > fb.t.exponent)
|
|
/* fa has a larger exponent, return true */
|
|
result.q = IEEE_TRUE;
|
|
else if (fa.t.exponent == fb.t.exponent
|
|
&& fa.t.fraction > fb.t.fraction)
|
|
/* compare fractions */
|
|
result.q = IEEE_TRUE;
|
|
else
|
|
result.q = IEEE_FALSE;
|
|
} else {
|
|
/* fa is positive */
|
|
if (fb.t.sign)
|
|
/* fb is negative, return false */
|
|
result.q = IEEE_FALSE;
|
|
else if (fb.t.exponent > fa.t.exponent)
|
|
/* fb has a larger exponent, return true */
|
|
result.q = IEEE_TRUE;
|
|
else if (fa.t.exponent == fb.t.exponent
|
|
&& fa.t.fraction < fb.t.fraction)
|
|
/* compare fractions */
|
|
result.q = IEEE_TRUE;
|
|
else
|
|
result.q = IEEE_FALSE;
|
|
}
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
fp_register_t
|
|
ieee_cmple(fp_register_t fa, fp_register_t fb, u_int64_t *status)
|
|
{
|
|
fp_register_t result;
|
|
if (isNaN(fa) || isNaN(fb)) {
|
|
if (isSNaN(fa) || isSNaN(fb))
|
|
*status |= FPCR_INV;
|
|
result.q = IEEE_FALSE;
|
|
} else {
|
|
if (isZERO(fa) && isZERO(fb))
|
|
result.q = IEEE_TRUE;
|
|
else if (fa.t.sign) {
|
|
/* fa is negative */
|
|
if (!fb.t.sign)
|
|
/* fb is positive, return true */
|
|
result.q = IEEE_TRUE;
|
|
else if (fa.t.exponent > fb.t.exponent)
|
|
/* fa has a larger exponent, return true */
|
|
result.q = IEEE_TRUE;
|
|
else if (fa.t.exponent == fb.t.exponent
|
|
&& fa.t.fraction >= fb.t.fraction)
|
|
/* compare fractions */
|
|
result.q = IEEE_TRUE;
|
|
else
|
|
result.q = IEEE_FALSE;
|
|
} else {
|
|
/* fa is positive */
|
|
if (fb.t.sign)
|
|
/* fb is negative, return false */
|
|
result.q = IEEE_FALSE;
|
|
else if (fb.t.exponent > fa.t.exponent)
|
|
/* fb has a larger exponent, return true */
|
|
result.q = IEEE_TRUE;
|
|
else if (fa.t.exponent == fb.t.exponent
|
|
&& fa.t.fraction <= fb.t.fraction)
|
|
/* compare fractions */
|
|
result.q = IEEE_TRUE;
|
|
else
|
|
result.q = IEEE_FALSE;
|
|
}
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
fp_register_t
|
|
ieee_convert_S_T(fp_register_t f, int rnd,
|
|
u_int64_t control, u_int64_t *status)
|
|
{
|
|
/*
|
|
* Handle exceptional values.
|
|
*/
|
|
if (isNaN(f)) {
|
|
/* Instructions Descriptions (I) section 4.7.10.1 */
|
|
f.t.fraction |= QNAN_BIT;
|
|
*status |= FPCR_INV;
|
|
}
|
|
if (isQNaN(f) || isINF(f))
|
|
return f;
|
|
|
|
/*
|
|
* If the number is a denormalised float, renormalise.
|
|
*/
|
|
if (isDENORM(f))
|
|
return makefloat(f.t.sign,
|
|
getexp(f, S_FORMAT),
|
|
getfrac(f),
|
|
T_FORMAT, rnd, control, status);
|
|
else
|
|
return f;
|
|
}
|
|
|
|
fp_register_t
|
|
ieee_convert_T_S(fp_register_t f, int rnd,
|
|
u_int64_t control, u_int64_t *status)
|
|
{
|
|
/*
|
|
* Handle exceptional values.
|
|
*/
|
|
if (isNaN(f)) {
|
|
/* Instructions Descriptions (I) section 4.7.10.1 */
|
|
f.t.fraction |= QNAN_BIT;
|
|
f.t.fraction &= ~S_FRACMASK;
|
|
*status |= FPCR_INV;
|
|
}
|
|
if (isQNaN(f) || isINF(f))
|
|
return f;
|
|
|
|
return makefloat(f.t.sign,
|
|
getexp(f, T_FORMAT),
|
|
getfrac(f),
|
|
S_FORMAT, rnd, control, status);
|
|
}
|
|
|
|
fp_register_t
|
|
ieee_convert_Q_S(fp_register_t f, int rnd,
|
|
u_int64_t control, u_int64_t *status)
|
|
{
|
|
u_int64_t frac = f.q;
|
|
int sign, exponent;
|
|
|
|
if (frac >> 63) {
|
|
sign = 1;
|
|
frac = -frac;
|
|
} else
|
|
sign = 0;
|
|
|
|
/*
|
|
* We shift up one bit to leave the sticky bit clear. This is
|
|
* possible unless frac == (1<<63), in which case the sticky
|
|
* bit is already clear.
|
|
*/
|
|
exponent = T_FRACBITS + FRAC_SHIFT;
|
|
if (frac < (1ULL << 63)) {
|
|
frac <<= 1;
|
|
exponent--;
|
|
}
|
|
|
|
return makefloat(sign, exponent, frac, S_FORMAT, rnd,
|
|
control, status);
|
|
}
|
|
|
|
fp_register_t
|
|
ieee_convert_Q_T(fp_register_t f, int rnd,
|
|
u_int64_t control, u_int64_t *status)
|
|
{
|
|
u_int64_t frac = f.q;
|
|
int sign, exponent;
|
|
|
|
if (frac >> 63) {
|
|
sign = 1;
|
|
frac = -frac;
|
|
} else
|
|
sign = 0;
|
|
|
|
/*
|
|
* We shift up one bit to leave the sticky bit clear. This is
|
|
* possible unless frac == (1<<63), in which case the sticky
|
|
* bit is already clear.
|
|
*/
|
|
exponent = T_FRACBITS + FRAC_SHIFT;
|
|
if (frac < (1ULL << 63)) {
|
|
frac <<= 1;
|
|
exponent--;
|
|
}
|
|
|
|
return makefloat(sign, exponent, frac, T_FORMAT, rnd,
|
|
control, status);
|
|
}
|
|
|
|
fp_register_t
|
|
ieee_convert_T_Q(fp_register_t f, int rnd,
|
|
u_int64_t control, u_int64_t *status)
|
|
{
|
|
u_int64_t frac;
|
|
int exp;
|
|
|
|
/*
|
|
* Handle exceptional values.
|
|
*/
|
|
if (isNaN(f)) {
|
|
/* Instructions Descriptions (I) section 4.7.10.1 */
|
|
if (isSNaN(f))
|
|
*status |= FPCR_INV;
|
|
f.q = 0;
|
|
return f;
|
|
}
|
|
if (isINF(f)) {
|
|
/* Instructions Descriptions (I) section 4.7.10.1 */
|
|
*status |= FPCR_INV;
|
|
f.q = 0;
|
|
return f;
|
|
}
|
|
|
|
exp = getexp(f, T_FORMAT) - (T_FRACBITS + FRAC_SHIFT);
|
|
frac = getfrac(f);
|
|
|
|
if (exp > 0) {
|
|
if (exp > 64 || frac >= (1 << (64 - exp)))
|
|
*status |= FPCR_IOV | FPCR_INE;
|
|
if (exp < 64)
|
|
frac <<= exp;
|
|
else
|
|
frac = 0;
|
|
} else if (exp < 0) {
|
|
u_int64_t mask;
|
|
u_int64_t fraclo, frachi;
|
|
u_int64_t diffhi, difflo;
|
|
exp = -exp;
|
|
if (exp > 64) {
|
|
fraclo = 0;
|
|
diffhi = 0;
|
|
difflo = 0;
|
|
if (frac) {
|
|
frachi = 1;
|
|
*status |= FPCR_INE;
|
|
} else
|
|
frachi = 0;
|
|
} else if (exp == 64) {
|
|
fraclo = 0;
|
|
if (frac) {
|
|
frachi = 1;
|
|
difflo = frac;
|
|
diffhi = -frac;
|
|
*status |= FPCR_INE;
|
|
} else {
|
|
frachi = 0;
|
|
difflo = 0;
|
|
diffhi = 0;
|
|
}
|
|
} else {
|
|
mask = (1 << exp) - 1;
|
|
fraclo = frac >> exp;
|
|
if (frac & mask) {
|
|
frachi = fraclo + 1;
|
|
difflo = frac - (fraclo << exp);
|
|
diffhi = (frachi << exp) - frac;
|
|
*status |= FPCR_INE;
|
|
} else {
|
|
frachi = fraclo;
|
|
difflo = 0;
|
|
diffhi = 0;
|
|
}
|
|
}
|
|
switch (rnd) {
|
|
case ROUND_CHOP:
|
|
frac = fraclo;
|
|
break;
|
|
case ROUND_MINUS_INF:
|
|
if (f.t.sign)
|
|
frac = frachi;
|
|
else
|
|
frac = fraclo;
|
|
break;
|
|
case ROUND_NORMAL:
|
|
#if 0
|
|
/*
|
|
* Round to nearest.
|
|
*/
|
|
if (difflo < diffhi)
|
|
frac = fraclo;
|
|
else if (diffhi > difflo)
|
|
frac = frachi;
|
|
else if (fraclo & 1)
|
|
frac = frachi;
|
|
else
|
|
frac = fraclo;
|
|
#else
|
|
/*
|
|
* Round to zero.
|
|
*/
|
|
frac = fraclo;
|
|
#endif
|
|
break;
|
|
case ROUND_PLUS_INF:
|
|
if (f.t.sign)
|
|
frac = fraclo;
|
|
else
|
|
frac = frachi;
|
|
break;
|
|
}
|
|
}
|
|
|
|
if (f.t.sign) {
|
|
if (frac > (1ULL << 63))
|
|
*status |= FPCR_IOV | FPCR_INE;
|
|
frac = -frac;
|
|
} else {
|
|
if (frac > (1ULL << 63) - 1)
|
|
*status |= FPCR_IOV | FPCR_INE;
|
|
}
|
|
|
|
f.q = frac;
|
|
return f;
|
|
}
|
|
|
|
fp_register_t
|
|
ieee_convert_S_Q(fp_register_t f, int rnd,
|
|
u_int64_t control, u_int64_t *status)
|
|
{
|
|
f = ieee_convert_S_T(f, rnd, control, status);
|
|
return ieee_convert_T_Q(f, rnd, control, status);
|
|
}
|
|
|
|
#ifndef _KERNEL
|
|
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
#include <stdlib.h>
|
|
|
|
union value {
|
|
double d;
|
|
fp_register_t r;
|
|
};
|
|
|
|
|
|
static double
|
|
random_double()
|
|
{
|
|
union value a;
|
|
int exp;
|
|
|
|
a.r.t.fraction = ((long long)random() & (1ULL << 20) - 1) << 32
|
|
| random();
|
|
exp = random() & 0x7ff;
|
|
#if 1
|
|
if (exp == 0)
|
|
exp = 1; /* no denorms */
|
|
else if (exp == 0x7ff)
|
|
exp = 0x7fe; /* no NaNs and INFs */
|
|
#endif
|
|
|
|
a.r.t.exponent = exp;
|
|
a.r.t.sign = random() & 1;
|
|
return a.d;
|
|
}
|
|
|
|
static float
|
|
random_float()
|
|
{
|
|
union value a;
|
|
int exp;
|
|
|
|
a.r.t.fraction = ((long)random() & (1ULL << 23) - 1) << 29;
|
|
exp = random() & 0xff;
|
|
#if 1
|
|
if (exp == 0)
|
|
exp = 1; /* no denorms */
|
|
else if (exp == 0xff)
|
|
exp = 0xfe; /* no NaNs and INFs */
|
|
#endif
|
|
|
|
/* map exponent from S to T format */
|
|
if (exp == 255)
|
|
a.r.t.exponent = 0x7ff;
|
|
else if (exp & 0x80)
|
|
a.r.t.exponent = 0x400 + (exp & 0x7f);
|
|
else if (exp)
|
|
a.r.t.exponent = 0x380 + exp;
|
|
else
|
|
a.r.t.exponent = 0;
|
|
a.r.t.sign = random() & 1;
|
|
|
|
return a.d;
|
|
}
|
|
|
|
/*
|
|
* Ignore epsilon errors
|
|
*/
|
|
int
|
|
equal_T(union value a, union value b)
|
|
{
|
|
if (isZERO(a.r) && isZERO(b.r))
|
|
return 1;
|
|
if (a.r.t.sign != b.r.t.sign)
|
|
return 0;
|
|
if (a.r.t.exponent != b.r.t.exponent)
|
|
return 0;
|
|
|
|
return a.r.t.fraction == b.r.t.fraction;
|
|
}
|
|
|
|
int
|
|
equal_S(union value a, union value b)
|
|
{
|
|
int64_t epsilon = 1ULL << 29;
|
|
|
|
if (isZERO(a.r) && isZERO(b.r))
|
|
return 1;
|
|
if (a.r.t.sign != b.r.t.sign)
|
|
return 0;
|
|
if (a.r.t.exponent != b.r.t.exponent)
|
|
return 0;
|
|
|
|
return ((a.r.t.fraction & ~(epsilon-1))
|
|
== (b.r.t.fraction & ~(epsilon-1)));
|
|
}
|
|
|
|
#define ITER 1000000
|
|
|
|
static void
|
|
test_double_add()
|
|
{
|
|
union value a, b, c, x;
|
|
u_int64_t status = 0;
|
|
int i;
|
|
|
|
for (i = 0; i < ITER; i++) {
|
|
a.d = random_double();
|
|
b.d = random_double();
|
|
status = 0;
|
|
c.r = ieee_add(a.r, b.r, T_FORMAT, ROUND_NORMAL,
|
|
0, &status);
|
|
/* ignore NaN and INF */
|
|
if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
|
|
continue;
|
|
x.d = a.d + b.d;
|
|
if (!equal_T(c, x)) {
|
|
printf("bad double add, %g + %g = %g (should be %g)\n",
|
|
a.d, b.d, c.d, x.d);
|
|
c.r = ieee_add(a.r, b.r, T_FORMAT, ROUND_NORMAL,
|
|
0, &status);
|
|
}
|
|
}
|
|
}
|
|
|
|
static void
|
|
test_single_add()
|
|
{
|
|
union value a, b, c, x, t;
|
|
float xf;
|
|
u_int64_t status = 0;
|
|
int i;
|
|
|
|
for (i = 0; i < ITER; i++) {
|
|
#if 0
|
|
if (i == 0) {
|
|
a.r.q = 0xb33acf292ca49700ULL;
|
|
b.r.q = 0xcad3191058a693aeULL;
|
|
}
|
|
#endif
|
|
a.d = random_float();
|
|
b.d = random_float();
|
|
status = 0;
|
|
c.r = ieee_add(a.r, b.r, S_FORMAT, ROUND_NORMAL,
|
|
0, &status);
|
|
/* ignore NaN and INF */
|
|
if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
|
|
continue;
|
|
xf = a.d + b.d;
|
|
x.d = xf;
|
|
t.r = ieee_convert_S_T(c.r, ROUND_NORMAL, 0, &status);
|
|
if (!equal_S(t, x)) {
|
|
printf("bad single add, %g + %g = %g (should be %g)\n",
|
|
a.d, b.d, t.d, x.d);
|
|
c.r = ieee_add(a.r, b.r, S_FORMAT, ROUND_NORMAL,
|
|
0, &status);
|
|
}
|
|
}
|
|
}
|
|
|
|
static void
|
|
test_double_mul()
|
|
{
|
|
union value a, b, c, x;
|
|
u_int64_t status = 0;
|
|
int i;
|
|
|
|
for (i = 0; i < ITER; i++) {
|
|
a.d = random_double();
|
|
b.d = random_double();
|
|
status = 0;
|
|
c.r = ieee_mul(a.r, b.r, T_FORMAT, ROUND_NORMAL,
|
|
0, &status);
|
|
/* ignore NaN and INF */
|
|
if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
|
|
continue;
|
|
x.d = a.d * b.d;
|
|
if (!equal_T(c, x)) {
|
|
printf("bad double mul, %g * %g = %g (should be %g)\n",
|
|
a.d, b.d, c.d, x.d);
|
|
c.r = ieee_mul(a.r, b.r, T_FORMAT, ROUND_NORMAL,
|
|
0, &status);
|
|
}
|
|
}
|
|
}
|
|
|
|
static void
|
|
test_single_mul()
|
|
{
|
|
union value a, b, c, x, t;
|
|
float xf;
|
|
u_int64_t status = 0;
|
|
int i;
|
|
|
|
for (i = 0; i < ITER; i++) {
|
|
a.d = random_double();
|
|
b.d = random_double();
|
|
status = 0;
|
|
c.r = ieee_mul(a.r, b.r, S_FORMAT, ROUND_NORMAL,
|
|
0, &status);
|
|
/* ignore NaN and INF */
|
|
if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
|
|
continue;
|
|
xf = a.d * b.d;
|
|
x.d = xf;
|
|
t.r = ieee_convert_S_T(c.r, ROUND_NORMAL, 0, &status);
|
|
if (!equal_S(t, x)) {
|
|
printf("bad single mul, %g * %g = %g (should be %g)\n",
|
|
a.d, b.d, t.d, x.d);
|
|
c.r = ieee_mul(a.r, b.r, T_FORMAT, ROUND_NORMAL,
|
|
0, &status);
|
|
}
|
|
}
|
|
}
|
|
|
|
static void
|
|
test_double_div()
|
|
{
|
|
union value a, b, c, x;
|
|
u_int64_t status = 0;
|
|
int i;
|
|
|
|
for (i = 0; i < ITER; i++) {
|
|
a.d = random_double();
|
|
b.d = random_double();
|
|
status = 0;
|
|
c.r = ieee_div(a.r, b.r, T_FORMAT, ROUND_NORMAL,
|
|
0, &status);
|
|
/* ignore NaN and INF */
|
|
if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
|
|
continue;
|
|
x.d = a.d / b.d;
|
|
if (!equal_T(c, x) && !isZERO(x.r)) {
|
|
printf("bad double div, %g / %g = %g (should be %g)\n",
|
|
a.d, b.d, c.d, x.d);
|
|
c.r = ieee_div(a.r, b.r, T_FORMAT, ROUND_NORMAL,
|
|
0, &status);
|
|
}
|
|
}
|
|
}
|
|
|
|
static void
|
|
test_single_div()
|
|
{
|
|
union value a, b, c, x, t;
|
|
float xf;
|
|
u_int64_t status = 0;
|
|
int i;
|
|
|
|
for (i = 0; i < ITER; i++) {
|
|
a.d = random_double();
|
|
b.d = random_double();
|
|
status = 0;
|
|
c.r = ieee_div(a.r, b.r, S_FORMAT, ROUND_NORMAL,
|
|
0, &status);
|
|
/* ignore NaN and INF */
|
|
if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
|
|
continue;
|
|
xf = a.d / b.d;
|
|
x.d = xf;
|
|
t.r = ieee_convert_S_T(c.r, ROUND_NORMAL, 0, &status);
|
|
if (!equal_S(t, x)) {
|
|
printf("bad single div, %g / %g = %g (should be %g)\n",
|
|
a.d, b.d, t.d, x.d);
|
|
c.r = ieee_mul(a.r, b.r, T_FORMAT, ROUND_NORMAL,
|
|
0, &status);
|
|
}
|
|
}
|
|
}
|
|
|
|
static void
|
|
test_convert_int_to_double()
|
|
{
|
|
union value a, c, x;
|
|
u_int64_t status = 0;
|
|
int i;
|
|
|
|
for (i = 0; i < ITER; i++) {
|
|
a.r.q = (u_int64_t)random() << 32
|
|
| random();
|
|
status = 0;
|
|
c.r = ieee_convert_Q_T(a.r, ROUND_NORMAL, 0, &status);
|
|
/* ignore NaN and INF */
|
|
if (isNaN(c.r) || isINF(c.r))
|
|
continue;
|
|
x.d = (double) a.r.q;
|
|
if (c.d != x.d) {
|
|
printf("bad convert double, (double)%qx = %g (should be %g)\n",
|
|
a.r.q, c.d, x.d);
|
|
c.r = ieee_convert_Q_T(a.r, ROUND_NORMAL, 0, &status);
|
|
}
|
|
}
|
|
}
|
|
|
|
static void
|
|
test_convert_int_to_single()
|
|
{
|
|
union value a, c, x, t;
|
|
float xf;
|
|
u_int64_t status = 0;
|
|
int i;
|
|
|
|
for (i = 0; i < ITER; i++) {
|
|
a.r.q = (unsigned long long)random() << 32
|
|
| random();
|
|
status = 0;
|
|
c.r = ieee_convert_Q_S(a.r, ROUND_NORMAL, 0, &status);
|
|
/* ignore NaN and INF */
|
|
if (isNaN(c.r) || isINF(c.r))
|
|
continue;
|
|
xf = (float) a.r.q;
|
|
x.d = xf;
|
|
t.r = ieee_convert_S_T(c.r, ROUND_NORMAL, 0, &status);
|
|
if (t.d != x.d) {
|
|
printf("bad convert single, (double)%qx = %g (should be %g)\n",
|
|
a.r.q, c.d, x.d);
|
|
c.r = ieee_convert_Q_S(a.r, ROUND_NORMAL, 0, &status);
|
|
}
|
|
}
|
|
}
|
|
|
|
static void
|
|
test_convert_double_to_int()
|
|
{
|
|
union value a, c;
|
|
u_int64_t status = 0;
|
|
int i;
|
|
|
|
for (i = 0; i < ITER; i++) {
|
|
a.d = random_double();
|
|
status = 0;
|
|
c.r = ieee_convert_T_Q(a.r, ROUND_NORMAL, 0, &status);
|
|
if ((int)c.r.q != (int)a.d) {
|
|
printf("bad convert double, (int)%g = %d (should be %d)\n",
|
|
a.d, (int)c.r.q, (int)a.d);
|
|
c.r = ieee_convert_T_Q(a.r, ROUND_NORMAL, 0, &status);
|
|
}
|
|
}
|
|
}
|
|
|
|
int
|
|
main(int argc, char* argv[])
|
|
{
|
|
srandom(0);
|
|
|
|
test_double_div();
|
|
test_single_div();
|
|
test_double_add();
|
|
test_single_add();
|
|
test_double_mul();
|
|
test_single_mul();
|
|
test_convert_int_to_double();
|
|
test_convert_int_to_single();
|
|
#if 0
|
|
/* x86 generates SIGFPE on overflows. */
|
|
test_convert_double_to_int();
|
|
#endif
|
|
|
|
return 0;
|
|
}
|
|
|
|
#endif
|