Julia: How to handle the precision of `BigFloat`s

[dpsanders]

Currently, the precision of BigFloats is determined by a global variable stored in an array DEFAULT_PRECISION, which is manipulated by set_bigfloat_precision.

This is not very Julian. I have been thinking about two possibilities:

1. The precision is given explicitly, e.g. as a second argument to the various BigFloat constructors and to the big function, e.g.

a = BigFloat(3.1, 100)
b = big(3.1, 200)

Here, though, the precision is still hidden inside the object.

1. An arguably more Julian approach, which I would favour, is that the precision be a parameter of the BigFloat{prec} type, so that we could write

a = BigFloat{100}(3.1)
b = BigFloat{200}(3.1)

This version would have the advantage that operations on BigFloats could explicitly track the precision of the objects being operated on (which MPFR explicitly states that it does not do). E.g., a + b in this example should have only 100 bits. (If I understand correctly, with MPFR there is the possibility to specify any precision for the result, but bits beyond 100 will be incorrect.)

If there is a consensus about whether either of these would be useful, or am I missing a reason why this would not be possible?

c.c. @simonbyrne, @andrioni, @lbenet

[stevengj]

Don't most computer-algebra systems with arbitrary-precision arithmetic use a global-like concept to control default precision? Maple has Digits, Mathematica has N[...] (sort of analogous to with_precision), Maxima has fpprec, mpmath has a global precision mp ...

I think you are confusing precision with accuracy here. Just because a number has n bits of precision doesn't mean that it has n accurate bits—it could have more, or it could have less. (e.g. if I compute 3! in floating-point arithmetic I will get an exact result, 6.0.) And just because the operands have n bits doesn't mean that you might not want more (or fewer) digits for the result of floating-point operations. (In fact, it is quite common to use arbitrary precision in problems that have numerical-stability problems, so that intermediate calculations need much more precision than the inputs in order to get any accurate digits in the result.) So, I'm suspicious (in an arbitrary-precision context) of using the precision of the operands to control the precision of the results.

cc: @jiahao

[StefanKarpinski]

I've considered this before for a somewhat different reason: we could potentially arrange for BigFloats to be immediate objects if their size is part of their type.

[JeffBezanson]

My guess is that the cost of recompiling everything for every different precision is not worth it.

[jiahao]

I vaguely recall from my first exposure to Julia that having the precision set globally was a conscious decision to discourage mixing BigFloats with different precisions in the same code.

