compute_shortest question!

[m3g4z0rdEX]

In the paper, in compute_shortest algorithm, can we send to input a, b, c, such that length(b) < length(c)? If yes, can we get bi(output final value) equal to zero, or zero + 1(by rounding)?

=====================================
For example in vacuum: a = 6, b = 8, c = 11, and as output we can get 1e1, which is not closest to input: right answer is just 8e0 ?

Where am i wrong?

[ulfjack]

You are correct that this is a corner case.

The way I defined "shortest" in the paper is this: if you compare two numbers, the one with the larger exponent is considered "shorter". In this case, 1e1 has a larger exponent than 8e0, so this is consistent with the definition.

This matches what other implementations are doing, as far as I can tell, and it only affects a very small set of numbers, if any. I tried a couple of values, but I don't immediately see any for which 9eX is closer and it goes up to 1e(X+1) due to this rule.

[StephanTLavavej]

If I understand the issue correctly - would the only possibly-affected numbers be those surrounding the powers of 10? I verified for double that Ryu emits different strings for all of the powers of 10 and their immediately preceding and following values. Using Ryu-powered charconv:

C:\Temp>type ryu.cpp
#include <charconv>
#include <limits>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <string>
#include <system_error>
using namespace std;

string print_ryu(const double val) {
    char buf[100];
    const auto to_result = to_chars(buf, end(buf), val, chars_format::scientific);
    if (to_result.ec != errc{}) {
        puts("FAIL");
        exit(EXIT_FAILURE);
    } else {
        return string(buf, to_result.ptr - buf);
    }
}

int main() {
    constexpr auto infinity = numeric_limits<double>::infinity();
    for (int i = -323; i <= 308; ++i) {
        char buf[100];
        sprintf(buf, "1e%+03d", i);
        double d2;
        const auto from_result = from_chars(buf, buf + strlen(buf), d2, chars_format::scientific);
        double d1 = nextafter(d2, -infinity);
        double d3 = nextafter(d2, infinity);
        const string s1 = print_ryu(d1);
        const string s2 = print_ryu(d2);
        const string s3 = print_ryu(d3);
        if (s1 == s2 || s1 == s3 || s2 == s3 || s2 != buf) {
            puts("FAIL");
            exit(EXIT_FAILURE);
        }
        if (i == 137) {
            puts(s1.c_str());
            puts(s2.c_str());
            puts(s3.c_str());
        }
    }
    puts("DONE");
}

C:\Temp>cl /EHsc /nologo /W4 /std:c++17 /D_CRT_SECURE_NO_WARNINGS ryu.cpp && ryu
ryu.cpp
9.999999999999998e+136
1e+137
1.0000000000000002e+137
DONE

[ulfjack]

Let the number we search for be f. f must be slightly less than a power of 10, let's say 1eX, i.e., f < 1eX. It also needs to be closer to 9e(X-1) than to 1eX, i.e., f - 9e(X-1) < 1eX - f. Finally, the next larger floating point number f' must be larger than 1eX, or more precisely, the half-way point between f and f' must be larger than 1eX, i.e., (f+f')/2 > 1eX.

This gives us a constraint for f'-f:
(f+f+f'-f)/2 > 1eX <=> (f'-f)/2 > 1eX - f > f - 9e(X-1)

[ulfjack]

Let me try again. That isn't quite accurate. f can be less than 9e(X-1).

If we look at the number line, we can look at the interval that rounds to f, which is from half-way to the next smaller number to half-way to the next larger number. This interval must contain both 1eX and 9e(X-1) in order to have a conflicting choice here. In addition, f must be less than the half-way point between 1eX and 9e(X-1) or 1eX would be closer to f than any small integer multiple of e(X-1) (such as 9e(X-1)).

The half-way point to the next larger number is (f + f')/2 as before. Therefore, the first constraint is (f+f')/2 > 1eX. Our second constraint is f < (1eX+9e(X-1))/2 (note that 1eX+9e(X-1))/2 = 9.5e(X-1).

