jk-jeon / dragonbox

Reference implementation of Dragonbox in C++
Apache License 2.0
607 stars 39 forks source link

Paper/implementation mismatch? #31

Closed vyivel closed 2 years ago

vyivel commented 2 years ago

In practice, computing the parity of x(i) is still a heavy operation compared to others. Hence, it can be beneficial to delay it until it is absolutely necessary. In Algorithm 5.2, we need to compute the parity of x(i) when r and δ(i) turn out to be the same. However, in fact if either the interval I does not include the left endpoint or x is not an integer, we always conclude I ∩ 10^(−k0+1)Z is empty regardless of the parity of x(i). Therefore, if the given rounding mode indicates that the left endpoint is not in I, we can skip computing the parity of x(i).

(section 5.1.6, page 21)

If I understand this correctly, the condition of the left endpoint being excluded or x not being an integer is enough to say that the intersection in question is empty, and the parity doesn't need to be checked.

However, the algorithm implementation seems to disagree:

                    if (!interval_type.include_left_endpoint() ||
                        exponent < case_fc_pm_half_lower_threshold ||
                        exponent > divisibility_check_by_5_threshold) {
                        // If the left endpoint is not included, the condition for
                        // success is z^(f) < delta^(f) (odd parity).
                        // Otherwise, the inequalities on exponent ensure that
                        // x is not an integer, so if z^(f) >= delta^(f) (even parity), we in fact
                        // have strict inequality.
                        if (!compute_mul_parity(two_fl, cache, beta).parity) {
                            goto small_divisor_case_label;
                        }
                    }
                    else {
                        auto const [xi_parity, x_is_integer] =
                            compute_mul_parity(two_fl, cache, beta);
                        if (!xi_parity && !x_is_integer) {
                            goto small_divisor_case_label;
                        }
                    }

In fact, the comments explicitly state that checking parity is required. From quick testing it seems that the code is right here.

vyivel commented 2 years ago

Side note: I think that part can be simplified to something like the following pseudocode.

xi_parity = compute_mul_parity().parity
if (!xi_parity) {
    if (!include_left_endpoint) {
        goto small_divisor_case_label
    }
    // x(i) is an even number and the left endpoint is included, check if x is an integer
    if (exponent < case_fc_pm_half_lower_threshold || exponent > case_fc_pm_half_lower_threshold) {
        goto small_divisor_case_label
    }
    x_is_integer = compute_mul_parity().is_integer
    if (!x_is_integer) {
        goto small_divisor_case_label
    }
}

This matches Algorithm 5.2 (step 7) more closely.

jk-jeon commented 2 years ago

Hmm. This is indeed weird. I don't know why I wrote it like that 😂 Probably I came up with some wrong idea, tested it, figured out I was wrong, so corrected the code but in a stupid way, and also forgot to rewrite the paper. Thanks a lot for catching this!! I'll correct it right away.

jk-jeon commented 2 years ago

Thanks again, it's fixed now.