Skip to content

Commit

Permalink
squash some precision issues
Browse files Browse the repository at this point in the history
  • Loading branch information
@ellie-idb
ellie-idb committed Jan 28, 2022
1 parent e3f6743 commit a16f586
Showing 1 changed file with 116 additions and 62 deletions.
178 changes: 116 additions & 62 deletions source/mir/bignum/low_level_view.d
View file
Original file line number Diff line number Diff line change
Expand Up @@ -269,33 +269,6 @@ void divm(scope ref BigUIntView!uint _u, scope BigUIntView!uint _v, scope ref Bi
_q.coefficients = _q.normalized.coefficients;
}

void divrem(scope BigUIntView!uint u, scope BigUIntView!uint v, scope ref BigUIntView!uint q, scope ref BigUIntView!uint r)
{
r.coefficients = u.coefficients.dup;
divm(r, v, q);
}
///
version(mir_bignum_test)
unittest
{
import mir.bignum.low_level_view : BigUIntView;
{
// Test division by a single-digit divisor here.
auto dividend = BigUIntView!uint.fromHexString("000000005");
auto divisor = BigUIntView!uint.fromHexString("5");
auto quotient = BigUIntView!uint.fromHexString("0");
auto remainder = BigUIntView!uint.fromHexString("0");

divrem(dividend, divisor, quotient, remainder);

assert(dividend == BigUIntView!uint.fromHexString("000000005"));
assert(divisor == BigUIntView!uint.fromHexString("5"));
assert(quotient == BigUIntView!uint.fromHexString("1"));
assert(remainder.coefficients.length == 0);

}
}

///
version(mir_bignum_test)
unittest
Expand Down Expand Up @@ -1669,6 +1642,8 @@ struct BigUIntView(W, WordEndian endian = TargetEndian)
return str.length - i;
}