It holds that (f+f')/2 = f+(f'-f)/2, therefore (f'-f)/2 > 1eX - f. Making f larger doesn't change the inequality, so we replace it using the second constraint to obtain (f'-f)/2 > 1e(X-1)/2 (because 1eX-(1eX+9e(X-1))/2 = (1eX-9e(X-1))/2) and f'-f > 1e(X-1).

We now deconstruct f as f=m*2^e:
f = m*2^e and f' = (m+1)*2^e. Then f'-f = 2^e.

We obtain 2^e > 1e(X-1) and m*2^e < 9.5e(X-1).

Now we obtain a constraint on m by taking m*2^e from the second inequality, and dividing by the value for 2^e from the first inequality. This gives m < 9.5e(X-1) / 1e(X-1) < 10.

Next, you have to remember that the IEEE format has an implicit leading 1 in the mantissa. The only case where the mantissa can be less than 10 is for denormal numbers, which do not have an implicit leading 1, that is, the numbers where the stored IEEE exponent is 0 corresponding to an exponent of 1-bias. We therefore only need to check the first ten floating-point numbers.

Let's look at 64-bit numbers (I'm printing with Java, which always prints 2 digits):
0.0
4.9E-324
1.0E-323
1.5E-323
2.0E-323
2.5E-323
3.0E-323
3.5E-323
4.0E-323
4.4E-323

The only candidate here is 1.0E-323.

The exact value (same exponent for exact value, lower and upper half-way points) is 988131291682493088353137585736442744730119605228649528851171365001351014540417503730599672723271984759593129390891435461853313420711879592797549592021563756252601426380622809055691634335697964207377437272113997461446100012774818307129968774624946794546339230280063430770796148252477131182342053317113373536374079120621249863890543182984910658610913088802254960259419999083863978818160833126649049514295738029453560318710477223100269607052986944038758053621421498340666445368950667144166486387218476578691673612021202301233961950615668455463665849580996504946155275185449574931216955640746893939906729403594535543517025132110239826300978220290207572547633450191167477946719798732961988232841140527418055848553508913045817507736501283943653106689453125000.

The lower half-way point is 741098468761869816264853189302332058547589703921487146638378523751013260905313127797949754542453988569694847043168576596389985065533909694598162194016172817189451069785467106791768725751773473155533077954085498096084575009581113730347476580968710095909754422710047573078097111189357848386756539987835030152280559340465937397917907387238682993958184816601691220194564999312897984113620624844986787135721803522090170239032857917325202205289740208029068540216066123755499834026713000358124864790413857434018755209015901725925471462961751341597749387185747378709616456389087181198412716730560170454930047052695901657637768849082679869725733665217655679410725087643375608460039849049721491174630855395563541886415131684784363130802375962957739830017089843750.

The upper half-way point is 1235164114603116360441421982170553430912649506535811911063964206251688768175521879663249590904089980949491411738614294327316641775889849490996936990026954695315751782975778511319614542919622455259221796590142496826807625015968522883912460968281183493182924037850079288463495185315596413977927566646391716920467598900776562329863178978731138323263641361002818700324274998854829973522701041408311311892869672536816950398388096528875337008816233680048447567026776872925833056711188333930208107984023095723364592015026502876542452438269585569329582311976245631182694093981811968664021194550933617424883411754493169429396281415137799782876222775362759465684541812738959347433399748416202485291051425659272569810691886141307271884670626604929566383361816406250.

After cutting off most of the digits, we get:

d-= 74
d  = 98
d+=123

Here, 1e-323 is closer to the exact value than 9e-324.

For 32-bit float, the only candidate is m=7, which prints as 1e-44, and, again, this is closer to the exact value than 9e-45. The exact value is 98090892502737194964661070830294129189618335931356104022994779872285375788010242104064673185348510742187500.

[ulfjack]

It's possible that there are other formats where this makes a difference, e.g., for 16, 80, or 128 bits. We have shown that there are at most 10 potential candidates in any case. I think I can additionally show that there's at most one candidate that's close to a power of 10 in the first ten floating-point numbers.