`pi < pi` gives error

[dpsanders]

This seems to me to be a bug.

julia> pi < pi
ERROR: < not defined for MathConst{:π}
 in < at promotion.jl:199

It can be fixed as follows, but this does not seem like a very good idea:

julia> <(a::MathConst, b::MathConst) = float(a) < float(b)
< (generic function with 42 methods)

julia> pi < pi
false

[JeffBezanson]

There is quite a cottage industry of filing issues about MathConst. I think the reason for this is that the only way to fix all such issues is to implement a full symbolic system, which we so far do not intend to do.

[simonbyrne]

Another open question is whether float(pi) == pi (since we have 1//10 != 0.1). But I agree that this is hardly a priority.

[iamed2]

A weirder thing is this:

julia> pi == pi
true

julia> pi <= pi
ERROR: stack overflow
 in <= at promotion.jl:170 (repeats 79996 times)

[StefanKarpinski]

I don't think it's going to be reasonable to do general operations with math consts. We should certainly make the error messages better but there's just no end to this and there's no good reason to want to do any of this.

[pao]

...there's no good reason to want to do any of this.

There's no reason to want to do any of this intentionally, but library code and user code could certainly come together to cause some unintended interactions of this kind.

MathConst seems like it has something in common with C++14 variable templates. What behavior does that specify, and is it possible to emulate if that makes sense?

[cdsousa]

I imagine that this can happen in the rare case that some function is defined like

function f(x)
    clamp(x, 0, π)
end

then

f(π)
ERROR: < not defined for MathConst{:π}

[cdsousa]

Is there any issue about adding

(<){T}(a::MathConst{T}, b::MathConst{T}) = false
(>){T}(a::MathConst{T}, b::MathConst{T}) = false
(<=){T}(a::MathConst{T}, b::MathConst{T}) = true
(>=){T}(a::MathConst{T}, b::MathConst{T}) = true

to Base?

[StefanKarpinski]

No, that seems ok. Then the next question becomes how to compare different MathConsts – assuming that they differ as Float64s seems viable. What about comparing Float64 and MathConst? How would one do that? Widen the Float64 to BigFloat and then compare? What about comparing a BigFloat and a MathConst?

[StefanKarpinski]

One idea is defining nextfloat and prevfloat for MathConst where the exact value is between those two. You would need a type argument, of course.

[jiahao]

@cdsousa's suggestion seems to work fine for the MathConsts in Base:

c = MathConst[π, e, γ, catalan, φ]
for x in c, y in c, ∘ in [<, ≤, == ,> ,≥ ]
    fx, fy = float(x), float(y)
    @assert (x ∘ y) == (fx ∘ fy)
end

[cdsousa]

Mixing Float64 and MathConst also works out of the box

for x in c, y in c, ∘ in [<, ≤, == ,> ,≥ ]
    fx, fy = float(x), float(y)
    @assert (x ∘ y) == (fx ∘ fy) == (x ∘ fy) == (fx ∘ y)
end

[jiahao]

The issue of comparing MathConsts sounds very much like the corresponding issue for Rationals. #2960

[StefanKarpinski]

MathConsts are generally irrational whereas Rationals aren't. Then again, overflow is less of an issue. There are many problems with the current state of affairs:

julia> pi == float(pi)
true

julia> pi == big(pi)
true

julia> float(pi) == big(pi)
false

julia> float(pi) < pi
false

julia> float(pi) < big(pi)
true

I'm starting to think we should have just gone with making pi a function that takes a type argument.

[pao]

The float-to-big comparisons make sense--you're clearly working with different quantities there. In the perspective that MathConst is like a templated variable, the other three also make sense, since the MathConst just takes whichever type is inferred. So I'm not sure what's wrong with any of the six examples you've posed.

[StefanKarpinski]

The first three results show that pi violates transitivity of ==. The last two show that pi is ambivalent about whether it is bigger than float(pi) or not – if you ask it without converting to big it's not, but if you ask for more digits via big then it's bigger. I'm sure that worse violations of inequality can be found.

[StefanKarpinski]

Here's a fun example with the current comparison behavior:

julia> float(e) < big(e)
true

julia> float(e) <= e <= big(e)
true

julia> float(e) < e
false

julia> e < big(e)
false

[pao]

But that assumes that pi is some kind of single object. But this isn't a symbolic system; pi is one of any of a number of floating-point representations. You're not "converting" pi to anything, you're selecting which representation of pi you get, and it's not terribly reasonable to expect it to act any other way given our overall approach to machine representations.

I would accept lack of transparency as to which representation you're going to get is a fair critique--we call them MathConsts, but they're not mathy things.

[StefanKarpinski]

If they're not single values then all of these operations should be errors.

[pao]

And this is why I brought up C++14 variable templates. I'd be interested how many of these kinds of expressions are compilation errors, and how many get inferred instantiations in that environment.

[StefanKarpinski]

I don't know anything about C++14 variable templates so that particular reference was entirely lost on me. The thing is we can make all of these thing behave sensibly, it's just a question of how.

[pao]

A pi example appears in the proposal, though (being C++) all examples in the proposal use explicit template parameters: http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2013/n3651.pdf

[nalimilan]

When comparing a float with a BitFloat, shouldn't the comparison stop after the precision of the float is exhausted? When using a Float64, we don't know whether 1.0 + 0.00000000000000001 < 1.0 + 0.00000000000000002, because they are stored as the same value. So why should 1.0 + 0.00000000000000001 < BigFloat("1.0") + BigFloat("0.000000000000002") return true? It should raise an an InexactError, or false.

Then the same principle could be applied when comparing MathConst with floats.

[simonbyrne]

Here's my proposal for handling comparisons:

	Int	Float	BigInt	BigFl.	Rat{Int}	Rat{BI}	MathC.
1 Int/Uint	intrinsic
2 Float32/64	intrinsic	intrinsic
3 BigInt	promote 3,3	promote 4,3	libcall
4 BigFloat	promote 4,4	promote 4,4	libcall	libcall
5 Rational{Int}	divrem	div, widen cp	divrem	div, cp	widen cp
6 Rational{BigInt}	divrem	div, cp	divrem	div, cp	cp	cp
7 MathConst	prev/next	prev/next	prev/next	cmp. with higher prec	cont. frac.		float

A quick explanation:

• intrinsic: use llvm intrinsic
• libcall: use gmp/mpfr library function
• promote: promote arguments to use other method
• divrem: divide rational to get integer value, break ties with remainder.
• cp: cross product (a/b < c/d <=> a_d < b_c).
• widen: use widen to so multiplication is exact (Int64 -> Int128, etc.)
• prev/next: compare to float below or above as appropriate (this assumes that MathConsts represent irrational numbers)
• cmp. with higher prec.: calculate MathConst in higher precision than current value: technically this won't always be an exact comparison (if the extra bits are all 0), but if we use say 32 extra bits, the chance of this should be small.
• cont. frac.: sequentially compare elements of continued fraction until unequal. Continued fractions of MathConsts can be precomputed to desired precision (enough for Rational{Int128}).
• float: convert to float as proposed at start by @dpsanders

The only one that that this doesn't cover is Rational{BigInt} vs. MathConst. For some MathConsts the continued fraction is a known sequence (e.g. e), so could be used. This doesn't work for all (e.g. pi).

[StefanKarpinski]

@simonbyrne: That looks like a really good proposal (and impressively comprehensive).

@nalimilan: Floating-point numbers aren't intervals – they have precise rational values. The fact that 0.00000000000000001 and 0.00000000000000002 are two different input forms for the same Float64 doesn't mean that either of those decimal values is the rational Float64 value they both map to.

[dpsanders]

The pi < pi error occurred naturally in the interval arithmetic package I am writing when I passed pi as both arguments to a function that checks a < b (i.e. to create a thin interval around pi); I was not just trying to break or explore corner cases of MathConst, as was implied.

I like the MathConst type, but it should certainly be possible in a language designed for numerical computation to reason about pi like this! @simonbyrne 's proposal indeed seems (without having checked the details) a good approach.

[StefanKarpinski]

Sorry, @dpsanders, I didn't mean to imply that you were doing anything stupid. It's just a question of how far down this rabbit hole we want to go. Do we throw an error message when this happens and force the user to pick a type for pi up front, or do we try to make more numeric operations work for math consts? I like @simonbyrne's proposal and support trying this on the premise that it may be enough that people don't end up needing an unbounded amount of numeric operations on math consts.

[kmsquire]

It might be worth considering how FixedPointNumbers.jl and Nemo/Flint fit in with any scheme used for inter-type comparisons.

Of course, details would probably have to be taken care of in the packages themselves, but how easy would it be?

[simonbyrne]

Does anyone have any suggestions for implementing prevfloat/nextfloat for MathConsts? As the macro doesn't have access to the full decimal expansion, the only method seems to be using BigFloats:

prevfloat(x::MathConst) = with_rounding(()->float64(big(x)),BigFloat,RoundDown)
nextfloat(x::MathConst) = with_rounding(()->float64(big(x)),BigFloat,RoundUp)

but it doesn't make sense to make it a function call, as the answer will always be constant. However we can't farm it out to the @math_const macro as big(x) won't have been defined in the macro scope. Any ideas?

[StefanKarpinski]

Staged functions seem like they would help with this.

[simonbyrne]

This should now be fixed.