Julia: inv(::Irrational) errors

[MasonProtter]

Current behaviour:

julia> inv(π)
ERROR: MethodError: no method matching Irrational{:π}(::Int64)
Closest candidates are:
  Irrational{:π}(::T<:Number) where T<:Number at boot.jl:741
  Irrational{:π}() where sym at irrationals.jl:18
  Irrational{:π}(::Complex) where T<:Real at complex.jl:37
  ...
Stacktrace:
 [1] convert(::Type{Irrational{:π}}, ::Int64) at ./number.jl:7
 [2] one(::Type{Irrational{:π}}) at ./number.jl:297
 [3] one(::Irrational{:π}) at ./number.jl:298
 [4] inv(::Irrational{:π}) at ./number.jl:220
 [5] top-level scope at none:0

Probably the ideal behaviour should be

inv(x::Irrational) = 1/x

[StefanKarpinski]

I suppose there might be cases where 1/x is less accurate than doing the inversion at infinite precision and then rounding to floating-point. We could do a compile-time BigFloat computation to try to avoid that. But maybe that's just too annoying.

[MasonProtter]

Looks like there is no accuracy advantage to doing a compile time big float inverse of at least the irrationals in base other than for the golden ratio.

using Base.MathConstants, BenchmarkTools

macro compiletimeinv(T)
    Float64(1/(big(eval(T))))
end

Base.inv(::Irrational{:π}) = @compiletimeinv(π)
Base.inv(::Irrational{:ℯ}) = @compiletimeinv(ℯ)
Base.inv(::Irrational{:φ}) = @compiletimeinv(φ)
Base.inv(::Irrational{:γ}) = @compiletimeinv(γ)
Base.inv(::Irrational{:catalan}) = @compiletimeinv(catalan)


julia> for x in (π, ℯ, γ, φ, catalan)
           println("\n #----------------")   
           println(x)
           @show inv(x)
           @show 1/x
           @show inv(x) - 1/x
       end

 #----------------
π = 3.1415926535897...
inv(x) = 0.3183098861837907
1 / x = 0.3183098861837907
inv(x) - 1 / x = 0.0

 #----------------
ℯ = 2.7182818284590...
inv(x) = 0.36787944117144233
1 / x = 0.36787944117144233
inv(x) - 1 / x = 0.0

 #----------------
γ = 0.5772156649015...
inv(x) = 1.7324547146006335
1 / x = 1.7324547146006335
inv(x) - 1 / x = 0.0

 #----------------
φ = 1.6180339887498...
inv(x) = 0.6180339887498949
1 / x = 0.6180339887498948
inv(x) - 1 / x = 1.1102230246251565e-16

 #----------------
catalan = 0.9159655941772...
inv(x) = 1.0917440637039062
1 / x = 1.0917440637039062
inv(x) - 1 / x = 0.0

[stevengj]

I agree that inv(x::AbstractIrrational) = 1/x is a reasonable fallback definition. We could also default to inv(Float64(x)), but that seems to be basically equivalent.

Individual Irrational constants can always override this if needed. (I don't think its worth it to get one more ulp for the golden ratio — why single out inv to get extra precision?)

[StefanKarpinski]

why single out inv to get extra precision?

It's not about singling out inv it's just that we try to be as exact as we can.

[stevengj]

It's not about singling out inv it's just that we try to be as exact as we can.

We could also define exactly rounded versions of exp, sqrt, and many other math functions, but we don't. So in that sense we are singling out inv here if we try to make it exactly rounded.

I guess your argument is that if we have to explicitly define a method for irrational types, then we usually try to be as accurate as we can. But in fact we don't always do that. For example, we don't guarantee exact rounding for pi + e or similar, even though we have a specialized + method for two irrationals.

I think inv(x::AbstractIrrational) = 1/x is about as accurate as any other operation on irrationals, and it is not worth it to try to guarantee an exactly rounded result in general.

[arnavs]

I opened PR https://github.com/JuliaLang/julia/pull/31068 to add the following to base/irrationals.jl:

# inv
inv(x::AbstractIrrational) = 1/x
@test inv(π) ≈ 0.3183098861837907
@inferred inv(π)

In terms of accuracy, I agree that it would be ideal for us to implement this in a full-precision way. But, I think defaulting to 1/x is perfectly reasonable/intuitive behavior. I also think that if people need extremely high accuracy, they'll probably be taking their own precautions for it.