Julia: typeof(2*pi*1f0)

[goszlanyi]

I was recently bitten by this behaviour:

julia> typeof(2*1f0)
Float32

julia> typeof(pi*1f0)
Float32

julia> typeof(2*pi*1f0)
Float64

julia> typeof(2*1f0*pi)
Float32

I wonder if Float64 is according to design or is it a bug?

[yuyichao]

It's by design and there doesn't seem to be anything we can do about it.

[goszlanyi]

OK, but I would also appreciate a bit more explanation.

Should I always read the types from left to right,
so that typeof(2*pi) is already Float64 before multiplying it by 1f0?

[yuyichao]

Yes

[goszlanyi]

Although I understand the current design,
please consider the following logic:

julia> typeof(pi)
Irrational{:π}

and it behaves in a special/chameleon way in type conversions.

julia> typeof(pi*1.0)
Float64

julia> typeof(pi*1f0)
Float32

julia> typeof(pi*1)
Float64

Then 2*pi could be just an other irrational,
that would also behave in a special/chamelon way in type conversions.

Then exp(2pi*x) would not cause any surprises
for those who payed attention for defining the type of x,
but did not perform the less readable reordering exp(2*x*pi).

[stevengj]

In principle, we could make a special type struct Pi{T} <: Real; x::T; end that represents x*pi (analogous to the UniformScaling), so that expressions like 2*pi have an exact representation.

Working with integer (and rational) multiples of pi is common enough that this might be worthwhile.

pi*Int could return a Pi{Int}, and pi/Int could return a Pi{Rational{Int}}, whereas Pi*float could just return a float.

(Note that this could also be done in a package: Julia is powerful enough that this doesn't require any alterations to the core language. But it might arguably be more useful in Base, since it's the sort of thing where people wouldn't even think to go looking for a package.)

[GunnarFarneback]

Pet peeve: promote_rule(::Type{Float64}, ::Type{Float32}) = Float32 would also solve this problem.

[stevengj]

Note: it looks like it would be a bit easier to define a Pi type if we defined an AbstractIrrational type, since many of the methods are the same as for Irrational.

[stevengj]

Here is a partial implementation just to give the flavor. As I said above, it would be a lot easier if we defined AbstractIrrational in Base:

import Base: promote_rule, float, convert, ==, hash, big, show, *, /, \
struct Pi{T<:Union{Integer,Rational{<:Integer}}} <: Real
    x::T
end
show(io::IO, p::Pi) = print(io, "$(p.x)π = ", string(p.x*float(π))[1:15], "...")
promote_rule(::Type{<:Pi}, ::Type{Float16}) = Float16
promote_rule(::Type{<:Pi}, ::Type{Float32}) = Float32
promote_rule(::Type{<:Pi}, ::Type{T}) where {T<:Number} = promote_type(Float64, T)
float(::Type{<:Pi}) = Float64
convert(::Type{AbstractFloat}, x::Pi) = Float64(x.x)*Float64(π)
convert(::Type{T}, x::Pi) where {T<:AbstractFloat} = T(x.x*T(π))
convert(::Type{Float16}, x::Pi) = Float16(Float32(x))
convert(::Type{Complex{T}}, x::Pi) where {T<:Real} = convert(Complex{T}, convert(T,x))
convert(::Type{Rational{T}}, x::Pi) where T<:Integer = Rational{T}(x.x) * convert(Rational{T}, π)
big(x::Pi) = convert(BigFloat, x)
big(::Type{<:Pi}) = BigFloat
==(x::Pi, y::Pi) = x.x == y.x
==(x::Pi, ::Irrational{:π}) = x.x == 1
==(::Irrational{:π}, x::Pi) = x.x == 1
==(x::Pi, y::Real) = false
==(x::Real, y::Pi) = false
hash(x::Pi, h::UInt) = hash(x.x, hash(pi, h))
*(x::Pi, y::Union{Integer,Rational{<:Integer}}) = Pi(x.x * y)
*(y::Union{Integer,Rational{<:Integer}}, x::Pi) = x * y
*(::Irrational{:π}, x::Union{Integer,Rational{<:Integer}}) = Pi(x)
*(x::Union{Integer,Rational{<:Integer}}, ::Irrational{:π}) = Pi(x)
/(x::Pi, y::Union{Integer,Rational{<:Integer}}) = Pi(x.x // y)
\(y::Union{Integer,Rational{<:Integer}}, x::Pi) = x / y
/(::Irrational{:π}, x::Union{Integer,Rational{<:Integer}}) = Pi(1//x)
\(x::Union{Integer,Rational{<:Integer}}, ::Irrational{:π}) = Pi(1//x)

With these definitions, you can do e.g.:

julia> 2pi
2π = 6.2831853071795...

julia> 2pi * 0.5
3.141592653589793

julia> 2pi * big(0.5) - big(pi)
0.000000000000000000000000000000000000000000000000000000000000000000000000000000

[StefanKarpinski]

I've occasionally thought about adding a scalar multiple field for Irrational. Limiting it to integer multiples makes sense but rationals could be useful too, as in (2//3)π. The question is how to preserve full precision without doing BigFloat computations at runtime, which seems like it would require generated functions.

[stevengj]

@StefanKarpinski, I don't see the problem with

convert(::Type{T}, x::Pi) where {T<:AbstractFloat} = T(x.x*T(π))

This won't be exactly rounded, it's true, but that's not the point — it's close enough. The goal is to be able to represent 2π (and 2π/3 etc.) in a way that expands to (nearly) the precision of its operands, not to do "exact" computations on π.

(It might also be worthwhile to implement e.g sincos(π * rational) more accurately — this is possible to do quickly, without going to BigFloat, and can be quite useful. e.g. we did this inside FFTW. And of course sincos(π * integer) can be computed exactly.)

[stevengj]

Rather than adding a "multiple" field to Irrational, I think it would make more sense to define

struct IrrationalMultiple{T<:Union{Integer,Rational{<:Integer}}, S<:AbstractIrrational} <: AbstractIrrational
    multiple::T
    irrational::S
end

that wraps around another irrational type. (Note that for a singleton type like Irrational{sym}, the irrational field will take no space.)