static if (W.sizeof == size_t.sizeof && endian == TargetEndian)
///
version(mir_bignum_test)
Expand Down Expand Up @@ -3151,6 +3126,7 @@ struct DecimalView(W, WordEndian endian = TargetEndian, Exp = sizediff_t)
@optStrategy("minsize")
package T algoM(T)(T ret, scope BigUIntView!(const(uint)) coeff, int exponent)
{
import mir.bignum.fixed : UInt;
import mir.math.ieee: ldexp, frexp, nextDown, nextUp;

BigUIntAccumulator!uint u, v;
Expand All @@ -3159,52 +3135,125 @@ package T algoM(T)(T ret, scope BigUIntView!(const(uint)) coeff, int exponent)

v.coefficients.length++;
v.put(1);

if (exponent < 0)
{
v.mulPow5(-exponent);
v.mulPow2(-exponent);
auto degree = -exponent;

if (v.approxCanMulPow5(degree)) {
v.mulPow5(degree);
}
else
{
enum n = MaxWordPow5!uint;
auto growFactor = degree / n + (degree % n != 0);
v.coefficients.length += growFactor;
v.mulPow5(degree);
}

if (v.canMulPow2(degree))
{
v.mulPow2(degree);
}
else
{
import mir.bitop: ctlz;
enum n = uint.sizeof * 8;
auto growFactor = degree / n + (degree % n > ctlz(v.view.mostSignificant));
v.coefficients.length += growFactor;
v.mulPow2(degree);
}

}
else
{
auto degree = exponent;

if (u.approxCanMulPow5(degree))
{
u.mulPow5(degree);
}
else
{
enum n = MaxWordPow5!uint;
auto growFactor = degree / n + (degree % n != 0);

}
u.mulPow5(exponent);
u.mulPow2(exponent);
}

// We need to introduce a new digit position at the MSB -- do that here
u.coefficients.length++;
u.length = coeff.coefficients.length + 1;
// XXX: is this needed?
// v.coefficients = v.view.normalized.coefficients;

// Adopted from Rust's rawfp implementation (TBD: do we have these already?)
enum EXP_BITS = (T.sizeof * 8) - T.mant_dig;
enum MAX_EXP = (1 << (EXP_BITS - 1)) - 1;
enum MIN_EXP = -MAX_EXP + 1;
enum MAX_EXP_INT = MAX_EXP - (T.mant_dig - 1);
enum MIN_EXP_INT = MIN_EXP - (T.mant_dig - 1);

// MAX_SIG
enum betaN = (ulong(1) << T.mant_dig) - 1;
// MIN_SIG
enum betaNMinus1 = ulong(1) << (T.mant_dig - 1);

int k = 0;

// Ugly ugly ugly, but we don't want to cast to ulong to compare
// as there's a potential for truncation
auto min_sig = cast(BigUIntView!uint)UInt!64(betaNMinus1).view;
auto max_sig = cast(BigUIntView!uint)UInt!64(betaN).view;

BigUIntView!uint x, rem;
while (true)
{
// dynamically resize the quotient / remainder buffers
auto m = u.view.normalized.coefficients.length;
auto n = v.view.normalized.coefficients.length;
x.coefficients.length = (m - n) + 1;
rem.coefficients.length = n;
divrem(u.view, v.view, x, rem);
// As part of "scaling", we'll initially start with m < n as true,
// so we need to have this fall back here
if (m < n) {
x.coefficients.length = m;
}
else {
x.coefficients.length = (m - n) + 1;
}

// Division requires us to introduce
// a new digit in the MSB-place -- we do that here
rem.coefficients.length = u.view.coefficients.length + 1;
rem.coefficients[0 .. $ - 1] = u.view.coefficients;
rem.coefficients[$ - 1] = 0;
divm(rem, v.view, x);

// Check if we're about to underflow
if (k == MAX_EXP_INT)
{
if (x >= min_sig && x <= max_sig)
{
break;
}
// XXX: deal with underflow
assert(0, "underflow");
}

if (betaNMinus1 <= cast(ulong)x && cast(ulong)x < betaN)
// Check if our exponent is out of range
if (k > MAX_EXP_INT)
{
break;
return T.infinity;
}

if (cast(ulong)x < betaNMinus1)
if (x < min_sig)
{
if (!u.canMulPow2(1))
{
u.coefficients.length++;
}

u.mulPow2(1);
u.mulPow2(1);
k -= 1;
}
else if (cast(ulong)x >= betaN)
else if (max_sig <= x)
{
if (!v.canMulPow2(1))
{
Expand All @@ -3214,21 +3263,23 @@ package T algoM(T)(T ret, scope BigUIntView!(const(uint)) coeff, int exponent)
v.mulPow2(1);
k += 1;
}
else
{
break;
}
}

// This part is specified as ratio->float within the paper
ulong q = cast(ulong)x;
debug {
import std.stdio;
writefln("q: %s (x: %s)", q, x);
writefln("q: %s", cast(T)q);
}
T z = ldexp(cast(T)q, k);
v.view -= rem;
if (rem < v.view)
auto vmr = v.view;
vmr -= rem;

if (rem < vmr)
{
ret = z;
}
else if (rem > v.view)
else if (rem > vmr)
{
ret = nextUp(z);
}
Expand All @@ -3248,28 +3299,31 @@ package T algoM(T)(T ret, scope BigUIntView!(const(uint)) coeff, int exponent)
version(mir_bignum_test)
unittest
{
alias AliasSeq(T...) = T;

{
auto view = DecimalView!ulong(false, -8, BigUIntView!ulong.fromHexString("BEBC2000000011E1A3"));
auto coeff = (cast(BigUIntView!uint)view.coefficient).lightConst;
assert (algoM!double(0.0, coeff, cast(int)view.exponent) == 3.518437208883201171875E+013);
}

// TBD: triggers underflow
// {
// auto view = DecimalView!ulong(false, 0, BigUIntView!ulong.fromHexString("88BF4748507FB9900ADB624CCFF8D78897DC900FB0460327D4D86D327219"));
// auto coeff = (cast(BigUIntView!uint)view.coefficient).lightConst;
// debug {
// import std.stdio;
// writefln("%s", algoM!float(0.0, coeff, cast(int)view.exponent));
// writefln("%s", algoM!double(0.0, coeff, cast(int)view.exponent));
// }
// assert (algoM!float(0.0, coeff, cast(int)view.exponent) == float.infinity);
// assert (algoM!double(0.0, coeff, cast(int)view.exponent) == 0x1.117e8e90a0ff7p+239);
// }

{
auto view = DecimalView!ulong(false, 0, BigUIntView!ulong.fromHexString("88BF4748507FB9900ADB624CCFF8D78897DC900FB0460327D4D86D327219"));
auto view = DecimalView!ulong(false, -324, BigUIntView!ulong.fromHexString("4F0CEDC95A718E"));
auto coeff = (cast(BigUIntView!uint)view.coefficient).lightConst;
debug {
import std.stdio;
writefln("%s", algoM!float(0.0, coeff, cast(int)view.exponent));
writefln("%s", algoM!double(0.0, coeff, cast(int)view.exponent));
}
assert (algoM!float(0.0, coeff, cast(int)view.exponent) == float.infinity);
assert (algoM!double(0.0, coeff, cast(int)view.exponent) == 0x1.117e8e90a0ff7p+239);
assert (algoM!float(0.0, coeff, cast(int)view.exponent) == 0);
assert (algoM!double(0.0, coeff, cast(int)view.exponent) == 2.2250738585072014e-308);
}



}

@optStrategy("minsize")
Expand Down

0 comments on commit a16f586

Please sign in to comment.