Re: Julia and Ocaml

[Andreas Enge]

TLDR: Solution for julia in the last paragraph.


Am Wed, Apr 12, 2023 at 05:39:23PM +0200 schrieb Simon Tournier:
> About Julia [1], it comes from the update of mpfr and I have no clue how
> to fix that.  Any idea?

Well, mixing base 10 and base 2 calls for trouble... I will try to go
through this, maybe we can get an idea of the cause. But it depends on
the Julia functions parse and string, so it could mean quite some work
to go through their source code.

> Error in testset mpfr:
> Test Failed at /tmp/guix-build-julia-1.8.3.drv-0/julia-1.8.3/test/mpfr.jl:667
>   Expression: string(parse(BigFloat, "0.1")) == "0.10000002"
>    Evaluated: "0.1" == "0.10000002"

This comes from the following code:
    setprecision(21) do
        @test string(parse(BigFloat, "0.1")) == "0.10000002"
...
I assume that means that we are working with 21 bit floating point
precision.
0.1 lies between 1677721*2^-24 and 1677722*2^-24, but is closer to the
latter. So with rounding to nearest, I would expect the result of
parse(BigFloat, "0.1")) to be 1677722*2^-24. Now one would need to know
what the string function does. MPFR, for instance, has a function that
either allows the user to fix the number of digits, or (by passing 0)
outputs a few extra digits so that, when reading the number again into
a variable of 21 bits, one is sure to recover the same number 1677722*2^-24.
Of course here "0.1" would be a good output; but the semantic does not
specify the output, anything sufficiently close to 0.1 would also be
considered correct. Here 0.10000002*2^24 = 1677721.93554432, and it
would also be read in again as 1677722*2^-24.
(Otherwise said, the equivalent functions of parse and string from mpfr
satisfy parse o string == id; instead of checking whether string o parse
outputs a specific value, I think it would have been better to test
whether parse o string o parse == parse on some values).

I looked at the source code of julia, but am not sufficiently familiar with
it to understand where and how parse and string are defined. I think they
might use mpfr_set_str for parse and mpfr_asprintf for string. At least
I see this in base/mpfr.jl:
function string_mpfr(x::BigFloat, fmt::String)
    pc = Ref{Ptr{UInt8}}()
    n = ccall((:mpfr_asprintf,libmpfr), Cint,
              (Ptr{Ptr{UInt8}}, Ptr{UInt8}, Ref{BigFloat}...),
              pc, fmt, x)
...
function _string(x::BigFloat, fmt::String)::String
    isfinite(x) || return string(Float64(x))
    _prettify_bigfloat(string_mpfr(x, fmt))
end
_string(x::BigFloat) = _string(x, "%Re")
_string(x::BigFloat, k::Integer) = _string(x, "%.$(k)Re")

The format string "%Re" means output an mpfr ("R") and use the precision
convention as detailed above ("e").

I tried to mimick this with pure mpfr code; both libraries, the 4.1.0 we
have in Guix master and 4.2.0, read in 0.1 as 1677722*2^-24, then output
it again as 1.0000002e-01. The julia function _prettify_bigfloat from
mpfr.jl probably transforms this into 0.10000002.

So I have no idea why in the test apparently string(parse(BigFloat, "0.1"))
yields "0.1".

Oh, I do; here is a change:
-_string(x::BigFloat) = _string(x, "%.Re")
+_string(x::BigFloat) = _string(x, "%Re")
With the . in the format string, mpfr@4.2 indeed prints out 1e-01,
whereas mpfr@4.1 used to print 1.0000002e-01. (News item from MPFR:
"Bug fixes: In particular, for the formatted output functions (mpfr_printf,
etc.), the case where the precision consists only of a period has been
fixed to be like .0 as specified in the ISO C standard, and the manual has
been corrected and clarified.")

This is part of a larger commit
commit 1e5fdb29f8858f3244f6aff116ee12e4c8247f3a
Author: Simon Byrne <simon.byrne@gmail.com>
Date:   Tue Jan 10 14:52:36 2023 -0800
    update MPFR to 4.2.0 (#48165)
    Co-authored-by: Mosè Giordano <giordano@users.noreply.github.com>

This should be part of the upcoming 1.9.0 release of julia.
You could try to apply this commit in isolation, or update to 1.8.5 and
apply the commit then if it does not apply otherwise.

Good luck!

Andreas