Add BigFloat constructors with precision and rounding mode

[dpsanders]

This adds new constructors for BigFloat of the form

BigFloat(x, precision)
BigFloat(x, precision, rounding_mode)

It also changes the behaviour of eps, prevfloat and nextfloat acting on BigFloats so that they are calculated using the precision of x itself. Note that this is one of the ambiguities forced on us by the fact that the precision is not a parameter of the BigFloat type; c.f. #10040

[dpsanders]

I forgot to say that it also changes convert for x that is a BigFloat to return a version of x with the current (global) precision, something that was previously missing.

This was implemented by calling an MPFR, but could alternatively implemented by multiplying x by BigFloat(1).

[dpsanders]

Bump. This is ready for comments. Cc @simonbyrne @pabloferz

[simonbyrne]

cc: @stevengj

• I'm in favour of adding precision and rounding mode to the BigFloat constructor method.
• I'm not in favour of adding it to big (especially since there isn't a big(x::AbstractString) method). I would be okay with adding them to the parse(BigFloat,...) method though.
• Good call on the change to eps, but I think we can improve it further:
  • it should avoid allocating any interim variables
  • it should give correct values for prevfloat(big(Inf))
  • throw an error when the result underflows (e.g. for nextfloat(big(0.0))).

[stevengj]

I'm inclined to agree with @simonbyrne on adding this to BigFloat but not big.

[dpsanders]

@simonbyrne and @stevengj Thanks for the comments.

The idea is to provide a less wordy version via big:

big("3.1", 100, RoundUp) instead of parse(BigFloat, "3.1", 100, RoundUp) is a win in my book, which allows for a nicer interactive experience.

I would also be happy to add big("3.1"). Since big"3" gives a BigInt, the function big(x::AbstractString) can be reserved for producing BigFloats.

(An alternative that I would not like so much would be to add a new function bigf for this purpose.)

[dpsanders]

@simonbyrne:

• eps does not allocate any interim variables as far as I can see:

function eps(x::BigFloat)
    with_bigfloat_precision(precision(x)) do
    nextfloat(x) - x
    end
end

eps(::Type{BigFloat}) = eps(BigFloat(1))

• Currently with this branch:

julia> eps(prevfloat(BigFloat(Inf)))
inf

This seems like a reasonable answer to me. On master it gives nan. This would apparently require a special case in the function for a rather unlikely use case.

• I guess the underflow occurs because mpfr doesn't handle subnormal numbers. I'll add a method.

By the way, I often find myself wanting to throw a DomainError, but also wanting to give an error message. Could we change DomainError to accept a string argument like ArgumentError? Where in the code base is this defined?

[dpsanders]

@simonbyrne I have implemented the suggestion that eps should throw an error when nextfloat(x)-x is zero. EDIT: It now apparently needs to define an interim variable...

In the process, I realised that in fact part of the problem was in nextfloat and prevfloat themselves, which should respect the precision of x.

For these kinds of things, we are in need of a BigFloat constructor that takes an optional precision. I tried to add a keyword argument prec to the first constructor, but this leads to some inference issues apparently. For the moment I am just wrapping a lot of stuff in with_bigfloat_precision, which is ugly.

[simonbyrne]

• I think you should make the BigFloat constructor an inner constructor: that way you can specify the precision directly, avoiding using with_bigfloat_precision.
• The big"" macro can construct a BigInt and BigFloat, depending on the format of the number. If we were to make big(x::AbstractString) have the same behaviour, then it wouldn't be type-stable.
• Having prevfloat(big(Inf)) == Inf doesn't match the behaviour of Float32/Float64
• I don't know enough about the internals, but the reason adding DomainError doesn't have a message is that it would require allocating memory (which would be problematic, as DomainError is used in otherwise allocation-free functions).
• For eps, I was thinking it might be possible to calculate it by calling precision and exponent directly.

[dpsanders]

• Re. the inner constructor: that's what I had in mind, but the problem is that I couldn't find a suitable syntax, since

BigFloat(100)

is already taken for converting 100 to a BigFloat.

I was trying via a keyword argument with

BigFloat(; prec=get_bigfloat_precision())

but so far this has led to some kind of inference recursion and stack overflows that I haven't been able to get rid of.

• My suggestion was that big(x::AbstractString) should always create a BigFloat. An alternative would be to bring back the BigFloat(x::AbstractString) syntax -- it seems to me that there should definitely be some nice, easy syntax for creating a BigFloat from a string; cf. Maybe deprecate BigFloat(f::FloatingPoint) constructor #4278

In any case, I'll move the big stuff to a separate PR.

• For Float64 we have

julia> eps(prevfloat(Inf))
1.99584030953472e292

which seems to me to just be wrong -- Inf seems a better answer, since that's the "distance to the next Float64".

• Thanks, I made an issue about DomainError but indeed there is some subtle problem.
• For calculating eps directly, that's an interesting idea, I'll have a think about it.

[simonbyrne]

Before we merge this, it would be useful to hash out the rules for BigFloat precision: when do we use the precision of the current value, and when do we use the global precision? (and document it all, cf #7535).

[dpsanders]

I agree that this needs thinking about.

For example, to me it seems natural that one-input functions like sin should return a result with the same precision as the input x; however, this may be controversial since this is not how MPFR does it. (MPFR, if I understand correctly, returns a result with the current global precision, as Julia does currently.)

However, for two-input functions like + the answer is less clear, and the current solution is perhaps already simple and clear.

[simonbyrne]

MPFR requires that you specify the precision of every output.

[dpsanders]

Yes, you're right of course, sorry.

[simonbyrne]

Perhaps a useful distinction is whether the function is about the number (e.g. sin) or the format (eps/nextfloat)? In other words, whether or not the function would make sense if applied to a real number as opposed to a floating point number.

[dpsanders]

That's certainly a useful distinction to have in mind.
What is your feeling about what that implies for sin, though?

[simonbyrne]

sin is clearly a "number" function, so should use the global rounding mode.

[simonbyrne]

Some edge cases that might be worth thinking about:

• significand
• frexp
• ldexp

[andrioni]

IIRC, as @simonbyrne said, precision is mostly kept global because MPFR requires the precision of the output in order to do almost any operation. I once tried using the precision of the arguments whenever possible, but it got out of control pretty fast, and ended up being harder to understand than using with_bigfloat_precision every other line.

In this other comment I also talk a bit about related issues.

[andrioni]

You might be doing double rounding in this case, no? The result of the subtraction is calculated using the global precision, and then the result is converted down to precision of x.

[StefanKarpinski]

@andrioni, nice to have you back around here!

[tkelman]

Bump. @dpsanders any interest in revisiting this?

[dpsanders]

@tkelman Yes, I certainly think this should be revisited, thanks for the bump. I won't be able to get to it for the next couple of weeks, though.

[tkelman]

what part of this was not covered by #17217 ?

[dpsanders]

The part where nextfloat(x), prevfloat(x) and eps(x) are changed to work with the precision of x, rather than the current global precision. I still think this is important and should be done, and I'll try to get round to it.

[dpsanders]

All content here is now superseded.