This version would have the advantage that operations on BigFloats could explicitly track the precision of the objects being operated on (which MPFR explicitly states that it does not do

I have not used MPFR directly myself, but the sample program suggests otherwise. We see explicitly that each call to mpfr_init2 sets a working precision for each given variable.

If we're going to do this, we would want the precision as a type parameter. Otherwise any sort of matrix operation on BigFloats is going to have to worry about the possibility of mixing precisions of various elements at runtime. Ick.

And just because the operands have n bits doesn't mean that you might not want more (or fewer) digits for the result of floating-point operations. (In fact, it is quite common to use arbitrary precision in problems that have numerical-stability problems, so that intermediate calculations need much more precision than the inputs in order to get any accurate digits in the result.)

I think MPFR does try to guarantee that the output of a floating-point computation is correct (roughly) to the working precision of the least precise variable. Their technical paper says: (doi:10.1145/1236463.1236468, Sec. 2.3)

The semantics chosen in MPFR is the following: for each assignment a ← b  c or a ← f (b, c), the variables a, b, c may all have different precisions; the inputs b, c are considered to their full precision, and a correct rounding to the full target precision of a is computed.

[ivarne]

You'd only have to recompile anything if you for some reason don't like the default number of bits. In the cases you do change something, we'll only recompile the code you actually use with the different precision.

[simonbyrne]

I don't have a problem with the first option: operations with different precisions already work correctly, this would just be equivalent to

x = with_bigfloat_precision(200) do
    BigFloat(3.1)
end

I think it's only worth including the precision in the type parameter if there is some sort of performance advantage, which at the moment there doesn't seem to be.

As @stevengj says, it is important not to conflate precision with accuracy (as Mathematica does): see, for example, slide 32 onwards of this presentation by Kahan.

[ivarne]

BigFloat(3.1) is just such a ugly expression, because too few programmers realize that 3.1 is a Float64, and gets rounded before the BigFloat conversion.

[JeffBezanson]

Given how slow BigFloats are likely to be anyway, I doubt there's much to gain by specializing their representation per precision. The underlying C code is already designed for run-time flexibility in precision, which we should be thankful for --- it could have been a C++ template library.

[dpsanders]

Yes I see that I have been confusing precision and accuracy, thanks.

But both @jiahao 's MPFR example, and @stevengj 's Mathematica example seem to show that there is not in fact always a global precision defined, and the precision is specified for each variable independently.

This could be dealt with via an _optional_ precision argument for the BigFloat constructors, which could take as a default value the globally-defined precision.

For the type of application that @stevengj discussed, I had envisaged explicitly changing the precision of variables, something like

x = BigFloat("3.1", 100)
...
x = BigFloat(x, 200)     

so that the mathematical functions would then automatically use the input precision as the output precision.
Though I admit that if there are many variables involved, that could be painful (though alleviated with a macro) and the global precision seems to win in terms of ease of use; I guess in the end it's rather a question of taste.

I see that for now there is no advantage, and apparently some significant disadvantage, in having the precision as a type parameter.

@StefanKarpinski I don't know what you mean by "immediate object"? I presume something to do with allocation?

@ivarne You are right about BigFloat(3.1); I had forgotten about that subtletly since I wrapped it all up in a macro in the ValidatedNumerics package, so that you can do @interval(3.1) and get the correct thing (though it is subtle, as pointed out in #4278).

[mmaechler]

Hmm, getting back to the beginning and the above from @dpsanders , I still think that for big(), it is very UN-natural to not provide an optional precision (in bits) as second argument; think e.g., of big(pi)... and then why not also for BigFloat(.) ?
I'm the author/maintainer of the Rmpfr package (the R <--> MPFR interface) and hence have some experience with these numbers and their application. In Julia you have the nicety to have BigFloat as part of the language, which I find very attractive.
I think you currently cripple functionality by not allowing different BigFloats to have their own precision each. Also, think of saving Julia objects and loading them. There, the bigfloats' precision has to exist independently of a global default precision.
Note that the underlying MPFR library does entirely support in all its functions to _mix_ bigfloats with different precisions... it's all well defined ... and Rmpfr does support it too...

Also, in my experience there are quite a bit of speed differences between, say using 128 bit or 1024 bit precision.... and both are useful/ important to have available.

[ScottPJones]

:+1: :100: to the @dpsanders' to the approach with the precision as part of the type.
This seems to be another place, like with unums, where Julia could support big floats better than many other languages.

[stevengj]

@mmaechler, I agree that it would be good to support a precision argument in the BigFloat constructor. However, I don't think that the precision of the operands should determine the precision of the results; the precision of the results should still be controlled by the global precision setting.

[mmaechler]

@stevengj : I think Julia should follow what MPFR does itself. MPFR is a very mature project with _many_ years of growth and consolidation, and all their arithmetic (and special) functions basically use the principle that the precision of the operands determines the precision of the result.

AFAIU, Julia strives to be a functional language (wherever sensible). Global settings such as the bigfloat precision are actually quite contrary to functional programming in the above sense.

[nalimilan]

Indeed, in general global options appear to be frowned upon by core Julia developers (and rightfully so).

[simonbyrne]

The precision of the result in MPFR is determined by the precision of the variable created to store the result: the precision of the operands is irrelevant.

[StefanKarpinski]

I think that @simonbyrne is correct: MPFR never decides how much precision a result should have – that's always determined by the precision of the first argument which is the return value, mpfr_t rop. So MPFR itself doesn't actually provide any guidance here. I'm not a huge fan of the global setting.

I do like the idea of having the precision as part of the type, a la BigFloat{256}(pi). I know this worries @JeffBezanson because it might lead to too many types, but I think that in practice, people will only use a few precisions for BigFloat. It's a little tricky to make the d field of BigFloat a tuple and pass the precision through, it seems to require a second type parameter that is always equal to ceil(Int,p/8).

[mmaechler]

@simonbyrne : You are right about MPFR, I was wrong; but I was lead by the very typical use; notably in cases as ours: You translate a BigFloat to an mpfr_t, say mx; now it is natural and well supported by MPFR, to use _the_ same argument, e.g., mpfr_sin (mx, mx, MPFR_RNDN) and then translate the result back to BigFloat in Julia. Of course, if your mx uses the default precision instead of the precision of the Julia BigFloat, the result will also be in the default precision;...
For two arguments, for similar reasons it makes sense to use e.g. mpfr_add(mx, mx, my, MPFR_RNDN) where you''d make sure that mx has the larger of the two precisions (if they differ).

[dpsanders]

From http://www.mpfr.org/mpfr-current/mpfr.html:

MPFR does not keep track of the accuracy of a computation. This is left to the user or to a higher layer

This is what I had in mind in one of my previous comments; indeed I was getting confused between accuracy and precision.

Maybe Julia could be the layer in which the _accuracy_ (number of _correct_ bits) is tracked: for example, multiplying a number in which 100 bits are guaranteed correct (known from previous calculations) with one with 200 correct bits should (presumably) result in a number with 100 bits correct. However, this sounds like it could be difficult.

[ivarne]

Sub ussues:

1. Parameterize BigFloat with regard to precision.
2. Make the precision of the results of of BigFloat operations depend on the operands instead of a global variable.
3. Add an optional argument to specify the result precision of some BigFloat operations. The pattern should probably extended beyond constructors.
4. (Create a Julia layer to keep track of accuracy)

Initially I loved the idea of 1. but it has some major drawbacks.

1. It creates a new performance trap with the new possibilities for writing type unstable code.
2. We will waste time recompiling MPFR wrappers, while the the backend calculations will be done in non-specialized C code. (If some functions can be translated to Julia, we could exploit recompilation for efficiency though)
3. Jeff worries that some programs will use so many different precision that it will cause strain on the type system.

The benefit of 1. would of course be that Julia could allocate storage for BigFloat objects inline on the stack, in arrays and structs.

[ViralBShah]

I feel that the burden on the type system and compiler is unlikely to be too much for most use cases, where one may have at best, a few different precisions.

[andrioni]

I once tried to parametrize the precision and avoid using a global variable for this, but I never managed to solve some issues in a satisfactory manner, like the precision of the returned value.

A small example is trying to do arithmetic when you want the result with a higher precision than either the arguments, say:

set_bigfloat_precision(2)
x = big(1.)
y = big(3.)
with_bigfloat_precision(()-> x / y, 1024)

The pre-merge version of BigFloats had the precision as a type parameter, and one of the first commits tried to allow each function to set the precision of the results, while using either the input or a generic value as the default in other cases, but it made very hard to reason about any code.

[StefanKarpinski]

I really feel that we should experiment with having the precision of a BigFloat be part of its type. In practice, programs don't use many different precisions. This could make reasoning about how big the result of an operation should be easy. It would also potentially allow stack allocation of BigFloats.

[andrioni]

The areas I'm most familiar with, computational number theory and interval arithmetic, relied heavily on changing the precision dynamically to obtain sufficient bounds (this paper, for example, dynamically adjusts the precision to isolate polynomial roots efficiently).

That said, I like the idea!

[stevengj]

Making the precision of the result depend on the precision of the operands seems wrong to me. Reasoning about how big the result of an operation from the operands is _not_ easy — arithmetic precision is a global (or dynamically scoped, ala with_precision) property of a computation, not a local property of individual primitive operations.

Each variable can (and already does) have its own precision, but that precision should _not_ determine the precision of the result of arithmetic with that variable.

See also my comments above — we should take the hint here from essentially _every other language_ that supports arbitrary precision well. It's not an accident that they all use dynamically scoped and/or global precision.

[stevengj]

(And if you're not using the precision of the operands to determine the precision of the result, what is to be gained by parameterizing the BigFloat type on the precision?)

[stevengj]

(Note that number theory is _not_ a good guide here. When you do BigInt arithmetic for number theory, you're trying to make every operation _exact_. That's not the case in floating-point arithmetic, even in arbitrary precision.)

[JeffreySarnoff]

_in contradistinction to the overwhelming tide of shared opinion on housebreaking BigFloat precision_

I have made and remade some precision imbued number types. The initial versions kept a current global precision setting and maintained the precisions of type realizations either directly as a field or indirectly as associated dict entries. As the rewrites made better type expression, the extra clutter of managing the precision as an internal externality nagged and distracted. I still keep a current precision for each type that may take different precisions.

Now type definitions include a precision parameter. This requires

convert{P,Q}(::Type{Numelle{P}}, x::Numelle{Q}) #  that handles P>Q and Q>P correctly  
convert{P}(::Type{Numellle{P}}, x::Numelle{P}) = x # to avoid calling convert{P,P}
promote_rule{P,Q}(::Type{Numelle{P}, ::Type{Numelle{Q}) = P>Q ? Numelle{P} : Numelle{Q}

• Asking for x = Numelle{128}(2.0) sets x to 2.0 with 128 bits of precision.
• Asking for x = Numelle(2.0) sets x to 2.0 with precision(Numelle) bits of precision
  _the type's current precision_ bitsofprecision = 64; setprecision(Numelle, bitsofprecision)
• Asking for y = x sets y to 2.0 with precision(x) bits of precision
• Asking for y = Numelle{256}(x) sets y to expand(2.0 from 128 to 256 bits) with 256 bits of precision
• Asking for y = Numelle{32}(x) sets y to round(2.0 from 256 bits to 32 bits) with 32 bits of precision

Parameterizing the precision has led to shorter, crisper and clearer source code. Even where there is cause to set the precision explicitly. Another benefit is the ease of doubling the working precision:

import Base: *
(*){P, Q}(a::Numelle{P}, b::Numelle{Q}) = (*)(promote(a, b)...)

function (*){P}(a::Numelle{P}, b::Numelle{P}, c::Numelle{P})
     z = Numelle{2*P}(a)
     z = z * b
     z = z * c
    return Numellle{P}(z)
end

[mmaechler]

As someone not routinely working with Julia (but quite some experience in the subject),
I congratulate on your decision and user interface implementation... I think this is -- in spirit, as I'm not fluent enough in Julia -- a very nice interface and indeed very close to what I was proposing / thinking about..

[JeffreySarnoff]

@mmaechler good of you to say so, appreciated

[StefanKarpinski]

I like the API. Should we explore this approach in a package and see how it pans out?

[JeffreySarnoff]

@StefanKarpinski Are you using the indirect singular form of the communal we?

[JeffreySarnoff]

What is the reasoning behind the decision that arithmetic with BigFloats of different precision results in a BigFloat with the smaller of the two precisions? While I see a reasoned argument that one might prefer addition and subtraction to resolve to the lesser precision, the more compelling argument (imo) is that multiplication and division (and sqrt) resolve to the greater precision. And it seems reasonable that all arithmetic resolve mixed precisions the same way.

[StefanKarpinski]

Are you using the indirect singular form of the communal we?

Yes, I probably won't be involved except as a commentator. Too much to do for 1.0...

[JeffreySarnoff]

point of clarification:

# v0.4

type BigFloat <: AbstractFloat
    prec::Clong
    sign::Cint
    exp::Clong
    d::Ptr{Limb}
    ...
end

# proposed
type BigFloat{P} <: AbstractFloat
    sign::Cint
    exp::Clong
    d::Ptr{Limb}
    ...
end

precision{P}(::Type{:BigFloat{P}}) = P
precision{P}(x::BigFloat{P}) = P

# with the proposed type def,
#   what is the distinction, apart from brevity, between
function add{P}(a::BigFloat{P}, b::BigFloat{P}) p=P; ... end
# and
function add{T<:BigFloat}(a::T, b::T) p=precision(T); ... end

[JeffreySarnoff]

To try BigFloat{Precision} visit https://github.com/JuliaNumberTypes/MPFR3.jl

If you run into an issue -- please note it (best with the solution).
If you want to see this become into Julia, gotta help some.

[ettersi]

Another issue with a global precision parameter is that numerical functions like sqrt, exp and friends are no longer pure:

julia> x = big(π)/3;
       setprecision(()->sin(x), 2) - setprecision(()->sin(x),3)
-1.250[...]e-01

That is actually quite a bummer. I would like to write formulae like in a textbook and trust on the compiler to eliminate common subexpressions for me.

[StefanKarpinski]

Came up again recently: https://discourse.julialang.org/t/parametric-bits/26117/10. Would changing operations on BigFloats to produce results with the max precision of their arguments be too breaking? It seems quite bad that the size of results is controlled by a global setting. It's not even dynamically scoped or per-thread, just straight up global.

[JeffreySarnoff]

^I agree, and do not find it breaking .. rather fixing.

[perrutquist]

If the result is the max of precision of the arguments _and_ the setprecision global state, then it would only be "breaking" in the sense of sometimes making things slower and more exact.

Something that might otherwise break is re-running a calculation with a higher setprecision in order to compute the truncation error. For example:

setprecision(BigFloat, 256)
a = rand(BigFloat)
b = rand(BigFloat)
c = a*b
setprecision(BigFloat, 512)
d = a*b - c

Here it would be very confusing if d suddenly became zero.

[StefanKarpinski]

Including the global precision in the decision of the result precision feels like a half-measure 😐

[JeffreySarnoff]

At 256 bits, the global precision will dominate the precision of the arg[s] too often to provide a more meaningful computational strategy. From the abacus' perspective, there is at least one (unary ops) and often at least two (binary ops) contributing precisions .. and another precision may have an important role. Of these, none is the global precision setting. There are the precision[s] of the arg[s] and, sometimes another precision relative to those for increasing/decreasing the operational precision at any given, opportune place in a calculation's flow.

[perrutquist]

The user could opt in to only using the precision of the args by doing setprecision(0) so that it never affects the maximum. (This could be done from startup.jl by those who so desire.)

One could also make it an error to setprecision to a nonzero number after it has been set to zero, to prevent accidental re-activation of the old behavior, and to make it easier to find and correct code that relies on it. (Introducing this error is unlikely to break anything, as there is no reason to set the precision to zero other than to opt in to the new behavior.)

It should still be possible to use setprecision(::Function, ...) (typically a do construct) when a different precision is desired for a specific calculation. (Ideally, this should be implemented in a thread-safe way, of course.) Probably, this should override the precision of the args.

[dpsanders]

Bump.

This has now become pressing since BigFloats are not thread-safe due to their global precision state (as far as I can see).

We definitely need a version of BigFloat parametrized by the precision.

[stevengj]

Whether the type is parameterized by precision is orthogonal to this issue. Even for an unparameterized BigFloat type, you could use the precision of the operands to set the precision of the result. However, that would be a hugely breaking change to apply by default and I don't see any way to do it in 1.x.

To make it less breaking, we would need any new precision behavior to be opt-in, e.g. by defining a new BigFloat-like type.

[simonbyrne]

This has now become pressing since _BigFloats are not thread-safe_ due to their global precision state (as far as I can see).

Presumably we can just change DEFAULT_PRECISION and ROUNDING_MODE to be Threads.Atomic values?

[simonbyrne]

Note that it looks like MPFR is not compiled to be thread-safe:

julia> ccall((:mpfr_buildopt_tls_p, :libmpfr), Cint, ())
0

See https://www.mpfr.org/mpfr-current/mpfr.html#Memory-Handling for how their threading works.

[JeffreySarnoff]

from the MPFR doc "Note: when MPFR is built with the ‘--enable-thread-safe’ configure option, the default precision is local to each thread." So it could be (and that change would be considered nonbreaking, I assume).

[ViralBShah]

from the MPFR doc "Note: when MPFR is built with the ‘--enable-thread-safe’ configure option, the default precision is local to each thread." So it could be (and that change would be considered nonbreaking, I assume).

Should probably open a new issue for this. Will have to be done in Yggdrasil and in julia deps. Perhaps we should also consider enabling --enable-gmp-internals. Are the float128 and decimal float capabilities also worth enabling in the build?

[stevengj]

Decimal float support in MPFR is only _Decimal64, right? That seems superseded by the DecFP.jl package.

[JeffreySarnoff]

Yes. DecFP.jl is more robust, offering 32, 64 and 128 bit decimal floats.

[simonbyrne]

from the MPFR doc "Note: when MPFR is built with the ‘--enable-thread-safe’ configure option, the default precision is local to each thread." So it could be (and that change would be considered nonbreaking, I assume).

I don't think we use the MPFR default precision or rounding modes, as we use the custom interface:
https://github.com/JuliaLang/julia/blob/ddf7ce9a595b0c84fbed1a42e8c987a9fdcddaac/base/mpfr.jl#L113
so I don't think this will be breaking. We store the default precisions and rounding modes on the Julia side.

Note that using atomics + the above fix won't fix https://github.com/JuliaIntervals/IntervalArithmetic.jl/issues/376: that would require much larger breaking changes (e.g. using task local storage, or some other sort of change)

[JeffreySarnoff]

Is Compile MPFR to be thread safe relevant to IntervalArithmetic.jl#376?

MPFR internal data such as flags, the exponent range, the default precision and rounding mode,and caches (i.e., data that are not accessed via parameters) are either global (if MPFR has not been compiled as thread safe) or per-thread (thread local storage, TLS).

[ViralBShah]

Certainly sounds like it.