Subversion Repositories HelenOS

Compare Revisions

Ignore whitespace Rev 827 → Rev 828

//uspace/trunk/softfloat/include/div.h
30,6 → 30,8
#define __DIV_H__
 
float32 divFloat32(float32 a, float32 b);
float64 divFloat64(float64 a, float64 b);
 
__u64 divFloat64estim(__u64 a, __u64 b);
#endif
 
//uspace/trunk/softfloat/generic/softfloat.c
120,6 → 120,14
return divFloat32(fa, fb).f;
}
 
double __divdf3(double a, double b)
{
float64 da, db;
da.d = a;
db.d = b;
return divFloat64(da, db).d;
}
 
float __negsf2(float a)
{
float32 fa;
//uspace/trunk/softfloat/generic/div.c
28,7 → 28,9
 
#include<sftypes.h>
#include<add.h>
#include<div.h>
#include<comparison.h>
#include<mul.h>
 
float32 divFloat32(float32 a, float32 b)
{
140,7 → 142,7
/* pack and round */
/* TODO: find first nonzero digit and shift result and detect possibly underflow */
/* find first nonzero digit and shift result and detect possibly underflow */
while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) {
cexp--;
cfrac <<= 1;
147,7 → 149,6
/* TODO: fix underflow */
};
cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/
if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
185,3 → 186,208
return result;
}
 
float64 divFloat64(float64 a, float64 b)
{
float64 result;
__s32 aexp, bexp, cexp;
__u64 afrac, bfrac, cfrac;
__u64 remlo, remhi;
result.parts.sign = a.parts.sign ^ b.parts.sign;
if (isFloat64NaN(a)) {
if (isFloat64SigNaN(a)) {
/*FIXME: SigNaN*/
}
/*NaN*/
return a;
}
if (isFloat64NaN(b)) {
if (isFloat64SigNaN(b)) {
/*FIXME: SigNaN*/
}
/*NaN*/
return b;
}
if (isFloat64Infinity(a)) {
if (isFloat64Infinity(b)) {
/*FIXME: inf / inf */
result.binary = FLOAT64_NAN;
return result;
}
/* inf / num */
result.parts.exp = a.parts.exp;
result.parts.fraction = a.parts.fraction;
return result;
}
 
if (isFloat64Infinity(b)) {
if (isFloat64Zero(a)) {
/* FIXME 0 / inf */
result.parts.exp = 0;
result.parts.fraction = 0;
return result;
}
/* FIXME: num / inf*/
result.parts.exp = 0;
result.parts.fraction = 0;
return result;
}
if (isFloat64Zero(b)) {
if (isFloat64Zero(a)) {
/*FIXME: 0 / 0*/
result.binary = FLOAT64_NAN;
return result;
}
/* FIXME: division by zero */
result.parts.exp = 0;
result.parts.fraction = 0;
return result;
}
 
afrac = a.parts.fraction;
aexp = a.parts.exp;
bfrac = b.parts.fraction;
bexp = b.parts.exp;
/* denormalized numbers */
if (aexp == 0) {
if (afrac == 0) {
result.parts.exp = 0;
result.parts.fraction = 0;
return result;
}
/* normalize it*/
afrac <<= 1;
/* afrac is nonzero => it must stop */
while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) {
afrac <<= 1;
aexp--;
}
}
 
if (bexp == 0) {
bfrac <<= 1;
/* bfrac is nonzero => it must stop */
while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) {
bfrac <<= 1;
bexp--;
}
}
 
afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 );
bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1);
 
if ( bfrac <= (afrac << 1) ) {
afrac >>= 1;
aexp++;
}
cexp = aexp - bexp + FLOAT64_BIAS - 2;
cfrac = divFloat64estim(afrac, bfrac);
if (( cfrac & 0x1FF ) <= 2) { /*FIXME:?? */
mul64integers( bfrac, cfrac, &remlo, &remhi);
/* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/
remhi = afrac - remhi - ( remlo > 0);
remlo = - remlo;
while ((__s64) remhi < 0) {
cfrac--;
remlo += bfrac;
remhi += ( remlo < bfrac );
}
cfrac |= ( remlo != 0 );
}
/* pack and round */
/* find first nonzero digit and shift result and detect possibly underflow */
while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1 ) )))) {
cexp--;
cfrac <<= 1;
/* TODO: fix underflow */
};
cfrac >>= 1;
++cexp;
cfrac += (0x1 << (64 - FLOAT64_FRACTION_SIZE - 3));
 
if (cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1 ))) {
++cexp;
cfrac >>= 1;
}
 
/* check overflow */
if (cexp >= FLOAT64_MAX_EXPONENT ) {
/* FIXME: overflow, return infinity */
result.parts.exp = FLOAT64_MAX_EXPONENT;
result.parts.fraction = 0;
return result;
}
 
if (cexp < 0) {
/* FIXME: underflow */
result.parts.exp = 0;
if ((cexp + FLOAT64_FRACTION_SIZE) < 0) {
result.parts.fraction = 0;
return result;
}
cfrac >>= 1;
while (cexp < 0) {
cexp ++;
cfrac >>= 1;
}
return result;
} else {
cexp ++; /*normalized*/
result.parts.exp = (__u32)cexp;
}
result.parts.fraction = ((cfrac >>(64 - FLOAT64_FRACTION_SIZE - 2 ) ) & (~FLOAT64_HIDDEN_BIT_MASK));
return result;
}
 
__u64 divFloat64estim(__u64 a, __u64 b)
{
__u64 bhi;
__u64 remhi, remlo;
__u64 result;
if ( b <= a ) {
return 0xFFFFFFFFFFFFFFFFull;
}
bhi = b >> 32;
result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32;
mul64integers(b, result, &remlo, &remhi);
remhi = a - remhi - (remlo > 0);
remlo = - remlo;
 
b <<= 32;
while ( (__s64) remhi < 0 ) {
result -= 0x1ll << 32;
remlo += b;
remhi += bhi + ( remlo < b );
}
remhi = (remhi << 32) | (remlo >> 32);
if (( bhi << 32) <= remhi) {
result |= 0xFFFFFFFF;
} else {
result |= remhi / bhi;
}
return result;
}
 
//uspace/trunk/softfloat/generic/add.c
89,7 → 89,7
 
if (exp2 == 0) {
/* second operand is denormalized */
--expdiff;
--expdiff;
} else {
/* add hidden bit to second operand */
frac2 |= FLOAT32_HIDDEN_BIT_MASK;
118,7 → 118,8
frac1 >>= 1;
};
if ((a.parts.exp == FLOAT32_MAX_EXPONENT ) || (a.parts.exp < exp1)) {
if ((exp1 == FLOAT32_MAX_EXPONENT ) || (exp2 > exp1)) {
/* overflow - set infinity as result */
a.parts.exp = FLOAT32_MAX_EXPONENT;
a.parts.fraction = 0;
//uspace/trunk/softfloat/generic/mul.c
329,21 → 329,22
alow = a & 0xFFFFFFFF;
blow = b & 0xFFFFFFFF;
a <<= 32;
b <<= 32;
a >>= 32;
b >>= 32;
low = (__u64)alow * blow;
low = ((__u64)alow) * blow;
middle1 = a * blow;
middle2 = alow * b;
high = a * b;
 
middle1 += middle2;
high += ((__u64)(middle1 < middle2) << 32) + (middle1 >> 32);
high += (((__u64)(middle1 < middle2)) << 32) + (middle1 >> 32);
middle1 <<= 32;
low += middle1;
high += (low < middle1);
*lo = low;
*hi = high;
return;
}