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]

👍 💯 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).