Reasoning about decimal digits of floats

[triska]

Dear all,

trying to address #2771, I came up with the following code to relate a float to decimal digits represented as a list of characters:

:- use_module(library(dcgs)).

float_with_n_decimal_digits(F, N0, Cs) :-
        I is truncate(F),
        number_chars(I, Is),
        N1 is N0 + 1,
        phrase((seq(Is),".",float_decimal_digits(F,1,N1)), Cs).

float_decimal_digits(_, N, N) --> !.
float_decimal_digits(F, N0, N) -->
        { D is abs(truncate(F*(10^N0))) mod 10,
          number_chars(D, Cs),
          N1 is N0 + 1 },
        seq(Cs),
        float_decimal_digits(F, N1, N).

Thanks to the great recent improvements to truncate/1 and related features by @adri326 in #2777, this seems to solve the concrete issue nicely:

?- float_with_n_decimal_digits(1.0e-6, 3, Cs).
   Cs = "0.000". % expected

It also works nicely for several other cases, such as:

?- float_with_n_decimal_digits(1.1, 100, Cs).
   Cs = "1.10000000000000006 ...". % expected
?- float_with_n_decimal_digits(-1.1, 100, Cs).
   Cs = "-1.1000000000000000 ...". % (probably?) expected
?- float_with_n_decimal_digits(5.9999, 60, Ds).
    Ds = "5.99990000000000000 ...". % expected

However, the logic unfortunately falls short in cases like the following:

?- float_with_n_decimal_digits(4.1, 10, Cs).
   Cs = "4.1900990990", unexpected.

I would greatly appreciate any help with this!

Thank you a lot,
Markus

[hurufu]

