Skip to content
Permalink
Browse files

Fix drift in Num->Str roundtrip; Make nqp::div_In 14x slower

Fixes R#1651 rakudo/rakudo#1651

The original code converted bigints to doubles right away, before dividing them,
which caused some precision loss that was responsible for repeated Num->Str
conversion to drift the value.

Fix by scaling up the divident by the bits of the divisor + additional
64 bits (that's needed to get correct answers for stuff like 1e-300). Then
performing integer division, and then converting the result to double, and
reversing the original scaling inside the to-double conversion routine.

Since we combine the scaling and denormals handling, we need to break up
the shift by power of two we do in the to-double conversion routine to do
it in 2**1023 steps (otherwise, the 2 raised to a larger power becomes Inf
(or denormal/-Inf for negative powers).

This work makes nqp::div_In op 14x slower, but I'm hopeful that after the
CaR grant, to have extra knowledge to speed it up.
  • Loading branch information...
zoffixznet committed Apr 4, 2018
1 parent 02dd490 commit b735866ddee9bd719440e5c822045430b335f27b
Showing with 42 additions and 9 deletions.
  1. +42 −9 src/math/bigintops.c
@@ -63,7 +63,7 @@ static const double mp_get_double_multiplier = (double)(MP_MASK + 1);

static MVMnum64 mp_get_double(mp_int *a, int shift) {
MVMnum64 d;
int i, limit;
int i, limit, final_shift;
d = 0.0;

mp_clamp(a);
@@ -78,8 +78,21 @@ static MVMnum64 mp_get_double(mp_int *a, int shift) {

if (a->sign == MP_NEG)
d *= -1.0;

return d * pow(2.0, i * DIGIT_BIT - shift);
final_shift = i * DIGIT_BIT - shift;
if (final_shift < 0) {
while (final_shift < -1023) {
d *= pow(2.0, -1023);
final_shift += 1023;
}
}
else {
while (final_shift > 1023) {
d *= pow(2.0, 1023);
final_shift -= 1023;
}
}
d *= pow(2.0, final_shift);
return d;
}

static void from_num(MVMnum64 d, mp_int *a) {
@@ -935,12 +948,32 @@ MVMnum64 MVM_bigint_div_num(MVMThreadContext *tc, MVMObject *a, MVMObject *b) {
mp_int *ia = force_bigint(ba, tmp);
mp_int *ib = force_bigint(bb, tmp);

int max_bits = MAX(mp_count_bits(ia), mp_count_bits(ib));
if (max_bits > MAX_BIGINT_BITS_IN_DOUBLE) {
c = mp_get_double(ia, max_bits - MAX_BIGINT_BITS_IN_DOUBLE)
/ mp_get_double(ib, max_bits - MAX_BIGINT_BITS_IN_DOUBLE);
} else {
c = mp_get_double(ia, 0) / mp_get_double(ib, 0);
mp_clamp(ib);
if (ib->used == 0) { /* zero-denominator special case */
if (ia->sign == MP_NEG)
c = -1e0/0e0;
else
c = 1e0/0e0;
/*
* we won't have NaN case here, since the branch requires at
* least one bigint to be big
*/
}
else {
mp_int scaled;
int bbits = mp_count_bits(ib)+64;

if (mp_init(&scaled) != MP_OKAY)
MVM_exception_throw_adhoc(tc,
"Failed to initialize bigint for scaled divident");
if (mp_mul_2d(ia, bbits, &scaled) != MP_OKAY)
MVM_exception_throw_adhoc(tc, "Failed to scale divident");
// simply re-use &scaled for result
if (mp_div(&scaled, ib, &scaled, NULL) != MP_OKAY)
MVM_exception_throw_adhoc(tc,
"Failed to preform bigint division");
c = mp_get_double(&scaled, bbits);
mp_clear(&scaled);
}
clear_temp_bigints(tmp, 2);
} else {

0 comments on commit b735866

Please sign in to comment.
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.