664a31e496
is an application space macro and the applications are supposed to be free to use it as they please (but cannot). This is consistant with the other BSD's who made this change quite some time ago. More commits to come.
1521 lines
31 KiB
C
1521 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.
|
|
*
|
|
* $FreeBSD$
|
|
*/
|
|
|
|
/*
|
|
* 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/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
|