This happens because 4.1 is represented as 4.099... (as per https://www.h-schmidt.net/FloatConverter/IEEE754.html) and it messes-up with an algorithm.

[hurufu]

But floats are beast to work with, I have no idea if it can introduce yet another errors – so the "ideal" algorithm should operate directly on binary representation of a given float – but we don't have such functionality in Scryer.

[triska]

There is a related issue and discussion about this feature and its relation to formatting: #1343

[hurufu]

When I was talking about "introducing new errors" I have meant that any arithmetic operation will inevitably introduce some errors and if you will do repeated computation (like 10^X or X+Y and so on) errors will just add-up and it can lead to incorrect decimal expansion when dealing with high precision formatting. Also you need to make sure that you never work with floats less then epsilon.

[hurufu]

offtopic: Floats are evil, for example if square operation is implemented naively then you can find such small positive a that this expression will be true: 0.01 + a < 0.1² 🤯

[hurufu]

offtopic part 2: Just don't mess with floats :) Reasoning about them is hard and nobody needs to deal with both very small and very large numbers in a same computation – using integers and multiplier is a solutions, and is much easier to deal with.

[UWN]

float_with_n_decimal_digits2(F, N) -->
   {  Fr is abs(float_fractional_part(F)),
      FrR0 is round(Fr*10^N),
      I0 is truncate(F),
      (  FrR0 >= 10^N
      -> I is I0+truncate(sign(F)), FrR = FrR0
      ;  I = I0, FrR is FrR0+10^N
      )
   },
   ( { I=0, F<0, FrR>10^N } -> "-0" ; { number_chars(I, Is) }, seq(Is) ),
   ".",
   ( { FrR = 1 } -> "0" ; { number_chars(FrR, ['1'|FrRs]) }, seq(FrRs) ).

[UWN]

This works precisely only somewhere below 23 digits:

?- length(_,N), D1 is round(10**N), D1 =\= 10^N.
   N = 23, D1 = 100000000000000008388608
;  N = 24, D1 = 999999999999999983222784
;  ... .

[triska]

Thank you a lot! The following surprises me, since 0.5 can be exactly represented in base 2 via 2^(-1):

?- phrase(float_with_n_decimal_digits(0.5, 100), Cs) ; false.
   Cs = ...
   Cs = "0.5000000000000000079514455548799590234180404281972640694890663778873919386085190530406734992928407552"
;  false.

Where do this additional values greater than 7.0e-18 come from?

[notoria]

There are floats for which it doesn't work:

?- repeat,
   random(F0),
   phrase(float_with_n_decimal_digits2(F0, 24), Cs),
   number_chars(F1, Cs),
   F0 =\= F1.
   F0 = 0.30884494382772054, Cs = "0.30884494382772050 ...", F1 = 0.3088449438277205
   F0 = 0.30884494382772054, Cs = "0.308844943827720504082432", F1 = 0.3088449438277205

[adri326]

0.5 * 10 ^ 100 > 2 ^ 53, so the floating point representation rounds the number to some power of 2 (here 2 ^ 99)

[UWN]

This works precisely only somewhere below 23 digits:

?- length(_,N), D1 is round(10**N), D1 =\= 10^N.
   N = 23, D1 = 100000000000000008388608
;  N = 24, D1 = 999999999999999983222784
;  ... .

Please reread this! All your examples are with N larger than 23.

[notoria]

The conversion from E notation could be done in Rust.

:- use_module(library(format)).
:- use_module(library(dif)).
:- use_module(library(reif)).
:- use_module(library(random)).
:- use_module(library(iso_ext)).
:- use_module(library(dcgs)).
:- use_module(library(lists)).

float_with_n_decimal_digits(F, N, Cs) :-
    number_chars(F, Cs0),
    if_(memberd_t(e, Cs0), from_e_notation(Cs0, Cs1), Cs1 = Cs0),
    phrase((seq(Ds0), ".", seq(Ds1)), Cs1),
    length(Ds1, N1),
    compare(R, N1, N),
    last_digits(R, N1, N, Ds0, Ds1, Cs1, Cs).

last_digits(<, N0, N, _, _, Cs0, Cs) :-
    N1 is N-N0,
    length(Ds, N1),
    append(Ds, _, ['0'|Ds]),
    append(Cs0, Ds, Cs).
last_digits(=, _, _, _, _, Cs, Cs).
last_digits(>, _, N, Ds0, Ds1, _, Cs) :-
    length(Ds2, N),
    append(Ds2, _, Ds1),
    zero(Ds2, Ds),
    append(Ds0, ['.'|Ds], Cs).

zero([], "0").
zero([C|Cs], [C|Cs]).

from_e_notation(Cs0, Cs) :-
    phrase((seq(Ds), ".", seq(Ds1), ("e";"E"), ("+";[]),seq(Es)), Cs0),
    [E0|_] = Es, dif(E0, '+'),
    number_chars(_, Ds1),
    number_chars(E, Es),
    [D0|_] = Ds,
    if_(D0 = (-), [_|Ds0] = Ds, Ds0 = Ds),
    number_chars(_, Ds0),
    reverse(Ds0, Rs0),
    digits(E, Rs0, Ds1, Cs1),
    if_(D0 = (-), Cs = [D0|Cs1], Cs = Cs1).

digits(E, Ds0, Ds1, Cs) :-
    compare(R, 0, E),
    '@digits'(R, E, Ds0, Ds1, Cs).


digit('0', [], []).
digit(C, [C|Cs], Cs).

'@digits'(<, E0, Ds0, Ds1, Cs) :-
    E is E0-1,
    digit(D, Ds1, Ds),
    digits(E, [D|Ds0], Ds, Cs).
'@digits'(=, _, Ds0, Ds1, Cs) :-
    digit(D0, Ds0, Ds2),
    reverse([D0|Ds2], Rs0),
    digit(D1, Ds1, Ds3),
    Rs1 = [D1|Ds3],
    append(Rs0, ['.'|Rs1], Cs).
'@digits'(>, E0, Ds0, Ds1, Cs) :-
    E is E0+1,
    digit(D, Ds0, Ds),
    digits(E, Ds, [D|Ds1], Cs).

[UWN]

Reconsider rounding:

?- between(0,200,B),F is 2/3,phrase(float_with_n_decimal_digits2(F,B),Cs),float_with_n_decimal_digits(F,B,Cs2), Cs\=Cs2.
   B = 0, F = 0.6666666666666666, Cs = "1.0", Cs2 = "0.0"
;  B = 1, F = 0.6666666666666666, Cs = "0.7", Cs2 = "0.6"
;  B = 2, F = 0.6666666666666666, Cs = "0.67", Cs2 = "0.66"
;  B = 3, F = 0.6666666666666666, Cs = "0.667", Cs2 = "0.666"
;  B = 4, F = 0.6666666666666666, Cs = "0.6667", Cs2 = "0.6666"
;  ... .

[notoria]

With round/1:

float_with_n_decimal_digits(F0, N, Cs) :-
    round_float(N, F0, F),
    number_chars(F, Cs0),
    if_(memberd_t(e, Cs0), from_e_notation(Cs0, Cs1), Cs1 = Cs0),
    phrase((seq(Ds0), ".", seq(Ds1)), Cs1),
    length(Ds1, N1),
    compare(R, N1, N),
    last_digits(R, N1, N, Ds0, Ds1, Cs1, Cs).

round_float(_, F0, F) :-
    abs(F0) > 2**53,
    float(F0),
    !,
    F = F0.
round_float(N, F0, F) :-
    N > 15,
    integer(N),
    !,
    F = F0.
round_float(N, F0, F) :-
    N >= 0,
    integer(N),
    !,
    F is round(F0*10**N)/10**N.

Maybe there is an input for which it's incorrect. An alternative rounding method is to round in the string format.

[adri326]

Here's my attempt at it:

get_two_exponent(N10, N2) :-
    %    10 ^ (3*N10) ~= 2 ^ (10*N2)
    % -> 10 ^ N10     ~= 2 ^ (10*N2/3)
    % The + 5 is for good measure, and there's otherwise no need to go above 53 bits.
    N2 is min(ceiling(N10) * 10 // 3 + 5, 53).

format_float3(F, N) -->
    {
        Int is truncate(F),
        Frac is abs(float_fractional_part(F)),
        get_two_exponent(N, N2),
        TwoPow is 2 ^ N2,
        FracInt is floor(Frac * TwoPow),
        N1 is N + 1,
        fractional_part_to_int(FracInt, TwoPow, N1, FracInt10),
        FracIntRounded is round(FracInt10 rdiv 10)
    },
    ( { Int = 0, F < 0, FracInt > TwoPow } -> "-0" ; { number_chars(Int, Is) }, seq(Is) ),
    ".",
    format_float_fractional_part(FracIntRounded, N).

fractional_part_to_int(F, B, N, Out) :-
    N > 0,
    N1 is N - 1,
    F10 is F * 10,
    NextF is F10 mod B,
    Digit is F10 // B,
    fractional_part_to_int(NextF, B, N1, Prev),
    Out is Prev + Digit * 10 ^ N1.
fractional_part_to_int(_, _, 0, 0).

format_float_fractional_part(I, N) -->
    {
        N > 0,
        N1 is N - 1,
        Digit is I // 10 ^ N1,
        NextI is I mod 10 ^ N1,
        number_chars(Digit, Ds)
    },
    seq(Ds),
    format_float_fractional_part(NextI, N1).
format_float_fractional_part(_, 0) --> [].

[hurufu]

Also we need to decide first what do we mean by "correct" string representation. I would argue that the important property is reading it back should yield exactly the same number up to rounding factor. Another also "correct" representation would be actual value that is stored.

[hurufu]

I've actually googled it a little bit an as I initially thought printing floats correctly is very nontrivial task, there is an algorithm "Dragon4" and its variants that is proven to be correct though.

Handling NaN will be hard, and also there are different registers in FPU of you processors that can control how floats behave, so there is no even a universal agreement on what is float: https://www.intel.com/content/dam/develop/external/us/en/documents/fp-control-2013-03.pdf

[UWN]

Before going too much into such details, first check the value 5.0e-323 with 1080 digits which in a precise implementation should result in this value: #2779 (comment)