Norm of Number type I define

General Usage
,
1

I have implemented my own Type of Real (let’s call it My)

Julia’s boast holds true: having defined “all” the operations for My, Julia is “unreasonably” effective at giving me all sorts of matrix algebra operations at zero additional work on my side.

But then I want to compute the LinearAlgebra.norm, and things go sour, although I have implemented all the maths I expect to be needed for an Euclidean norm.

What I do not understand is what do I need to provide to get unreasonnable success.

My attempt below is problematic because

  1. I don’t get a norm
  2. Error messages have pushed me into implementing conversion to Float64. In the real world (as opposed to this MnonWE) this conversion would be lossy, and my norm would be wrong
  3. Error messages have pushed me to implement LinearAlgebra.norm. Here I have lost the"unreasonnable effectivity" of not having to think what packages My is going to compile together with.

Of course I could write my own norm function, but then I loose any LinearAlgebra functionality that needs a norm…

So here are my attempts.

struct My <: Real
    x::Float64
end
Base.:(+)(a::My,b::My)   = My(a.x+b.x)
Base.:(-)(a::My,b::My)   = My(a.x-b.x)
Base.:(*)(a::My,b::My)   = My(a.x*b.x)
Base.:(/)(a::My,b::My)   = My(a.x/b.x)
Base.:(^)(a::My,b::My)   = My(a.x^b.x)
Base.:(^)(a::My,b::Real) = My(a.x^b)
Base.abs(a::My)          = My(abs(a.x))
Base.sqrt(a::My)         = My(sqrt(a.x))
Base.isnan(a::My)        = isnan(a.x)
Base.:(<)(a::My,b::My)   = a.x<b.x
Base.AbstractFloat(a::My) = a.x
Base.convert(::Type{Float64},a::My) = a.x

using LinearAlgebra
LinearAlgebra.norm(a::My) = abs(a)
v = My.([1,2,3])
@show norm(v)
2

You need to define promotion rules. Right now, +(::My, ::My) works, but +(::My, ::AnotherNumberType) and +(::AnotherNumberType, ::My) do not. In your example, norm fails because it calls +(x::Float64, y::My), which has a fallback that promotes its arguments to a common type. To handle this, you need to define what common type that should be. This is done by adding promote_rule methods. One example could be

Base.promote_rule(::Type{My}, ::Type{T}) where {T<:Number} = promote_type(Float64, T)

That way, the type My behaves like a Float64 when promoted.

3

You don’t need to implement conversion to Float64, you only need to implement promotion with Float64. e.g.

Base.promote_rule(::Type{My}, ::Type{Float64}) = My

(The reason for this is that norm needs to promote any type to a floating-point type for accumulation, and for accuracy/overflow reasons asks you to promote to at least double precision.)

You’ll also want:

Base.float(x::My) = x

(My is already a floating-point type … more generally, if you had My{T} then you want to convert to My{float(T)}).

2 Likes
4

Yes, I have a MWE now, thanks to your advice. The norm is of type My, so that resolves my concern about a lossy conversion. I had to implement Base.isnan, Base.inf and as you said Base.float, and no need for overloading LinearAlgebra.norm (good!!!)

I had to make another change though:

struct My <: Number

instead of

struct My <: Real

Now I wonder how my real code is going to like that! :roll_eyes: - let’s find out! … YES :grin:

5

In practice you’d probably want a My type out of promotion. e.g. suppose My is a quaternion type.

More generally, you probably want something like My to be a parameterized type:

struct My{T<:Real} <: Real
    x::T
end

Base.:(+)(a::My,b::My)   = My(a.x+b.x)
Base.:(-)(a::My,b::My)   = My(a.x-b.x)
Base.:(*)(a::My,b::My)   = My(a.x*b.x)
Base.:(/)(a::My,b::My)   = My(a.x/b.x)
Base.:(^)(a::My,b::My)   = My(a.x^b.x)
Base.:(^)(a::My,b::Real) = My(a.x^b)
Base.abs(a::My)          = My(abs(a.x))
Base.sqrt(a::My)         = My(sqrt(a.x))
Base.isnan(a::My)        = isnan(a.x)
Base.isinf(a::My)        = isinf(a.x)
Base.isfinite(a::My)        = isfinite(a.x)
Base.:(<)(a::My,b::My)   = a.x<b.x

Base.promote_rule(::Type{My{T}}, ::Type{S}) where {T,S<:Real} = My{promote_type(Float64, T)}
Base.float(x::My) = My(float(x.x))
Base.convert(::Type{My{T}}, x::Real) where {T} = My(convert(T, x))
Base.convert(::Type{My{T}}, x::My) where {T} = My(convert(T, x.x))

and then you get

julia> using LinearAlgebra

julia> norm(My.([1,2,3]))
My{Float64}(3.7416573867739413)
1 Like
6

Note also that it’s easier to write this sort of code as:

for f in (:abs, :sqrt)
    @eval Base.$f(a::My) = My($f(a))
end
for f in (:isnan, :isinf, :isfinite)
    @eval Base.$f(a::My) = $f(a)
end
1 Like
7

Yes!

My real code now works, thanks to you, and it has all the features you describe: parametric type, eval and convert.

But interestingly, you guys seem to get it to work with My <: Real while I had to move to My<:Number (well, if it works…)

8

To get everything to work, inclusive LinearAlgebra.vecnormalise I ended up overloading Base.isinf, isginite, typemax, typemin, floatmax, floatmin, eps and float.